From e0ab2c006ef6375c5a92108fba2094b0554406e0 Mon Sep 17 00:00:00 2001 From: Jonas Schwab Date: Mon, 27 Jul 2020 17:23:16 +0200 Subject: [PATCH 1/2] Introduce Subroutine Obser_vec_measure and make subroutines for type obse_vec typebound. --- Prog/observables_mod.F90 | 131 ++++++++++++++++++++++++--------------- 1 file changed, 82 insertions(+), 49 deletions(-) diff --git a/Prog/observables_mod.F90 b/Prog/observables_mod.F90 index 74a773af..82424eb2 100644 --- a/Prog/observables_mod.F90 +++ b/Prog/observables_mod.F90 @@ -1,31 +1,31 @@ ! Copyright (C) 2016 - 2020 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 +! +! - 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. @@ -33,40 +33,47 @@ !> @author !> ALF-project ! -!> @brief -!> This module defines the Obser_Vec and Obser_Latt types and provides +!> @brief +!> This module defines the Obser_Vec and Obser_Latt types and provides !> routine to initialize them and to print out the bins ! !-------------------------------------------------------------------- Module Observables Use Files_mod - - Type Obser_Vec -!> Data structure for -!> < O_n > n : =1, size(Obs,1) - Integer :: N ! Number of measurements - real (Kind=Kind(0.d0)) :: Ave_Sign ! Averarge sign + + Type :: Obser_Vec +!> Data structure for +!> < O_n > n : =1, size(Obs,1) + !private + Integer :: N ! Number of measurements + real (Kind=Kind(0.d0)) :: Ave_Sign ! Averarge sign + Character (len=64) :: name ! Name of file in which the bins will be written out complex (Kind=Kind(0.d0)), pointer :: Obs_vec(:) ! Vector of observables - Character (len=64) :: File_Vec ! Name of file in which the bins will be written out + + contains + procedure :: make => Obser_vec_make + procedure :: init => Obser_vec_init + procedure :: print_bin => print_bin_vec + procedure :: measure => Obser_vec_measure end type Obser_Vec - + Type Obser_Latt -!> Data structure for -!> < O^{dagger}(i,tau)_n O(j,0)_m> - < O_n> -!> where it is assumed that translation symmetry as specified by the lattice Latt is present. -!> Obs_Latt(i-j,tau,n,m) = < O^{dagger}(i,tau)_n O(j,0)_m> +!> Data structure for +!> < O^{dagger}(i,tau)_n O(j,0)_m> - < O_n> +!> where it is assumed that translation symmetry as specified by the lattice Latt is present. +!> Obs_Latt(i-j,tau,n,m) = < O^{dagger}(i,tau)_n O(j,0)_m> !> Obs_Latt0(n) = < O_n> -!> For equal time correlation functions, tau runs from 1,1 -!> For unequal time correlation functions, tau runs from 1,Ltrot+1 +!> For equal time correlation functions, tau runs from 1,1 +!> For unequal time correlation functions, tau runs from 1,Ltrot+1 Integer :: N ! Number of measurements Real (Kind=Kind(0.d0)) :: Ave_Sign ! Averarge sign - complex (Kind=Kind(0.d0)), pointer :: Obs_Latt (:,:,:,:) ! i-j, tau, norb, norb - complex (Kind=Kind(0.d0)), pointer :: Obs_Latt0(:) ! norb + complex (Kind=Kind(0.d0)), pointer :: Obs_Latt (:,:,:,:) ! i-j, tau, norb, norb + complex (Kind=Kind(0.d0)), pointer :: Obs_Latt0(:) ! norb Character (len=64) :: File_Latt ! Name of file in which the bins will be written out end type Obser_Latt - + Contains @@ -75,7 +82,7 @@ Implicit none Type (Obser_Latt), intent(INOUT) :: Obs Integer, Intent(IN) :: Ns,Nt,No - Character (len=64), Intent(IN) :: Filename + Character (len=64), Intent(IN) :: Filename Allocate (Obs%Obs_Latt (Ns,Nt,No,No)) Allocate (Obs%Obs_Latt0(No) ) Obs%File_Latt = Filename @@ -90,29 +97,55 @@ Obs%N = 0 Obs%Ave_Sign = 0.d0 end subroutine Obser_Latt_Init - + !-------------------------------------------------------------------- Subroutine Obser_Vec_make(Obs,N,Filename) Implicit none - Type (Obser_vec), intent(INOUT) :: Obs + class(Obser_vec), intent(INOUT) :: Obs Integer, Intent(IN) :: N - Character (len=64), Intent(IN) :: Filename + Character (len=64), Intent(IN) :: Filename Allocate (Obs%Obs_vec(N)) - Obs%File_Vec = Filename + Obs%name = Filename end subroutine Obser_Vec_make !-------------------------------------------------------------------- - + Subroutine Obser_Vec_Init(Obs) Implicit none - Type (Obser_vec), intent(INOUT) :: Obs + class(Obser_vec), intent(INOUT) :: Obs Obs%Obs_vec = cmplx(0.d0,0.d0,kind(0.d0)) Obs%N = 0 Obs%Ave_Sign= 0.d0 end subroutine Obser_Vec_Init !-------------------------------------------------------------------- - + + Subroutine Obser_vec_measure(obs, value, Phase) + Implicit none + + class(Obser_vec), Intent(Inout) :: Obs + complex(Kind=Kind(0.d0)), Intent(In) :: value(:) ! Vector of observables + complex(Kind=Kind(0.d0)), Intent(IN), optional :: Phase + !Local + Complex (Kind=Kind(0.d0)) :: ZP, ZS + + obs%N = obs%N + 1 + + if ( present(Phase) ) then + ZP = PHASE/Real(Phase, kind(0.D0)) + ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) + + obs%Ave_sign = obs%Ave_sign + real(ZS, kind(0.D0)) + obs%obs_vec = obs%obs_vec + value *ZS*ZP + else + obs%Ave_sign = obs%Ave_sign + 1.d0 + obs%obs_vec = obs%obs_vec + value + endif + + end Subroutine Obser_vec_measure + +!-------------------------------------------------------------------- + Subroutine Print_bin_Latt(Obs,Latt,dtau,Group_Comm) Use Lattices_v3 #ifdef MPI @@ -128,7 +161,7 @@ ! Local Integer :: Ns,Nt, Norb, no, no1, I , Ntau Complex (Kind=Kind(0.d0)), allocatable :: Tmp(:,:,:,:) - Real (Kind=Kind(0.d0)) :: x_p(2) + Real (Kind=Kind(0.d0)) :: x_p(2) Complex (Kind=Kind(0.d0)) :: Sign_bin Character (len=64) :: File_pr, File_suff #ifdef MPI @@ -148,8 +181,8 @@ Ns = Size(Obs%Obs_Latt,1) Ntau = Size(Obs%Obs_Latt,2) Norb = Size(Obs%Obs_Latt,3) - if ( .not. (Latt%N == Ns ) ) then - Write(6,*) 'Error in Print_bin' + if ( .not. (Latt%N == Ns ) ) then + Write(6,*) 'Error in Print_bin' Stop endif If (Ntau == 1) then @@ -163,7 +196,7 @@ Obs%Obs_Latt0 = Obs%Obs_Latt0/dble(Obs%N*Ns*Ntau) Obs%Ave_sign = Obs%Ave_Sign /dble(Obs%N ) -#if defined(MPI) +#if defined(MPI) I = Ns*Ntau*Norb*Norb Tmp = cmplx(0.d0, 0.d0, kind(0.D0)) CALL MPI_REDUCE(Obs%Obs_Latt,Tmp,I,MPI_COMPLEX16,MPI_SUM, 0,Group_Comm,IERR) @@ -183,7 +216,7 @@ If (Irank_g == 0 ) then #endif -#if defined(TEMPERING) +#if defined(TEMPERING) write(File_pr ,'(A,I0,A,A,A)') "Temp_",igroup,"/",trim(Obs%File_Latt),trim(File_suff ) #endif @@ -204,7 +237,7 @@ Write(10,*) Obs%Obs_Latt0(no) enddo do I = 1,Latt%N - x_p = dble(Latt%listk(i,1))*Latt%b1_p + dble(Latt%listk(i,2))*Latt%b2_p + x_p = dble(Latt%listk(i,1))*Latt%b1_p + dble(Latt%listk(i,2))*Latt%b2_p Write(10,*) X_p(1), X_p(2) Do nt = 1,Ntau do no = 1,Norb @@ -215,7 +248,7 @@ enddo enddo close(10) -#if defined(MPI) +#if defined(MPI) Endif #endif @@ -231,7 +264,7 @@ #endif Implicit none - Type (Obser_vec), intent(Inout) :: Obs + class(Obser_vec), intent(Inout) :: Obs Integer, INTENT(IN) :: Group_Comm ! Local @@ -252,10 +285,10 @@ Obs%Obs_vec = Obs%Obs_vec /dble(Obs%N) Obs%Ave_sign = Obs%Ave_sign/dble(Obs%N) File_suff ="_scal" - File_pr = file_add(Obs%File_Vec,File_suff) + File_pr = file_add(Obs%name,File_suff) #if defined(MPI) - No = size(Obs%Obs_vec, 1) + No = size(Obs%Obs_vec, 1) Allocate (Tmp(No) ) Tmp = cmplx(0.d0,0.d0,kind(0.d0)) CALL MPI_REDUCE(Obs%Obs_vec,Tmp,No,MPI_COMPLEX16,MPI_SUM, 0,Group_Comm,IERR) @@ -269,13 +302,13 @@ if (Irank_g == 0 ) then #endif -#if defined(TEMPERING) - write(File_pr,'(A,I0,A,A,A)') "Temp_",igroup,"/",trim(Obs%File_Vec),trim(File_suff) +#if defined(TEMPERING) + write(File_pr,'(A,I0,A,A,A)') "Temp_",igroup,"/",trim(Obs%name),trim(File_suff) #endif Open (Unit=10,File=File_pr, status="unknown", position="append") WRITE(10,*) size(Obs%Obs_vec,1)+1, (Obs%Obs_vec(I), I=1,size(Obs%Obs_vec,1)), Obs%Ave_sign close(10) -#if defined(MPI) +#if defined(MPI) endif #endif -- GitLab From 638c9b8ae548d0c8ecc914df08b86a258aea59c1 Mon Sep 17 00:00:00 2001 From: Jonas Schwab Date: Mon, 27 Jul 2020 17:54:33 +0200 Subject: [PATCH 2/2] Apply scalar observables typebound procedures in Hamiltonian_Hubbard --- Prog/Global_mod.F90 | 286 +++++++++--------- .../Hamiltonian_Hubbard_include.h | 79 +++-- Prog/Hamiltonians/Hamiltonian_Hubbard_mod.F90 | 234 +++++++------- 3 files changed, 297 insertions(+), 302 deletions(-) diff --git a/Prog/Global_mod.F90 b/Prog/Global_mod.F90 index b61f240b..1afdd6df 100644 --- a/Prog/Global_mod.F90 +++ b/Prog/Global_mod.F90 @@ -1,67 +1,67 @@ ! Copyright (C) 2016 - 2018 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 +! +! - 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 +!> @author !> ALF-project ! -!> @brief +!> @brief !> Handles global updates and parallel tempering ! !-------------------------------------------------------------------- Module Global_mod Use Hamiltonian - Use MyMats + Use MyMats Use Operator_mod Use Control Use Observables Use Fields_mod - + Implicit none Type (Obser_Vec ), private :: Tempering_acceptance - + Contains #if defined(TEMPERING) !-------------------------------------------------------------------- -!> @author +!> @author !> ALF-project ! -!> @brief +!> @brief !> Handles parrallel tempering \n !> This subroutine is called only if the tempering flag is switched on \n !> Requires MPI @@ -99,7 +99,7 @@ Module Global_mod !> \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 @@ -124,7 +124,7 @@ Module Global_mod INTEGER :: NVAR END SUBROUTINE CGR end Interface - + ! Arguments COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT) :: Phase CLASS(UDV_State), intent(inout), allocatable, Dimension(:) :: udvl, udvr @@ -133,8 +133,8 @@ Module Global_mod 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, NT, NT1, NVAR,N, N1,N2, I, NC, I_Partner, n_step, N_count, N_part Type (Fields) :: nsigma_old @@ -156,7 +156,7 @@ Module Global_mod Integer :: Isize, Irank, Ierr, irank_g, isize_g, igroup Integer :: STATUS(MPI_STATUS_SIZE) - Character (Len=64) :: storage + Character (Len=64) :: storage CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) @@ -171,18 +171,18 @@ Module Global_mod 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 ( 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. + + ! Compute for each core the old weights. L_test = .false. if (Tempering_calc_det) then - ! Set old weight. + ! Set old weight. storage = "Full" Call Compute_Fermion_Det(Phase_det_old,Det_Vec_old, udvl, udvst, Stab_nt, storage) Phase_old = cmplx(1.d0,0.d0,kind=kind(0.d0)) @@ -225,32 +225,32 @@ Module Global_mod 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 + 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 + ! 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_REAL8, List_partner(IRANK), 0, & & nsigma%f , n, MPI_REAL8, 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 + ! 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) @@ -260,11 +260,11 @@ Module Global_mod Phase_new = Phase_new*Phase_det_new(nf) Enddo Call Op_phase(Phase_new,OP_V,Nsigma,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) - + & 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 @@ -272,7 +272,7 @@ Module Global_mod 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) @@ -282,12 +282,12 @@ Module Global_mod 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. + 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) @@ -299,10 +299,10 @@ Module Global_mod 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) + Call Control_upgrade_Temp (Toggle) If (toggle) then ! Move has been accepted if (Tempering_calc_det) then @@ -319,7 +319,7 @@ Module Global_mod nsigma_irank = nsigma_old_irank endif enddo - + ! Finalize if (Tempering_calc_det) then ! If move has been accepted, no use to recomute storage @@ -350,9 +350,9 @@ Module Global_mod CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf), udvl(nf)) Phase = Phase*Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + call Op_phase(Phase,OP_V,Nsigma,N_SUN) else - ! Send >> Phase, GR, udvr, udvl, udvst << to new node + ! 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 @@ -377,11 +377,11 @@ Module Global_mod 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 = 1,N_Fl CALL udvr(nf)%MPI_Sendrecv(nsigma_old_irank, 0, nsigma_irank, 0, STATUS, IERR) enddo @@ -398,13 +398,13 @@ Module Global_mod call nsigma_old%clear - + if (Tempering_calc_det) then - Deallocate ( Det_vec_old, Det_vec_new ) + 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 @@ -445,13 +445,13 @@ Module Global_mod !> \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 Implicit none - + Interface SUBROUTINE WRAPUL(NTAU1, NTAU, udvl) Use Hamiltonian @@ -470,7 +470,7 @@ Module Global_mod INTEGER :: NVAR END SUBROUTINE CGR end Interface - + ! Arguments COMPLEX (Kind=Kind(0.d0)), INTENT(INOUT) :: Phase CLASS (UDV_State), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: udvl, udvr @@ -478,8 +478,8 @@ Module Global_mod 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 Real (Kind=Kind(0.d0)) :: T0_Proposal_ratio, Weight @@ -490,22 +490,22 @@ Module Global_mod Complex (Kind=Kind(0.d0)) :: Ratio(2) Logical :: TOGGLE, L_Test - Real (Kind=Kind(0.d0)) :: size_clust + Real (Kind=Kind(0.d0)) :: size_clust Real (Kind=Kind(0.d0)) :: ratio_2_test - Character (Len=64) :: storage - - - + 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), Det_vec_new(NDIM,N_FL), Det_vec_test(NDIM,N_FL) ) + Allocate ( Det_vec_old(NDIM,N_FL), Det_vec_new(NDIM,N_FL), Det_vec_test(NDIM,N_FL) ) Allocate ( Phase_Det_new(N_FL), Phase_Det_old(N_FL) ) - + L_test = .false. - ! Set old weight. + ! Set old weight. storage = "Full" Call Compute_Fermion_Det(Phase_det_old,Det_Vec_old, udvl, udvst, Stab_nt, storage) Phase_old = cmplx(1.d0,0.d0,kind=kind(0.d0)) @@ -513,9 +513,9 @@ Module Global_mod Phase_old = Phase_old*Phase_det_old(nf) Enddo Call Op_phase(Phase_old,OP_V,Nsigma,N_SUN) - - If (L_test) then - ! Testing + + If (L_test) then + ! Testing Do nf = 1,N_FL if (Projector) then CALL udvr(nf)%reset('r',WF_R(nf)%P) @@ -529,7 +529,7 @@ Module Global_mod CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf), udvl(nf)) Phase = Phase*Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + call Op_phase(Phase,OP_V,Nsigma,N_SUN) Do Nf = 1,N_FL Call DET_C_LU(GR(:,:,nf),Det_vec_test(:,nf),Ndim) Z = Phase_det_old(nf) @@ -543,11 +543,11 @@ Module Global_mod 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. + ! 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 @@ -564,20 +564,20 @@ Module Global_mod Phase_new = Phase_new*Phase_det_new(nf) Enddo Call Op_phase(Phase_new,OP_V,Nsigma,N_SUN) - + Ratiotot = Compute_Ratio_Global(Phase_Det_old, Phase_Det_new, & - & Det_vec_old, Det_vec_new, nsigma_old, T0_Proposal_ratio, Ratio) - + & 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. + + + TOGGLE = .false. if ( Weight > ranf_wrap() ) Then TOGGLE = .true. Phase_old = Phase_new @@ -592,7 +592,7 @@ Module Global_mod Call Control_upgrade_Glob(TOGGLE,size_clust) endif Enddo - + If (NC > 0 ) then If (.not.TOGGLE) then DO nf = 1,N_FL @@ -621,17 +621,17 @@ Module Global_mod CALL CGR(Z, NVAR, GR(:,:,nf), udvr(nf), udvl(nf)) Phase = Phase*Z Enddo - call Op_phase(Phase,OP_V,Nsigma,N_SUN) + call Op_phase(Phase,OP_V,Nsigma,N_SUN) endif - - + + call nsigma_old%clear - Deallocate ( Det_vec_old , Det_vec_new, Det_vec_test ) + Deallocate ( Det_vec_old , Det_vec_new, Det_vec_test ) Deallocate ( Phase_Det_new, Phase_Det_old ) - - + + End Subroutine Global_Updates - + !-------------------------------------------------------------------- !> @author @@ -642,7 +642,7 @@ Module Global_mod !> !> @details !> \verbatim -!> Ratio_Global = T0(nsigma--> nsigma_old) W(nsigma)/T0(nsigma_old--> nsigma) W(nsigma_old) = +!> 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. @@ -663,21 +663,21 @@ Module Global_mod 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 + 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 + + ! Local Integer :: Nf, i, nt Complex (Kind=Kind(0.d0)) :: Z, Z1 Real (Kind=Kind(0.d0)) :: X, Ratio_2 - + Ratio = cmplx(0.d0,0.d0,kind(0.d0)) Ratio_2 = 0.d0 !X = 1.d0 @@ -693,14 +693,14 @@ Module Global_mod !Z = Z*Phase_Det_new(nf)/Phase_Det_old(nf) Ratio(1) = Ratio(1) * Phase_Det_new(nf)/Phase_Det_old(nf) enddo - !Z = Z**N_SUN + !Z = Z**N_SUN Ratio(1) = Ratio(1)**N_SUN Ratio_2 = real(N_SUN,kind(0.d0))*Ratio_2 - + 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) ) + !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) Enddo @@ -712,17 +712,17 @@ Module Global_mod !Z = Z * cmplx( Delta_S0_global(Nsigma_old),0.d0,kind(0.d0) ) !Z = Z * cmplx( T0_Proposal_ratio, 0.d0,kind(0.d0)) Ratio_2 = Ratio_2 + log(Delta_S0_global(Nsigma_old)) + log(T0_Proposal_ratio) - + Ratio(2) = Ratio_2 Compute_Ratio_Global = Ratio(1)*exp(Ratio(2)) - - + + end Function Compute_Ratio_Global - + !-------------------------------------------------------------------- -!> @author -!> -!> @brief +!> @author +!> +!> @brief !> Computes the fermion determinant from scratch or from the storage. !> !> @details @@ -742,7 +742,7 @@ Module Global_mod !> \endverbatim !> @param[out] Phase_det Complex, Dimension(N_FL) !> \verbatim -!> The phase, per flavor +!> The phase, per flavor !> \endverbatim !> @param[out] Det_vec Real, Dimension(:,N_FL) !> @param[in] Stab_nt Integer, Dimension(:) @@ -754,14 +754,14 @@ Module Global_mod !> 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 - + Implicit none Interface @@ -773,15 +773,15 @@ Module Global_mod Integer :: NTAU1, NTAU END SUBROUTINE WRAPUL end Interface - + 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 + Character (Len=64), Intent(IN) :: storage + - !> Local variables Integer :: N_size, NCON, I, J, N_part, info, NSTM, N, nf, nst, nt, nt1 Integer, allocatable :: ipiv(:) @@ -790,7 +790,7 @@ Module Global_mod COMPLEX (Kind=Kind(0.d0)), Dimension(:,:), Allocatable :: TP!, U, V COMPLEX (Kind=Kind(0.d0)), Dimension(:), Allocatable :: D - + if(udvl(1)%side .ne. "L" .and. udvl(1)%side .ne. "l" ) then write(*,*) "calling wrong decompose" Stop @@ -815,9 +815,9 @@ Module Global_mod ENDDO ENDDO NT1 = stab_nt(1) - CALL WRAPUL(NT1,0, udvl) + CALL WRAPUL(NT1,0, udvl) Endif - + if (Projector) then N_part=udvst(1,1)%N_part N_size=udvl(1)%ndim @@ -860,14 +860,14 @@ Module Global_mod 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) + CALL udvlocal%alloc(N_size) Do nf = 1,N_FL TP = udvl(nf)%U !udvl stores U^dag instead of U !CT(udvl%U) #if !defined(LOG) @@ -899,7 +899,7 @@ Module Global_mod 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)%U(1,1), N_size, udvlocal%U(1,1), N_size, beta, TP, N_size) - Z1 = Det_C(TP, N_size) + Z1 = Det_C(TP, N_size) Phase_det(nf) = Z*Z1/ABS(Z*Z1) #if !defined(LOG) #if !defined(STAB3) @@ -931,33 +931,33 @@ Module Global_mod Enddo Deallocate (TP) CALL udvlocal%dealloc - + end subroutine Compute_Fermion_Det - - + + !-------------------------------------------------------------------- -!> @author +!> @author !> The ALF Project contributors !> -!> @brief +!> @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 +!> @author !> The ALF Project contributors !> -!> @brief +!> @brief !> Allocates the Tempering_acceptance (Type Obser_vec) variable that monitors the exchange acceptance !-------------------------------------------------------------------- Subroutine Global_Tempering_setup @@ -965,26 +965,26 @@ Module Global_mod Character (len=64) :: Filename N = 1 Filename = "Acc_Temp" - Call Obser_Vec_make(Tempering_acceptance,N,Filename) + Call Tempering_acceptance%make(N,Filename) End Subroutine Global_Tempering_setup !-------------------------------------------------------------------- -!> @author +!> @author !> The ALF Project contributors !> -!> @brief +!> @brief !> Initializes Tempering_acceptance !-------------------------------------------------------------------- Subroutine Global_Tempering_init_obs - Call Obser_vec_Init( Tempering_acceptance ) + Call Tempering_acceptance%Init() end Subroutine Global_Tempering_init_obs !-------------------------------------------------------------------- -!> @author +!> @author !> The ALF Project contributors !> -!> @brief +!> @brief !> Measures Tempering_acceptance !-------------------------------------------------------------------- @@ -992,22 +992,24 @@ Module Global_mod 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)) - + if (toggle) then + call Tempering_acceptance%measure( [cmplx(1.d0,0.d0,kind(0.d0))] ) + else + call Tempering_acceptance%measure( [cmplx(0.d0,0.d0,kind(0.d0))] ) + endif + end Subroutine Global_Tempering_obser !-------------------------------------------------------------------- -!> @author +!> @author !> The ALF Project contributors !> -!> @brief +!> @brief !> Prints Tempering_acceptance !-------------------------------------------------------------------- Subroutine Global_Tempering_Pr Implicit none - Call Print_bin_Vec(Tempering_acceptance,Group_Comm) + Call Tempering_acceptance%print_bin(Group_Comm) end Subroutine Global_Tempering_Pr - - + + end Module Global_mod diff --git a/Prog/Hamiltonians/Hamiltonian_Hubbard_include.h b/Prog/Hamiltonians/Hamiltonian_Hubbard_include.h index 5a86b262..51e1cc49 100644 --- a/Prog/Hamiltonians/Hamiltonian_Hubbard_include.h +++ b/Prog/Hamiltonians/Hamiltonian_Hubbard_include.h @@ -1,9 +1,9 @@ - + !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> -!> @brief +!> @brief !> Prints out the bins. No need to change this routine. !------------------------------------------------------------------- Subroutine Pr_obs(LTAU) @@ -11,13 +11,13 @@ Implicit none Integer, Intent(In) :: Ltau - - !Local + + !Local Integer :: I Do I = 1,Size(Obs_scal,1) - Call Print_bin_Vec(Obs_scal(I),Group_Comm) + Call Obs_scal(I)%print_bin(Group_Comm) enddo Do I = 1,Size(Obs_eq,1) Call Print_bin_Latt(Obs_eq(I),Latt,dtau,Group_Comm) @@ -31,23 +31,23 @@ end Subroutine Pr_obs !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> -!> @brief +!> @brief !> Initializes observables to zero before each bins. No need to change !> this routine. !------------------------------------------------------------------- - Subroutine Init_obs(Ltau) + Subroutine Init_obs(Ltau) Implicit none Integer, Intent(In) :: Ltau - - ! Local + + ! Local Integer :: I Do I = 1,Size(Obs_scal,1) - Call Obser_vec_Init(Obs_scal(I)) + Call Obs_scal(I)%init() Enddo Do I = 1,Size(Obs_eq,1) @@ -63,7 +63,7 @@ end Subroutine Init_obs !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief @@ -101,8 +101,8 @@ Subroutine Global_move_tau(T0_Proposal_ratio, S0_ratio, & & Flip_list, Flip_length,Flip_value,ntau) - - Implicit none + + Implicit none Real (Kind = Kind(0.d0)),INTENT(OUT) :: T0_Proposal_ratio, S0_ratio Integer , INTENT(OUT) :: Flip_list(:) Real (Kind = Kind(0.d0)),INTENT(OUT) :: Flip_value(:) @@ -115,11 +115,11 @@ Real (Kind=Kind(0.d0)) :: T0_proposal Flip_length = nranf(4) - do n = 1,flip_length + do n = 1,flip_length n_op = nranf(size(OP_V,1)) Flip_list(n) = n_op Flip_value(n) = nsigma%flip(n_op,ntau) - If ( OP_V(n_op,1)%type == 1 ) then + If ( OP_V(n_op,1)%type == 1 ) then S0_ratio = S0(n_op,ntau,Flip_value(n)) T0_Proposal = 1.d0 - 1.d0/(1.d0+S0_ratio) ! No move prob If ( T0_Proposal > Ranf_wrap() ) then @@ -132,11 +132,11 @@ S0_ratio = 1.d0 endif Enddo - + end Subroutine Global_move_tau !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief @@ -149,22 +149,22 @@ !> the initial field !> \endverbatim !-------------------------------------------------------------------- - Subroutine Hamiltonian_set_nsigma(Initial_field) + Subroutine Hamiltonian_set_nsigma(Initial_field) Implicit none Real (Kind=Kind(0.d0)), allocatable, dimension(:,:), Intent(OUT) :: Initial_field - + end Subroutine Hamiltonian_set_nsigma !-------------------------------------------------------------------- -!> @author +!> @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 !-------------------------------------------------------------------- @@ -175,7 +175,7 @@ end Subroutine Overide_global_tau_sampling_parameters !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief @@ -185,7 +185,7 @@ !> a spin flip of Operator n on time slice nt !> @details !-------------------------------------------------------------------- - Real (Kind=Kind(0.d0)) function S0(n,nt,Hs_new) + Real (Kind=Kind(0.d0)) function S0(n,nt,Hs_new) Implicit none !> Operator index Integer, Intent(IN) :: n @@ -193,22 +193,22 @@ Integer, Intent(IN) :: nt !> New local field on time slice nt and operator index n Real (Kind=Kind(0.d0)), Intent(In) :: Hs_new - + Integer :: nt1,I !Write(6,*) "Hi1" - + S0 = 1.d0 - + end function S0 !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief !> Global moves -!> +!> !> @details -!> This routine generates a +!> 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 @@ -216,7 +216,7 @@ !> \endverbatim !> @param [OUT] T0_Proposal_ratio Real !> \verbatimam -!> T0_Proposal_ratio = T0( sigma_new -> sigma_old ) / T0( sigma_old -> sigma_new) +!> T0_Proposal_ratio = T0( sigma_new -> sigma_old ) / T0( sigma_old -> sigma_new) !> \endverbatim !> @param [OUT] Size_clust Real !> \verbatim @@ -225,7 +225,7 @@ !------------------------------------------------------------------- ! Functions for Global moves. These move are not implemented in this example. 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 @@ -245,15 +245,15 @@ n2 = nranf(N_tau) nsigma%f(n1,n2) = nsigma_old%flip(n1,n2) enddo - + End Subroutine Global_move !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief !> Computes the ratio exp(S0(new))/exp(S0(old)) -!> +!> !> @details !> This function computes the ratio \verbatim e^{-S0(nsigma)}/e^{-S0(nsigma_old)} \endverbatim !> @param [IN] nsigma_old, Type(Fields) @@ -264,17 +264,16 @@ Real (Kind=kind(0.d0)) Function Delta_S0_global(Nsigma_old) ! This function computes the ratio: e^{-S0(nsigma)}/e^{-S0(nsigma_old)} - Implicit none - + Implicit none + ! Arguments Type (Fields), INTENT(IN) :: nsigma_old ! Local Integer :: I,n,n1,n2,n3,n4,nt,nt1, nc_F, nc_J, nc_h_p, nc_h_m - + Delta_S0_global = 1.d0 end Function Delta_S0_global - diff --git a/Prog/Hamiltonians/Hamiltonian_Hubbard_mod.F90 b/Prog/Hamiltonians/Hamiltonian_Hubbard_mod.F90 index 8b1d26a0..539a8afd 100644 --- a/Prog/Hamiltonians/Hamiltonian_Hubbard_mod.F90 +++ b/Prog/Hamiltonians/Hamiltonian_Hubbard_mod.F90 @@ -1,5 +1,5 @@ ! Copyright (C) 2016 - 2020 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 @@ -19,22 +19,22 @@ ! 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 +! 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 +! - 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 +!> @author !> ALF-project !> -!> @brief +!> @brief !> This module defines the Hamiltonian and observables. Here, we have included a !> set of predefined Hamiltonians. They include the Hubbard and SU(N) tV models !> on honeycomb, pi-flux and square lattices. @@ -42,16 +42,16 @@ !> @details !> The public variables of this module are the following !> -!> +!> !> @param [public] OP_V !> \verbatim -!> Type (Operator), dimension(:,:), allocatable +!> Type (Operator), dimension(:,:), allocatable !> List of operators of type=1,2 and 3 describing the sequence of interactions on a time slice. !> The first index runs over this sequence. The second corresponds to the flavor index. \endverbatim -!> +!> !> @param [public] OP_T !> \verbatim -!> Type (Operator), dimension(:,:), allocatable +!> Type (Operator), dimension(:,:), allocatable !> Sequence of operators accounting for the hopping on a time slice. This can include various !> checkerboard decompositions. The first index runs over this sequence. The second corresponds to !> the flavor index. \endverbatim @@ -59,7 +59,7 @@ !> \f$ \prod_{\tau} \; \; \prod_{n=1}^{N_V}e^{V_n(\tau)} \prod_{n=1}^{N_T}e^{T_n} \f$. That is !> first the hopping and then the potential energy. !> -!>@param [public] WF_L +!>@param [public] WF_L !> \verbatim Type (WaveFunction), dimension(:), allocatable !> Left trial wave function. \endverbatim !> @@ -83,35 +83,35 @@ !> @param [public] N_SUN !> \verbatim Integer !> # of colors. Propagation is color independent. \endverbatim -!> +!> !> @param [public] Ltrot !> \verbatim Integer !> Available measurment interval in units of Delta Tau. \endverbatim !> -!> @param [public] Thtrot +!> @param [public] Thtrot !> \verbatim Integer !> Effective projection parameter in units of Delta Tau. (Only relevant if projective option is turned on) \endverbatim !> !> @param [public] Projector !> \verbatim Logical !> Flag for projector. If true then the total number of time slices will correspond to Ltrot + 2*Thtrot \endverbatim -!> -!> @param [public] Group_Comm +!> +!> @param [public] Group_Comm !> \verbatim Integer !> Defines MPI communicator \endverbatim ! !> @param [public] Symm !> \verbatim Logical \endverbatim !> If set to true then the green functions will be symmetrized -!> before being sent to the Obser, ObserT subroutines. +!> before being sent to the Obser, ObserT subroutines. !> In particular, the transformation, \f$ \tilde{G} = e^{-\Delta \tau T /2 } G e^{\Delta \tau T /2 } \f$ !> will be carried out and \f$ \tilde{G} \f$ will be sent to the Obser and ObserT subroutines. Note that !> if you want to use this feature, then you have to be sure the hopping and interaction terms are decomposed !> symmetrically. If Symm is true, the propagation reads: !> \f$ \prod_{\tau} \; \; \prod_{n=N_T}^{1}e^{T_n/2} \prod_{n=1}^{N_V}e^{V_n(\tau)} \prod_{n=1}^{N_T}e^{T_n/2} \f$ !> -!> -!> You still have to add some docu for the other private variables in this module. +!> +!> You still have to add some docu for the other private variables in this module. !> !-------------------------------------------------------------------- @@ -119,8 +119,8 @@ Use Operator_mod Use WaveFunction_mod - Use Lattices_v3 - Use MyMats + Use Lattices_v3 + Use MyMats Use Random_Wrap Use Files_mod Use Matrix @@ -129,11 +129,11 @@ Use Predefined_Hoppings Use LRC_Mod - + Implicit none - - Type (Operator), dimension(:,:), allocatable :: Op_V + + Type (Operator), dimension(:,:), allocatable :: Op_V Type (Operator), dimension(:,:), allocatable :: Op_T Type (WaveFunction), dimension(:), allocatable :: WF_L Type (WaveFunction), dimension(:), allocatable :: WF_R @@ -142,7 +142,7 @@ Integer :: N_FL Integer :: N_SUN Integer :: Ltrot - Integer :: Thtrot + Integer :: Thtrot Logical :: Projector Integer :: Group_Comm Logical :: Symm @@ -153,7 +153,7 @@ Integer, private :: L1, L2 Type (Hopping_Matrix_type), Allocatable, private :: Hopping_Matrix(:) real (Kind=Kind(0.d0)), private :: ham_T , ham_U, Ham_chem - real (Kind=Kind(0.d0)), private :: ham_T2, ham_U2, ham_Tperp ! For Bilayers + real (Kind=Kind(0.d0)), private :: ham_T2, ham_U2, ham_Tperp ! For Bilayers real (Kind=Kind(0.d0)), private :: Phi_Y, Phi_X Integer , private :: N_Phi real (Kind=Kind(0.d0)), private :: Dtau, Beta, Theta @@ -166,19 +166,19 @@ Type (Obser_Vec ), private, dimension(:), allocatable :: Obs_scal Type (Obser_Latt), private, dimension(:), allocatable :: Obs_eq Type (Obser_Latt), private, dimension(:), allocatable :: Obs_tau - - contains + + contains !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief !> Sets the Hamiltonian !-------------------------------------------------------------------- Subroutine Ham_Set - + #if defined (MPI) || defined(TEMPERING) Use mpi #endif @@ -186,19 +186,19 @@ integer :: ierr, N_part, nf Character (len=64) :: file_info, file_para - - + + ! L1, L2, Lattice_type, List(:,:), Invlist(:,:) --> Lattice information ! Ham_T, Chem, Phi_X, XB_B, Checkerboard, Symm --> Hopping ! Interaction --> Model - ! Simulation type --> Finite T or Projection Symmetrize Trotter. - + ! Simulation type --> Finite T or Projection Symmetrize Trotter. + NAMELIST /VAR_Lattice/ L1, L2, Lattice_type, Model NAMELIST /VAR_Model_Generic/ Checkerboard, N_SUN, N_FL, Phi_X, Phi_Y, Symm, Bulk, N_Phi, Dtau, Beta, Theta, Projector NAMELIST /VAR_Hubbard/ ham_T, ham_chem, ham_U, ham_T2, ham_U2, ham_Tperp, Mz - + #ifdef MPI Integer :: Isize, Irank, irank_g, isize_g, igroup @@ -209,7 +209,7 @@ Checkerboard = .false. Symm = .false. Projector = .false. - Bulk = .true. + Bulk = .true. Phi_X = 0.d0 Phi_Y = 0.d0 N_Phi = 0 @@ -217,7 +217,7 @@ Ham_Tperp = 0.d0 Ham_U2 = 0.d0 - + #ifdef MPI CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR) CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR) @@ -227,9 +227,9 @@ #endif File_Para = "parameters" File_info = "info" -#if defined(TEMPERING) +#if defined(TEMPERING) write(File_para,'(A,I0,A)') "Temp_",igroup,"/parameters" - write(File_info,'(A,I0,A)') "Temp_",igroup,"/info" + write(File_info,'(A,I0,A)') "Temp_",igroup,"/info" #endif #ifdef MPI @@ -258,7 +258,7 @@ else N_FL = 1 endif - + #ifdef MPI Endif CALL MPI_BCAST(L1 ,1 ,MPI_INTEGER, 0,Group_Comm,ierr) @@ -289,11 +289,11 @@ ! Setup the Bravais lattice Call Ham_Latt - + ! Setup the hopping / single-particle part Call Ham_Hop - - + + ! Setup the interaction. call Ham_V @@ -302,12 +302,12 @@ #endif OPEN(Unit = 50,file=file_info,status="unknown",position="append") Write(50,*) '=====================================' - Write(50,*) 'Model is : ', Model + Write(50,*) 'Model is : ', Model Write(50,*) 'Lattice is : ', Lattice_type Write(50,*) '# of orbitals : ', Ndim Write(50,*) 'Flux_1 : ', Phi_X Write(50,*) 'Flux_2 : ', Phi_Y - If (Bulk) then + If (Bulk) then Write(50,*) 'Twist as phase factor in bulk' Else Write(50,*) 'Twist as boundary condition' @@ -337,18 +337,18 @@ Write(50,*) 't : ', Ham_T Write(50,*) 'Ham_U : ', Ham_U Write(50,*) 'Ham_chem : ', Ham_chem - Close(50) + Close(50) #ifdef MPI Endif #endif ! Setup the trival wave function, in case of a projector approach if (Projector) Call Ham_Trial(File_info) - + end Subroutine Ham_Set - + !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief @@ -357,21 +357,21 @@ Subroutine Ham_Latt Use Predefined_Lattices - + Implicit none ! Use predefined stuctures or set your own lattice. Call Predefined_Latt(Lattice_type, L1,L2,Ndim, List,Invlist,Latt,Latt_Unit) - + end Subroutine Ham_Latt !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief !> Sets the Hopping !-------------------------------------------------------------------- Subroutine Ham_Hop - + Implicit none Real (Kind=Kind(0.d0) ) :: Ham_Lambda = 0.d0 @@ -395,7 +395,7 @@ Ham_T2_vec = Ham_T2 Ham_Lambda_vec = Ham_Lambda N_Phi_vec = N_Phi - + Select case (Lattice_type) Case ("Square") Call Set_Default_hopping_parameters_square(Hopping_Matrix,Ham_T_vec, Ham_Chem_vec, Phi_X_vec, Phi_Y_vec, & @@ -418,15 +418,15 @@ & List, Invlist, Latt, Latt_unit ) end Select - + Call Predefined_Hoppings_set_OPT(Hopping_Matrix,List,Invlist,Latt, Latt_unit, Dtau, Checkerboard, Symm, OP_T ) - + Deallocate (Ham_T_vec, Ham_T2_vec, Ham_Tperp_vec, Ham_Chem_vec, Phi_X_vec, Phi_Y_vec, & & N_Phi_vec, Ham_Lambda_vec ) - + end Subroutine Ham_Hop !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief @@ -440,10 +440,10 @@ #endif Use Predefined_Trial - Implicit none + Implicit none Character (len=64), intent(in) :: file_info - + Integer :: N_part, nf #ifdef MPI Integer :: Isize, Irank, irank_g, isize_g, igroup, ierr @@ -460,7 +460,7 @@ Call Predefined_TrialWaveFunction(Lattice_type ,Ndim, List,Invlist,Latt, Latt_unit, & & N_part, N_FL, WF_L, WF_R) - + #ifdef MPI If (Irank_g == 0) then #endif @@ -477,7 +477,7 @@ end Subroutine Ham_Trial !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief @@ -486,12 +486,12 @@ Subroutine Ham_V Use Predefined_Int - Implicit none - + Implicit none + Integer :: nf, I, I1, I2, nc, nc1, J, no Real (Kind=Kind(0.d0)) :: X Real (Kind=Kind(0.d0)), allocatable :: Ham_U_vec(:) - + Allocate (Ham_U_vec(Latt_unit%Norb)) Ham_U_vec(:) = Ham_U @@ -505,7 +505,7 @@ Do I1 = 1,Latt%N do no = 1, Latt_unit%Norb I = invlist(I1,no) - Call Predefined_Int_U_MZ ( OP_V(I,1), OP_V(I,2), I, DTAU, Ham_U_vec(no) ) + Call Predefined_Int_U_MZ ( OP_V(I,1), OP_V(I,2), I, DTAU, Ham_U_vec(no) ) enddo enddo else @@ -517,21 +517,21 @@ Enddo Enddo Endif - + Deallocate (Ham_U_vec) - + end Subroutine Ham_V !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> !> @brief !> Specifiy the equal time and time displaced observables !> @details !-------------------------------------------------------------------- - Subroutine Alloc_obs(Ltau) + Subroutine Alloc_obs(Ltau) Implicit none !> Ltau=1 if time displaced correlations are considered. @@ -554,11 +554,11 @@ case (4) N = 1; Filename ="Ener" case default - Write(6,*) ' Error in Alloc_obs ' + Write(6,*) ' Error in Alloc_obs ' end select - Call Obser_Vec_make(Obs_scal(I),N,Filename) + Call Obs_scal(I)%make(N,Filename) enddo - + ! Equal time correlators If ( Mz ) Then Allocate ( Obs_eq(5) ) @@ -575,13 +575,13 @@ case (5) Ns = Latt%N; No = Norb; Filename ="Den" case default - Write(6,*) ' Error in Alloc_obs ' + Write(6,*) ' Error in Alloc_obs ' end select Nt = 1 Call Obser_Latt_make(Obs_eq(I),Ns,Nt,No,Filename) enddo - - If (Ltau == 1) then + + If (Ltau == 1) then ! Equal time correlators Allocate ( Obs_tau(5) ) Do I = 1,Size(Obs_tau,1) @@ -597,7 +597,7 @@ case (5) Ns = Latt%N; No = Norb; Filename ="Den" case default - Write(6,*) ' Error in Alloc_obs ' + Write(6,*) ' Error in Alloc_obs ' end select Nt = Ltrot+1-2*Thtrot Call Obser_Latt_make(Obs_tau(I),Ns,Nt,No,Filename) @@ -615,13 +615,13 @@ case (3) Ns = Latt%N; No = Norb; Filename ="Den" case default - Write(6,*) ' Error in Alloc_obs ' + Write(6,*) ' Error in Alloc_obs ' end select Nt = 1 Call Obser_Latt_make(Obs_eq(I),Ns,Nt,No,Filename) enddo - - If (Ltau == 1) then + + If (Ltau == 1) then ! Equal time correlators Allocate ( Obs_tau(3) ) Do I = 1,Size(Obs_tau,1) @@ -633,56 +633,56 @@ case (3) Ns = Latt%N; No = Norb; Filename ="Den" case default - Write(6,*) ' Error in Alloc_obs ' + Write(6,*) ' Error in Alloc_obs ' end select Nt = Ltrot+1-2*Thtrot Call Obser_Latt_make(Obs_tau(I),Ns,Nt,No,Filename) enddo endif endif - + End Subroutine Alloc_obs !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> -!> @brief +!> @brief !> Computes equal time observables !> @details -!> @param [IN] Gr Complex(:,:,:) +!> @param [IN] Gr Complex(:,:,:) !> \verbatim !> Green function: Gr(I,J,nf) = on time slice ntau !> \endverbatim !> @param [IN] Phase Complex !> \verbatim -!> Phase +!> Phase !> \endverbatim !> @param [IN] Ntau Integer !> \verbatim -!> Time slice +!> Time slice !> \endverbatim !------------------------------------------------------------------- subroutine Obser(GR,Phase,Ntau) Use Predefined_Obs - + Implicit none - + Complex (Kind=Kind(0.d0)), INTENT(IN) :: GR(Ndim,Ndim,N_FL) Complex (Kind=Kind(0.d0)), Intent(IN) :: PHASE Integer, INTENT(IN) :: Ntau - - !Local + + !Local Complex (Kind=Kind(0.d0)) :: GRC(Ndim,Ndim,N_FL), ZK Complex (Kind=Kind(0.d0)) :: Zrho, Zkin, ZPot, Z, ZP,ZS, ZZ, ZXY Integer :: I,J, imj, nf, dec, I1, J1, no_I, no_J,n Real (Kind=Kind(0.d0)) :: X - + ZP = PHASE/Real(Phase, kind(0.D0)) ZS = Real(Phase, kind(0.D0))/Abs(Real(Phase, kind(0.D0))) - - + + Do nf = 1,N_FL Do I = 1,Ndim Do J = 1,Ndim @@ -693,17 +693,12 @@ Enddo ! GRC(i,j,nf) = < c^{dagger}_{j,nf } c_{j,nf } > - ! Compute scalar observables. - Do I = 1,Size(Obs_scal,1) - Obs_scal(I)%N = Obs_scal(I)%N + 1 - Obs_scal(I)%Ave_sign = Obs_scal(I)%Ave_sign + Real(ZS,kind(0.d0)) - Enddo - + Zkin = cmplx(0.d0, 0.d0, kind(0.D0)) - Call Predefined_Hoppings_Compute_Kin(Hopping_Matrix,List,Invlist, Latt, Latt_unit, GRC, ZKin) + Call Predefined_Hoppings_Compute_Kin(Hopping_Matrix,List,Invlist, Latt, Latt_unit, GRC, ZKin) Zkin = Zkin* dble(N_SUN) - Obs_scal(1)%Obs_vec(1) = Obs_scal(1)%Obs_vec(1) + Zkin *ZP* ZS + call Obs_scal(1)%measure([Zkin], Phase) ZPot = cmplx(0.d0, 0.d0, kind(0.D0)) @@ -720,19 +715,18 @@ Enddo Zpot = Zpot*ham_U Endif - Obs_scal(2)%Obs_vec(1) = Obs_scal(2)%Obs_vec(1) + Zpot * ZP*ZS + call Obs_scal(2)%measure([Zpot], Phase) Zrho = cmplx(0.d0,0.d0, kind(0.D0)) Do nf = 1,N_FL Do I = 1,Ndim - Zrho = Zrho + Grc(i,i,nf) + Zrho = Zrho + Grc(i,i,nf) enddo enddo Zrho = Zrho* dble(N_SUN) - Obs_scal(3)%Obs_vec(1) = Obs_scal(3)%Obs_vec(1) + Zrho * ZP*ZS - - Obs_scal(4)%Obs_vec(1) = Obs_scal(4)%Obs_vec(1) + (Zkin + Zpot)*ZP*ZS + call Obs_scal(3)%measure([Zrho], Phase) + call Obs_scal(4)%measure([Zkin + Zpot], Phase) ! Standard two-point correlations If ( Mz ) then @@ -744,31 +738,31 @@ Call Predefined_Obs_eq_SpinSUN_measure( Latt, Latt_unit, List, GR, GRC, N_SUN, ZS, ZP, Obs_eq(2) ) Call Predefined_Obs_eq_Den_measure ( Latt, Latt_unit, List, GR, GRC, N_SUN, ZS, ZP, Obs_eq(3) ) endif - + end Subroutine Obser !-------------------------------------------------------------------- -!> @author +!> @author !> ALF Collaboration !> -!> @brief +!> @brief !> Computes time displaced observables !> @details !> @param [IN] NT, Integer !> \verbatim !> Imaginary time !> \endverbatim -!> @param [IN] GT0, GTT, G00, GTT, Complex(:,:,:) +!> @param [IN] GT0, GTT, G00, GTT, Complex(:,:,:) !> \verbatim !> Green functions: -!> GT0(I,J,nf) = -!> G0T(I,J,nf) = -!> G00(I,J,nf) = -!> GTT(I,J,nf) = +!> GT0(I,J,nf) = +!> G0T(I,J,nf) = +!> G00(I,J,nf) = +!> GTT(I,J,nf) = !> \endverbatim !> @param [IN] Phase Complex !> \verbatim -!> Phase +!> Phase !> \endverbatim !------------------------------------------------------------------- Subroutine ObserT(NT, GT0,G0T,G00,GTT, PHASE) @@ -776,11 +770,11 @@ Use Predefined_Obs Implicit none - + Integer , INTENT(IN) :: NT Complex (Kind=Kind(0.d0)), INTENT(IN) :: GT0(Ndim,Ndim,N_FL),G0T(Ndim,Ndim,N_FL),G00(Ndim,Ndim,N_FL),GTT(Ndim,Ndim,N_FL) Complex (Kind=Kind(0.d0)), INTENT(IN) :: Phase - + !Locals Complex (Kind=Kind(0.d0)) :: Z, ZP, ZS, ZZ, ZXY Real (Kind=Kind(0.d0)) :: X @@ -801,10 +795,10 @@ Call Predefined_Obs_tau_SpinSUN_measure( Latt, Latt_unit, List, NT, GT0,G0T,G00,GTT, N_SUN, ZS, ZP, Obs_tau(2) ) Call Predefined_Obs_tau_Den_measure ( Latt, Latt_unit, List, NT, GT0,G0T,G00,GTT, N_SUN, ZS, ZP, Obs_tau(3) ) endif - + end Subroutine OBSERT -#include "Hamiltonian_Hubbard_include.h" +#include "Hamiltonian_Hubbard_include.h" + - end Module Hamiltonian -- GitLab