TABLE OF CONTENTS
ABINIT/entropyrec [ Functions ]
NAME
entropyrec
FUNCTION
This routine computes the local part of the entropy at a point using a path integral, in the recursion method.
COPYRIGHT
Copyright (C) 2008-2018 ABINIT group ( ). 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
an, bn2 : coefficient given by the recursion. nrec=order of recursion trotter=trotter parameter multce=a multiplicator for computing entropy ; 2 for non-spin-polarized system debug_rec=debug variable n_pt_integ=number of points of integration for the path integral xmax =max point of integration on the real axis
OUTPUT
ent_out=entropy at the point ent_out1,ent_out2,ent_out3,ent_out4=debug entropy at the point
PARENTS
vtorhorec
CHILDREN
timab,wrtout
NOTES
at this time :
- multce should be not used
- the routine should be integraly rewrited and use the routine recursion.
- only modified for p /= 0
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 51 subroutine entropyrec(an,bn2,nrec,trotter,ent_out,multce,debug_rec, & 52 & n_pt_integ,xmax,& 53 & ent_out1,ent_out2,ent_out3,ent_out4) 54 55 use m_profiling_abi 56 57 use defs_basis 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'entropyrec' 63 use interfaces_14_hidewrite 64 use interfaces_18_timing 65 !End of the abilint section 66 67 implicit none 68 69 !Arguments ------------------------------- 70 !scalars 71 integer,intent(in) :: n_pt_integ,nrec,trotter 72 logical,intent(in) :: debug_rec 73 real(dp), intent(in) :: multce,xmax 74 real(dp),intent(out) :: ent_out,ent_out1,ent_out2,ent_out3,ent_out4 75 !arrays 76 real(dp),intent(in) :: an(0:nrec),bn2(0:nrec) 77 78 !Local variables------------------------------- 79 !scalars 80 integer, parameter :: level = 7 81 integer, save :: first_en = 1 82 integer :: ii,kk,n_pt_integ_path2,n_pt_integ_path3 83 real(dp) :: arg,epsilo,step,twotrotter,xmin,dr_step 84 complex(dpc) :: D,Dnew,Dold,N,Nnew,Nold,dz_path,ent_acc,ent_acc1,ent_acc2 85 complex(dpc) :: ent_acc3,ent_acc4 86 complex(dpc) :: funczero,z_path,zj 87 complex(dpc) ::delta_calc 88 character(len=500) :: msg 89 !arrays 90 real(dp) :: tsec(2) 91 real(dp) :: iif,factor 92 93 94 ! ************************************************************************* 95 96 97 call timab(610,1,tsec) 98 99 !structured debugging if debug_rec=T : print detailled result the first time we enter entropyrec 100 101 if(debug_rec .and. first_en==1)then 102 write(msg,'(a)')' ' 103 call wrtout(std_out,msg,'PERS') 104 write(msg,'(a)')' entropyrec : enter ' 105 call wrtout(std_out,msg,'PERS') 106 write(msg,'(a,i6)')'n_pt_integ ' , n_pt_integ 107 call wrtout(std_out,msg,'COLL') 108 end if 109 110 ent_out = zero 111 ent_out1 = zero 112 ent_out2 = zero 113 ent_out3 = zero 114 ent_out4 = zero 115 ent_acc = czero 116 ent_acc1 = czero 117 ent_acc2 = czero 118 ent_acc3 = czero 119 ent_acc4 = czero 120 121 !path parameters 122 twotrotter = max(two*real(trotter,dp),one) 123 if(trotter==0)then 124 factor = tol5 125 arg =pi*three_quarters 126 zj = cmplx(-one,one-sin(arg),dp) 127 else 128 factor = xmax/ten 129 arg = pi/twotrotter 130 zj = cmplx( cos(arg) , sin(arg),dp ) 131 end if 132 133 epsilo = factor*sin( arg ) 134 xmin = factor*cos( arg ) 135 step = (xmax-xmin)/real(n_pt_integ,dp) 136 137 !#################################################################### 138 ![xmax + i*epsilo,xmin + i*epsilo] 139 dr_step = one/real(n_pt_integ,dp) 140 path1: do ii = 0,n_pt_integ 141 z_path = cmplx(xmin+real(ii,dp)*(xmax-xmin)*dr_step,epsilo,dp) 142 dz_path = -cmplx((xmax-xmin)*dr_step,zero,dp) 143 144 Nold = czero 145 Dold = cone 146 N = cone 147 D = z_path - cmplx(an(0),zero,dp) 148 149 do kk=1,nrec 150 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 151 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 152 153 Nold = N 154 Dold = D 155 N = Nnew 156 D = Dnew 157 158 if(kk/=nrec)then 159 if((bn2(kk+1)<tol14))exit 160 end if 161 end do 162 163 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 164 delta_calc = func1_rec(z_path**twotrotter)*(N/D)*dz_path 165 if(ii==0.or.ii==n_pt_integ)then 166 ent_acc = ent_acc + half*delta_calc 167 ent_acc1 = ent_acc1 + half*delta_calc 168 else 169 ent_acc = ent_acc + delta_calc 170 ent_acc1 = ent_acc1 + delta_calc 171 end if 172 end do path1 173 174 175 !#################################################################### 176 ![1/2zj,0] 177 if(epsilo/step>100.d0)then 178 n_pt_integ_path2 = int((factor*abs(zj))/step)+1 179 else 180 n_pt_integ_path2 = 100 181 end if 182 183 if(trotter/=0)then 184 n_pt_integ_path3 = 0 185 dr_step = one/real(n_pt_integ_path2,dp) 186 dz_path = -cmplx(xmin,epsilo,dp)*dr_step 187 path5: do ii = 0,n_pt_integ_path2 188 z_path = cmplx(real(ii,dp)*xmin,real(ii,dp)*epsilo,dp)*dr_step 189 if(abs(z_path)>tol14)then 190 Nold = czero 191 Dold = cone 192 N = cone 193 D = z_path - cmplx(an(0),zero,dp) 194 do kk=1,nrec 195 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 196 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 197 Nold = N 198 Dold = D 199 N = Nnew 200 D = Dnew 201 if(kk/=nrec)then 202 if((bn2(kk+1)<tol14))exit 203 end if 204 end do 205 206 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 207 if(abs(z_path)**twotrotter>tiny(one)) then 208 funczero = func1_rec(z_path**twotrotter) 209 else 210 funczero = czero 211 end if 212 delta_calc = funczero*N/D*dz_path 213 if(ii==0.or.ii==n_pt_integ_path2)then 214 ent_acc = ent_acc + half*delta_calc 215 if(debug_rec) ent_acc3 = ent_acc3 + half*delta_calc 216 else 217 ent_acc = ent_acc + funczero*delta_calc 218 if(debug_rec) ent_acc3 = ent_acc3 + funczero*delta_calc 219 end if 220 end if 221 end do path5 222 223 else ! trotter==0 224 225 n_pt_integ_path3 = max(100,int((epsilo*half*pi)/real(step,dp))+1) 226 dr_step = one/real(n_pt_integ_path3,dp) 227 path6: do ii = 0,n_pt_integ_path3 228 iif=half*pi*real(ii,dp)*dr_step 229 z_path = epsilo*cmplx(-cos(iif),1-sin(iif),dp) 230 dz_path = epsilo*cmplx(sin(iif),-cos(iif),dp)*half*pi*dr_step 231 if(abs(z_path)**twotrotter>tol14)then 232 Nold = czero 233 Dold = cone 234 N = cone 235 D = z_path - cmplx(an(0),zero,dp) 236 do kk=1,nrec 237 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 238 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 239 Nold = N 240 Dold = D 241 N = Nnew 242 D = Dnew 243 if(kk/=nrec .and. bn2(kk+1)<tol14) exit !-EXIT 244 end do 245 246 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 247 delta_calc = func1_rec(z_path**twotrotter) * N/D * dz_path 248 if(ii==0.or.ii==n_pt_integ_path3)then 249 ent_acc = ent_acc + half*delta_calc 250 if(debug_rec) ent_acc3 = ent_acc3 + half*delta_calc 251 else 252 ent_acc = ent_acc + delta_calc !<r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 253 if(debug_rec) ent_acc3 = ent_acc3 + delta_calc !<r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 254 end if 255 end if 256 end do path6 257 258 end if 259 260 if(first_en==1 .and. debug_rec) then 261 write(msg,'(a,i5,2a,i5,2a,i5,2a,es11.4,2a,es11.4,2a,es11.4)')& 262 & 'n_pt_path =',n_pt_integ,ch10,& 263 & 'n_pt_path2 =',n_pt_integ_path2,ch10,& 264 & 'n_pt_path3 =',n_pt_integ_path3,ch10,& 265 & 'xmin =',xmin,ch10,& 266 & 'xmax =',xmax,ch10,& 267 & 'epsilon =',epsilo 268 call wrtout(std_out,msg,'COLL') 269 first_en = 0 270 end if 271 272 !#################################################################### 273 ![xmax,xmax+i*epsilo] 274 dr_step = one/real(n_pt_integ_path2,dp) 275 dz_path = cmplx(zero,epsilo*dr_step,dp) 276 path4: do ii = 0,n_pt_integ_path2 277 z_path = cmplx(xmax,real(ii,dp)*epsilo*dr_step,dp) 278 279 Nold = czero 280 Dold = cone 281 N = cone 282 D = z_path - cmplx(an(0),zero,dp) 283 284 do kk=1,nrec 285 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 286 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 287 288 Nold = N 289 Dold = D 290 N = Nnew 291 D = Dnew 292 293 if(kk/=nrec)then 294 if((bn2(kk+1)<tol14))exit 295 end if 296 end do 297 298 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 299 delta_calc = func1_rec(z_path**twotrotter)*N/D*dz_path 300 if(ii==0.or.ii==n_pt_integ_path2)then 301 302 ent_acc = ent_acc + half*delta_calc 303 if(debug_rec) ent_acc2 = ent_acc2 + half*delta_calc 304 else 305 ent_acc = ent_acc + delta_calc 306 if(debug_rec) ent_acc2 = ent_acc2 + delta_calc 307 end if 308 end do path4 309 310 311 ent_out = multce*real(ent_acc*cmplx(zero,-piinv,dp),dp) 312 if(debug_rec) then 313 ent_out1 = multce*real(ent_acc1*cmplx(zero,-piinv,dp),dp) 314 ent_out2 = multce*real(ent_acc2*cmplx(zero,-piinv,dp),dp) 315 ent_out3 = multce*real(ent_acc3*cmplx(zero,-piinv,dp),dp) 316 ent_out4 = multce*real(ent_acc4*cmplx(zero,-piinv,dp),dp) 317 end if 318 319 call timab(610,2,tsec) 320 321 contains 322 323 !function to integrate over the path 324 !func1_rec(z_path,twotrotter) = ( z_path**twotrotter/(1+z_path**twotrotter)*log(1+1/z_path**twotrotter)+& !- f*ln(f) 325 !&1/(1+z_path**twotrotter)*log(1+z_path**twotrotter)) !- (1-f)*ln(1-f) 326 327 !func1_rec(z_path_pow) = z_path_pow/(cone+z_path_pow)*log(cone+cone/z_path_pow)+& !- f*ln(f) 328 !&cone/(cone+z_path_pow)*log(cone+z_path_pow) !- (1-f)*ln(1-f) 329 330 !other expression of func for a path like ro(t)*exp(2*i*pi/(2*p)*(j+1/2)) 331 332 function func1_rec(z) 333 334 335 !This section has been created automatically by the script Abilint (TD). 336 !Do not modify the following lines by hand. 337 #undef ABI_FUNC 338 #define ABI_FUNC 'func1_rec' 339 !End of the abilint section 340 341 implicit none 342 343 complex(dpc) :: func1_rec 344 complex(dpc),intent(in) :: z 345 346 func1_rec = z/(cone+z)*log(cone+cone/z)+ cone/(cone+z)*log(cone+z) 347 348 end function func1_rec 349 350 end subroutine entropyrec