TABLE OF CONTENTS
ABINIT/gran_potrec [ Functions ]
NAME
gran_potrec
FUNCTION
This routine computes the local part of the grand-potential 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 mult=a multiplicator for computing grand-potential (2 for non-spin-polarized system) debug_rec=debugging variable n_pt_integ=points for computation of path integral xmax= maximum point on the x-axis for integration
OUTPUT
ene_out=grand-potential at the point if debug_rec=T then ene_out1,ene_out2,ene_out3,ene_out4 are the different path branch contriubutions to the grand-potential. In reality it is not the gren potential but the grand-potential (omega=-PV) divided by -T
PARENTS
vtorhorec
CHILDREN
timab,wrtout
NOTES
in reality it is not the gren potential but the grand-potential (omega=-PV) divided by -T
at this time :
- mult should be not used
- the routine should be integraly rewrited and use the routine recursion.
- only modified for p /= 0
SOURCE
49 #if defined HAVE_CONFIG_H 50 #include "config.h" 51 #endif 52 53 #include "abi_common.h" 54 55 56 subroutine gran_potrec(an,bn2,nrec,trotter,ene_out, mult, & 57 & debug_rec,n_pt_integ,xmax,& 58 & ene_out1,ene_out2,ene_out3,ene_out4) 59 60 use defs_basis 61 use m_profiling_abi 62 63 !This section has been created automatically by the script Abilint (TD). 64 !Do not modify the following lines by hand. 65 #undef ABI_FUNC 66 #define ABI_FUNC 'gran_potrec' 67 use interfaces_14_hidewrite 68 use interfaces_18_timing 69 !End of the abilint section 70 71 implicit none 72 73 !Arguments ------------------------------- 74 !scalars 75 integer,intent(in) :: n_pt_integ,nrec,trotter 76 logical,intent(in) :: debug_rec 77 real(dp), intent(in) :: mult,xmax 78 real(dp),intent(inout) :: ene_out,ene_out1,ene_out2,ene_out3,ene_out4 !vz_i 79 !arrays 80 real(dp), intent(in) :: an(0:nrec),bn2(0:nrec) 81 82 !Local variables------------------------------- 83 !scalars 84 integer, parameter :: level = 7 85 integer, save :: first = 1 86 integer :: ii,kk,n_pt_integ_path2 87 real(dp) :: epsilon,step,twotrotter,xmin,dr_step 88 complex(dpc) :: D,Dnew,Dold,N,Nnew,Nold,dz_path,ene_acc,ene_acc1,ene_acc2 89 complex(dpc) :: ene_acc3,ene_acc4 90 complex(dpc) :: z_path,delta_calc 91 character(len=500) :: message 92 !arrays 93 real(dp) :: tsec(2) 94 95 ! ************************************************************************* 96 97 98 call timab(611,1,tsec) 99 100 !structured debugging if debug_rec=T : print detailled result the first time we enter gran_potrec 101 if(debug_rec .and. first==1)then 102 write(message,'(a)')' ' 103 call wrtout(std_out,message,'PERS') 104 write(message,'(a)')' gran_potrec : enter ' 105 call wrtout(std_out,message,'PERS') 106 write(message,'(a,i8)')'n_pt_integ ' , n_pt_integ 107 call wrtout(std_out,message,'COLL') 108 first=0 109 end if 110 111 ene_out = zero 112 ene_acc = czero 113 ene_acc1 = czero 114 ene_acc2 = czero 115 ene_acc3 = czero 116 ene_acc4 = czero 117 118 119 !path parameters 120 !n_pt_integ = 2500 121 xmin = -half 122 step = (xmax-xmin)/real(n_pt_integ,dp) 123 if(trotter==0)then 124 twotrotter = one 125 epsilon = .5d-1 126 else 127 twotrotter = two*real(trotter,dp) 128 epsilon = half*sin( pi/twotrotter) 129 end if 130 131 !xmin = -abs(xmin)**(1.d0/twotrotter) 132 133 !#################################################################### 134 ![xmax + i*epsilon,xmin + i*epsilon] 135 dr_step = one/real(n_pt_integ,dp) 136 dz_path = -cmplx((xmax-xmin)*dr_step,zero,dp) 137 path1: do ii = 0,n_pt_integ 138 ! z_path = cmplx(xmin + real(ii,dp)*(xmax-xmin)*dr_step,epsilon,dp) 139 z_path = cmplx(xmin,epsilon,dp) - real(ii,dp)*dz_path 140 Nold = czero 141 Dold = cone 142 N = cone 143 D = z_path - cmplx(an(0),zero,dp) 144 145 do kk=1,nrec 146 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 147 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 148 149 Nold = N 150 Dold = D 151 N = Nnew 152 D = Dnew 153 154 if(kk/=nrec)then 155 if((bn2(kk+1)<tol14))exit 156 end if 157 158 end do 159 160 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 161 delta_calc = func_rec(z_path,twotrotter)* N/D *dz_path 162 if(ii==0.or.ii==n_pt_integ)then 163 ene_acc = ene_acc + half*delta_calc 164 if(debug_rec) ene_acc1 = ene_acc1 + half*delta_calc 165 else 166 ene_acc = ene_acc + delta_calc 167 if(debug_rec) ene_acc1 = ene_acc1 + delta_calc 168 end if 169 end do path1 170 171 !#################################################################### 172 ![xmin + i*epsilon,xmin] 173 if(epsilon/step>4.d0)then 174 n_pt_integ_path2 = int(epsilon/step)+1 175 else 176 n_pt_integ_path2 = 5 177 end if 178 n_pt_integ_path2 = n_pt_integ 179 dr_step = one/real(n_pt_integ_path2,dp) 180 dz_path = -cmplx(zero,epsilon*dr_step,dp) 181 path2: do ii = 0,n_pt_integ_path2 182 ! z_path = cmplx(xmin,real(ii,dp)*epsilon*dr_step,dp) 183 z_path = cmplx(xmin,zero,dp)-dz_path*real(ii,dp) 184 Nold = czero 185 Dold = cone 186 N = cone 187 D = z_path - cmplx(an(0),zero,dp) 188 189 do kk=1,nrec 190 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 191 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 192 193 Nold = N 194 Dold = D 195 N = Nnew 196 D = Dnew 197 198 if(kk/=nrec)then 199 if((bn2(kk+1)<tol14))exit 200 end if 201 202 end do 203 204 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 205 delta_calc = func_rec(z_path,twotrotter)* N/D *dz_path 206 if(ii==0.or.ii==n_pt_integ_path2)then 207 ene_acc = ene_acc + half*delta_calc 208 if(debug_rec) ene_acc3 = ene_acc3 + half*delta_calc 209 else 210 ene_acc = ene_acc + delta_calc 211 if(debug_rec) ene_acc3 = ene_acc3 + delta_calc 212 end if 213 end do path2 214 215 216 217 !#################################################################### 218 ![xmin,0] 219 if(xmin/=czero)then 220 dr_step = one/real(n_pt_integ,dp) 221 dz_path = cmplx(xmin*dr_step,zero,dp) 222 path3: do ii = 1,n_pt_integ !the integrand is 0 at 0 223 ! z_path = cmplx(real(ii,dp)*xmin*dr_step,zero,dp) 224 z_path = real(ii,dp)*dz_path 225 226 Nold = czero 227 Dold = cone 228 N = cone 229 D = z_path - cmplx(an(0),zero,dp) 230 231 do kk=1,nrec 232 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 233 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 234 235 Nold = N 236 Dold = D 237 N = Nnew 238 D = Dnew 239 240 if(kk/=nrec)then 241 if((bn2(kk+1)<tol14))exit 242 end if 243 end do 244 245 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 246 delta_calc = func_rec(z_path,twotrotter) * N/D *dz_path 247 if(ii==n_pt_integ)then 248 ene_acc = ene_acc +half*delta_calc 249 if(debug_rec) ene_acc4 = ene_acc4 + half*delta_calc 250 else 251 ene_acc = ene_acc + delta_calc 252 if(debug_rec) ene_acc4 = ene_acc4 +delta_calc 253 end if 254 end do path3 255 end if 256 257 !#################################################################### 258 ![xmax,xmax+i*epsilon] 259 dr_step = one/real(n_pt_integ_path2,dp) 260 dz_path = cmplx(zero,epsilon*dr_step,dp) 261 path4: do ii = 0,n_pt_integ_path2 262 ! z_path = cmplx(xmax,real(ii,dp)*epsilon*dr_step,dp) 263 z_path = cmplx(xmax,0,dp)+real(ii,dp)*dz_path 264 265 Nold = czero 266 Dold = cone 267 N = cone 268 D = z_path - cmplx(an(0),zero,dp) 269 270 do kk=1,nrec 271 Nnew = (z_path - cmplx(an(kk),zero,dp))*N - cmplx(bn2(kk),zero,dp)*Nold 272 Dnew = (z_path - cmplx(an(kk),zero,dp))*D - cmplx(bn2(kk),zero,dp)*Dold 273 274 Nold = N 275 Dold = D 276 N = Nnew 277 D = Dnew 278 279 if(kk/=nrec)then 280 if((bn2(kk+1)<tol14))exit 281 end if 282 283 end do 284 285 ! <r|1/(z-e**(-beta/(2p)*(H-mu)))|r> dz 286 delta_calc = func_rec(z_path,twotrotter) * N/D *dz_path 287 if(ii==0.or.ii==n_pt_integ_path2)then 288 ene_acc = ene_acc + half*delta_calc 289 if(debug_rec) ene_acc2 = ene_acc2 + half*delta_calc 290 else 291 ene_acc = ene_acc + delta_calc 292 if(debug_rec) ene_acc2 = ene_acc2 + delta_calc 293 end if 294 end do path4 295 296 ene_out = mult*real(ene_acc*cmplx(zero,-piinv,dp),dp) 297 if(debug_rec) then 298 ene_out1 = mult*real(ene_acc1*cmplx(zero,-piinv,dp),dp) 299 ene_out2 = mult*real(ene_acc2*cmplx(zero,-piinv,dp),dp) 300 ene_out3 = mult*real(ene_acc3*cmplx(zero,-piinv,dp),dp) 301 ene_out4 = mult*real(ene_acc4*cmplx(zero,-piinv,dp),dp) 302 end if 303 304 call timab(611,2,tsec) 305 306 contains 307 308 !func_rec(z_path,twotrotter) = log(cone+z_path**twotrotter) 309 310 function func_rec(z,x) 311 312 313 !This section has been created automatically by the script Abilint (TD). 314 !Do not modify the following lines by hand. 315 #undef ABI_FUNC 316 #define ABI_FUNC 'func_rec' 317 !End of the abilint section 318 319 implicit none 320 321 complex(dpc) :: func_rec 322 complex(dpc),intent(in) :: z 323 real(dp),intent(in) :: x 324 325 func_rec = log(cone+z**x) 326 327 end function func_rec 328 329 end subroutine gran_potrec