diff --git a/Documentation/Doc.idx b/Documentation/Doc.idx new file mode 100644 index 000000000..0ef194095 --- /dev/null +++ b/Documentation/Doc.idx @@ -0,0 +1,14 @@ +\indexentry{\texttt{Checkerboard}}{28} +\indexentry{\texttt{Symm}}{30} +\indexentry{\texttt{Symm}}{30} +\indexentry{\texttt{Symm}}{31} +\indexentry{\texttt{Precision Green} }{33} +\indexentry{\texttt{Precision Phase}}{33} +\indexentry{\path{NWrap}}{33} +\indexentry{\texttt{Square} }{83} +\indexentry{\path{Bilayer_square}}{84} +\indexentry{\path{N_leg_ladder}}{84} +\indexentry{\path{Honeycomb}}{84} +\indexentry{\texttt{Bilayer\_honeycomb}}{84} +\indexentry{\texttt{Triangular}}{85} +\indexentry{\texttt{Kagome}}{85} diff --git a/Documentation/doc.pdf b/Documentation/doc.pdf new file mode 100644 index 000000000..3ec88fd31 Binary files /dev/null and b/Documentation/doc.pdf differ diff --git a/Documentation/hmc.tex b/Documentation/hmc.tex index a689fd0af..ba8457844 100644 --- a/Documentation/hmc.tex +++ b/Documentation/hmc.tex @@ -95,3 +95,18 @@ \subsubsection{Hybrid Monte dynamics} \label{sec:hmc} To test the code, we have carried out high precision calculations of the 6-sites Hubbard chain with open boundary conditions (see Sec.~\ref{sec:langevin}) at $\beta t = 4$. Using \textit{only} hybrid molecular dynamics updates we obtain: $ \langle \hat{H} \rangle = -2.81715 \pm 0.00023 $. This compares very well with the discrete field sequential updating run: $ \langle \hat{H} \rangle = -2.81706 \pm 0.00013 $. +Specifically to start the above HMC run, the user again has to define the subroutine \texttt{Ham\_Langevin\_HMC\_S0} in the base Hamiltonian module that computes the derivative of the bosonic part action with respect to the fields +and in the parameter file chose the option: +\begin{lstlisting}[style=fortran] +&VAR_QMC + +...... +Sequential = .F. ! Sequential moves +...... +Delta_t_Langevin_HMC = 0.5 ! Default time step for Langevin and HMC updates. +N_HMC_sweeps = 1 ! # of HMC updates between sequential ones +Leapfrog_steps = 100 ! Number of leapfrog iterations +...... +/ +\end{lstlisting} +For this parameter choice, since \texttt{Sequential = .F.}, the code will carry out only hybrid molecular dynamics updates. The user can of course choose to carry out both sequential and hybrid molecular dynamics updates by setting \texttt{Sequential = .T.} and choosing a finite value for \texttt{N\_HMC\_sweeps}. It is also interesting to note that HMC time step can be chosen much larger than for the Langevin dynamics. The Langevin integrator suffers from systematic errors in the time step whereas the HMC integrator does not. diff --git a/Documentation/langevin.tex b/Documentation/langevin.tex index b6b9517e4..9ce3f9083 100644 --- a/Documentation/langevin.tex +++ b/Documentation/langevin.tex @@ -173,6 +173,25 @@ \subsubsection*{Example - Hubbard chain at half-filling} Controlling Langevin dynamics when the action has logarithmic divergences is a challenge, and it is not a given that the results are satisfactory. For our specific problem we can solve this issue by considering open boundary conditions. Following an argument put forward in \cite{Assaad07}, we can show, using world lines, that the determinant is always positive. In this case the action does not have logarithmic divergences and the Langevin dynamics works beautifully well, see Fig.~\ref{Langevin.fig}. +Practically to start the Langevin dynamic run the user will have to provide the function +\begin{lstlisting}[style=fortran] +Subroutine Ham_Langevin_HMC_S0(Forces_0) + Real (Kind=Kind(0.d0)), Intent(inout), allocatable :: Forces_0(:,:) + ! Forces_0(n,nt) = \frac{\partial S_0}{\partial s_{n,\tau}}. + ! Here n is the index of the field in the operator list and nt is the time slice. +End Subroutine Ham_Langevin_HMC_S0 +\end{lstlisting} +and in the parameter file chose the option: +\begin{lstlisting}[style=fortran] +&VAR_QMC +...... +Langevin = .T. ! Langevin update +Delta_t_Langevin_HMC = 0.01 ! Default time step for Langevin and HMC updates. +Max_Force = 5.0 ! Max Force for Langevin +...... +/ +\end{lstlisting} + \begin{figure}[H] \begin{center} \includegraphics[scale=0.9,clip]{Langevin.pdf} diff --git a/Documentation/updating.tex b/Documentation/updating.tex index 323e5c183..f9d6cd53d 100644 --- a/Documentation/updating.tex +++ b/Documentation/updating.tex @@ -6,41 +6,88 @@ % !TEX root = doc.tex %------------------------------------------------------------ -\subsection{Updating schemes}\label{sec:updating} +\subsection{Sequential updating schemes}\label{sec:updating} %------------------------------------------------------------ % -The program allows for different types of updating schemes, which are described below and summarized in Tab.~\ref{table:Updating_schemes}. With the exception of Langevin dynamics, for a given configuration $C$, we propose a new one, $C'$, with a given probability $T_0(C \rightarrow C')$ and accept it according to the Metropolis-Hastings acceptance-rejection probability, + +The program allows for different types of updating schemes, which are described below and summarized in Tab.~\ref{table:Updating_schemes}. Except for Langevin dynamics, for a given configuration $C$, we propose a new one, $C'$, with a given probability $T_0(C \rightarrow C')$ and accept it according to the Metropolis-Hastings acceptance-rejection probability, \begin{equation} P(C \rightarrow C') = \text{min} \left( 1, \frac{T_0(C' \rightarrow C) W(C')}{T_0(C \rightarrow C') W(C)} \right), \end{equation} -so as to guarantee the stationarity condition. Here, $ W(C) = \left| \Re \left[ e^{-S(C)} \right] \right| $. +to guarantee the stationarity condition. Here, $ W(C) = \left| \Re \left[ e^{-S(C)} \right] \right| $. Predicting how efficient a certain Monte Carlo update scheme will turn out to be for a given simulation is very hard, so one must typically resort to testing to find out which option produces best results. Methods to optimize the acceptance of global moves include Hybrid Monte Carlo \cite{Duane85} as well as self-learning techniques \cite{LiuJ17,Xu17a}. Langevin dynamics stands apart, and as we will see does not depend on the Metropolis-Hastings acceptance-rejection scheme. + + +% \begin{table}[h] +% \begin{tabular}{@{} p{0.22\columnwidth} p{0.11\columnwidth} p{0.61\columnwidth} @{}} +% \toprule +% Updating schemes & Type & Description \\ +% \midrule +% \texttt{Sequential} & \texttt{logical} & (internal variable) If true, the configurations moves through sequential, single spin flips \\ +% \texttt{Propose\_S0} & \texttt{logical} & If true, proposes sequential local moves according to the probability $e^{-S_0}$, where $S_0$ is the free Ising action. This option only works for +% \texttt{type=1} operator where the field corresponds to an Ising variable \\ +% \texttt{Global\_tau\_moves} & \texttt{logical} & Whether to carry out global moves on a single time slice. +% For a given time slice the user can define which part of the operator string is to be computed sequentially. This is specified by the variable \texttt{N\_sequential\_start} and \texttt{N\_sequential\_end}. A number of \texttt{N\_tau\_Global} user-defined global moves on the given time slice will then be carried out \\ +% \texttt{Global\_moves} & \texttt{logical} & If true, allows for global moves in space and time. A user-defined number \texttt{N\_Global} of global moves in space and time will be carried out at the end of each sweep \\ +% \texttt{Langevin} & \texttt{logical} & If true, Langevin dynamics is used exclusively (i.e., can only be used in association with tempering) \\ +% \texttt{Tempering} & Compiling option & Requires MPI and runs the code in a parallel tempering mode, also see Sec.~\ref{Parallel_tempering.sec},~\ref{sec:compilation} \\ +% \bottomrule +% \end{tabular} +% \caption{Variables required to control the updating scheme. Per default the program carries out sequential, single spin-flip sweeps, and logical variables are set to \texttt{.false.}} +% \label{table:Updating_schemes} +% \end{table} + \begin{table}[h] - \begin{tabular}{@{} p{0.22\columnwidth} p{0.11\columnwidth} p{0.61\columnwidth} @{}} + \begin{tabular}{@{} p{0.45\columnwidth} p{0.35\columnwidth} p{0.15\columnwidth} @{}} \toprule - Updating schemes & Type & Description \\ + Updating schemes & Method & Section \\ \midrule - \texttt{Sequential} & \texttt{logical} & (internal variable) If true, the configurations moves through sequential, single spin flips \\ - \texttt{Propose\_S0} & \texttt{logical} & If true, proposes sequential local moves according to the probability $e^{-S_0}$, where $S_0$ is the free Ising action. This option only works for - \texttt{type=1} operator where the field corresponds to an Ising variable \\ - \texttt{Global\_tau\_moves} & \texttt{logical} & Whether to carry out global moves on a single time slice. - For a given time slice the user can define which part of the operator string is to be computed sequentially. This is specified by the variable \texttt{N\_sequential\_start} and \texttt{N\_sequential\_end}. A number of \texttt{N\_tau\_Global} user-defined global moves on the given time slice will then be carried out \\ - \texttt{Global\_moves} & \texttt{logical} & If true, allows for global moves in space and time. A user-defined number \texttt{N\_Global} of global moves in space and time will be carried out at the end of each sweep \\ - \texttt{Langevin} & \texttt{logical} & If true, Langevin dynamics is used exclusively (i.e., can only be used in association with tempering) \\ - \texttt{Tempering} & Compiling option & Requires MPI and runs the code in a parallel tempering mode, also see Sec.~\ref{Parallel_tempering.sec},~\ref{sec:compilation} \\ + Sequential space-time discrete fields & Metropolis, spin flip & \ref{sec:sequential},\ref{sec:S0} \\ + Sequential space-time continuous fields & Metropolis, box move & \ref{sec:sequential_box} \\ + Sequential space-time continuous fields & Metropolis assisted Langevin & \ref{sec:sequential_mala} \\ + Sequential time, global space & Metropolis, User defined & \ref{sec:global_space} \\ + Sequential time, global space & Metropolis assisted Langevin & \ref{sec:global_space_mala} \\ + Sequential space-time, global space & Combinations of moves & \ref{sec:sequential_combination} \\ + Global space-time, discrete & User defined & \ref{sec:global_space_time} \\ + Global space-time, continuous & Parallel tempering & \ref{Parallel_tempering.sec} \\ + Global space-time, continuous & Langevin dynamics & \ref{sec:langevin} \\ + Global space-time, continuous & Hybrid Monte Carlo & \ref{sec:hmc} \\ \bottomrule \end{tabular} - \caption{Variables required to control the updating scheme. Per default the program carries out sequential, single spin-flip sweeps, and logical variables are set to \texttt{.false.}} + \caption{Summary of possible updating schemes. Aside form the parallel tempering schemes, the user can test and compare all these schemes using our implementation of the Hubbard model. Here we + recommend using a one-dimensional lattice with open boundary conditions and a Hubbard-Stratonovitch decoupling that couples to the magnetization. As mention in Sec.~\ref{sec:langevin}, this choice + guarantees the absence of zero in the fermion determinant such that all updating schemes will work out of the box and can be compared. } \label{table:Updating_schemes} \end{table} + +% \begin{table}[h] +% \begin{tabular}{@{} p{0.22\columnwidth} p{0.35\columnwidth} p{0.35\columnwidth} @{}} +% \toprule +% Updating schemes & Variables & Subroutines \\ +% \midrule +% \texttt{Sequential} & & S0 \\ +% \texttt{Propose\_S0} & & S0 \\ +% \texttt{Propose\_MALA} & Delta\_t\_MALA\_sequential, Max\_Force\_MALA\_sequential & S0, Ham\_Langevin\_HMC\_S0\_single \\ +% \texttt{Global\_tau\_moves} & N\_Global\_tau, Nt\_sequential\_start, Nt\_sequential\_end & Global\_move\_tau, Overide\_global\_tau\_sampling\_parameters \\ +% \texttt{Global\_tau\_MALA\_moves} & N\_Global\_tau\_MALA, Nt\_sequential\_start, Nt\_sequential\_end, Delta\_t\_MALA\_global\_tau, MAX\_Force\_MALA\_global\_tau & Global\_MALA\_move\_tau, Overide\_global\_tau\_sampling\_parameters, Ham\_Langevin\_HMC\_S0\_single \\ +% \texttt{Global\_moves} & N\_Global & Global\_move, Delta\_S0\_global \\ +% \texttt{Langevin} & Delta\_t\_Langevin\_HMC, Max\_Force & Ham\_Langevin\_HMC\_S0 \\ +% \texttt{HMC} & N\_HMC\_sweeps, Delta\_t\_Langevin\_HMC, Leapfrog\_steps & Delta\_S0\_global, Ham\_Langevin\_HMC\_S0 \\ +% \texttt{MALA} & N\_MALA\_sweeps, Delta\_t\_MALA\_global, Max\_Force\_MALA\_global & Global\_MALA\_move, Delta\_S0\_global, Ham\_Langevin\_HMC\_S0 \\ +% \texttt{Tempering} & & Delta\_S0\_global \\ +% \bottomrule +% \end{tabular} +% \caption{Variables required to control the updating scheme. Per default the program carries out sequential, single spin-flip sweeps, and logical variables are set to \texttt{.false.}} +% \label{table:Updating_schemes} +% \end{table} % %------------------------------------------------------------ -\subsubsection{Sequential single spin flips} +\subsubsection{Space-time sequential single spin flips: discrete fields} %------------------------------------------------------------ -% +%d \label{sec:sequential} %If the boolean variable \texttt{Sequential\_Sweeps} is set to true, then the The program adopts per default a sequential, single spin-flip strategy. It will visit sequentially each HS field in the space-time operator list and propose a spin flip. Consider the Ising spin $s_{i,\tau}$. By default (\texttt{Propose\_S0=.false.}), we will flip it with probability $1$, such that for this local move the proposal matrix is symmetric. If we are considering the HS field $l_{i,\tau}$ we will propose with probability $1/3$ one of the other three possible fields. For a continuous field, we modify it with a box distribution of width \texttt{Amplitude} centered around the origin. The default value of \texttt{Amplitude} is set to unity. These updating rules are defined in the \texttt{Fields\_mod.F90} module (see Sec.~\ref{sec:fields}). Again, for these local moves, the proposal matrix is symmetric. Hence in all cases we will accept or reject the move according to @@ -54,11 +101,13 @@ \subsubsection{Sequential single spin flips} % %------------------------------------------------------------ -\subsubsection{Sampling of $e^{-S_0}$} +\paragraph{Sampling of $e^{-S_0}$} \label{sec:S0} %------------------------------------------------------------ % -The package can also propose single spin-flip updates according to a non-vanishing free bosonic action $S_0(C)$. This sampling scheme is used if the logical variable \texttt{Propose\_S0} is set to \texttt{.true.}. As mentioned previously, this option only holds for Ising variables. +The package can also propose single spin-flip updates according to a non-vanishing free bosonic +action $S_0(C)$. This sampling scheme is used if the logical variable \texttt{Propose\_S0} is set to \texttt{.true.}. +As mentioned previously, this option only holds for Ising variables. Consider an Ising spin at space-time $i,\tau$ in the configuration $C$. Flipping this spin generates the configuration $C'$ and we propose this move according to \begin{equation} @@ -74,8 +123,72 @@ \subsubsection{Sampling of $e^{-S_0}$} % % +As mentioned above, this sequential updating scheme is space-time is the default and requires the following settings in the + +\begin{lstlisting}[style=fortran,escapechar=\#,breaklines=true] +&VAR_QMC +...... +Sequential = .T. ! Flag for sequential moves +Propose_S0 = .F. ! Whether to propose moves according to e^{-S0} +...... +/ +\end{lstlisting} +%------------------------------------------------------------ + +%------------------------------------------------------------ +\subsubsection{Space-time sequential single spin flips: Continuous fields } +%------------------------------------------------------------ +% +\label{sec:sequential_continuous} +If the field is a continuous real variable, \texttt{type=3}, (See Sec.~\ref{sec:fields}) then the program provides two types of updating schemes. + +\paragraph{Box distribution} +\label{sec:sequential_box} +The first one is the default and consists in proposing a new value for the field according to a box distribution of width \texttt{Amplitude} centered around the origin. This updating scheme is symmetric, and the move is accepted with probability $P(C \rightarrow C') = \text{min} \left( 1, \frac{ W(C')}{W(C)} \right)$. This is the default updating scheme for continuous fields. + +\paragraph{Metropolis assisted Langevin algorithm (MALA)} \label{sec:sequential_mala} +This updating scheme is based on the Langevin equation that we will discuss in Sec.~\ref{sec:langevin}. For a given field $s_{n,\tau}$ (See. Eq.~\ref{Btau.eq}), we propose a new value of the field according to a single step of the discretized Langevin equation just for this field: +\begin{equation} + s_{n,\tau}' = s_{n,\tau} - \frac{\partial S}{\partial s_{n,\tau}} \delta t + \sqrt{ 2 \delta t} \eta +\end{equation} +where $\eta$ is a Gaussian random variable with zero mean and unit variance, and $\delta t$ is the time step of the Langevin dynamics. This updating scheme is not symmetric, and the move is accepted with probability the Metropolis-Hastings acceptance-rejection probability, +\begin{equation}P(C \rightarrow C') = \text{min} \left( 1, \frac{T_0(C' \rightarrow C) W(C')}{T_0(C \rightarrow C') W(C)} \right),\end{equation} +where $T_0(C \rightarrow C')$ is the probability to propose the move $C \rightarrow C'$ and reads: +\begin{equation} + T_0(C \rightarrow C') = \frac{1}{\sqrt{2 \pi \delta t}} e^{-\frac{(s_{n,\tau}' - s_{n,\tau} + \frac{\partial S}{\partial s_{n,\tau}} \delta t)^2}{2 \delta t}}. +\end{equation} +The user can control the time step $\delta t$ by setting the variable \texttt{Delta\_t\_MALA\_sequential} in the input file. The maximum value of the force $\frac{\partial S}{\partial s_{n,\tau}}$ is controlled by the variable \texttt{Max\_Force\_MALA\_sequential}, which is used to avoid proposing moves with very low acceptance. The force is then defined as $\text{sign}(\frac{\partial S}{\partial s_{n,\tau}}) \text{min}(\left|\frac{\partial S}{\partial s_{n,\tau}}\right|, \texttt{Max\_Force\_MALA\_sequential})$. + + +\begin{lstlisting}[style=fortran] +&VAR_QMC +...... +Sequential = .T. ! Flag for sequential moves +Propose_S0 = .F. ! Whether to propose moves according to e^{-S0} +Amplitude = 0.5D0 ! Amplitude of the box distribution +Propose_MALA = .true. ! MALA sequential moves for continuous fields +Max_Force_MALA_sequential = 1.0 ! Maximum value of the force +Delta_t_MALA_sequential = 0.1 ! Time step for the MALA +...... +/ +\end{lstlisting} + +We note that we have two routines to compute the force of the bosonic part of the action, \texttt{S0}. For the Langevin dynamics described in Sec.~\ref{sec:langevin} the force for each field has to be computed. This is done in the routine \texttt{Ham\_Langevin\_HMC\_S0} that returns the two-dimension array +of forces. For the MALA sequential moves, we only need the force for a single field, which is computed in the routine \texttt{Ham\_Langevin\_HMC\_S0\_single}. If this routine is not provided by the user, the code will call \texttt{Ham\_Langevin\_HMC\_S0} and extract the force for the given field. This is computationally expensive, such that a warning is printed on the standard output at execution. The following routine must be provided by the user before this option can be used: +\begin{lstlisting}[style=fortran] +Subroutine Ham_Langevin_HMC_S0(Forces_0) + Real (Kind=Kind(0.d0)), Intent(inout), allocatable :: Forces_0(:,:) + ! Forces_0(n,nt) = \frac{\partial S_0}{\partial s_{n,\tau}}. + ! Here n is the index of the field in the operator list and nt is the time slice. +End Subroutine Ham_Langevin_HMC_S0 +\end{lstlisting} + %------------------------------------------------------------ -\subsubsection{Global updates in space} + +%------------------------------------------------------------ + +\subsubsection{Time sequential, global in space updates.} +\paragraph{User-defined global moves on a single time slice} \label{sec:global_space} This option allows one to carry out user-defined global moves on a single time slice. This option is enabled by setting the logical variable \texttt{Global\_tau\_moves} to \texttt{.true.}. Recall that the propagation over a time step $\Delta \tau$ (see Eq. \ref{Btau.eq}) can be written as: \begin{equation} @@ -84,7 +197,7 @@ \subsubsection{Global updates in space} where $e^{-V_{n}(s_{n})}$ denotes one element of the operator list containing the HS fields. One can provide an interval of indices, \texttt{[Nt\_sequential\_start, Nt\_sequential\_end]}, in which the operators will be updated sequentially. Setting \texttt{Nt\_sequential\_start}$\; = 1$ and \texttt{Nt\_sequential\_end}$\; = M_I+M_V$ reproduces the sequential single spin flip strategy of the above section. -The variable \texttt{N\_tau\_Global} sets the number of global moves carried out on each time slice \texttt{ntau}. Each global move is generated in the routine \texttt{Global\_move\_tau}, which is provided by the user in the Hamiltonian file. In order to define this move, one specifies the following variables: +The variable \texttt{N\_Global\_tau} sets the number of global moves carried out on each time slice \texttt{ntau}. Each global move is generated in the routine \texttt{Global\_move\_tau}, which is provided by the user in the Hamiltonian file. In order to define this move, one specifies the following variables: \begin{itemize} \item \texttt{Flip\_length}: An integer stipulating the number of spins to be flipped. \item \texttt{Flip\_list(1:Flip\_length)}: Integer array containing the indices of the operators to be flipped. @@ -97,13 +210,102 @@ \subsubsection{Global updates in space} Since we allow for a stochastic generation of the global move, it may very well be that no change is proposed. In this case, \texttt{T0\_Proposal\_ratio} takes the value 0 upon exit of the routine \texttt{Global\_move\_tau} and no update is carried out. \item \texttt{S0\_ratio}: Real number containing the ratio $e^{-S_0(C')}/e^{-S_0(C)}$. \end{itemize} +Specifically, the routine \texttt{Global\_move\_tau} must be provided by the user and has the following structure: +\begin{lstlisting}[style=fortran] +Subroutine Global_move_tau(T0_Proposal_ratio, S0_ratio, & + & Flip_list, Flip_length,Flip_value,ntau) + + Real (Kind = Kind(0.d0)), INTENT(OUT) :: T0_Proposal_ratio, S0_ratio + Integer , INTENT(OUT) :: Flip_list(:) + Complex (Kind = Kind(0.d0)),INTENT(OUT) :: Flip_value(:) + Integer, INTENT(OUT) :: Flip_length + Integer, INTENT(IN) :: ntau + +end Subroutine Global_move_tau +\end{lstlisting} +The subroutine \texttt{Overide\_global\_tau\_sampling\_parameters} can be used to overrule the default values of the parameters controlling +the global moves on a single time slice, which are set in the input file. +\begin{lstlisting}[style=fortran] +Subroutine Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end,& + & N_Global_tau, N_Global_tau_MALA) + Integer, Intent(OUT) :: N_Global_tau, Nt_sequential_start, Nt_sequential_end +end Subroutine Overide_global_tau_sampling_parameters +\end{lstlisting} +%------------------------------------------------------------ +As an example, one can run the Hubbard Hamiltonian with the parameters: +\begin{lstlisting}[style=fortran] +&VAR_QMC +...... +Sequential = .T. ! Flag for sequential moves +Global_tau_moves = .T. ! Global moves on a time slice +N_Global_tau = 1 ! Number of global moves will be set + ! to the length of OP_V by the override routine +...... +\end{lstlisting} +Here we randomly chose a vertex and update the field. + %------------------------------------------------------------ +\paragraph{Metropolis assisted Langevin algorithm (MALA) on a single time slice} +\label{sec:global_space_mala} +This option will carry out a global move on a single time slice using the Metropolis assisted Langevin algorithm. To use this option, the + +\begin{lstlisting}[style=fortran] +&VAR_QMC +...... +Sequential = .T. ! Flag for sequential moves +Global_tau_MALA_moves = .T. ! Global MALA moves on a single time slice +N_Global_tau_MALA = 1 ! Number of global MALA moves will be set +Delta_t_MALA_global_tau = 0.1 ! Time step for the MALA global-tau move +Max_Force_MALA_global_tau = 1.0 ! Maximum force for the MALA global-tau move +...... +\end{lstlisting} +The user can override the number of \texttt{N\_Global\_tau\_MALA} by using the aforementioned +\texttt{Overide\_global\_tau\_sampling\_parameters} routine. The routine: + +\begin{lstlisting}[style=fortran] +subroutine Global_MALA_move_tau(Flip_list, Flip_length,ntau) + Integer , INTENT(OUT) :: Flip_list(:) + Integer, INTENT(OUT) :: Flip_length + Integer, INTENT(IN) :: ntau + +end Subroutine Global_MALA_move_tau +\end{lstlisting} +Allows to specify the list of field to be updated in the MALA move. By default, all the field on the time slice are visited and +a warning is issued on the standard output. + %------------------------------------------------------------ -\subsubsection{Global updates in time and space} +\subsubsection{Combining various sequential updating schemes.} +\label{sec:sequential_combination} +The user can combine the various sequential updating schemes described above. For example one can carry out a number of sequential single spin-flip sweeps, followed +by a number of global ones. The range of the sequential sweeps are specified by the variables \texttt{Nt\_sequential\_start} and \texttt{Nt\_sequential\_end}, while the number of global moves is set by \texttt{N\_Global\_tau}. +For example, the parameter file: +\begin{lstlisting}[style=fortran] +&VAR_QMC +...... +Sequential = .T. ! Sequential moves +Propose_S0 = .F. ! Proposes single spin flip moves with probability exp(-S0) +Nt_sequential_start = 1 ! Starting index of the sequential moves on a time slice +Nt_sequential_end = 3 ! Ending index of the sequential moves on a time slice +Propose_MALA = .T. ! MALA sequential moves for continuous fields +Max_Force_MALA_sequential = 1.0 ! Max force for the MALA sequential moves +Delta_t_MALA_sequential = 0.1 ! Time step for the MALA sequential moves +Global_tau_moves = .true. ! Global moves on a time slice +N_Global_tau = 6 ! Number of global moves on a time slice +...... +\end{lstlisting} +Will carry out sequential single spin-flip moves on the first three operators of the operator list using the MALA updating scheme \ref{sec:sequential_mala}, followed by \texttt{N\_Global\_tau=6} global moves on the time slice. + + +\subsection{Global updates in time and space} + + \label{sec:global_slice} %------------------------------------------------------------ % -The code allows for global updates as well. The user must then provide two additional functions (see \texttt{Hamiltonian\_Hubbard\_include.h}): \texttt{Global\_move} and \texttt{Get\_\-Delta\_\-S0\_\-global( Nsigma\_old )}. + +\subsubsection{User specified global moves in space and time} +\label{sec:global_space_time} +The code allows for global updates as well. The user must then provide two additional functions \texttt{Global\_move} and \texttt{Get\_\-Delta\_\-S0\_\-global( Nsigma\_old )}. The subroutine \texttt{Global\_move(T0\_Proposal\_ratio,nsigma\_old,size\_clust)} proposes a global move. Its single input is the variable \texttt{nsigma\_old} of type \texttt{Field} (see Section~\ref{sec:fields}) that contains the full configuration $C$ stored in \texttt{nsigma\_old\%f(M\_V + M\_I, Ltrot)}. On output, the new configuration $C'$, determined by the user, is stored in the two-dimensional array \texttt{nsigma}, which is a global variable declared in the Hamiltonian module. @@ -115,6 +317,32 @@ \subsubsection{Global updates in time and space} \texttt{Get\_\-Delta\_\-S0}\texttt{\_\-global(nsigma\_old)} that computes the difference $\Delta S_0=S_0(C)-S_0(C')$. Again, the configuration $C'$ is given by the field \texttt{nsigma}. The variable \texttt{N\_Global} determines the number of global updates performed per sweep. Note that global updates are expensive, since they require a complete recalculation of the weight. + +Practically one will combine global update with local ones. In the implementation of the Hubbard model, the parameter file: +\begin{lstlisting}[style=fortran] +&VAR_QMC +...... +Sequential = .T. ! Flag for sequential moves +Global_moves = .T. ! Flag for global moves in space and time +N_Global = 6 ! Number of global moves in space and time +...... +/ +\end{lstlisting} +will carry out 6 global moves in space and time at the end of each sweep. A possible choice for the global move is to randomly choose a +lattice point and change the sign of the field at this point for all time slices. The routine that will do this global update reads: +\begin{lstlisting}[style=fortran] +Subroutine Global_move(T0_Proposal_ratio, nsigma_old, size_clust) + Real (Kind = Kind(0.d0)), INTENT(OUT) :: T0_Proposal_ratio, size_clust + Type(Field), INTENT(IN) :: nsigma_old + + size_clust = Ltrot + n = nranf(Size(Op_V,1)) + nsigma%f(n,:) = - nsigma_old%f(n,:) + T0_Proposal_ratio = 1 +end Subroutine Global_move +\end{lstlisting} + + % %------------------------------------------------------------ diff --git a/Prog/Hamiltonian_main_mod.F90 b/Prog/Hamiltonian_main_mod.F90 index 13ec11746..50f5a024f 100644 --- a/Prog/Hamiltonian_main_mod.F90 +++ b/Prog/Hamiltonian_main_mod.F90 @@ -149,13 +149,16 @@ Module Hamiltonian_main procedure, nopass :: Pr_obs => Pr_obs_base procedure, nopass :: Init_obs => Init_obs_base procedure, nopass :: Global_move_tau => Global_move_tau_base + procedure, nopass :: Global_MALA_move_tau => Global_MALA_move_tau_base procedure, nopass :: Hamiltonian_set_nsigma => Hamiltonian_set_nsigma_base procedure, nopass :: Overide_global_tau_sampling_parameters => Overide_global_tau_sampling_parameters_base procedure, nopass :: Global_move => Global_move_base + procedure, nopass :: Global_MALA_move => Global_MALA_move_base procedure, nopass :: Delta_S0_global => Delta_S0_global_base procedure, nopass :: Get_Delta_S0_global => Get_Delta_S0_global_base procedure, nopass :: S0 => S0_base procedure, nopass :: Ham_Langevin_HMC_S0 => Ham_Langevin_HMC_S0_base + procedure, nopass :: Ham_Langevin_HMC_S0_single => Ham_Langevin_HMC_S0_single_base procedure, nopass :: weight_reconstruction => weight_reconstruction_base procedure, nopass :: GR_reconstruction => GR_reconstruction_base procedure, nopass :: GRT_reconstruction => GRT_reconstruction_base @@ -626,6 +629,94 @@ Subroutine Global_move_tau_base(T0_Proposal_ratio, S0_ratio, & CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) end Subroutine Global_move_tau_base +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Specify a global langevin move on a given time slice tau. +!> +!> @details +!> @param[in] ntau Integer +!> \verbatim +!> Time slice +!> \endverbatim +!> @param[out] Flip_length Integer +!> \verbatim +!> Number of flips stored in the first Flip_length entries of the array Flip_values. +!> Has to be smaller than NDIM +!> \endverbatim +!> @param[out] Flip_list Integer(Ndim) +!> \verbatim +!> List of spins to be flipped: nsigma%f(Flip_list(1),ntau) ... nsigma%f(Flip_list(Flip_Length),ntau) +!> Note that Ndim = size(Op_V,1) +!> \endverbatim +!-------------------------------------------------------------------- + Subroutine Global_MALA_move_tau_base(Flip_list, Flip_length,ntau) + + Implicit none + Integer , INTENT(OUT) :: Flip_list(:) + Integer, INTENT(OUT) :: Flip_length + Integer, INTENT(IN) :: ntau + + Logical, save :: first_call=.True. + integer :: n + + Flip_length = size(flip_list,1) + do n = 1, Flip_length + flip_list(n) = n + enddo + + If (first_call) then + write(output_unit,*) + write(output_unit,*) "ATTENTION: Base implementation of Global_MALA_move_tau is being called!" + write(output_unit,*) "All fields of type 3 will be updated in the global tau MALA moves." + write(output_unit,*) "Suppressing further printouts of this message." + write(output_unit,*) + first_call=.false. + endif + + end Subroutine Global_MALA_move_tau_base + +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Specify a global Metropolis-adjusted langevin move. +!> +!> @details +!> @param[out] Flip_length Integer +!> \verbatim +!> Number of flips stored in the first Flip_length entries of the array Flip_values. +!> Has to be smaller than NDIM*Ltrot +!> \endverbatim +!> @param[out] Flip_list Integer(Ndim,Ltrot) +!> \verbatim +!> List of spins to be flipped: nsigma%f(Flip_list(1,1),Flip_list(1,2)) ... nsigma%f(Flip_list(Flip_Length,1),Flip_list(Flip_Length,2)) +!> Note that Ndim = size(Op_V,1) +!> \endverbatim +!-------------------------------------------------------------------- + Subroutine Global_MALA_move_base(Flip_list) + + Implicit none + Integer , INTENT(OUT) :: Flip_list(:,:) + + Logical, save :: first_call=.True. + + If (first_call) then + write(output_unit,*) + write(output_unit,*) "ATTENTION: Base implementation of Global_MALA_move is being called!" + write(output_unit,*) "All fields of type 3 will be updated in the global MALA moves." + write(output_unit,*) "Suppressing further printouts of this message." + write(output_unit,*) + first_call=.false. + endif + + Flip_list = 1 + + end Subroutine Global_MALA_move_base + !-------------------------------------------------------------------- !> @author !> ALF Collaboration @@ -665,10 +756,11 @@ end Subroutine Hamiltonian_set_nsigma_base !> @details !> \endverbatim !-------------------------------------------------------------------- - Subroutine Overide_global_tau_sampling_parameters_base(Nt_sequential_start,Nt_sequential_end,N_Global_tau) + Subroutine Overide_global_tau_sampling_parameters_base(Nt_sequential_start,Nt_sequential_end, & + & N_Global_tau, N_Global_tau_MALA) Implicit none - Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau + Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau, N_Global_tau_MALA !!$ write(output_unit,*) !!$ write(output_unit,*) "ATTENTION: Base implementation of Overide_global_tau_sampling_parameters is getting calling!" @@ -698,9 +790,9 @@ Subroutine Ham_Langevin_HMC_S0_base(Forces_0) if (first_call) then write(output_unit,*) - write(output_unit,*) "ATTENTION: Base implementation of Ham_Langevin_HMC_S0 is getting calling!" + write(output_unit,*) "ATTENTION: Base implementation of Ham_Langevin_HMC_S0 is being called!" write(output_unit,*) "This assumes trivial S0 action and is likely incorrect!" - write(output_unit,*) "Consider overwritting this routine according to the model in your Hamiltonian." + write(output_unit,*) "Consider overwriting this routine according to the model in your Hamiltonian." write(output_unit,*) "Suppressing further printouts of this message." write(output_unit,*) first_call=.False. @@ -709,6 +801,43 @@ Subroutine Ham_Langevin_HMC_S0_base(Forces_0) end Subroutine Ham_Langevin_HMC_S0_base +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Forces_0(n,nt) = \partial S_0 / \partial s(n,nt) are calculated and returned to main program. +!> +!------------------------------------------------------------------- + Subroutine Ham_Langevin_HMC_S0_single_base(Force_0,n,nt) + + Implicit none + + Real (Kind=Kind(0.d0)), Intent(inout) :: Force_0 + integer , intent(in) :: n, nt + + logical, save :: first_call=.True. + real (kind=kind(0.d0)), allocatable :: forces_0(:,:) + + allocate(forces_0(size(nsigma%f,1),size(nsigma%f,2))) + call ham%Ham_Langevin_HMC_S0(Forces_0) + force_0 = forces_0(n,nt) + + if (first_call) then + write(output_unit,*) + write(output_unit,*) "ATTENTION: Base implementation of Ham_Langevin_HMC_S0_single is being called!" + write(output_unit,*) "This uses the subroutine Ham_Langevin_HMC_S0, which provides all bosonic forces." + write(output_unit,*) "Since Ham_Langevin_HMC_S0_single only requires one force this might be inefficient." + write(output_unit,*) "Consider overwriting this routine according to the model in your Hamiltonian." + write(output_unit,*) "Suppressing further printouts of this message." + write(output_unit,*) + first_call=.False. + endif + + deallocate(forces_0) + + end Subroutine Ham_Langevin_HMC_S0_single_base + !-------------------------------------------------------------------- !> @author diff --git a/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 index 0c24334ed..5b3e1b32a 100644 --- a/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Hubbard_smod.F90 @@ -142,6 +142,9 @@ procedure, nopass :: S0 procedure, nopass :: Ham_Langevin_HMC_S0 procedure, nopass :: Get_Delta_S0_global + procedure, nopass :: Global_move_tau + procedure, nopass :: Overide_global_tau_sampling_parameters + procedure, nopass :: Global_move #ifdef HDF5 procedure, nopass :: write_parameters_hdf5 #endif @@ -960,5 +963,130 @@ Real (Kind=kind(0.d0)) Function Get_Delta_S0_global(Nsigma_old) ! S0 = exp( (-Hs_new**2 + nsigma%f(n,nt)**2 ) /2.d0 ) end Function Get_Delta_S0_global - - end submodule ham_Hubbard_smod + + !-------------------------------------------------------------------- + !> @author + !> ALF Collaboration + !> + !> @brief + !> Specify a global move on a given time slice tau. + !> + !> @details + !> @param[in] ntau Integer + !> \verbatim + !> Time slice + !> \endverbatim + !> @param[out] T0_Proposal_ratio, Real + !> \verbatim + !> T0_Proposal_ratio = T0( sigma_new -> sigma ) / T0( sigma -> sigma_new) + !> \endverbatim + !> @param[out] S0_ratio, Real + !> \verbatim + !> S0_ratio = e^( S_0(sigma_new) ) / e^( S_0(sigma) ) + !> \endverbatim + !> @param[out] Flip_length Integer + !> \verbatim + !> Number of flips stored in the first Flip_length entries of the array Flip_values. + !> Has to be smaller than NDIM + !> \endverbatim + !> @param[out] Flip_list Integer(Ndim) + !> \verbatim + !> List of spins to be flipped: nsigma%f(Flip_list(1),ntau) ... nsigma%f(Flip_list(Flip_Length),ntau) + !> Note that Ndim = size(Op_V,1) + !> \endverbatim + !> @param[out] Flip_value Real(Ndim) + !> \verbatim + !> Flip_value(:)= nsigma%flip(Flip_list(:),ntau) + !> Note that Ndim = size(Op_V,1) + !> \endverbatim + !-------------------------------------------------------------------- + Subroutine Global_move_tau(T0_Proposal_ratio, S0_ratio, & + & Flip_list, Flip_length,Flip_value,ntau) + + Implicit none + Real (Kind = Kind(0.d0)), INTENT(OUT) :: T0_Proposal_ratio, S0_ratio + Integer , INTENT(OUT) :: Flip_list(:) + Complex (Kind = Kind(0.d0)),INTENT(OUT) :: Flip_value(:) + Integer, INTENT(OUT) :: Flip_length + Integer, INTENT(IN) :: ntau + + If (.not.Continuous) then + Write(6,*) "Error: Global_move_tau_base is implemented only continuous HS fields. Please implement it or set Continuous = False in the input file. " + CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) + endif + Flip_length = 1 + Flip_list(1) = nranf(Size(Op_V,1)) + Flip_value(1) = nsigma%flip(Flip_list(1),ntau) + T0_Proposal_ratio = 1.d0 + S0_ratio = S0(Flip_list(1),ntau,Flip_value(1)) + + + end Subroutine Global_move_tau + +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> This routine allows to user to determine the global_tau sampling parameters at run time +!> It is especially usefull if these parameters are dependent on other parameters. +!> +!> @details +!> \endverbatim +!-------------------------------------------------------------------- + Subroutine Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end, & + & N_Global_tau, N_Global_tau_MALA) + + Implicit none + Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau, N_Global_tau_MALA + + N_Global_tau = Size(Op_V,1) + + end Subroutine Overide_global_tau_sampling_parameters + +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Global moves +!> +!> @details +!> This routine generates a +!> global update and returns the propability T0_Proposal_ratio = T0( sigma_out-> sigma_in ) / T0( sigma_in -> sigma_out) +!> @param [IN] nsigma_old, Type(Fields) +!> \verbatim +!> Old configuration. The new configuration is stored in nsigma. +!> \endverbatim +!> @param [OUT] T0_Proposal_ratio Real +!> \verbatimam +!> T0_Proposal_ratio = T0( sigma_new -> sigma_old ) / T0( sigma_old -> sigma_new) +!> \endverbatim +!> @param [OUT] Size_clust Real +!> \verbatim +!> Size of cluster that will be flipped. +!> \endverbatim +!------------------------------------------------------------------- + Subroutine Global_move(T0_Proposal_ratio, nsigma_old, size_clust) + + Implicit none + Real (Kind=Kind(0.d0)), intent(out) :: T0_Proposal_ratio, size_clust + Type (Fields), Intent(IN) :: nsigma_old + + Integer :: nt, n + + If (.not.Continuous) then + Write(6,*) "Error: Global_move_tau_base is implemented only continuous HS fields. Please implement it or set Continuous = False in the input file. " + CALL Terminate_on_error(ERROR_HAMILTONIAN,__FILE__,__LINE__) + endif + size_clust = Ltrot + n = nranf(Size(Op_V,1)) + do nt = 1,Ltrot + nsigma%f(n,nt) = -nsigma_old%f(n,nt) + enddo + T0_Proposal_ratio = 1 + + + End Subroutine Global_move + + end submodule ham_Hubbard_smod diff --git a/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 b/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 index c60b8cfc5..0a926cce0 100644 --- a/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_LRC_smod.F90 @@ -823,10 +823,10 @@ end Subroutine Global_move_tau !> @details !> \endverbatim !-------------------------------------------------------------------- - Subroutine Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end,N_Global_tau) + Subroutine Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end,N_Global_tau, N_Global_tau_MALA) Implicit none - Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau + Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau, N_Global_tau_MALA If ( str_to_upper(Model) == "LRC" ) then Nt_sequential_start = 1 diff --git a/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 b/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 index 05a099c26..fa5b653f9 100644 --- a/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Z2_Matter_smod.F90 @@ -1331,10 +1331,10 @@ end Subroutine Hamiltonian_set_Z2_matter !> @details !> \endverbatim !-------------------------------------------------------------------- - Subroutine Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end,N_Global_tau) + Subroutine Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end,N_Global_tau, N_Global_tau_MALA) Implicit none - Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau + Integer, Intent(INOUT) :: Nt_sequential_start,Nt_sequential_end, N_Global_tau, N_Global_tau_MALA Nt_sequential_start = 1 diff --git a/Prog/Langevin_HMC_mod.F90 b/Prog/Langevin_HMC_mod.F90 index e866737f5..fd69fc9f0 100644 --- a/Prog/Langevin_HMC_mod.F90 +++ b/Prog/Langevin_HMC_mod.F90 @@ -53,12 +53,13 @@ Module Langevin_HMC_mod Private - Public :: Langevin_HMC, Langevin_HMC_type, Langevin_HMC_Reset_storage + Public :: Langevin_HMC, Metropolis_Langevin, Langevin_HMC_type, Langevin_HMC_Reset_storage, calculate_Force enum, bind(c) enumerator :: Scheme_none = 0 enumerator :: Scheme_Langevin = 1 enumerator :: Scheme_HMC = 2 + enumerator :: Scheme_MALA = 3 end enum Type Langevin_HMC_type private @@ -85,7 +86,7 @@ Module Langevin_HMC_mod procedure :: calc_Forces => Langevin_HMC_Forces end type Langevin_HMC_type - Type (Langevin_HMC_type) :: Langevin_HMC + Type (Langevin_HMC_type) :: Langevin_HMC, Metropolis_Langevin Contains @@ -115,7 +116,7 @@ SUBROUTINE Langevin_HMC_Forces(this, Phase, GR, GR_Tilde, Test, udvr, udvl, Sta COMPLEX (Kind=Kind(0.d0)), intent(inout), allocatable, dimension(:,:,:) :: GR, GR_Tilde Integer, intent(in), dimension(:), allocatable :: Stab_nt Integer, intent(in) :: LOBS_ST, LOBS_EN - Logical, intent(in) :: Calc_Obser_eq + Logical, intent(in) :: Calc_Obser_eq !Local @@ -354,6 +355,7 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab !Local Integer :: N_op, n, nt, n1, n2, i, j, t_leap, nf, nf_eff Real (Kind=Kind(0.d0)) :: X, Xmax,E_kin_old, E_kin_new,T0_Proposal_ratio, weight, cluster_size + Real (kind=kind(0.d0)) :: Delta_t_running_old, Delta_t_running_new Logical :: Calc_Obser_eq, toggle Real (Kind=Kind(0.d0)), allocatable :: Det_vec_old(:,:), Det_vec_new(:,:) Complex (Kind=Kind(0.d0)), allocatable :: Phase_Det_new(:), Phase_Det_old(:) @@ -361,6 +363,9 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab Type (Fields) :: nsigma_old Character (Len=64) :: storage Complex (Kind=Kind(0.d0)) :: Ratio(2), Phase_old, Ratiotot,Phase_new, Z + Complex (Kind=Kind(0.d0)), allocatable :: Forces_old (:,:) + Real (Kind=Kind(0.d0)), allocatable :: Forces_0_old(:,:) + Integer, allocatable :: Flip_list(:,:) Complex (Kind=Kind(0.d0)) :: Phase_array(N_FL) select case (this%scheme) !(trim(this%Update_scheme)) @@ -576,6 +581,132 @@ SUBROUTINE Langevin_HMC_update(this,Phase, GR, GR_Tilde, Test, udvr, udvl, Stab Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst) !if accepted ! LATER (optimization idea) restore Phase, GR, udvr, udvl and don't reset storage + case(Scheme_MALA) !("MALA") + + n1 = size(nsigma%f,1) + n2 = size(nsigma%f,2) + call nsigma_old%make(n1, n2) + nsigma_old%f = nsigma%f + nsigma_old%t = nsigma%t + Phase_old = Phase + + Allocate ( Phase_Det_new(N_FL_eff), Det_vec_new(NDIM,N_FL_eff)) + Allocate ( Phase_Det_old(N_FL_eff), Det_vec_old(NDIM,N_FL_eff)) + Allocate ( Forces_old(n1,n2) , Forces_0_old(n1,n2) ) + Allocate ( Flip_list(Size(nsigma%f,1),size(nsigma%f,2)) ) + + storage = "Full" + If ( .not. this%L_Forces) then + Call Compute_Fermion_Det(Phase_det_old,Det_Vec_old, udvl, udvst, Stab_nt, storage) + else + Det_vec_old = this%Det_vec_old + Phase_det_old = this%Phase_det_old + endif + + Calc_Obser_eq = .false. + If ( .not. this%L_Forces) then + Call this%calc_Forces(Phase, GR, GR_Tilde, Test, udvr, udvl, Stab_nt, udvst,& + & LOBS_ST, LOBS_EN, Calc_Obser_eq ) + endif + Call ham%Ham_Langevin_HMC_S0( this%Forces_0) + forces_old = this%Forces + forces_0_old = this%Forces_0 + + call ham%Global_MALA_move(Flip_list) + + call Control_MALA_Global(forces_old, forces_0_old, flip_list) + + Xmax = 0.d0 + do n = 1,n1 + do nt = 1,n2 + if ( flip_list(n,nt) == 1 ) then + X = abs(Real(this%Forces (n,nt), Kind(0.d0))) + if (X > Xmax) Xmax = X + X = abs(Real(this%Forces_0(n,nt), Kind(0.d0))) + if (X > Xmax) Xmax = X + endif + enddo + enddo + Delta_t_running_old = this%Delta_t_Langevin_HMC + If (Xmax > this%Max_Force) Delta_t_running_old = this%Max_Force*this%Delta_t_Langevin_HMC/Xmax + + do n = 1, n1 + if (OP_V(n,1)%type == 3 ) then + do nt = 1, n2 + if ( flip_list(n,nt) == 1 ) then + nsigma%f(n,nt) = nsigma%f(n,nt) - ( this%Forces_0(n,nt) + & + & real( Phase*this%Forces(n,nt),kind(0.d0)) / Real(Phase,kind(0.d0)) ) * Delta_t_running_old + & + & sqrt( 2.d0 * Delta_t_running_old) * rang_wrap() + endif + enddo + endif + enddo + + Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst) + Call Compute_Fermion_Det(Phase_det_new,Det_Vec_new, udvl, udvst, Stab_nt, storage) + Call this%calc_Forces(Phase, GR, GR_Tilde, Test, udvr, udvl, Stab_nt, udvst,& + & LOBS_ST, LOBS_EN, Calc_Obser_eq ) + Call ham%Ham_Langevin_HMC_S0( this%Forces_0) + call Control_MALA_Global(this%forces, this%forces_0, flip_list) + + Xmax = 0.d0 + do n = 1,n1 + do nt = 1,n2 + if ( flip_list(n,nt) == 1 ) then + X = abs(Real(this%Forces (n,nt), Kind(0.d0))) + if (X > Xmax) Xmax = X + X = abs(Real(this%Forces_0(n,nt), Kind(0.d0))) + if (X > Xmax) Xmax = X + endif + enddo + enddo + Delta_t_running_new = this%Delta_t_Langevin_HMC + If (Xmax > this%Max_Force) Delta_t_running_new = this%Max_Force*this%Delta_t_Langevin_HMC/Xmax + + t0_proposal_ratio = 1.d0 + do n = 1, n1 + if (OP_V(n,1)%type == 3 ) then + do nt = 1, n2 + if ( flip_list(n,nt) == 1 ) then + t0_proposal_ratio = t0_proposal_ratio * sqrt(Delta_t_running_old/Delta_t_running_new) * & + & exp(-0.25d0/Delta_t_running_new * (Abs(nsigma_old%f(n,nt) - nsigma%f(n,nt) + & + & Delta_t_running_new*(this%forces_0(n,nt) + real( Phase*this%forces(n,nt),kind(0.d0)) & + & / Real(Phase,kind(0.d0))) )**2 ) + 0.25d0/Delta_t_running_old * ( & + & Abs(nsigma%f(n,nt) - nsigma_old%f(n,nt) + Delta_t_running_old*(forces_0_old(n,nt) + & + & real( phase_old*forces_old(n,nt),kind(0.d0)) / Real(phase_old,kind(0.d0))) )**2 ) ) + endif + enddo + endif + enddo + + Ratiotot = Compute_Ratio_Global(Phase_Det_old, Phase_Det_new, & + & Det_vec_old, Det_vec_new, nsigma_old, T0_Proposal_ratio, Ratio) + Weight = abs( real( Phase_old * Ratiotot, kind=Kind(0.d0))/real(Phase_old,kind=Kind(0.d0)) ) + + + TOGGLE = .false. + if ( Weight > ranf_wrap() ) Then + TOGGLE = .true. + this%Det_vec_old = Det_vec_new + this%Phase_Det_old = Phase_det_new + else + this%Det_vec_old = Det_vec_old + this%Phase_Det_old = Phase_det_old + nsigma%t = nsigma_old%t + nsigma%f = nsigma_old%f + endif + + Z = Phase_old * Ratiotot/ABS(Ratiotot) + Call Control_PrecisionP_MALA(Z,Phase) + Call Control_upgrade_MALA(TOGGLE) + + Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst) + + deallocate ( Phase_Det_new, Det_vec_new ) + deallocate ( Phase_Det_old, Det_vec_old ) + deallocate ( Forces_old , Forces_0_old) + deallocate ( Flip_list ) + case default WRITE(error_unit,*) 'Unknown Global_update_scheme ', trim(this%Update_scheme) WRITE(error_unit,*) 'Global_update_scheme is Langevin or HMC' @@ -595,7 +726,7 @@ end SUBROUTINE Langevin_HMC_update !-------------------------------------------------------------------- - SUBROUTINE Langevin_HMC_setup(this,Langevin,HMC, Delta_t_Langevin_HMC, Max_Force, Leapfrog_steps ) + SUBROUTINE Langevin_HMC_setup(this,Langevin,HMC, MALA, Delta_t_Langevin_HMC, Max_Force, Leapfrog_steps ) Implicit none @@ -604,7 +735,7 @@ SUBROUTINE Langevin_HMC_setup(this,Langevin,HMC, Delta_t_Langevin_HMC, Max_Forc class (Langevin_HMC_type) :: this - Logical , Intent(in) :: Langevin, HMC + Logical , Intent(in) :: Langevin, HMC, MALA Integer , Intent(in) :: Leapfrog_steps Real (Kind=Kind(0.d0)), Intent(in) :: Delta_t_Langevin_HMC, Max_Force @@ -681,6 +812,26 @@ SUBROUTINE Langevin_HMC_setup(this,Langevin,HMC, Delta_t_Langevin_HMC, Max_Forc this%Delta_t_running = 1.0d0 ! WRITE(error_unit,*) 'HMC step is not yet implemented' ! CALL Terminate_on_error(ERROR_GENERIC,__FILE__,__LINE__) + elseif (MALA) then + Nr = size(nsigma%f,1) + Nt = size(nsigma%f,2) + Do i = 1, Nr + if ( nsigma%t(i) /= 3 ) then + write(output_unit,*) + WRITE(output_unit,*) 'Warning: Not all fields are of type 3.' + WRITE(output_unit,*) 'Fields that are not of type 3 will not be updated in MALA updates.' + write(output_unit,*) + exit + endif + enddo + Allocate ( this%Forces(Nr,Nt), this%Forces_0(Nr,Nt) ) + Allocate ( this%Det_vec_old(NDIM,N_FL), this%Phase_Det_old(N_FL) ) + this%Update_scheme = "MALA" + this%scheme = Scheme_MALA + this%Delta_t_Langevin_HMC = Delta_t_Langevin_HMC + this%Max_Force = Max_Force + this%L_Forces = .False. + this%Delta_t_running = 1.0d0 else this%Update_scheme = "None" endif @@ -738,6 +889,9 @@ SUBROUTINE Langevin_HMC_clear(this) Deallocate ( this%Det_vec_old, this%Phase_Det_old ) ! WRITE(error_unit,*) 'HMC step is not yet implemented' ! CALL Terminate_on_error(ERROR_GENERIC,__FILE__,__LINE__) + case(Scheme_MALA) + Deallocate ( this%Forces, this%Forces_0 ) + Deallocate ( this%Det_vec_old, this%Phase_Det_old ) case default end select end SUBROUTINE Langevin_HMC_clear @@ -800,12 +954,12 @@ end function Langevin_HMC_get_Update_scheme !> @brief !> Sets the update_scheme !-------------------------------------------------------------------- - subroutine Langevin_HMC_set_Update_scheme(this, Langevin, HMC ) + subroutine Langevin_HMC_set_Update_scheme(this, Langevin, HMC, MALA ) Implicit none class (Langevin_HMC_type) :: this - Logical, intent(in) :: Langevin, HMC + Logical, intent(in) :: Langevin, HMC, MALA If (Langevin) then this%Update_scheme = "Langevin" @@ -813,6 +967,9 @@ subroutine Langevin_HMC_set_Update_scheme(this, Langevin, HMC ) elseif (HMC) then this%Update_scheme = "HMC" this%scheme = Scheme_HMC + elseif (MALA) then + this%Update_scheme = "MALA" + this%scheme = Scheme_MALA else this%Update_scheme = "None" this%scheme = Scheme_none @@ -836,4 +993,47 @@ function Langevin_HMC_get_Delta_t_running(this) Langevin_HMC_get_Delta_t_running = this%Delta_t_running end function Langevin_HMC_get_Delta_t_running +!-------------------------------------------------------------------- +!> @author +!> ALF Collaboration +!> +!> @brief +!> Returns Force for a given operator n and time slice ntau +!-------------------------------------------------------------------- + function calculate_Force(n,ntau,gr) + Implicit none + + complex(Kind=Kind(0.d0)) :: calculate_Force + integer, intent(in) :: n, ntau + Complex (Kind=Kind(0.d0)), intent(in), dimension(:,:,:) :: Gr + + integer :: nf_eff, nf, i, j + Complex (Kind=Kind(0.d0)) :: Z(N_FL), Z1, g_loc, force + + force = cmplx(0.d0, 0.d0, kind(0.d0)) + If (Op_V(n,1)%type == 3) then + Z = cmplx(0.d0,0.d0,Kind(0.d0)) + Do nf_eff = 1, N_Fl_eff + nf=Calc_FL_map(nf_eff) + do I = 1,size(OP_V(n,nf)%P,1) + do J = 1,size(OP_V(n,nf)%P,1) + Z1 = cmplx(0.d0,0.d0,Kind(0.d0)) + if ( I == J ) Z1 = cmplx(1.d0,0.d0,Kind(0.d0)) + Z(nf) = Z(nf) + Op_V(n,nf)%O(I,J) * ( Z1 - Gr(Op_V(n,nf)%P(J),Op_V(n,nf)%P(I), nf) ) + Enddo + Enddo + Z(nf) = Z(nf) + Op_V(n,nf)%alpha + Enddo + if (reconstruction_needed) call ham%weight_reconstruction(Z) + Do nf = 1, N_Fl + g_loc = Op_V(n,nf)%g + if (Op_V(n,nf)%get_g_t_alloc() ) g_loc = Op_V(n,nf)%g_t(ntau) + force = force - g_loc * Z(nf) * cmplx(real(N_SUN,Kind(0.d0)), 0.d0, Kind(0.d0)) + Enddo + endif + + calculate_Force = force + + end function calculate_Force + end Module Langevin_HMC_mod diff --git a/Prog/Wrapgr_mod.F90 b/Prog/Wrapgr_mod.F90 index 5466219a1..2b358fb64 100644 --- a/Prog/Wrapgr_mod.F90 +++ b/Prog/Wrapgr_mod.F90 @@ -52,9 +52,13 @@ Module Wrapgr_mod Use Hamiltonian_main Use Hop_mod use upgrade_mod + use Langevin_HMC_mod Implicit none + INTERFACE wrapgr_sort + MODULE PROCEDURE wrapgr_sort_value, wrapgr_sort_langevin + END INTERFACE !> Privat Complex (Kind=Kind(0.d0)), private, allocatable :: GR_ST(:,:,:) @@ -78,13 +82,15 @@ Subroutine Wrapgr_dealloc end Subroutine Wrapgr_dealloc !-------------------------------------------------------------------- - SUBROUTINE WRAPGRUP(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_end, N_Global_tau) + SUBROUTINE WRAPGRUP(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_end, N_Global_tau, & + & Propose_MALA, Delta_t_MALA_sequential, Max_Force_MALA_sequential, & + & N_Global_tau_MALA, Delta_t_MALA_global_tau, Max_Force_MALA_global_tau) !-------------------------------------------------------------------- !> @author !> ALF-project ! !> @brief -!> Given the green function matrix GR at time NTAU the routine propagates +!> Given the green function matrix GR at time NTAU at m = 0, the routine propagates !> it to time NTAU + 1 and carries out an update of the fields at time NTAU+1 !> NTAU: [0:LTROT-1] ! @@ -95,17 +101,22 @@ SUBROUTINE WRAPGRUP(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT), allocatable :: GR(:,:,:) COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT) :: PHASE INTEGER, INTENT(IN) :: NTAU - LOGICAL, INTENT(IN) :: Propose_S0 - INTEGER, INTENT(IN) :: Nt_sequential_start, Nt_sequential_end, N_Global_tau + LOGICAL, INTENT(IN) :: Propose_S0, Propose_MALA + INTEGER, INTENT(IN) :: Nt_sequential_start, Nt_sequential_end, N_Global_tau, N_Global_tau_MALA + Real (kind=Kind(0.d0)), intent(in) :: Delta_t_MALA_sequential, Delta_t_MALA_global_tau + real (kind=kind(0.d0)), intent(in) :: Max_Force_MALA_sequential, Max_Force_MALA_global_tau !Local Integer :: nf, nf_eff, N_Type, NTAU1,n, m Complex (Kind=Kind(0.d0)) :: Prev_Ratiotot, HS_Field, HS_New + Complex (Kind=Kind(0.d0)) :: force_old, force_new, phase_st, nsigma_st Real (Kind=Kind(0.d0)) :: T0_proposal, T0_Proposal_ratio, S0_ratio + real (kind=kind(0.d0)) :: force_0_old, force_0_new, weight + real (kind=kind(0.d0)) :: delta_t_running_old, Xmax, Delta_t_running_new Character (Len=64) :: Mode Logical :: Acc, toggle1 - - ! Wrap up, upgrade ntau1. with B^{1}(tau1) + + ! Wrap up, upgrade ntau1. with B^{1}(tau1) NTAU1 = NTAU + 1 Do nf_eff = 1,N_FL_eff nf=Calc_Fl_map(nf_eff) @@ -115,14 +126,35 @@ SUBROUTINE WRAPGRUP(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ Do n = Nt_sequential_start,Nt_sequential_end Do nf_eff = 1, N_FL_eff nf=Calc_Fl_map(nf_eff) - HS_Field = nsigma%f(n,ntau1) + HS_Field = nsigma%f(n,ntau1) N_type = 1 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),HS_Field,Ndim,N_Type,ntau1) enddo nf = 1 T0_proposal = 1.5D0 T0_Proposal_ratio = 1.D0 - Hs_New = nsigma%flip(n,ntau1) + if (Propose_MALA .and. Op_V(n,nf)%type == 3) then + nsigma_st = nsigma%f(n,ntau1) + phase_st = Phase + Gr_st = Gr + N_type = 2 + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) + Call Op_Wrapup(Gr_st(:,:,nf),Op_V(n,nf),HS_Field,Ndim,N_Type,ntau1) + enddo + force_old = calculate_force(n,ntau1,Gr_st) + Call ham%Ham_Langevin_HMC_S0_single( force_0_old,n, ntau1) + call Control_MALA_sequential(force_old, force_0_old) + Xmax = abs(dble(force_old)) + if (abs(force_0_old) > Xmax) Xmax = abs(force_0_old) + delta_t_running_old = Delta_t_MALA_sequential + if (Xmax > Max_Force_MALA_sequential) delta_t_running_old = Max_Force_MALA_sequential*Delta_t_MALA_sequential/Xmax + HS_New = nsigma%f(n,ntau1) - ( force_0_old + & + & real( Phase*force_old,kind(0.d0)) / Real(Phase,kind(0.d0)) ) * delta_t_running_old + & + & sqrt( 2.d0 * delta_t_running_old) * rang_wrap() + else + Hs_New = nsigma%flip(n,ntau1) + Endif S0_ratio = ham%S0(n,ntau1, Hs_New ) if ( Propose_S0 ) then If ( Op_V(n,nf)%type == 1) then @@ -131,10 +163,14 @@ SUBROUTINE WRAPGRUP(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ endif Endif If ( T0_proposal > ranf_wrap() ) Then - !Write(6,*) 'Hi', n, Op_V(n,nf)%type, T0_Proposal_ratio, S0_ratio - mode = "Final" + !Write(6,*) 'Hi', n, Op_V(n,nf)%type, T0_Proposal_ratio, S0_ratio + if (Propose_MALA .and. Op_V(n,1)%type == 3) then + mode = "Intermediate" + else + mode = "Final" + endif Prev_Ratiotot = cmplx(1.d0,0.d0,kind(0.d0)) - Call Upgrade2(GR,n,ntau1,PHASE,HS_new, Prev_Ratiotot, S0_ratio,T0_Proposal_ratio, Acc, mode ) + Call Upgrade2(GR,n,ntau1,PHASE,HS_new, Prev_Ratiotot, S0_ratio,T0_Proposal_ratio, Acc, mode ) else toggle1 = .false. Call Control_upgrade_eff(toggle1) @@ -144,26 +180,66 @@ SUBROUTINE WRAPGRUP(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ N_type = 2 Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),HS_Field,Ndim,N_Type,ntau1) enddo + + if (Propose_MALA .and. Op_V(n,1)%type == 3) then + phase = Phase * Prev_Ratiotot/sqrt(Prev_Ratiotot*conjg(Prev_Ratiotot)) + Call ham%Ham_Langevin_HMC_S0_single( force_0_new,n,ntau1) + force_new = calculate_force(n,ntau1,GR) + call Control_MALA_sequential(force_new, force_0_new) + Xmax = abs(dble(force_new)) + if (abs(force_0_new) > Xmax) Xmax = abs(force_0_new) + Delta_t_running_new = Delta_t_MALA_sequential + if (Xmax > Max_Force_MALA_sequential) Delta_t_running_new = Max_Force_MALA_sequential*Delta_t_MALA_sequential/Xmax + + t0_Proposal_ratio = sqrt(delta_t_running_old/Delta_t_running_new) * exp(-0.25d0/Delta_t_running_new * (Abs(nsigma_st - hs_new + & + & Delta_t_running_new*(force_0_new + real( Phase*force_new,kind(0.d0)) / Real(Phase,kind(0.d0))) )**2) + 0.25d0/delta_t_running_old * ( & + & Abs(hs_new - nsigma_st + delta_t_running_old*(force_0_old + real( phase_st*force_old,kind(0.d0)) / Real(phase_st,kind(0.d0))) )**2 ) ) + + weight = S0_ratio * T0_proposal_ratio * abs( real(Phase_st * Prev_Ratiotot, kind=Kind(0.d0))/real(Phase_st,kind=Kind(0.d0)) ) + + if (weight > ranf_wrap()) then + acc = .true. + else + acc = .false. + nsigma%f(n,ntau1) = nsigma_st + Gr = Gr_st + phase = phase_st + endif + + Call Control_upgrade(acc) + Call Control_upgrade_eff(acc) + + endif Enddo - If ( N_Global_tau > 0 ) then - m = Nt_sequential_end + if (Nt_sequential_end >= Nt_sequential_start) then + m = Nt_sequential_end + else + m = 0 + endif + If ( N_Global_tau > 0 ) then !if ( Nt_sequential_start > Nt_sequential_end ) m = Nt_sequential_start Call Wrapgr_Random_update(GR,m,ntau1, PHASE, N_Global_tau ) - Call Wrapgr_PlaceGR(GR,m, Size(OP_V,1), ntau1) Endif + If ( N_Global_tau_MALA > 0 ) then + !if ( Nt_sequential_start > Nt_sequential_end ) m = Nt_sequential_start + Call Wrapgr_Langevin_update(GR,m,ntau1, PHASE, N_Global_tau_MALA, Delta_t_MALA_global_tau, Max_Force_MALA_global_tau ) + Endif + if ( N_Global_tau > 0 .or. N_Global_tau_MALA > 0 ) Call Wrapgr_PlaceGR(GR,m, Size(OP_V,1), ntau1) END SUBROUTINE WRAPGRUP !-------------------------------------------------------------------- - SUBROUTINE WRAPGRDO(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_end, N_Global_tau) + SUBROUTINE WRAPGRDO(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_end, N_Global_tau, & + & Propose_MALA, Delta_t_MALA_sequential, Max_Force_MALA_sequential, & + & N_Global_tau_MALA, Delta_t_MALA_global_tau, Max_Force_MALA_global_tau) !-------------------------------------------------------------------- !> @author !> ALF-project ! !> @brief -!> Given the green function matrix GR at time NTAU the routine carries out an +!> Given the green function matrix GR at time NTAU at m = size(OP_V,1), the routine carries out an !> update of the fields at time NTAU and propagates it to time NTAU-1 !> NTAU: [LTROT:1] ! @@ -176,26 +252,47 @@ SUBROUTINE WRAPGRDO(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT), allocatable :: GR(:,:,:) COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT) :: PHASE Integer, INTENT(IN) :: NTAU - LOGICAL, INTENT(IN) :: Propose_S0 - INTEGER, INTENT(IN) :: Nt_sequential_start, Nt_sequential_end, N_Global_tau + LOGICAL, INTENT(IN) :: Propose_S0, Propose_MALA + INTEGER, INTENT(IN) :: Nt_sequential_start, Nt_sequential_end, N_Global_tau, N_Global_tau_MALA + real (kind=kind(0.d0)), intent(in) :: Delta_t_MALA_sequential, Delta_t_MALA_global_tau + real (kind=kind(0.d0)), intent(in) :: Max_Force_MALA_sequential, Max_Force_MALA_global_tau ! Local - Integer :: nf, nf_eff, N_Type, n, m + Integer :: nf, nf_eff, N_Type, n, m, m1 Complex (Kind=Kind(0.d0)) :: Prev_Ratiotot, HS_Field, HS_New + Complex (Kind=Kind(0.d0)) :: force_old, force_new, phase_st, nsigma_st Real (Kind=Kind(0.d0)) :: T0_proposal, T0_Proposal_ratio, S0_ratio + real (kind=kind(0.d0)) :: force_0_old, force_0_new, weight + real (kind=kind(0.d0)) :: delta_t_running_old, Xmax, Delta_t_running_new Character (Len=64) :: Mode Logical :: Acc, toggle1 - + m = Size(OP_V,1) If ( N_Global_tau > 0 ) then - m = Size(OP_V,1) !Write(6,*) 'Call Ran_up ', m,ntau Call Wrapgr_Random_update(GR,m,ntau, PHASE, N_Global_tau ) - Call Wrapgr_PlaceGR(GR,m, Nt_sequential_end, ntau) Endif + If ( N_Global_tau_MALA > 0 ) then + !Write(6,*) 'Call Ran_up ', m,ntau + Call Wrapgr_Langevin_update(GR,m,ntau, PHASE, N_Global_tau_MALA, Delta_t_MALA_global_tau, Max_Force_MALA_global_tau ) + Endif + + !HERE + If (Nt_sequential_end >= Nt_sequential_start) then + m1 = Nt_sequential_end + else + m1 = 0 + endif + if ( N_Global_tau > 0 .or. N_Global_tau_MALA > 0 ) Call Wrapgr_PlaceGR(GR,m, m1, ntau) + Do n = Nt_sequential_end, Nt_sequential_start, -1 + if (Propose_MALA .and. Op_V(n,1)%type == 3) then + nsigma_st = nsigma%f(n,ntau) + phase_st = phase + force_old = calculate_force(n,ntau,Gr) + endif N_type = 2 nf = 1 HS_Field = nsigma%f(n,ntau) @@ -207,7 +304,20 @@ SUBROUTINE WRAPGRDO(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ nf = 1 T0_proposal = 1.5D0 T0_Proposal_ratio = 1.D0 - HS_new = nsigma%flip(n,ntau) + if ( Propose_MALA .and. Op_V(n,nf)%type == 3) then + Gr_st = Gr + Call ham%Ham_Langevin_HMC_S0_single( force_0_old, n,ntau) + call Control_MALA_sequential(force_old, force_0_old) + Xmax = abs(dble(force_old)) + if (abs(force_0_old) > Xmax) Xmax = abs(force_0_old) + delta_t_running_old = Delta_t_MALA_sequential + if (Xmax > Max_Force_MALA_sequential) delta_t_running_old = Max_Force_MALA_sequential*Delta_t_MALA_sequential/Xmax + HS_New = nsigma%f(n,ntau) - ( force_0_old + & + & real( Phase*force_old,kind(0.d0)) / Real(Phase,kind(0.d0)) ) * delta_t_running_old + & + & sqrt( 2.d0 * delta_t_running_old) * rang_wrap() + else + HS_new = nsigma%flip(n,ntau) + endif S0_ratio = ham%S0(n,ntau,HS_new) if ( Propose_S0 ) then If ( Op_V(n,nf)%type == 1) then @@ -216,9 +326,53 @@ SUBROUTINE WRAPGRDO(GR,NTAU,PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_ endif Endif If ( T0_proposal > ranf_wrap() ) Then - mode = "Final" Prev_Ratiotot = cmplx(1.d0,0.d0,kind(0.d0)) - Call Upgrade2(GR,n,ntau,PHASE,HS_new, Prev_Ratiotot, S0_ratio,T0_Proposal_ratio, Acc, mode ) + if (Propose_MALA .and. Op_V(n,1)%type == 3) then + mode = "Intermediate" + Call Upgrade2(GR,n,ntau,PHASE,HS_new, Prev_Ratiotot, S0_ratio,T0_Proposal_ratio, Acc, mode ) + phase = Phase * Prev_Ratiotot/sqrt(Prev_Ratiotot*conjg(Prev_Ratiotot)) + + N_type = 2 + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) + Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),HS_Field,Ndim,N_Type,ntau) + enddo + + Call ham%Ham_Langevin_HMC_S0_single( force_0_new,n,ntau) + force_new = calculate_force(n,ntau,GR) + call Control_MALA_sequential(force_new, force_0_new) + Xmax = abs(dble(force_new)) + if (abs(force_0_new) > Xmax) Xmax = abs(force_0_new) + Delta_t_running_new = Delta_t_MALA_sequential + if (Xmax > Max_Force_MALA_sequential) Delta_t_running_new = Max_Force_MALA_sequential*Delta_t_MALA_sequential/Xmax + + t0_Proposal_ratio = sqrt(delta_t_running_old/Delta_t_running_new) * exp(-0.25d0/Delta_t_running_new * (Abs(nsigma_st - hs_new + & + & Delta_t_running_new*(force_0_new + real( Phase*force_new,kind(0.d0)) / Real(Phase,kind(0.d0))) )**2) + 0.25d0/delta_t_running_old *( & + & Abs(hs_new - nsigma_st + delta_t_running_old*(force_0_old + real( phase_st*force_old,kind(0.d0)) / Real(phase_st,kind(0.d0))) )**2 ) ) + + weight = S0_ratio * T0_proposal_ratio * abs( real(Phase_st * Prev_Ratiotot, kind=Kind(0.d0))/real(Phase_st,kind=Kind(0.d0)) ) + + if (weight > ranf_wrap()) then + acc = .true. + N_type = 2 + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) + Call Op_Wrapdo( Gr(:,:,nf), Op_V(n,nf), HS_Field, Ndim, N_Type,ntau) + enddo + else + acc = .false. + nsigma%f(n,ntau) = nsigma_st + Gr = Gr_st + phase = phase_st + endif + + Call Control_upgrade(acc) + Call Control_upgrade_eff(acc) + + else + mode = "Final" + Call Upgrade2(GR,n,ntau,PHASE,HS_new, Prev_Ratiotot, S0_ratio,T0_Proposal_ratio, Acc, mode ) + endif else toggle1 = .false. Call Control_upgrade_eff(toggle1) @@ -433,8 +587,189 @@ Subroutine Wrapgr_Random_update( GR,m,ntau, PHASE, N_Global_tau ) end Subroutine Wrapgr_Random_update + +!-------------------------------------------------------------------- + Subroutine Wrapgr_Langevin_update( GR,m,ntau, PHASE, N_Global_tau_MALA, Delta_t_MALA_global_tau, Max_Force ) +!-------------------------------------------------------------------- +!> @author +!> ALF-project +! +!> @brief +!> On input: +!> GR(tau,m) as defined in Global_tau_mod_PlaceGR and the direction of updating scheme +!> direction=u --> You are visiting the time slices from tau = 1 to tau =Ltrot +!> direction=d --> You are visiting the time slices from tau = Ltrot to tau = 1 +!> +!> The routine calls +!> Global_MALA_move_tau(Flip_list, Flip_length,ntau) +!> in the Hamiltonian module and then carries out the update +!> +!> On output +!> +!> Flip_length==1 +!> Green function is on GR(tau,Flip_list(1) +1 ) if direction = u +!> Green function is on GR(tau,Flip_list(1) -1 ) if direction = d +!> This is valid if the move has or has not been accepted. +!> +!> Flip_length > 1 +!> Let m_min = min(Flip_list), m_max = max(Flip_list) +!> direction = u --> On output Green on m_max is accepted. Green is on m_min if not accepted. +!> direction = d --> On output Green on m_min if accepted. Green is on m_max if not accepted. +!-------------------------------------------------------------------- + + Implicit none + + ! Arguments + COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), INTENT(INOUT), allocatable :: GR + Integer, INTENT(INOUT) :: m + Integer, INTENT(IN) :: ntau, N_Global_tau_MALA + Complex (Kind=Kind(0.d0)), INTENT(INOUT) :: PHASE + Real (kind=Kind(0.d0)), intent(in) :: Delta_t_MALA_global_tau + real (kind=kind(0.d0)), intent(in) :: Max_Force + + + + ! Space for local variables + Integer :: n, Flip_length, nf, nf_eff, N_Type, ng_c, Flip_count + Real (Kind=Kind(0.d0)) :: T0_Proposal_ratio, T0_proposal,S0_ratio, force_0 + COMPLEX (Kind=Kind(0.d0)) :: Prev_Ratiotot, HS_Field, HS_New, Phase_st + Logical :: Acc + Character (Len=64) :: Mode + Integer, allocatable :: Flip_list(:) + Complex (Kind=Kind(0.d0)), allocatable :: Flip_value(:), Flip_value_st(:), forces_old(:), forces_new(:) + Real (Kind=Kind(0.d0)) :: Zero = 10D-8, weight + real (kind=kind(0.d0)), allocatable :: Forces_0_old(:), Forces_0_new(:) + Real (Kind=Kind(0.d0)) :: X, Xmax, delta_t_running_old, Delta_t_running_new + + Allocate ( Flip_list(Size(Op_V,1)), Flip_value(Size(Op_V,1)), Flip_value_st(Size(Op_V,1)) ) + Allocate ( Forces_old(Size(Op_V,1)), Forces_new(Size(Op_V,1)), Forces_0_old(Size(Op_V,1)), Forces_0_new(Size(Op_V,1))) + + Do ng_c = 1,N_Global_tau_MALA + Phase_st = Phase + ! New configuration + Call ham%Global_MALA_move_tau(Flip_list, Flip_length,ntau ) + !Write(6,*) "Calling global move", m, Flip_list(1), nsigma(Flip_list(1),ntau) + ! Order the list + Call wrapgr_sort(Flip_length,Flip_list) + Do Flip_count = 1, Flip_length + Flip_value_st(Flip_count) = nsigma%f( Flip_list(Flip_count), ntau ) + Enddo + !Write(6,*) "-----", Flip_length + !Calculate forces with current nsigma + Xmax = 0.d0 + do Flip_count = 1,Flip_length + n = Flip_list(Flip_count) + !Write(6,*) "PlaceGR", m, n-1,ntau + Call Wrapgr_PlaceGR(GR,m, n, ntau) + !Write(6,*) "Back from PlaceGR", m, n-1,ntau + If ( Flip_count == 1 ) GR_st = Gr + forces_old(Flip_count) = calculate_force(n,ntau,GR) + call ham%Ham_Langevin_HMC_S0_single (force_0,n,ntau) + forces_0_old(Flip_count) = force_0 + X = abs(dble(forces_old(Flip_count))) + if (X > Xmax) Xmax = X + X = abs(force_0) + if (X > Xmax) Xmax = X + m = n + enddo + delta_t_running_old = Delta_t_MALA_global_tau + If (Xmax > Max_Force) delta_t_running_old = Max_Force*Delta_t_MALA_global_tau/Xmax + + call Control_MALA_Global_tau(Forces_old, forces_0_old, flip_length) + + do Flip_count = 1,Flip_length + n = Flip_list(Flip_count) + flip_value(Flip_count) = nsigma%f(n,ntau) - ( forces_0_old(flip_count) + & + & real( Phase*forces_old(flip_count),kind(0.d0)) / Real(Phase,kind(0.d0)) ) * delta_t_running_old + & + & sqrt( 2.d0 * delta_t_running_old) * rang_wrap() + enddo + + !Update Greens function + S0_ratio = 1.d0 + Prev_Ratiotot = cmplx(1.d0,0.d0,kind(0.d0)) + do Flip_count = 1,Flip_length + n = Flip_list(Flip_count) + !Write(6,*) "PlaceGR", m, n-1,ntau + Call Wrapgr_PlaceGR(GR,m, n-1, ntau) + !Write(6,*) "Back from PlaceGR", m, n-1,ntau + Do nf_eff = 1, N_FL_eff + nf=Calc_Fl_map(nf_eff) + HS_Field = nsigma%f(n,ntau) + N_type = 1 + Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),HS_Field,Ndim,N_Type,ntau) + enddo + nf = 1 + mode = "Intermediate" + HS_new = Flip_value(Flip_count) + S0_ratio = S0_ratio * ham%S0(n,ntau,hs_new) + Call Upgrade2(GR,n,ntau,PHASE, HS_new , & + & Prev_Ratiotot, S0_ratio, T0_Proposal_ratio, Acc, mode ) + do nf_eff = 1,N_FL_eff + nf=Calc_Fl_map(nf_eff) + N_type = 2 + Call Op_Wrapup(Gr(:,:,nf),Op_V(n,nf),HS_Field,Ndim,N_Type,ntau) + enddo + m = n + enddo + phase = Phase * Prev_Ratiotot/sqrt(Prev_Ratiotot*conjg(Prev_Ratiotot)) + + !Calculate forces with new nsigma + Xmax = 0.d0 + do Flip_count = 1, Flip_length + n = Flip_list(Flip_count) + !Write(6,*) "PlaceGR", m, n-1,ntau + Call Wrapgr_PlaceGR(GR,m, n, ntau) + !Write(6,*) "Back from PlaceGR", m, n-1,ntau + forces_new(Flip_count) = calculate_force(n,ntau,GR) + call ham%Ham_Langevin_HMC_S0_single (force_0,n,ntau) + forces_0_new(Flip_count) = force_0 + X = abs(dble(forces_new(Flip_count))) + if (X > Xmax) Xmax = X + X = abs(force_0) + if (X > Xmax) Xmax = X + m = n + enddo + Delta_t_running_new = Delta_t_MALA_global_tau + If (Xmax > Max_Force) Delta_t_running_new = Max_Force*Delta_t_MALA_global_tau/Xmax + + call Control_MALA_Global_tau(Forces_new, forces_0_new, flip_length) + + t0_proposal_ratio = 1.d0 + do Flip_count = 1, Flip_length + n = Flip_list(Flip_count) + t0_Proposal_ratio = t0_proposal_ratio * sqrt(delta_t_running_old/Delta_t_running_new) * exp(-0.25d0/Delta_t_running_new * & + & (Abs(Flip_value_st(Flip_count) - Flip_value(Flip_count) + Delta_t_running_new*(forces_0_new(Flip_count) + & + & real( Phase*forces_new(Flip_count),kind(0.d0)) / Real(Phase,kind(0.d0))) )**2 ) + 0.25d0/delta_t_running_old * & + & (Abs(Flip_value(Flip_count) - Flip_value_st(Flip_count) + delta_t_running_old*(forces_0_old(Flip_count) + & + & real( phase_st*forces_old(Flip_count),kind(0.d0)) / Real(phase_st,kind(0.d0))) )**2 ) ) + enddo + + weight = S0_ratio * T0_proposal_ratio * abs( real(Phase_st * Prev_Ratiotot, kind=Kind(0.d0))/real(Phase_st,kind=Kind(0.d0)) ) + + if (weight > ranf_wrap()) then + acc = .true. + else + acc = .false. + Gr = Gr_st + phase = phase_st + m = Flip_list(1) + Do Flip_count = 1, Flip_length + nsigma%f( Flip_list(Flip_count), ntau ) = Flip_value_st(Flip_count) + Enddo + endif + + Call Control_upgrade(acc) + Call Control_upgrade_eff(acc) + + Enddo + + Deallocate ( Flip_list, Flip_value, Flip_value_st ) + + + end Subroutine Wrapgr_Langevin_update + !---------------------------------------------------------------------------- - subroutine Wrapgr_sort(Flip_length,Flip_list,Flip_value) + subroutine Wrapgr_sort_value(Flip_length,Flip_list,Flip_value) ! Arguments Integer, INTENT(IN) :: Flip_length @@ -476,7 +811,57 @@ subroutine Wrapgr_sort(Flip_length,Flip_list,Flip_value) ! Write(6,*) Flip_list(nc), Flip_value(nc) !Enddo - end subroutine Wrapgr_sort + end subroutine Wrapgr_sort_value + +!---------------------------------------------------------------------------- + subroutine Wrapgr_sort_langevin(Flip_length,Flip_list) + + ! Arguments + Integer, INTENT(INOUT) :: Flip_length + Integer, INTENT(INOUT), allocatable :: Flip_list(:) + + ! Local + integer :: swaps ! number of swaps made in one pass + integer :: nc ! loop variable + integer :: temp, n ! temporary holder for making swap + integer :: length + + !Check if all fields are continuous and remove non-continuous fields from flip_list + length = 0 + do nc = 1, Flip_length + if (Op_V(flip_list(nc),1)%type == 3) then + length = length + 1 + flip_list(length) = flip_list(nc) + endif + enddo + flip_length = length + + if ( Flip_length == 1 ) return + + !Write(6,*) 'Before sort' + !DO nc = 1,Flip_length + ! Write(6,*) Flip_list(nc) + !Enddo + + do + swaps = 0 ! Initially, we've made no swaps + do nc = 1, (Flip_length - 1) + if ( Flip_list(nc) > Flip_list(nc+1) ) then + temp = Flip_list(nc ) + Flip_list(nc) = Flip_list(nc+1) + Flip_list(nc+1) = temp + swaps = swaps + 1 + end if + end do + if ( swaps == 0 ) exit ! do count swaps + end do + + !Write(6,*) 'After sort' + !DO nc = 1,Flip_length + ! Write(6,*) Flip_list(nc) + !Enddo + + end subroutine Wrapgr_sort_langevin !---------------------------------------------------------------------------- diff --git a/Prog/control_mod.F90 b/Prog/control_mod.F90 index d99e2b535..17142cbea 100644 --- a/Prog/control_mod.F90 +++ b/Prog/control_mod.F90 @@ -58,17 +58,29 @@ module Control Integer (Kind=kind(0.d0)), private, save :: NC_Glob_up, ACC_Glob_up Integer (Kind=kind(0.d0)), private, save :: NC_HMC_up, ACC_HMC_up Integer (Kind=kind(0.d0)), private, save :: NC_Temp_up, ACC_Temp_up + Integer (Kind=kind(0.d0)), private, save :: NC_MALA_up, ACC_MALA_up real (Kind=Kind(0.d0)), private, save :: XMAXP_Glob, XMEANP_Glob Integer (Kind=Kind(0.d0)), private, save :: NC_Phase_GLob real (Kind=Kind(0.d0)), private, save :: XMAXP_HMC, XMEANP_HMC Integer (Kind=Kind(0.d0)), private, save :: NC_Phase_HMC + real (Kind=Kind(0.d0)), private, save :: XMAXP_MALA, XMEANP_MALA + Integer (Kind=Kind(0.d0)), private, save :: NC_Phase_MALA real (Kind=Kind(0.d0)), private, save :: size_clust_Glob_up, size_clust_Glob_ACC_up real (Kind=Kind(0.d0)), private, save :: Force_max, Force_mean Integer, private, save :: Force_Count + real (Kind=Kind(0.d0)), private, save :: Force_max_MALA_seq , Force_mean_MALA_seq + real (Kind=Kind(0.d0)), private, save :: Force_0_max_MALA_seq, Force_0_mean_MALA_seq + Integer, private, save :: Force_Count_MALA_seq + real (Kind=Kind(0.d0)), private, save :: Force_max_MALA_gtau, Force_mean_MALA_gtau + real (Kind=Kind(0.d0)), private, save :: Force_0_max_MALA_gtau, Force_0_mean_MALA_gtau + Integer, private, save :: Force_Count_MALA_gtau + real (Kind=Kind(0.d0)), private, save :: Force_max_MALA_global, Force_mean_MALA_global + real (Kind=Kind(0.d0)), private, save :: Force_0_max_MALA_global, Force_0_mean_MALA_global + Integer, private, save :: Force_Count_MALA_global #ifdef MPI Integer , private, save :: Ierr, Isize, Irank, irank_g, isize_g, igroup #endif @@ -92,6 +104,8 @@ subroutine control_init(Group_Comm) XMAXP_Glob = 0.d0 XMEANP_HMC = 0.d0 XMAXP_HMC = 0.d0 + XMEANP_MALA= 0.d0 + XMAXP_MALA = 0.d0 NCG = 0 @@ -107,6 +121,10 @@ subroutine control_init(Group_Comm) NC_Phase_HMC = 0 NC_HMC_up = 0 ACC_HMC_up = 0 + + NC_Phase_MALA= 0 + NC_MALA_up = 0 + ACC_MALA_up = 0 NC_Temp_up = 0 ACC_Temp_up = 0 @@ -117,6 +135,24 @@ subroutine control_init(Group_Comm) Force_max = 0.d0 Force_mean = 0.d0 Force_count = 0 + + Force_max_MALA_seq = 0.d0 + Force_mean_MALA_seq = 0.d0 + Force_0_max_MALA_seq = 0.d0 + Force_0_mean_MALA_seq = 0.d0 + Force_Count_MALA_seq = 0 + + Force_max_MALA_gtau = 0.d0 + Force_mean_MALA_gtau = 0.d0 + Force_0_max_MALA_gtau = 0.d0 + Force_0_mean_MALA_gtau = 0.d0 + Force_Count_MALA_gtau = 0 + + Force_max_MALA_global = 0.d0 + Force_mean_MALA_global = 0.d0 + Force_0_max_MALA_global = 0.d0 + Force_0_mean_MALA_global = 0.d0 + Force_Count_MALA_global = 0 #ifdef MPI CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) @@ -160,6 +196,94 @@ Subroutine Control_Langevin(Forces, Group_Comm) end Subroutine Control_Langevin +!------------------------------------------------------------- + + Subroutine Control_MALA_sequential(Force, force_0) + + Implicit none + + Complex (Kind=Kind(0.d0)), Intent(In) :: Force + Real (Kind=Kind(0.d0)), Intent(In) :: Force_0 + + Force_Count_MALA_seq = Force_Count_MALA_seq + 1 + + If ( abs( Real(Force,kind(0.d0))) >= Force_max_MALA_seq ) & + & Force_max_MALA_seq = abs( Real(Force,kind(0.d0))) + If ( abs( Force_0) >= Force_0_max_MALA_seq ) & + & Force_0_max_MALA_seq = abs( Force_0) + Force_mean_MALA_seq = Force_mean_MALA_seq + abs( Real(Force,kind(0.d0))) + Force_0_mean_MALA_seq = Force_0_mean_MALA_seq + abs( Force_0) + + end Subroutine Control_MALA_sequential + +!------------------------------------------------------------- + + Subroutine Control_MALA_Global(Forces, forces_0, flip_list) + + + Implicit none + + Complex (Kind=Kind(0.d0)), Intent(In) :: Forces(:,:) + Real (Kind=Kind(0.d0)), Intent(In) :: Forces_0(:,:) + integer, intent(in) :: flip_list(:,:) + + Integer :: n1,n2, n, nt + Real (Kind = Kind(0.d0) ) :: X, Y + + n1 = size(Forces,1) + n2 = size(Forces,2) + + X = 0.d0 + Y = 0.d0 + do n = 1,n1 + do nt =1,n2 + if (flip_list(n,nt) == 1) then + Force_Count_MALA_global = Force_Count_MALA_global + 1 + If ( abs( Real(Forces(n,nt),kind(0.d0))) >= Force_max_MALA_global ) & + & Force_max_MALA_global = abs( Real(Forces(n,nt),kind(0.d0))) + X = X + abs( Real(Forces(n,nt),kind(0.d0)) ) + If ( abs( Real(Forces_0(n,nt),kind(0.d0))) >= Force_0_max_MALA_global ) & + & Force_0_max_MALA_global = abs( Real(Forces_0(n,nt),kind(0.d0))) + Y = Y + abs( Real(Forces_0(n,nt),kind(0.d0)) ) + endif + enddo + enddo + Force_mean_MALA_global = Force_mean_MALA_global + X + Force_0_mean_MALA_global = Force_0_mean_MALA_global + Y + + end Subroutine Control_MALA_Global + +!------------------------------------------------------------- + + Subroutine Control_MALA_Global_tau(Forces, forces_0, flip_length) + + + Implicit none + + Complex (Kind=Kind(0.d0)), Intent(In) :: Forces(:) + Real (Kind=Kind(0.d0)), Intent(In) :: Forces_0(:) + integer, intent(in) :: flip_length + + Integer :: n + Real (Kind = Kind(0.d0) ) :: X, Y + + X = 0.d0 + Y = 0.d0 + do n = 1,flip_length + Force_Count_MALA_gtau = Force_Count_MALA_gtau + 1 + If ( abs( Real(Forces(n),kind(0.d0))) >= Force_max_MALA_gtau ) & + & Force_max_MALA_gtau = abs( Real(Forces(n),kind(0.d0))) + X = X + abs( Real(Forces(n),kind(0.d0)) ) + If ( abs( Real(Forces_0(n),kind(0.d0))) >= Force_0_max_MALA_gtau ) & + & Force_0_max_MALA_gtau = abs( Real(Forces_0(n),kind(0.d0))) + Y = Y + abs( Real(Forces_0(n),kind(0.d0)) ) + enddo + Force_mean_MALA_gtau = Force_mean_MALA_gtau + X + Force_0_mean_MALA_gtau = Force_0_mean_MALA_gtau + Y + + end Subroutine Control_MALA_Global_tau + +!------------------------------------------------------------- Subroutine Control_upgrade(toggle) Implicit none @@ -203,6 +327,15 @@ Subroutine Control_upgrade_HMC(toggle) endif end Subroutine Control_upgrade_HMC + Subroutine Control_upgrade_MALA(toggle) + Implicit none + Logical :: toggle + NC_MALA_up = NC_MALA_up + 1 + if (toggle) then + ACC_MALA_up = ACC_MALA_up + 1 + endif + end Subroutine Control_upgrade_MALA + Subroutine Control_PrecisionG(A,B,Ndim) #ifdef MPI @@ -340,8 +473,18 @@ Subroutine Control_PrecisionP_HMC(Z,Z1) NC_Phase_HMC = NC_Phase_HMC + 1 End Subroutine Control_PrecisionP_HMC + Subroutine Control_PrecisionP_MALA(Z,Z1) + Implicit none + Complex (Kind=Kind(0.D0)), INTENT(IN) :: Z,Z1 + Real (Kind=Kind(0.D0)) :: X + X = ABS(Z-Z1) + if ( X > XMAXP_MALA ) XMAXP_MALA = X + XMEANP_MALA = XMEANP_MALA + X + NC_Phase_MALA = NC_Phase_MALA + 1 + End Subroutine Control_PrecisionP_MALA + - Subroutine Control_Print(Group_Comm, Global_update_scheme) + Subroutine Control_Print(Group_Comm, Global_update_scheme, MALA) #ifdef MPI Use mpi #endif @@ -349,10 +492,12 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) Integer, Intent(IN) :: Group_Comm Character (Len = 64), Intent(IN) :: Global_update_scheme + Logical, intent(in) :: MALA Character (len=64) :: file1 Real (Kind=Kind(0.d0)) :: Time, Acc, Acc_eff, Acc_Glob, Acc_Temp, size_clust_Glob, size_clust_Glob_ACC, Acc_HMC + Real (Kind=Kind(0.d0)) :: Acc_MALA #ifdef MPI REAL (Kind=Kind(0.d0)) :: X Integer :: Ierr, Isize, Irank, irank_g, isize_g, igroup @@ -380,6 +525,10 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) IF (NC_HMC_up > 0 ) then ACC_HMC = dble(ACC_HMC_up)/dble(NC_HMC_up) endif + ACC_MALA = 0.d0 + IF (NC_MALA_up > 0 ) then + ACC_MALA = dble(ACC_MALA_up)/dble(NC_MALA_up) + endif ACC_TEMP = 0.d0 @@ -389,7 +538,19 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) call system_clock(count_CPU_end) time = (count_CPU_end-count_CPU_start)/dble(count_rate) if (count_CPU_end .lt. count_CPU_start) time = (count_max+count_CPU_end-count_CPU_start)/dble(count_rate) - If (str_to_upper(Global_update_scheme) == "LANGEVIN") Force_mean = Force_mean/real(Force_count,kind(0.d0)) + If (str_to_upper(Global_update_scheme) == "LANGEVIN") Force_mean = Force_mean/real(Force_count,kind(0.d0)) + If (Force_Count_MALA_seq > 0) then + Force_mean_MALA_seq = Force_mean_MALA_seq/real(Force_Count_MALA_seq,kind(0.d0)) + Force_0_mean_MALA_seq = Force_0_mean_MALA_seq/real(Force_Count_MALA_seq,kind(0.d0)) + endif + If (Force_Count_MALA_gtau > 0) then + Force_mean_MALA_gtau = Force_mean_MALA_gtau/real(Force_Count_MALA_gtau,kind(0.d0)) + Force_0_mean_MALA_gtau = Force_0_mean_MALA_gtau/real(Force_Count_MALA_gtau,kind(0.d0)) + endif + If (MALA) then + Force_mean_MALA_global = Force_mean_MALA_global/real(Force_Count_MALA_global,kind(0.d0)) + Force_0_mean_MALA_global = Force_0_mean_MALA_global/real(Force_Count_MALA_global,kind(0.d0)) + endif #if defined(MPI) If (str_to_upper(Global_update_scheme) == "LANGEVIN") then @@ -399,6 +560,42 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) CALL MPI_REDUCE(Force_max,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) Force_max= X endif + If (Force_Count_MALA_seq > 0) then + X = 0.d0 + CALL MPI_REDUCE(Force_mean_MALA_seq,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + Force_mean_MALA_seq= X/dble(Isize_g) + CALL MPI_REDUCE(Force_max_MALA_seq,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + Force_max_MALA_seq= X + X = 0.d0 + CALL MPI_REDUCE(Force_0_mean_MALA_seq,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + Force_0_mean_MALA_seq= X/dble(Isize_g) + CALL MPI_REDUCE(Force_0_max_MALA_seq,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + Force_0_max_MALA_seq= X + endif + If (Force_Count_MALA_gtau > 0) then + X = 0.d0 + CALL MPI_REDUCE(Force_mean_MALA_gtau,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + Force_mean_MALA_gtau= X/dble(Isize_g) + CALL MPI_REDUCE(Force_max_MALA_gtau,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + Force_max_MALA_gtau= X + X = 0.d0 + CALL MPI_REDUCE(Force_0_mean_MALA_gtau,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + Force_0_mean_MALA_gtau= X/dble(Isize_g) + CALL MPI_REDUCE(Force_0_max_MALA_gtau,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + Force_0_max_MALA_gtau= X + endif + If (MALA) then + X = 0.d0 + CALL MPI_REDUCE(Force_mean_MALA_global,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + Force_mean_MALA_global= X/dble(Isize_g) + CALL MPI_REDUCE(Force_max_MALA_global,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + Force_max_MALA_global= X + X = 0.d0 + CALL MPI_REDUCE(Force_0_mean_MALA_global,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + Force_0_mean_MALA_global= X/dble(Isize_g) + CALL MPI_REDUCE(Force_0_max_MALA_global,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + Force_0_max_MALA_global= X + endif X = 0.d0 CALL MPI_REDUCE(ACC,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) ACC = X/dble(Isize_g) @@ -414,6 +611,9 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) X = 0.d0 CALL MPI_REDUCE(ACC_Temp ,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) ACC_Temp = X/dble(Isize_g) + X = 0.d0 + CALL MPI_REDUCE(ACC_MALA,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) + ACC_MALA = X/dble(Isize_g) X = 0.d0 CALL MPI_REDUCE(size_clust_Glob,X,1,MPI_REAL8,MPI_SUM, 0,Group_Comm,IERR) @@ -452,6 +652,9 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) CALL MPI_REDUCE(XMAXP_HMC,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) XMAXP_HMC = X + CALL MPI_REDUCE(XMAXP_MALA,X,1,MPI_REAL8,MPI_MAX, 0,Group_Comm,IERR) + XMAXP_MALA = X + #endif #if defined(TEMPERING) @@ -498,6 +701,24 @@ Subroutine Control_Print(Group_Comm, Global_update_scheme) Write(50,*) ' Mean Phase diff HMC : ', XMEANP_HMC Write(50,*) ' Max Phase diff HMC : ', XMAXP_HMC Endif + + if (Force_Count_MALA_seq > 0) Then + Write(50,*) ' Sequential MALA Force Mean, Max : ', Force_mean_MALA_seq, Force_max_MALA_seq + Write(50,*) ' Sequential MALA Force_0 Mean, Max : ', Force_0_mean_MALA_seq, Force_0_max_MALA_seq + Endif + + if (Force_Count_MALA_gtau > 0) Then + Write(50,*) ' Global tau MALA Force Mean, Max : ', Force_mean_MALA_gtau, Force_max_MALA_gtau + Write(50,*) ' Global tau MALA Force_0 Mean, Max : ', Force_0_mean_MALA_gtau, Force_0_max_MALA_gtau + Endif + + if (MALA) Then + Write(50,*) ' MALA Force Mean, Max : ', Force_mean_MALA_global, Force_max_MALA_global + Write(50,*) ' MALA Force_0 Mean, Max : ', Force_0_mean_MALA_global, Force_0_max_MALA_global + Write(50,*) ' Acceptance_MALA : ', ACC_MALA + Write(50,*) ' Mean Phase diff MALA : ', XMEANP_MALA + Write(50,*) ' Max Phase diff MALA : ', XMAXP_MALA + Endif Write(50,*) ' CPU Time : ', Time Close(50) diff --git a/Prog/main.F90 b/Prog/main.F90 index ca769d534..575decd6b 100644 --- a/Prog/main.F90 +++ b/Prog/main.F90 @@ -162,11 +162,11 @@ Program Main Complex (Kind=Kind(0.d0)) , allocatable, dimension(:,:) :: Initial_field ! Space for choosing sampling scheme - Logical :: Propose_S0, Tempering_calc_det - Logical :: Global_moves, Global_tau_moves + Logical :: Propose_S0, Tempering_calc_det, Propose_MALA + Logical :: Global_moves, Global_tau_moves, Global_tau_MALA_moves Integer :: N_Global Integer :: Nt_sequential_start, Nt_sequential_end, mpi_per_parameter_set - Integer :: N_Global_tau + Integer :: N_Global_tau, N_Global_tau_MALA Logical :: Sequential real (Kind=Kind(0.d0)) :: Amplitude ! Needed for update of type 3 and 4 fields. @@ -175,9 +175,11 @@ Program Main Logical :: file_exists #endif ! Space for reading in Langevin & HMC parameters - Logical :: Langevin, HMC - Integer :: Leapfrog_Steps, N_HMC_sweeps + Logical :: Langevin, HMC, MALA + Integer :: Leapfrog_Steps, N_HMC_sweeps, N_MALA_sweeps Real (Kind=Kind(0.d0)) :: Delta_t_Langevin_HMC, Max_Force + real (kind=kind(0.d0)) :: MAX_Force_MALA_global, Max_Force_MALA_global_tau, Max_Force_MALA_sequential + real (kind=kind(0.d0)) :: Delta_t_MALA_sequential, Delta_t_MALA_global_tau, Delta_t_MALA_global #if defined(TEMPERING) Integer :: N_exchange_steps, N_Tempering_frequency @@ -188,7 +190,12 @@ Program Main & Propose_S0,Global_moves, N_Global, Global_tau_moves, & & Nt_sequential_start, Nt_sequential_end, N_Global_tau, & & sequential, Langevin, HMC, Delta_t_Langevin_HMC, & - & Max_Force, Leapfrog_steps, N_HMC_sweeps, Amplitude + & Max_Force, Leapfrog_steps, N_HMC_sweeps, Amplitude, & + & Propose_MALA, Delta_t_MALA_sequential, & + & Max_Force_MALA_sequential, Global_tau_MALA_moves, & + & N_Global_tau_MALA, Delta_t_MALA_global_tau, & + & Max_Force_MALA_global_tau, MALA, N_MALA_sweeps, & + & MAX_Force_MALA_global, delta_t_MALA_global NAMELIST /VAR_HAM_NAME/ ham_name @@ -363,6 +370,10 @@ Program Main Global_tau_moves = .false.; sequential = .true.; Langevin = .false. ; HMC =.false. Delta_t_Langevin_HMC = 0.d0; Max_Force = 0.d0 ; Leapfrog_steps = 0; N_HMC_sweeps = 1 Nt_sequential_start = 1 ; Nt_sequential_end = 0; N_Global_tau = 0; Amplitude = 1.d0 + Propose_MALA = .false.; Delta_t_MALA_sequential = 0.d0; Max_Force_MALA_sequential = 0.d0 + Global_tau_MALA_moves = .false.; N_Global_tau_MALA = 0; Delta_t_MALA_global_tau = 0.d0 + Max_Force_MALA_global_tau = 0.d0; MALA = .false.; N_MALA_sweeps = 1 + delta_t_MALA_global = 0.d0; MAX_Force_MALA_global = 0.d0 OPEN(UNIT=5,FILE=file_para,STATUS='old',ACTION='read',IOSTAT=ierr) IF (ierr /= 0) THEN WRITE(error_unit,*) 'main: unable to open ', file_para, ierr @@ -397,6 +408,17 @@ Program Main CALL MPI_BCAST(Max_Force ,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) CALL MPI_BCAST(Delta_t_Langevin_HMC ,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) CALL MPI_BCAST(Amplitude ,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(Propose_MALA ,1 ,MPI_LOGICAL ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(Delta_t_MALA_sequential,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(Max_Force_MALA_sequential,1 ,MPI_REAL8,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(Global_tau_MALA_moves,1 ,MPI_LOGICAL ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(N_Global_tau_MALA ,1 ,MPI_Integer ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(delta_t_MALA_global_tau,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(MAX_Force_MALA_global_tau,1 ,MPI_REAL8,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(MALA ,1 ,MPI_LOGICAL ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(N_MALA_sweeps ,1 ,MPI_Integer ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(MAX_Force_MALA_global,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) + CALL MPI_BCAST(delta_t_MALA_global ,1 ,MPI_REAL8 ,0,MPI_COMM_i,ierr) CALL MPI_BCAST(ham_name ,64,MPI_CHARACTER,0,MPI_COMM_i,ierr) #endif @@ -489,14 +511,16 @@ Program Main LOBS_EN = Ltrot endif endif - If ( .not. Global_tau_moves ) then + if ( .not. Global_tau_moves ) N_Global_tau = 0 + if ( .not. Global_tau_MALA_moves ) N_Global_tau_MALA = 0 + if ( Global_tau_moves .or. Global_tau_MALA_moves ) then + ! Gives the possibility to set parameters in the Hamiltonian file + Call ham%Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end, & + & N_Global_tau,N_Global_tau_MALA) + else ! This corresponds to the default updating scheme Nt_sequential_start = 1 Nt_sequential_end = Size(OP_V,1) - N_Global_tau = 0 - else - ! Gives the possibility to set parameters in the Hamiltonian file - Call ham%Overide_global_tau_sampling_parameters(Nt_sequential_start,Nt_sequential_end,N_Global_tau) endif call nsigma%make(N_op, Ltrot) @@ -519,7 +543,7 @@ Program Main Call Hop_mod_init IF (ABS(CPU_MAX) > Zero ) NBIN = 10000000 - If (N_Global_tau > 0) then + If (N_Global_tau > 0 .or. N_Global_tau_MALA > 0 .or. Propose_MALA) then Call Wrapgr_alloc endif @@ -579,6 +603,10 @@ Program Main write(output_unit,*) "Langevin mode does not allow HMC updates." write(output_unit,*) "Overriding HMC=.True. from parameter files." endif + if (MALA) then + write(output_unit,*) "Langevin mode does not allow MALA updates." + write(output_unit,*) "Overriding MALA=.True. from parameter files." + endif if (Global_moves) then write(output_unit,*) "Langevin mode does not allow global updates." write(output_unit,*) "Overriding Global_moves=.True. from parameter files." @@ -587,6 +615,10 @@ Program Main write(output_unit,*) "Langevin mode does not allow global tau updates." write(output_unit,*) "Overriding Global_tau_moves=.True. from parameter files." endif + if (Global_tau_MALA_moves) then + write(output_unit,*) "Langevin mode does not allow global tau MALA updates." + write(output_unit,*) "Overriding Global_tau_MALA_moves=.True. from parameter files." + endif #if defined(TEMPERING) if ( N_exchange_steps > 0 ) then write(output_unit,*) "Langevin mode does not allow tempering updates." @@ -598,15 +630,22 @@ Program Main #endif Sequential = .False. HMC = .False. + MALA = .False. Global_moves = .False. Global_tau_moves = .False. + Global_tau_MALA_moves = .False. #if defined(TEMPERING) N_exchange_steps = 0 #endif endif - Call Langevin_HMC%make(Langevin, HMC , Delta_t_Langevin_HMC, Max_Force, Leapfrog_steps) + Call Langevin_HMC%make(Langevin, HMC , .False., Delta_t_Langevin_HMC, Max_Force, Leapfrog_steps) + else + Call Langevin_HMC%set_Update_scheme(Langevin, HMC, .False. ) + endif + if (MALA) then + Call Metropolis_Langevin%make(.False., .False., MALA, delta_t_MALA_global, MAX_Force_MALA_global, Leapfrog_steps) else - Call Langevin_HMC%set_Update_scheme(Langevin, HMC ) + Call Metropolis_Langevin%set_Update_scheme(.False., .False., MALA ) endif if ( .not. Sequential .and. Global_tau_moves) then @@ -615,8 +654,14 @@ Program Main write(output_unit,*) "Sequential is set to .False. ." endif - if ( .not. Sequential .and. .not. HMC .and. .not. Langevin .and. .not. Global_moves) then - write(output_unit,*) "Warning: no updates will occur as Sequential, HMC, Langevin, and" + if ( .not. Sequential .and. Global_tau_MALA_moves) then + write(output_unit,*) "Warning: Sequential = .False. and Global_tau_MALA_moves = .True." + write(output_unit,*) "in the parameter file. Global Langevin tau updates will not occur if" + write(output_unit,*) "Sequential is set to .False. ." + endif + + if ( .not. Sequential .and. .not. HMC .and. .not. Langevin .and. .not. Global_moves .and. .not. MALA) then + write(output_unit,*) "Warning: no updates will occur as Sequential, HMC, Langevin, MALA, and" write(output_unit,*) "Global_moves are all .False. in the parameter file." endif @@ -624,6 +669,19 @@ Program Main write(output_unit,*) "Warning: Nt_sequential_end is smaller than Nt_sequential_start" endif + if ( Propose_MALA .or. Global_tau_MALA_moves) then + Do n = 1,N_op + if ( nsigma%t(n) /= 3 ) then + write(output_unit,*) + WRITE(output_unit,*) 'Warning: Not all fields are of type 3.' + WRITE(output_unit,*) 'Fields that are not of type 3 will not be updated with MALA updates.' + write(output_unit,*) + exit + exit + endif + enddo + endif + #if defined(MPI) if ( Irank_g == 0 ) then #endif @@ -642,12 +700,23 @@ Program Main Write(50,*) '# of interacting Ops per time slice : ', Size(OP_V,1) If ( Propose_S0 ) & & Write(50,*) 'Propose Ising moves according to bare Ising action' + If ( Propose_MALA ) then + Write(50,*) 'Propose continuous moves according to Langevin equation' + Write(50,*) 'Langevin del_t: ', Delta_t_MALA_sequential + Write(50,*) 'Max Force MALA sequential : ', Max_Force_MALA_sequential + Endif If ( Global_moves ) Then Write(50,*) 'Global moves are enabled ' Write(50,*) '# of global moves / sweep :', N_Global Endif if ( sequential ) then - If ( Global_tau_moves ) Then + If ( Global_tau_MALA_moves ) Then + Write(50,*) 'Nt_sequential_start : ', Nt_sequential_start + Write(50,*) 'Nt_sequential_end : ', Nt_sequential_end + Write(50,*) 'N_Global_tau_MALA : ', N_Global_tau_MALA + Write(50,*) 'Langevin del_t : ', Delta_t_MALA_global_tau + Write(50,*) 'Max Force MALA global tau : ', MAX_Force_MALA_global_tau + elseif ( Global_tau_moves ) Then Write(50,*) 'Nt_sequential_start: ', Nt_sequential_start Write(50,*) 'Nt_sequential_end : ', Nt_sequential_end Write(50,*) 'N_Global_tau : ', N_Global_tau @@ -664,13 +733,18 @@ Program Main Write(50,*) 'Leapfrog_Steps: ', Leapfrog_Steps Write(50,*) 'HMC_Sweeps: ', N_HMC_sweeps endif + if ( MALA ) then + Write(50,*) 'MALA global del_t : ', delta_t_MALA_global + Write(50,*) 'MALA_Sweeps : ', N_MALA_sweeps + Write(50,*) 'Max Force MALA global : ', MAX_Force_MALA_global + endif !Write out info for amplitude and flip_protocol Toggle = .false. Do n = 1,N_op if (nsigma%t(n) == 3 .or. nsigma%t(n) == 4) Toggle = .true. Enddo - if ( Toggle ) then + if ( Toggle .and. (.not.Propose_MALA) ) then Write(50,*) 'Amplitude for t=3,4 vertices is set to: ', Amplitude endif Toggle = .false. @@ -848,6 +922,39 @@ Program Main endif endif + If ( str_to_upper(Metropolis_Langevin%get_Update_scheme()) == "MALA" ) then + if (Sequential .and. str_to_upper(Langevin_HMC%get_Update_scheme()) /= "HMC" ) then + call Metropolis_Langevin%set_L_Forces(.False.) + endif + Do n=1, N_MALA_sweeps + Call Metropolis_Langevin%update(Phase, GR, GR_Tilde, Test, udvr, udvl, Stab_nt, udvst, & + & LOBS_ST, LOBS_EN, LTAU) + if (n /= N_MALA_sweeps) then + Call Metropolis_Langevin%calc_Forces(Phase, GR, GR_Tilde, Test, udvr, udvl, Stab_nt, udvst,& + & LOBS_ST, LOBS_EN, .True. ) + Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst) + call Metropolis_Langevin%set_L_Forces(.true.) + endif + enddo + + !Do time-displaced measurements if needed, else set Calc_Obser_eq=.True. for the very first leapfrog ONLY + If ( .not. sequential) then + IF ( LTAU == 1 ) then + If (Projector) then + NST = 0 + Call Tau_p ( udvl, udvr, udvst, GR, PHASE, NSTM, STAB_NT, NST, LOBS_ST, LOBS_EN) + else + Call Tau_m( udvst, GR, PHASE, NSTM, NWRAP, STAB_NT, LOBS_ST, LOBS_EN ) + endif + else + Call Metropolis_Langevin%calc_Forces(Phase, GR, GR_Tilde, Test, udvr, udvl, Stab_nt, udvst,& + & LOBS_ST, LOBS_EN, .True. ) + Call Langevin_HMC_Reset_storage(Phase, GR, udvr, udvl, Stab_nt, udvst) + endif + call Metropolis_Langevin%set_L_Forces(.true.) + endif + endif + If (Sequential) then ! Propagation from 1 to Ltrot ! Set the right storage to 1 @@ -863,7 +970,9 @@ Program Main NST = 1 DO NTAU = 0, LTROT-1 NTAU1 = NTAU + 1 - CALL WRAPGRUP(GR,NTAU,PHASE,Propose_S0, Nt_sequential_start, Nt_sequential_end, N_Global_tau) + CALL WRAPGRUP(GR,NTAU,PHASE,Propose_S0, Nt_sequential_start, Nt_sequential_end, N_Global_tau, & + & Propose_MALA, Delta_t_MALA_sequential, Max_Force_MALA_sequential, & + & N_Global_tau_MALA, delta_t_MALA_global_tau, Max_Force_MALA_global_tau) If (NTAU1 == Stab_nt(NST) ) then NT1 = Stab_nt(NST-1) @@ -922,7 +1031,9 @@ Program Main NST = NSTM-1 DO NTAU = LTROT,1,-1 NTAU1 = NTAU - 1 - CALL WRAPGRDO(GR,NTAU, PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_end, N_Global_tau) + CALL WRAPGRDO(GR,NTAU, PHASE,Propose_S0,Nt_sequential_start, Nt_sequential_end, N_Global_tau, & + & Propose_MALA, Delta_t_MALA_sequential, Max_Force_MALA_sequential, & + & N_Global_tau_MALA, delta_t_MALA_global_tau, Max_Force_MALA_global_tau) IF (NTAU1.GE. LOBS_ST .AND. NTAU1.LE. LOBS_EN ) THEN !write(*,*) "GR before obser sum: ",sum(GR(:,:,1)) !write(*,*) "Phase before obser : ",phase @@ -1049,7 +1160,7 @@ Program Main DEALLOCATE(udvl, udvr, udvst) DEALLOCATE(GR, TEST, Stab_nt,GR_Tilde) if (Projector) DEALLOCATE(WF_R, WF_L) - If (N_Global_tau > 0) then + If (N_Global_tau > 0 .or. N_Global_tau_MALA > 0 .or. Propose_MALA) then Call Wrapgr_dealloc endif do nf = 1, N_FL @@ -1067,7 +1178,7 @@ Program Main call deallocate_all_shared_memory #endif - Call Control_Print(Group_Comm, Langevin_HMC%get_Update_scheme()) + Call Control_Print(Group_Comm, Langevin_HMC%get_Update_scheme(), MALA) #if defined(MPI) If (Irank_g == 0 ) then @@ -1083,6 +1194,7 @@ Program Main #endif Call Langevin_HMC%clean() + Call Metropolis_Langevin%clean() deallocate(Calc_Fl_map,Phase_array) ! Delete the file RUNNING since the simulation finished successfully