Draft: Resolve "Proposal of updates with Langevin equation"

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}

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.

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}

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 

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

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