Global_mod.F90

!  Copyright (C) 2016 - 2022 The ALF project
!
!  This file is part of the ALF project.
!
!     The ALF project is free software: you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation, either version 3 of the License, or
!     (at your option) any later version.
!
!     The ALF project is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with ALF.  If not, see http://www.gnu.org/licenses/.
!
!     Under Section 7 of GPL version 3 we require you to fulfill the following additional terms:
!
!     - It is our hope that this program makes a contribution to the scientific community. Being
!       part of that community we feel that it is reasonable to require you to give an attribution
!       back to the original authors if you have benefitted from this program.
!       Guidelines for a proper citation can be found on the project's homepage
!       http://alf.physik.uni-wuerzburg.de .
!
!     - We require the preservation of the above copyright notice and this license in all original files.
!
!     - We prohibit the misrepresentation of the origin of the original source files. To obtain
!       the original source files please visit the homepage http://alf.physik.uni-wuerzburg.de .
!
!     - If you make substantial changes to the program we require you to either consider contributing
!       to the ALF project or to mark your material in a reasonable way as different from the original version.
 
 
 
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> Handles global updates and parallel tempering
!
!--------------------------------------------------------------------
 
Module global_mod
 
      Use runtime_error_mod
      Use hamiltonian_main
      Use mymats
      Use operator_mod
      Use control
      Use observables
      Use fields_mod
      Use random_wrap
      use iso_fortran_env, only: output_unit, error_unit
 
      Implicit none
 
      Type (Obser_Vec ),  private  ::   Tempering_acceptance
 
 
    Contains
#if defined(TEMPERING)
!--------------------------------------------------------------------
!> @author
!> ALF-project
!
!> @brief
!> Handles parrallel tempering \n
!> This subroutine is called only if the tempering flag is switched on \n
!> Requires MPI
!>
!> @details
!> On entry and on exit the left storage is full, the Green function is on time slice 0 and the phase is set.
!> The same holds on exit but with updated quantities.
!>
!> @param[inout] Phase  Complex
!> \verbatim
!>  Is updated upon acceptance
!> \endverbatim
!> @param[inout] udvl, udvr Class(UDV_State)
!> \verbatim
!>  On entry udvl contains the last udv decomposition such that G=(1 + VL*DL*UL^{dag})^{-1}
!>  udvr  is used as working space
!> \endverbatim
!> @param[inout] GR Complex
!> \verbatim
!>  Green function. Is updated upon acceptance.
!> \endverbatim
!> @param[inout] udvst Class(UDV_State)
!> \verbatim
!>  Storage. Is updated with left propagation upon acceptance
!> \endverbatim
!> @param[in] Stab_nt Integer
!> \verbatim
!>  List of time slices for stabilization.
!> \endverbatim
!> @param[in] N_exchange_steps Integer
!> \verbatim
!> Number of exchange steps
!> \endverbatim
!> @param[in] Tempering_calc_det Logical
!> \verbatim
!>  If true then the fermion determinant is computed in the calulcation of the exchange weight.
!> \endverbatim
!>
!--------------------------------------------------------------------
      Subroutine exchange_step(Phase,GR, udvr, udvl, Stab_nt, udvst, N_exchange_steps, Tempering_calc_det)
        Use udv_state_mod
        Use mpi
        use wrapul_mod
        use cgr1_mod
        Implicit none
 
        !  Arguments
        COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT)                                :: Phase
        CLASS(udv_state), intent(inout), allocatable, Dimension(:)              :: udvl, udvr
        COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:), INTENT(INOUT), allocatable :: GR
        CLASS(udv_state), intent(inout), allocatable, Dimension(:, :) :: udvst
        INTEGER, dimension(:),     INTENT   (IN), allocatable         :: Stab_nt
        INTEGER, INTENT(IN) :: N_exchange_steps
        Logical, INTENT(IN) :: Tempering_calc_det
 
 
        !>  Local variables.
        Integer :: NST, NSTM, NF, nf_eff, NT, NT1, NVAR,N, N1,N2, I, NC, I_Partner, n_step,  N_count, N_part
        Type    (Fields)           :: nsigma_old
        Real    (Kind=kind(0.d0)) :: t0_proposal_ratio, weight, weight1, delta_s0_log, exp_delta_s0
        Complex (Kind=Kind(0.d0)) :: Z_ONE = cmplx(1.d0, 0.d0, kind(0.d0)), z, ratiotot, ratiotot_p, phase_old, phase_new
        Real    (Kind=kind(0.d0)), allocatable :: det_vec_old(:,:), det_vec_new(:,:)
        Complex (Kind=Kind(0.d0)), allocatable :: Phase_Det_new(:), Phase_Det_old(:)
        Complex (Kind=Kind(0.d0)) :: Ratio(2), Ratio_p(2),Phase_array(N_FL)
        Logical :: TOGGLE, L_Test
        Integer, allocatable :: List_partner(:), List_masters(:)
 
        ! Keep track of where the configuration originally came from
        Integer        :: nsigma_irank, nsigma_old_irank, nsigma_irank_temp
        Integer        :: n_GR
 
        !Integer, Dimension(:,:),  allocatable :: nsigma_orig, nsigma_test
        !Integer :: I1, I2
 
        Integer        :: Isize, Irank, Ierr, irank_g, isize_g, igroup
        Integer        :: STATUS(MPI_STATUS_SIZE)
 
        Character (Len=64)  :: storage
 
 
        CALL mpi_comm_size(mpi_comm_world,isize,ierr)
        CALL mpi_comm_rank(mpi_comm_world,irank,ierr)
        call mpi_comm_rank(group_comm, irank_g, ierr)
        call mpi_comm_size(group_comm, isize_g, ierr)
        igroup           = irank/isize_g
        nsigma_irank     = irank
 
        n1 = size(nsigma%f,1)
        n2 = size(nsigma%f,2)
        nstm = Size(udvst, 1)
        call nsigma_old%make(n1, n2)
        if (tempering_calc_det) then
           Allocate ( det_vec_old(ndim,n_fl), det_vec_new(ndim,n_fl) )
           Allocate ( phase_det_new(n_fl), phase_det_old(n_fl) )
        endif
        Allocate ( list_partner(0:isize-1), list_masters(isize/2) )
 
        !         Allocate ( nsigma_orig(n1,n2) )
        !         nsigma_orig = nsigma
 
        !  Compute for each core the old weights.
        l_test = .false.
        if (tempering_calc_det) then
           ! Set old weight.
           storage = "Full"
           Call compute_fermion_det(phase_det_old,det_vec_old, udvl, udvst, stab_nt, storage)
           phase_array=1.0d0
           Do nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              phase_array(nf) = phase_det_old(nf)
              Call op_phase(phase_array(nf),op_v,nsigma,nf)
           Enddo
           if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
           phase_old=product(phase_array)
           phase_old=phase_old**n_sun
        endif
        !> Store old configuration
        nsigma_old%f = nsigma%f
        nsigma_old%t = nsigma%t
        nsigma_old_irank = nsigma_irank
        ! Setup the list of masters
        nc = 0
        Do i  = 0,isize-1,2*isize_g
           do n = 0,isize_g-1
              nc = nc + 1
              list_masters(nc)  = i + n
           enddo
        Enddo
        if (l_test) then
           if (irank == 0)  then
              Write(6,*) 'List  of masters: '
              Do nc = 1,isize/2
                 Write(6,*) nc, list_masters(nc)
              Enddo
           endif
        endif
        DO n_count = 1,n_exchange_steps
           !  Set the partner rank on each core
           If (irank == 0 ) then
              n_step = isize_g
              if (  ranf_wrap() > 0.5d0 ) n_step = -isize_g
              Do i = 0,isize-1,2*isize_g
                 do n = 0,isize_g-1
                    list_partner(npbc_tempering(i  + n             ,isize)) =  npbc_tempering(i + n   + n_step ,isize)
                    list_partner(npbc_tempering(i  + n  + n_step  , isize)) =  npbc_tempering(i + n            ,isize)
                 enddo
              enddo
           endif
 
           CALL mpi_bcast(list_partner, isize  ,mpi_integer,   0,mpi_comm_world,ierr)
 
           If (l_test) then
              if (irank == 0 ) then
                 Write(6,*) 'Testing global : '
                 do  i = 0,isize -1
                    Write(6,*)   i, list_partner(i) !, Phase_old, Phase
                    Write(11,*)  i, list_partner(i) !, Phase_old, Phase
                 enddo
                 Write(11,*)
              endif
           Endif
 
 
           !  Exchange configurations
           !  The types do not change --> no need to exchange them
           n = size(nsigma_old%f,1)*size(nsigma_old%f,2)
           CALL mpi_sendrecv(nsigma_old%f    , n, mpi_complex16, list_partner(irank), 0, &
                    &        nsigma%f        , n, mpi_complex16, list_partner(irank), 0, mpi_comm_world,status,ierr)
           CALL mpi_sendrecv(nsigma_old_irank, 1, mpi_integer, list_partner(irank), 0, &
                    &        nsigma_irank    , 1, mpi_integer, list_partner(irank), 0, mpi_comm_world,status,ierr)
 
           !  Each node now has a new configuration nsigma
 
 
           if (tempering_calc_det) then
              !  HERE
              !  Compute ratio on weights one each rank
              storage = "Empty"
              Call compute_fermion_det(phase_det_new,det_vec_new, udvl, udvst, stab_nt, storage)
 
              phase_array = cmplx(1.d0,0.d0,kind=kind(0.d0))
              Do nf_eff = 1,n_fl_eff
                 nf=calc_fl_map(nf_eff)
                 phase_array(nf) = phase_det_new(nf)
                 Call op_phase(phase_array(nf),op_v,nsigma,nf)
              Enddo
              if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
              phase_new=product(phase_array)
              phase_new=phase_new**n_sun
 
              t0_proposal_ratio = 1.d0
              ratiotot = compute_ratio_global(phase_det_old, phase_det_new, &
                   &            det_vec_old, det_vec_new, nsigma_old, t0_proposal_ratio,ratio)
 
              If (l_test) Write(6,*) 'Ratio_global: Irank, Partner',irank,list_partner(irank), &
                   &                  ratiotot, ratio(1)*exp(ratio(2))
           else
              delta_s0_log = ham%Get_Delta_S0_global(nsigma_old)
              ratio(1) = 1.0d0 !!! @FAKHER @JONAS double check please
              ratio(2) = delta_s0_log
           endif
 
           !  Acceptace/rejection decision taken on master node after receiving information from slave
           Do nc = 1,isize/2 ! Loop over masters
              i = list_masters(nc)
              If (irank == i ) Then
                 CALL mpi_send(ratio   ,2, mpi_complex16  , list_partner(i), i+1024, mpi_comm_world,ierr)
              else if (irank == list_partner(i) ) Then
                 CALL mpi_recv(ratio_p    , 2, mpi_complex16,  i, i+1024, mpi_comm_world,status,ierr)
                 !Weight = abs(Ratiotot_p*Ratiotot)
                 weight= abs(ratio(1) * ratio_p(1) * exp( ratio_p(2) + ratio(2)  ) )
                 toggle = .false.
                 if ( weight > ranf_wrap() )  toggle =.true.
                 If (l_test) Write(6,*) 'Master : ', list_partner(i), i, weight, toggle
              endif
           enddo
 
           !>  Send result of acceptance/rejection decision form master to slave
           Do nc =  1,isize/2 ! Loop over masters
              i = list_masters(nc)
              If (irank == list_partner(i) ) Then
                 ! Write(6,*) 'Send from ', List_Partner(I), 'to, ', I, I + 512
                 CALL mpi_send(toggle, 1, mpi_logical, i, i+1024, mpi_comm_world,ierr)
              else if (irank == i ) Then
                 CALL mpi_recv(toggle , 1, mpi_logical,   list_partner(i), i+1024 ,mpi_comm_world,status,ierr)
                 If (l_test) Write(6,*) 'Slave : ', irank,  toggle
              endif
           enddo
 
           If (l_test) Write(6,*) 'Final: ',  irank, list_partner(irank), toggle
           Call global_tempering_obser(toggle)
           Call control_upgrade_temp  (toggle)
           If (toggle)  then
              !     Move has been accepted
              if (tempering_calc_det) then
                 phase_old     = phase_new
                 phase_det_old = phase_det_new
                 det_vec_old   = det_vec_new
              endif
              nsigma_old%f       = nsigma%f
              nsigma_old%t       = nsigma%t
              nsigma_old_irank = nsigma_irank
           else
              nsigma%f       = nsigma_old%f
              nsigma%t       = nsigma_old%t
              nsigma_irank = nsigma_old_irank
           endif
        enddo
 
        ! Finalize
        if (tempering_calc_det) then
           ! If move has been accepted, no use to recomute storage
           If (.not.toggle) then
              DO nf_eff = 1,n_fl_eff
                 nf=calc_fl_map(nf_eff)
                 if (projector) then
                    CALL udvl(nf_eff)%reset('l',wf_l(nf)%P)
                 else
                    CALL udvl(nf_eff)%reset('l')
                 endif
              ENDDO
              DO nst = nstm-1,1,-1
                 nt1 = stab_nt(nst+1)
                 nt  = stab_nt(nst  )
                 !Write(6,*) NT1,NT, NST
                 CALL wrapul(nt1,nt, udvl)
                 Do nf_eff = 1,n_fl_eff
                    udvst(nst, nf_eff) = udvl(nf_eff)
                 ENDDO
              ENDDO
              nt1 = stab_nt(1)
              CALL wrapul(nt1,0, udvl)
           Endif
           ! Compute the Green functions so as to provide correct starting point for the sequential updates.
           nvar  = 1
           phase_array = cmplx(1.d0,0.d0,kind(0.d0))
           do nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              CALL cgr(z, nvar, gr(:,:,nf), udvr(nf_eff),  udvl(nf_eff))
              call op_phase(z,op_v,nsigma,nf)
              phase_array(nf)=z
           Enddo
           if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
           phase=product(phase_array)
           phase=phase**n_sun
        else
           !  Send >> Phase, GR, udvr, udvl, udvst << to new node
           !  First step: Each node sends to IRANK=0 its value nsigma_irank,
           !  which is the node where its new Phase, GR, udvr, udvl, udvst is stored
           !  This node then tells each node where to send its now old Phase, GR, udvr, udvl, udvst
           !  Finally, the variables get submitted
           If (irank == 0) then
              Do i = 1,isize-1
                 CALL mpi_recv(nsigma_irank_temp , 1, mpi_integer, i, 0, mpi_comm_world,status,ierr)
                 If ( nsigma_irank_temp == 0) then
                    nsigma_old_irank = i
                 else
                    CALL mpi_send(i , 1, mpi_integer, nsigma_irank_temp, 0, mpi_comm_world,ierr)
                 endif
              enddo
              If ( nsigma_irank /= 0 ) then
                 CALL mpi_send(0 , 1, mpi_integer, nsigma_irank, 0, mpi_comm_world,ierr)
              endif
           else
              CALL mpi_send(nsigma_irank     , 1, mpi_integer, 0, 0, mpi_comm_world,ierr)
              CALL mpi_recv(nsigma_old_irank , 1, mpi_integer, 0, 0, mpi_comm_world,status,ierr)
           endif
 
           if ( nsigma_irank /= irank ) then
              CALL mpi_sendrecv_replace(phase, 1, mpi_complex16, nsigma_old_irank, 0, &
                   &        nsigma_irank, 0, mpi_comm_world, status, ierr)
 
              n_gr = size(gr,1)*size(gr,2)*size(gr,3)
              CALL mpi_sendrecv_replace(gr, n_gr, mpi_complex16, nsigma_old_irank, 0, &
                   &        nsigma_irank, 0, mpi_comm_world, status, ierr)
 
              do nf_eff = 1,n_fl_eff
                 CALL udvr(nf_eff)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, status, ierr)
              enddo
              do nf_eff = 1,n_fl_eff
                 CALL udvl(nf_eff)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, status, ierr)
              enddo
              do nst = 1, nstm
                 do nf_eff = 1,n_fl_eff
                    CALL udvst(nst, nf_eff)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, status, ierr)
                 enddo
              enddo
           endif
        endif
 
 
        call nsigma_old%clear
 
        if (tempering_calc_det) then
           Deallocate ( det_vec_old, det_vec_new )
           Deallocate ( phase_det_new, phase_det_old )
        endif
        Deallocate ( list_partner, list_masters )
 
      end Subroutine exchange_step
#endif
 
 
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> Carries out N_global global updates as defeined in the Global_move routine of the Hamiltonian module
!>
!> @details
!> On entry and on exit the left storage is full, the Green function is on time slice 0 and the phase is set.
!> The same holds on exit but with updated quantities.
!>
!> @param[inout] Phase  Complex
!> \verbatim
!>  Is updated upon acceptance
!> \endverbatim
!> @param[inout] udvl, udvr Class(UDV_State)
!> \verbatim
!>  On entry udvl contains the last udv decomposition such that G=(1 + VL*DL*UL^{dag})^{-1}
!>  udvr  is used as working space
!> \endverbatim
!> @param[inout] GR Complex
!> \verbatim
!>  Green function. Is updated upon acceptance.
!> \endverbatim
!> @param[inout] udvst Class(UDV_State)
!> \verbatim
!>  Storage. Is updated with left propagation upon acceptance
!> \endverbatim
!> @param[in] Stab_nt Integer
!> \verbatim
!>  List of time slices for stabilization.
!> \endverbatim
!> @param[in] N_Global Integer
!> \verbatim
!>  Number of global moves that will be caried out
!> \endverbatim
!>
!--------------------------------------------------------------------
      Subroutine global_updates(Phase,GR, udvr, udvl, Stab_nt, udvst,N_Global)
 
        Use udv_state_mod
        use wrapul_mod
        use cgr1_mod
        Implicit none
 
        !  Arguments
        COMPLEX (Kind=Kind(0.d0)),                                INTENT(INOUT) :: Phase
        class(udv_state), DIMENSION(:), ALLOCATABLE,           INTENT(INOUT) :: udvl, udvr
        COMPLEX (Kind=Kind(0.d0)), Dimension(:,:,:),  allocatable,INTENT(INOUT) :: GR
        CLASS(udv_state), Dimension(:,:), ALLOCATABLE,            INTENT(INOUT) :: udvst
        INTEGER, dimension(:),   allocatable,                     INTENT(IN)    :: Stab_nt
        Integer,                                                  INTENT(IN)    :: N_Global
 
 
        !  Local variables.
        Integer :: NST, NSTM, NF, NT, NT1, NVAR,N, N1,N2, I, NC, N_part,j, nf_eff
        Real    (Kind=kind(0.d0)) :: t0_proposal_ratio, weight
        Complex (Kind=Kind(0.d0)) :: Z_ONE = cmplx(1.d0, 0.d0, kind(0.d0)), z, ratiotot, phase_old, phase_new
        Complex (Kind=Kind(0.d0)), allocatable :: Det_vec_test(:,:), Phase_Det_new(:), Phase_Det_old(:)
        Real    (Kind=kind(0.d0)), allocatable :: det_vec_old(:,:), det_vec_new(:,:)
        Type   (Fields)   :: nsigma_old
 
        Complex (Kind=Kind(0.d0)) :: Ratio(2), Phase_array(N_Fl)
        Logical :: TOGGLE, L_Test
        Real    (Kind=kind(0.d0)) :: size_clust
        Real    (Kind=kind(0.d0)) :: ratio_2_test
        Character (Len=64)  :: storage
 
 
 
        n1 = size(nsigma%f,1)
        n2 = size(nsigma%f,2)
        nstm = Size(udvst, 1)
        call nsigma_old%make(n1, n2)
 
        Allocate ( det_vec_old(ndim,n_fl_eff), det_vec_new(ndim,n_fl_eff), det_vec_test(ndim,n_fl_eff) )
        Allocate ( phase_det_new(n_fl_eff), phase_det_old(n_fl_eff) )
 
        l_test = .false.
        ! Set old weight.
        storage = "Full"
        Call compute_fermion_det(phase_det_old,det_vec_old, udvl, udvst, stab_nt, storage)
        phase_array = cmplx(1.d0,0.d0,kind=kind(0.d0))
        Do nf_eff = 1,n_fl_eff
           nf=calc_fl_map(nf_eff)
           phase_array(nf) = phase_det_old(nf_eff)
           Call op_phase(phase_array(nf),op_v,nsigma,nf)
        Enddo
        if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
        phase_old=product(phase_array)
        phase_old=phase_old**n_sun
 
        If (l_test) then
           ! Testing
           Do nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              if (projector) then
                 CALL udvr(nf_eff)%reset('r',wf_r(nf)%P)
              else
                 CALL udvr(nf_eff)%reset('r')
              endif
           Enddo
           nvar = 1
           phase_array = cmplx(1.d0,0.d0,kind(0.d0))
           do nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              CALL cgr(z, nvar, gr(:,:,nf), udvr(nf_eff), udvl(nf_eff))
              call op_phase(z,op_v,nsigma,nf)
              phase_array(nf)=z
           Enddo
           if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
           phase=product(phase_array)
           phase=phase**n_sun
           Do nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              Call det_c_lu(gr(:,:,nf),det_vec_test(:,nf_eff),ndim)
              !!!ATTENTION HOW DOES BELOW DEPEND ON NF BLOCK symm
              z = phase_det_old(nf_eff)
              ratio_2_test=0.d0
              DO i = 1,ndim
                 z = z*det_vec_test(i,nf_eff)/abs(det_vec_test(i,nf_eff))
                 ratio_2_test=ratio_2_test+log(abs(det_vec_test(i,nf_eff)))+det_vec_old(i,nf_eff)
              Enddo
              z=z*cmplx(exp(ratio_2_test),0.d0,kind(0.d0))
              Write(6,*) 'Testing weight: ', z
              if(abs(ratio_2_test)>650) write(6,*) "Weight is about to reach double underflow!"
           Enddo
        Endif
 
        ! Store old configuration
        nsigma_old%f = nsigma%f
        nsigma_old%t = nsigma%t
        ! Phase_old, Phase_det_old and Det_vec_old  are all set.
        nc = 0
        Do n = 1,n_global
           ! Draw a new spin configuration. This is provided by the user in the Hamiltonian module
           ! Note that nsigma is a variable in the module Hamiltonian
           Call ham%Global_move(t0_proposal_ratio,nsigma_old,size_clust)
           If (t0_proposal_ratio > 0.d0) then
              nc = nc + 1
              ! Compute the new Green function
              storage = "Empty"
              Call compute_fermion_det(phase_det_new,det_vec_new, udvl, udvst, stab_nt, storage)
 
              phase_array = cmplx(1.d0,0.d0,kind=kind(0.d0))
              Do nf_eff = 1,n_fl_eff
                 nf=calc_fl_map(nf_eff)
                 phase_array(nf) = phase_det_new(nf)
                 Call op_phase(phase_array(nf),op_v,nsigma,nf)
              Enddo
              if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
              phase_new=product(phase_array)
              phase_new=phase_new**n_sun
 
              ratiotot = compute_ratio_global(phase_det_old, phase_det_new, &
                   &                          det_vec_old, det_vec_new, nsigma_old, t0_proposal_ratio, ratio)
 
              !Write(6,*) 'Ratio_global: ', Ratiotot
 
              weight = abs(  real( phase_old * ratiotot, kind=kind(0.d0))/real(phase_old,kind=kind(0.d0)) )
 
              z = phase_old * ratiotot/abs(ratiotot)
              Call control_precisionp_glob(z,phase_new)
              !Write(6,*) Z, Phase_new
 
 
              toggle = .false.
              if ( weight > ranf_wrap() )  Then
                 toggle = .true.
                 phase_old     = phase_new
                 phase_det_old = phase_det_new
                 nsigma_old%t    = nsigma%t
                 nsigma_old%f    = nsigma%f
                 det_vec_old     = det_vec_new
              else
                 nsigma%t = nsigma_old%t
                 nsigma%f = nsigma_old%f
              endif
              Call control_upgrade_glob(toggle,size_clust)
           else
              nsigma%t = nsigma_old%t
              nsigma%f = nsigma_old%f
           endif
        Enddo
 
        If (nc > 0 ) then
           If (.not.toggle) then
              DO nf_eff = 1,n_fl_eff
                nf=calc_fl_map(nf_eff)
                if (projector) then
                  CALL udvl(nf_eff)%reset('l',wf_l(nf)%P)
                else
                  CALL udvl(nf_eff)%reset('l')
                endif
              ENDDO
              DO nst = nstm-1,1,-1
                 nt1 = stab_nt(nst+1)
                 nt  = stab_nt(nst  )
                 !Write(6,*) NT1,NT, NST
                 CALL wrapul(nt1,nt, udvl)
                 Do nf_eff = 1,n_fl_eff
                    udvst(nst, nf_eff) = udvl(nf_eff)
                 ENDDO
              ENDDO
              nt1 = stab_nt(1)
              CALL wrapul(nt1,0, udvl)
           Endif
           !Compute the Green functions so as to provide correct starting point for the sequential updates.
           nvar  = 1
           phase_array = cmplx(1.d0,0.d0,kind(0.d0))
           do nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              CALL cgr(z, nvar, gr(:,:,nf), udvr(nf_eff), udvl(nf_eff))
              call op_phase(z,op_v,nsigma,nf)
              phase_array(nf)=z
           Enddo
           if (reconstruction_needed) call ham%weight_reconstruction(phase_array)
           phase=product(phase_array)
           phase=phase**n_sun
        endif
 
 
        call nsigma_old%clear
        Deallocate ( det_vec_old  , det_vec_new, det_vec_test  )
        Deallocate ( phase_det_new, phase_det_old )
 
 
      End Subroutine global_updates
 
 
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> This fucntion computes ratio of weights for global moves
!>
!> @details
!> \verbatim
!>  Ratio_Global = T0(nsigma--> nsigma_old) W(nsigma)/T0(nsigma_old--> nsigma) W(nsigma_old) =
!>                 T0_Proposal_ratio   W(nsigma)/  W(nsigma_old)
!>  On entry the old and new fermion determinants read: Phase_det*e^{ \sum_{n=1}^{ndim} Det_vec(n) }
!>  as obtained from the Compute_Fermion_det routine.
!>  On exit Ratio_Global = Ratio(1)*exp(Ratio(2))
!> \endverbatim
!>
!> @param[in]  Phase_det_new  Complex, Dimension(N_FL)
!> @param[in]  Phase_det_old  Complex, Dimension(N_FL)
!> @param[in]  Det_vec_new  Real, Dimension(:,N_FL)
!> @param[in]  Det_vec_old  Real, Dimension(:,N_FL)
!> @param[in]  T0_proposal_ratio   Real
!> @param[in]  nsigma_old Type(Fields)
!> \verbatim
!>  Old configuration. The new configuration is stored in nsigma. nsigma is a globale variable
!>  contained in the Hamiltonian module.
!> \endverbatim
!--------------------------------------------------------------------
      Complex (Kind=Kind(0.d0)) Function compute_ratio_global(Phase_Det_old, Phase_Det_new, &
           &                    Det_vec_old, Det_vec_new, nsigma_old, T0_Proposal_ratio,Ratio)
 
 
        Implicit none
 
        ! Arguments
        Complex (Kind=Kind(0.d0)), allocatable, INTENT(IN) :: phase_det_old(:), phase_det_new(:)
        REAL    (kind=kind(0.d0)), allocatable, INTENT(IN) :: det_vec_old(:,:), det_vec_new(:,:)
        Real    (kind=kind(0.d0)),    INTENT(IN)  :: t0_proposal_ratio
        Type    (fields),             INTENT(IN)  :: nsigma_old
        Complex (Kind=Kind(0.d0)),    INTENT(out) :: ratio(2)
 
        ! Local
        Integer  :: nf, i, nt, nf_eff
        Complex (Kind=Kind(0.d0)) :: z, z1, ratio_1_array(n_fl), ratio_2_array(n_fl), g_loc
        Real    (kind=kind(0.d0)) :: x, ratio_2, delta, log_delta
 
        ratio = cmplx(0.d0,0.d0,kind(0.d0))
        ratio_2_array = 0.d0
        !X = 1.d0
        Do nf_eff = 1,n_fl_eff
           nf=calc_fl_map(nf_eff)
           DO i = 1,ndim
              !X= X*real(Det_vec_new(I,nf),kind(0.d0)) / Real(Det_vec_old(I,nf),kind(0.d0) )
              ratio_2_array(nf) = ratio_2_array(nf) +  det_vec_new(i,nf_eff) - det_vec_old(i,nf_eff)
           enddo
           !Z = Z**N_SUN
           ratio_2_array(nf) = real(n_sun,kind(0.d0))*ratio_2_array(nf)
        enddo
        !Z = cmplx(X,0.d0,kind(0.d0))
        ratio_1_array = cmplx(1.d0,0.d0,kind(0.d0))
        Do nf_eff = 1,n_fl_eff
           nf=calc_fl_map(nf_eff)
           !Z = Z*Phase_Det_new(nf)/Phase_Det_old(nf)
           ratio_1_array(nf_eff) = ratio_1_array(nf) *  phase_det_new(nf_eff)/phase_det_old(nf_eff)
           !Z = Z**N_SUN
           ratio_1_array(nf) = ratio_1_array(nf)**n_sun
        enddo
 
        ratio(1)=1.d0
        Do i = 1,Size(op_v,1)
           x = 0.d0
           Do nt = 1,ltrot
              !Z = Z * cmplx( Gaml(nsigma(i,nt),2)/Gaml(nsigma_old(i,nt),2),0.d0,kind(0.d0) )
              ratio(1) = ratio(1) * cmplx( nsigma%Gama(i,nt)/nsigma_old%Gama(i,nt),0.d0,kind(0.d0) ) !You could put this in Ratio_2
              !X = X + nsigma%Phi(i,nt) - nsigma_old%Phi(i,nt)
              Do nf_eff = 1,n_fl_eff
                 nf=calc_fl_map(nf_eff)
                 z =  nsigma%Phi(i,nt) - nsigma_old%Phi(i,nt)
                 g_loc = op_v(i,nf)%g
                 if ( op_v(i,nf)%get_g_t_alloc()) g_loc = op_v(i,nf)%g_t(nt)
                 !Z = Z * exp(cmplx( X*Real(N_SUN,Kind(0.d0)), 0.d0,kind(0.d0)) * Op_V(i,nf)%g * Op_V(i,nf)%alpha )
                 ratio_1_array(nf) = ratio_1_array(nf) * exp(z*cmplx( real(n_sun,kind(0.d0)), 0.d0,kind(0.d0)) * g_loc * op_v(i,nf)%alpha )
              Enddo
           Enddo
        Enddo
        if (reconstruction_needed) then
            call ham%weight_reconstruction(ratio_1_array)
            call ham%weight_reconstruction(ratio_2_array)
        endif
        ratio(1)=ratio(1)*product(ratio_1_array)
        !Z =  Z * cmplx( ham%Delta_S0_global(Nsigma_old),0.d0,kind(0.d0) )
        !Z =  Z * cmplx( T0_Proposal_ratio, 0.d0,kind(0.d0))
        ratio(2) = sum(ratio_2_array)
        log_delta = ham%Get_Delta_S0_global(nsigma_old)
        ratio(2) = ratio(2) + log_delta + log(t0_proposal_ratio)
 
        compute_ratio_global = ratio(1)*exp(ratio(2))
 
 
      end Function compute_ratio_global
 
!--------------------------------------------------------------------
!> @author
!>
!> @brief
!> Computes the fermion determinant from scratch or from the storage.
!>
!> @details
!> \verbatim
!> Finite temperature: If the storage is full then computes 
!>                     det( 1 +  VL*DL*UL^{dag} )= Phase_det*e^{ \sum_{n=1}^{ndim} Det_vec(n)}
!>                     Note that  Udvl%U = UL, Udvl%D = DL and  Udvl%V = VL^{dag}
!>                     If storage is empty one first computes  Udvl and in doing so fills up the storage.
!> Zero temperature:   If storage is full then computes  det \Psi_L U(\theta_tot,0) \Psi_R  =  Phase_det*e^{ \sum_{n=1}^{N_part} Det_vec(n) }
!>                     For the projective code the required scales which one can throw away for the Green function are in udvst.
!>                     If the storage is not full, then things will be computed from scratch and intermediate results will be stored in udvst.
!> \endverbatim
!>
!> @param[inout] udvl, udvr Class(UDV_State)
!> \verbatim
!>  If storage=Full  : On entry and exit udvl contains the last udv decomposition such that G=(1 + VL*DL*UL^{dag})^{-1}
!>  If storage=Empty : On exit udv1 will contain the last udv decomposition such that G=(1 + VL*DL*UL^{dag})^{-1}
!> \endverbatim
!> @param[out] Phase_det  Complex, Dimension(N_FL)
!> \verbatim
!>  The phase, per flavor
!> \endverbatim
!> @param[out] Det_vec  Real, Dimension(:,N_FL)
!> @param[in] Stab_nt  Integer, Dimension(:)
!> \verbatim
!>  List of time slices for stabilization.
!> \endverbatim
!> @param[in] Storage Character
!> \verbatim
!>  storage = "Full".  Compute the fermion det  with the use of the storage.
!>  storage = "Empty". Compute the fermion det from scratch and  in doing so, fill the storage.
!> \endverbatim
!>
!--------------------------------------------------------------------
      Subroutine compute_fermion_det(Phase_det,Det_Vec, udvl, udvst, Stab_nt, storage)
!--------------------------------------------------------------------
 
        Use udv_wrap_mod
        Use udv_state_mod
        use wrapul_mod
 
        Implicit none
 
        REAL    (Kind=kind(0.d0)), Dimension(:,:), Intent(OUT)  ::  det_vec
        Complex (Kind=Kind(0.d0)), Dimension(:)  , Intent(OUT)  ::  Phase_det
        CLASS(udv_state), DIMENSION(:)  , ALLOCATABLE,  INTENT(INOUT) :: udvl
        CLASS(udv_state), Dimension(:,:), ALLOCATABLE,  INTENT(INOUT) :: udvst
        INTEGER,          dimension(:),   allocatable,    INTENT(IN)  :: Stab_nt
        Character (Len=64), Intent(IN) :: storage
 
 
        !> Local variables
        Integer ::  N_size, NCON, I,  J, N_part, info, NSTM, N, nf, nst, nt, nt1, nf_eff
        Integer, allocatable :: ipiv(:)
        COMPLEX (Kind=Kind(0.d0)) :: alpha,beta, Z, Z1
        TYPE(udv_state) :: udvlocal
        COMPLEX (Kind=Kind(0.d0)), Dimension(:,:), Allocatable ::  TP!, U, V
        COMPLEX (Kind=Kind(0.d0)), Dimension(:), Allocatable :: D
 
        !! TODO adapt to flavor symmetry
 
        if(udvl(1)%side .ne. "L" .and. udvl(1)%side .ne. "l" ) then
           write(error_unit,*) "Compute_Fermion_Det: calling wrong decompose"
           CALL terminate_on_error(error_global_updates,__file__,__line__)
        endif
 
        nstm = Size(udvst, 1)
        If (str_to_upper(storage) == "EMPTY" ) then
           DO nf_eff = 1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              if (projector) then
                 CALL udvl(nf_eff)%reset('l',wf_l(nf)%P)
              else
                 CALL udvl(nf_eff)%reset('l')
              endif
           ENDDO
           DO nst = nstm-1,1,-1
              nt1 = stab_nt(nst+1)
              nt  = stab_nt(nst  )
              !Write(6,*) 'Call wrapul', NT1,NT, udvl(1)%d(1), udvl(1)%N_part, udvl(1)%Ndim
              CALL wrapul(nt1,nt, udvl)
              Do nf_eff = 1,n_fl_eff
                 udvst(nst, nf_eff) = udvl(nf_eff)
              ENDDO
           ENDDO
           nt1 = stab_nt(1)
           CALL wrapul(nt1,0, udvl)
        Endif
 
        if (projector) then
           n_part=udvst(1,1)%N_part
           n_size=udvl(1)%ndim
           det_vec = 0.d0
           Allocate (tp(n_part,n_part), ipiv(n_part))
           do nf_eff=1,n_fl_eff
              nf=calc_fl_map(nf_eff)
              n_part=udvst(1,nf)%N_part
              do i=1,nstm-1
                 do n=1,n_part
#if !defined(STABLOG)
                    det_vec(n,nf_eff)=det_vec(n,nf_eff)+log(dble(udvst(i,nf)%D(n)))
#else
                    det_vec(n,nf_eff)=det_vec(n,nf_eff)+udvst(i,nf)%L(n)
#endif
                 enddo
              enddo
              do n=1,n_part
#if !defined(STABLOG)
                 det_vec(n,nf_eff)=det_vec(n,nf_eff)+log(dble(udvl(nf_eff)%D(n)))
#else
                 det_vec(n,nf_eff)=det_vec(n,nf_eff)+udvl(nf_eff)%L(n)
#endif
              enddo
              alpha=1.d0
              beta=0.d0
              CALL zgemm('C','N',n_part,n_part,n_size,alpha,udvl(nf_eff)%U(1,1),n_size,wf_r(nf)%P(1,1),n_size,beta,tp(1,1),n_part)
              ! ZGETRF computes an LU factorization of a general M-by-N matrix A
              ! using partial pivoting with row interchanges.
              call zgetrf(n_part, n_part, tp(1,1), n_part, ipiv, info)
              z=1.d0
              Do j=1,n_part
                 if (ipiv(j).ne.j) then
                    z = -z
                 endif
                 z =  z * tp(j,j)
              enddo
              phase_det(nf_eff) = z/abs(z)
              det_vec(1,nf) = det_vec(1,nf)+log(abs(z))
           enddo
           Deallocate(tp,ipiv)
           return
        endif
 
        !    N_size = SIZE(DL,1)
        n_size = udvl(1)%ndim
        ncon  = 0
        alpha = cmplx(1.d0,0.d0,kind(0.d0))
        beta  = cmplx(0.d0,0.d0,kind(0.d0))
        Allocate (tp(n_size,n_size),d(n_size))
        CALL udvlocal%alloc(n_size)
        Do nf_eff = 1,n_fl_eff
           nf=calc_fl_map(nf_eff)
           tp = udvl(nf_eff)%U !udvl stores U^dag instead of U !CT(udvl%U)
#if !defined(STABLOG)
#if !defined(STAB3)
           DO j = 1,n_size
              tp(:,j) = tp(:,j) +  udvl(nf_eff)%V(:,j)*udvl(nf_eff)%D(j)
           ENDDO
#else
           DO j = 1,n_size
              if ( dble(udvl(nf_eff)%D(j)) <= 1.d0 ) then
                 tp(:,j) = tp(:,j) +  udvl(nf_eff)%V(:,j)*udvl(nf_eff)%D(j)
              else
                 tp(:,j) = tp(:,j)/udvl(nf_eff)%D(j) +  udvl(nf_eff)%V(:,j)
              endif
           ENDDO
#endif
#else
           DO j = 1,n_size
              if ( udvl(nf_eff)%L(j) <= 0.d0 ) then
                 tp(:,j) = tp(:,j) +  udvl(nf_eff)%V(:,j)*cmplx(exp(udvl(nf_eff)%L(j)),0.d0,kind(0.d0))
              else
                 tp(:,j) = tp(:,j)*cmplx(exp(-udvl(nf_eff)%L(j)),0.d0,kind(0.d0)) +  udvl(nf_eff)%V(:,j)
              endif
           ENDDO
#endif
 
           !Call  UDV_WRAP_Pivot(TP,udvlocal%U, D, udvlocal%V, NCON,N_size,N_Size)
           Call  udv_wrap(tp,udvlocal%U, d, udvlocal%V, ncon )
           z  = det_c(udvlocal%V, n_size) ! Det destroys its argument
           !         Call MMULT(TP, udvl%U, udvlocal%U)
           CALL zgemm('C', 'N', n_size, n_size, n_size, alpha, udvl(nf_eff)%U(1,1), n_size, udvlocal%U(1,1), n_size, beta, tp, n_size)
           z1 = det_c(tp, n_size)
           phase_det(nf_eff)   = z*z1/abs(z*z1)
#if !defined(STABLOG)
#if !defined(STAB3)
           det_vec(:,nf_eff) = log(real(d(:)))
           det_vec(1,nf_eff) = log(real(d(1))) + log(abs(z*z1))
#else
           det_vec(1,nf_eff) = log(real(d(1))*abs(z*z1))
           if (dble(udvl(nf_eff)%D(1)) > 1.d0) det_vec(1,nf_eff)=det_vec(1,nf_eff)+log(dble(udvl(nf_eff)%D(1)))
           Do j=2,ndim
              if (dble(udvl(nf_eff)%D(j))<=1.d0) then
                 det_vec(j,nf_eff) = log(real(d(j)))
              else
                 det_vec(j,nf_eff) = log(real(d(j)))+log(dble(udvl(nf_eff)%D(j)))
              endif
           enddo
#endif
#else
           det_vec(:,nf_eff) = log(real(d(:)))
           det_vec(1,nf_eff) = log(real(d(1))*abs(z*z1))
           if (udvl(nf_eff)%L(1) > 0.d0) det_vec(1,nf_eff)=det_vec(1,nf_eff)+udvl(nf_eff)%L(1)
           Do j=2,ndim
              if (udvl(nf_eff)%L(j)<=0.d0) then
                 det_vec(j,nf_eff) = log(real(d(j)))
              else
                 det_vec(j,nf_eff) = log(real(d(j)))+udvl(nf_eff)%L(j)
              endif
           enddo
#endif
        Enddo
        Deallocate (tp)
        CALL udvlocal%dealloc
 
      end subroutine compute_fermion_det
 
 
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> Periodic boundary conditions required to defined master and slave for the tempering
!--------------------------------------------------------------------
      Integer function  npbc_tempering(n,Isize)
        implicit none
        Integer,  INTENT(IN)   :: isize,n
 
        npbc_tempering = n
        if (  npbc_tempering < 0       ) npbc_tempering = npbc_tempering +  isize
        if (  npbc_tempering > isize -1) npbc_tempering = npbc_tempering -  isize
 
      end function npbc_tempering
 
 
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> Allocates the Tempering_acceptance (Type Obser_vec) variable  that monitors the exchange acceptance
!--------------------------------------------------------------------
      Subroutine global_tempering_setup
        Integer    ::   N
        Character (len=64) ::  Filename
        n = 1
        filename = "Acc_Temp"
        Call obser_vec_make(tempering_acceptance,n,filename)
      End Subroutine global_tempering_setup
 
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> Initializes  Tempering_acceptance
!--------------------------------------------------------------------
 
      Subroutine global_tempering_init_obs
        Call obser_vec_init( tempering_acceptance )
      end Subroutine global_tempering_init_obs
 
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> Measures Tempering_acceptance
!--------------------------------------------------------------------
 
      Subroutine global_tempering_obser(toggle)
        Implicit none
 
        Logical, intent(in) :: toggle
        tempering_acceptance%N         =  tempering_acceptance%N + 1
        tempering_acceptance%Ave_sign  =  tempering_acceptance%Ave_sign + 1.d0
        if (toggle) tempering_acceptance%Obs_vec(1) = tempering_acceptance%Obs_vec(1) +  cmplx(1.d0,0.d0,kind(0.d0))
 
      end Subroutine global_tempering_obser
!--------------------------------------------------------------------
!> @author
!> The ALF Project contributors
!>
!> @brief
!> Prints Tempering_acceptance
!--------------------------------------------------------------------
      Subroutine global_tempering_pr
        Implicit none
        Call  print_bin_vec(tempering_acceptance,group_comm)
      end Subroutine global_tempering_pr
 
 
end Module global_mod