TABLE OF CONTENTS


ABINIT/testcode_ctqmc [ Functions ]

[ Top ] [ 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