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.



!--------------------------------------------------------------------
!
!
!--------------------------------------------------------------------

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)
!--------------------------------------------------------------------
!
!--------------------------------------------------------------------
      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


        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
        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
        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
              ratiotot = ham%Delta_S0_global(nsigma_old)
              ratio(1) = ratiotot
              ratio(2) = 0
           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

           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


!--------------------------------------------------------------------
!--------------------------------------------------------------------
      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


!--------------------------------------------------------------------
!--------------------------------------------------------------------
      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

        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)
        ratio(2) = ratio(2) + log(ham%Delta_S0_global(nsigma_old)) + log(t0_proposal_ratio)

        compute_ratio_global = ratio(1)*exp(ratio(2))


      end Function compute_ratio_global

!--------------------------------------------------------------------
!--------------------------------------------------------------------
      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


        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


!--------------------------------------------------------------------
!--------------------------------------------------------------------
      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


!--------------------------------------------------------------------
!--------------------------------------------------------------------
      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

!--------------------------------------------------------------------
!--------------------------------------------------------------------

      Subroutine global_tempering_init_obs
        Call obser_vec_init( tempering_acceptance )
      end Subroutine global_tempering_init_obs

!--------------------------------------------------------------------
!--------------------------------------------------------------------

      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
!--------------------------------------------------------------------
!--------------------------------------------------------------------
      Subroutine global_tempering_pr
        Implicit none
        Call  print_bin_vec(tempering_acceptance,group_comm)
      end Subroutine global_tempering_pr


end Module global_mod