TABLE OF CONTENTS
ABINIT/testcode_ctqmc [ Functions ]
NAME
testcode_ctqmc
FUNCTION
Setup ultra simple hybridization to test CTQMC in simple situations.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (BAmadon) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
gtmp_nd gw_tmp_nd temp = temperature dmftqmc_l = number of times slices nflavor = number of flavor testrot = 0/1 if rotation of hybridization is tested or not testcode = 1 if tests are activated. opt = 1/2 if pre or postprocessing of CTQMC data.
OUTPUT
fw1_nd = non diagonal hybridization fw1 = hybridization umod = value of U
SIDE EFFECTS
gtmp_nd gw_tmp_nd
NOTES
PARENTS
qmc_prep_ctqmc
CHILDREN
SOURCE
45 #if defined HAVE_CONFIG_H 46 #include "config.h" 47 #endif 48 49 #include "abi_common.h" 50 51 subroutine testcode_ctqmc(dmftqmc_l,fw1_nd,fw1,gtmp_nd,gw_tmp_nd,levels_ctqmc,hybri_limit,& 52 & nflavor,opt,temp,testrot,testcode,umod) 53 54 use defs_basis 55 use defs_datatypes 56 use defs_abitypes 57 use m_errors 58 use m_ctqmc 59 use m_CtqmcInterface 60 use m_greenhyb 61 !use m_self, only : self_type,initialize_self,destroy_self,print_self,rw_self 62 use m_io_tools, only : flush_unit 63 64 !This section has been created automatically by the script Abilint (TD). 65 !Do not modify the following lines by hand. 66 #undef ABI_FUNC 67 #define ABI_FUNC 'testcode_ctqmc' 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 !scalars 74 integer, intent(in) :: dmftqmc_l,nflavor,testrot,testcode,opt 75 real(dp), intent(in) :: temp 76 real(dp), intent(out) :: umod(2,2) 77 complex(dpc), intent(inout) :: gw_tmp_nd(:,:,:) 78 real(dp), intent(inout) :: gtmp_nd(:,:,:) 79 complex(dpc), intent(out) :: fw1(:,:) 80 complex(dpc), intent(out) :: fw1_nd(:,:,:) 81 real(dp), intent(inout) :: levels_ctqmc(:) 82 complex(dpc), intent(inout) :: hybri_limit(:,:) 83 84 !Local variables ------------------------------ 85 character(len=500) :: message 86 integer :: ifreq, itau,realrot,simplehyb 87 real(dp) :: omega 88 real(dp) :: tbi1,tbi2,e2,tbi3,tbi4,e3,e4,tbi21,tbi12,e3b,e4b,tbi21b,tbi12b 89 complex(dpc) :: e1 90 ! arrays 91 complex(dpc) :: RR(2,2) 92 complex(dpc) :: RR1(2,2) 93 complex(dpc) :: RRi(2,2) 94 complex(dpc) :: RRt(2,2) 95 ! ************************************************************************ 96 if (testcode==0) return 97 if (nflavor/=2) then 98 write(message,'(2a)') ch10,' testcode nflavor.ne.2' 99 MSG_ERROR(message) 100 end if 101 102 simplehyb=2 103 simplehyb=1 104 simplehyb=3 105 !========================= 106 ! Built rotation matrix 107 !========================= 108 realrot=0 109 realrot=2 110 if (realrot==1) then 111 ! Real rotation 112 !========================= 113 RR(1,1) = SQRT(3.d0)/2.d0 114 RR(1,2) = -1.d0/2.d0 115 RR(2,1) = 1.d0/2.d0 116 RR(2,2) = SQRT(3.d0)/2.d0 117 else if (realrot==2) then 118 ! Real rotation 119 !========================= 120 RR(1,1) = SQRT(1.d0/2.d0) 121 RR(1,2) = -SQRT(1.d0/2.d0) 122 RR(2,1) = SQRT(1.d0/2.d0) 123 RR(2,2) = SQRT(1.d0/2.d0) 124 else 125 ! Complex rotation 126 !========================= 127 RR(1,1) = CMPLX(one,two) 128 RR(1,2) = CMPLX(one,one) 129 RR(2,1) = CMPLX(one,-one) 130 RR(2,2) = CMPLX(-one,two) 131 RR=RR/sqrt(seven) 132 end if 133 ! Check rotation is unitary 134 !========================== 135 RRi(1,1) = conjg(RR(1,1)) 136 RRi(1,2) = conjg(RR(2,1)) 137 RRi(2,1) = conjg(RR(1,2)) 138 RRi(2,2) = conjg(RR(2,2)) 139 RR1(:,:) = MATMUL ( RR(:,:) , RRi(:,:) ) 140 !write(6,*) "RR1",RR1 141 if(abs(RR1(1,1)-one).gt.tol7.or.abs(RR1(1,2)).gt.tol7.or.abs(RR1(2,2)-one).gt.tol7.or.abs(RR1(2,1)).gt.tol7) then 142 write(message,'(2a)') ch10,' testcode error in rotation matrix' 143 MSG_ERROR(message) 144 end if 145 146 147 !================================= 148 ! Built hybridization for CTQMC 149 !================================= 150 if (opt==1) then 151 152 ! Parameters: tight-binding + U 153 ! firt test of the code try umod=0, and (tbi1,tbi2,e1,e2)=(2,1,0.5,0.0) testrot=1 154 ! second test of the code try umod=four, and (tbi1,tbi2,e1,e2)=(2,1,0.0,0.0) testrot=1 155 !======================================================================================= 156 fw1_nd(:,:,:)= czero 157 tbi1=2.0_dp 158 tbi2=1.0_dp 159 tbi3=1.0_dp 160 tbi4=1.0_dp 161 tbi12=2.5_dp 162 tbi12b=2.5_dp 163 tbi21=2.5_dp 164 tbi21b=2.5_dp 165 e1=cmplx(0.0,0.0,8) 166 e2=zero 167 e3=0.2 168 e4=0.3 169 e3b=0.3 170 e4b=-0.2 171 umod(:,:)=0.d0 172 173 if(testrot==1.and.(abs(tbi1-tbi2)<tol6)) then 174 write(message,'(3a)') ch10,' testrot=1 with tbi1=tbi2 is equivalent' & 175 ,'to testrot=0: change testrot' 176 MSG_WARNING(message) 177 end if 178 ! Built fw1_nd 179 !============== 180 do ifreq=1,dmftqmc_l 181 182 omega=pi*temp*(two*float(ifreq)-1) 183 184 if(simplehyb==1) then 185 fw1_nd(ifreq,1,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 186 fw1_nd(ifreq,2,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 187 fw1(ifreq,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 188 fw1(ifreq,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 189 hybri_limit(1,1)=tbi1**2 190 hybri_limit(2,2)=tbi2**2 191 hybri_limit(1,2)=0.d0 192 hybri_limit(2,1)=0.d0 193 else if(simplehyb==2) then 194 fw1_nd(ifreq,1,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1)+tbi3**2/(dcmplx(0.d0,omega)-e3) 195 fw1_nd(ifreq,2,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2)+tbi4**2/(dcmplx(0.d0,omega)-e4) 196 fw1(ifreq,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 197 fw1(ifreq,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 198 else if(simplehyb==3) then 199 fw1_nd(ifreq,1,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 200 fw1_nd(ifreq,2,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 201 fw1_nd(ifreq,1,2) = tbi12**2/(dcmplx(0.d0,omega)-e3)+tbi12b**2/(dcmplx(0.d0,omega)-e3b) 202 fw1_nd(ifreq,2,1) = tbi21**2/(dcmplx(0.d0,omega)-e4)+tbi21b**2/(dcmplx(0.d0,omega)-e4b) 203 fw1(ifreq,1) = -umod(1,1)/two+tbi1**2/(dcmplx(0.d0,omega)-e1) 204 fw1(ifreq,2) = -umod(1,1)/two+tbi2**2/(dcmplx(0.d0,omega)-e2) 205 hybri_limit(1,1)=tbi1**2 206 hybri_limit(2,2)=tbi2**2 207 hybri_limit(1,2)=tbi12**2+tbi12b**2 208 hybri_limit(2,1)=tbi21**2+tbi21b**2 209 end if 210 write(132,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1)) 211 write(133,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2)) 212 write(134,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1)) 213 write(135,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2)) 214 write(1234,*) omega, real(fw1(ifreq,1)),aimag(fw1(ifreq,1)) 215 end do 216 ! Built level and limit of hybridization 217 !======================================= 218 levels_ctqmc(1:nflavor)=-umod(1,1)/two 219 220 write(std_out,*) "fw1_nd" 221 write(std_out,*) fw1_nd(1,1,1), fw1_nd(1,1,2) 222 write(std_out,*) fw1_nd(1,2,1), fw1_nd(1,2,2) 223 write(std_out,*) "fw1" 224 write(std_out,*) fw1(1,1), fw1(1,2) 225 write(std_out,*) fw1(2,1), fw1(2,2) 226 227 ! Rotate hybridization if testrot=1 228 !================================== 229 if(testrot==1) then 230 231 do ifreq=1,dmftqmc_l 232 RRt(:,:) = MATMUL ( RR(:,:) , fw1_nd(ifreq,:,:) ) 233 !write(6,*) "RRt" 234 !write(6,*) RRt(1,1), RRt(1,2) 235 !write(6,*) RRt(2,1), RRt(2,2) 236 RR1(:,:) = MATMUL ( RRt(:,:) , RRi(:,:) ) 237 !write(6,*) "RR1" 238 !write(6,*) RR1(1,1), RR1(1,2) 239 !write(6,*) RR1(2,1), RR1(2,2) 240 fw1_nd(ifreq,:,:)=RR1(:,:) 241 omega=pi*temp*(two*float(ifreq)+1) 242 write(3322,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1)) 243 write(232,*) omega,real(fw1_nd(ifreq,1,1)),aimag(fw1_nd(ifreq,1,1)) 244 write(233,*) omega,real(fw1_nd(ifreq,1,2)),aimag(fw1_nd(ifreq,1,2)) 245 write(234,*) omega,real(fw1_nd(ifreq,2,1)),aimag(fw1_nd(ifreq,2,1)) 246 write(235,*) omega,real(fw1_nd(ifreq,2,2)),aimag(fw1_nd(ifreq,2,2)) 247 end do 248 249 ! Rotate limit of hybridization 250 !======================================= 251 RRt(:,:) = MATMUL ( RR(:,:) , hybri_limit(:,:) ) 252 RR1(:,:) = MATMUL ( RRt(:,:) , RRi(:,:) ) 253 hybri_limit(:,:)=RR1(:,:) 254 255 end if 256 ! rajouter test real(fw1_nd(1,:,:)) doit etre diagonale 257 258 !====================================== 259 ! Rotate Green's function from CTQMC 260 !====================================== 261 else if(opt==2) then 262 263 write(std_out,*) "gw_tmp_nd" 264 write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2) 265 write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2) 266 ! Rotate Green's function back 267 !============================== 268 if(testrot==1) then 269 do ifreq=1,dmftqmc_l 270 RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gw_tmp_nd(ifreq,1:nflavor,1:nflavor) ) 271 RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) ) 272 gw_tmp_nd(ifreq,1:nflavor,1:nflavor)=RR1(1:nflavor,1:nflavor) 273 end do 274 275 write(std_out,*) "gw_tmp_nd after rotation" 276 write(std_out,*) gw_tmp_nd(1,1,1), gw_tmp_nd(1,1,2) 277 write(std_out,*) gw_tmp_nd(1,2,1), gw_tmp_nd(1,2,2) 278 279 do itau=1,dmftqmc_l 280 RRt(1:nflavor,1:nflavor) = MATMUL ( RRi(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) ) 281 RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RR(1:nflavor,1:nflavor) ) 282 gtmp_nd(itau,1:nflavor,1:nflavor)=real(RR1(1:nflavor,1:nflavor)) 283 end do 284 285 ! Rotate Green's function for comparison with testrot=1 286 !====================================================== 287 else if (testrot==0) then ! produce rotated green's function to compare to testrot=1 case 288 289 do itau=1,dmftqmc_l 290 RRt(1:nflavor,1:nflavor) = MATMUL ( RR(1:nflavor,1:nflavor),gtmp_nd(itau,1:nflavor,1:nflavor) ) 291 RR1(1:nflavor,1:nflavor) = MATMUL ( RRt(1:nflavor,1:nflavor),RRi(1:nflavor,1:nflavor) ) 292 write(444,*) real(itau-1)/(temp*real(dmftqmc_l)),real(RR1(1,1)),real(RR1(2,2)),real(RR1(1,2)),real(RR1(2,1)) 293 end do 294 295 end if 296 297 ! Print out rotated Green's function 298 !===================================== 299 do itau=1,dmftqmc_l 300 write(555,'(e14.5,4(2e14.5,3x))') real(itau-1)/(temp*real(dmftqmc_l)),gtmp_nd(itau,1,1),& 301 & gtmp_nd(itau,2,2),gtmp_nd(itau,1,2),gtmp_nd(itau,2,1) 302 end do 303 304 write(message,'(2a)') ch10,' testcode end of test calculation' 305 MSG_ERROR(message) 306 307 end if 308 close(444) 309 close(555) 310 311 end subroutine testcode_ctqmc