TABLE OF CONTENTS
ABINIT/xcpositron [ Functions ]
NAME
xcpositron
FUNCTION
Compute electron-positron correlation potentials and energy density. Used electron-positron correlation functional is controlled by ixcpositron argument. Returns Fxc, Vxc_pos, Vxc_el from input rhor_pos and rhor_el for positron and electrons.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (GJ,MT) 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
grhoe2(ngr)=square of the gradient of electronic density rhoe (needed for GGA) ixcpositron=type of electron-positron correlation functional: 1: LDA zero positron density limit parametrized by Arponen & Pajanne and provided by Boronski & Nieminen [1,2] 11: LDA zero positron density limit parametrized by Arponen & Pajanne and fitted by Sterne & Kaiser [1,3] 2: LDA electron-positron correlation provided by Puska, Seitsonen, and Nieminen [1,4] 3: GGA zero positron density limit parametrized by Arponen & Pajanne and provided by Boronski & Nieminen [1,2,5] 31: GGA zero positron density limit parametrized by Arponen & Pajanne and fitted by Sterne & Kaiser [1,3,5] See references below ngr=size of grho2 array (0 if LDA, npt if GGA) npt=number of real space points on which density is provided posdensity0_limit=True if we are in the zero positron density limit rhoer(npt)=electron density (bohr^-3) rhopr(npt)=positron density (bohr^-3)
OUTPUT
fnxc(npt)=correlation energy per unit volume fxc vxce(npt)=correlation potential for electron dfxc/drhoe (hartree) vxcp(npt)=correlation potential for positron dfxc/drhop (hartree) vxcegr(ngr)= 1/|gradRhoe| dfxc/d|gradRhoe| (empty if LDA, i.e. ngr=0) Optional outputs: dvxce(npt)=partial second derivatives of the xc energy wr to the electronic density dvxce(:)=dVxce/dRhoe dvxcp(npt)=partial second derivatives of the xc energy wr to the positronic density dvxcp(:)=dVxcp/drhop
NOTES
References for electron-positron correlation functionals: [1] J. Arponen and E. Pajanne, Ann. Phys. (N.Y.) 121, 343 (1979). [2] Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986). [3] P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991). [4] M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994). [5] B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1994)
PARENTS
m_pawxc,rhohxcpositron,rhotoxc
CHILDREN
invcb
SOURCE
65 #if defined HAVE_CONFIG_H 66 #include "config.h" 67 #endif 68 69 #include "abi_common.h" 70 71 subroutine xcpositron(fnxc,grhoe2,ixcpositron,ngr,npt,posdensity0_limit,rhoer,rhopr,vxce,vxcegr,vxcp,& 72 & dvxce,dvxcp) ! optional arguments 73 74 use defs_basis 75 use m_errors 76 use m_profiling_abi 77 78 !This section has been created automatically by the script Abilint (TD). 79 !Do not modify the following lines by hand. 80 #undef ABI_FUNC 81 #define ABI_FUNC 'xcpositron' 82 use interfaces_41_xc_lowlevel, except_this_one => xcpositron 83 !End of the abilint section 84 85 implicit none 86 87 !Arguments ------------------------------------ 88 !scalars 89 integer,intent(in) :: ixcpositron,ngr,npt 90 logical,intent(in) :: posdensity0_limit 91 !arrays 92 real(dp),intent(in) :: grhoe2(ngr),rhoer(npt),rhopr(npt) 93 real(dp),intent(out) :: fnxc(npt),vxce(npt),vxcegr(ngr),vxcp(npt) 94 real(dp),intent(out),optional :: dvxce(npt),dvxcp(npt) 95 96 !Local variables------------------------------- 97 !scalars 98 integer,parameter :: idebug=0 99 integer :: ipt 100 logical :: gga,need_dvxce,need_dvxcp 101 real(dp),parameter :: alpha_gga=0.22_dp 102 real(dp),parameter :: ap_a1=-1.56_dp,ap_b1=0.051_dp,ap_c1=-0.081_dp,ap_d1=1.14_dp 103 real(dp),parameter :: ap_a2=-0.92305_dp,ap_b2=-0.05459_dp 104 real(dp),parameter :: ap_a3=-0.6298_dp,ap_b3=-13.15111_dp,ap_c3=2.8655_dp 105 real(dp),parameter :: ap_a4=-179856.2768_dp,ap_b4=186.4207_dp,ap_c4=-0.524_dp 106 real(dp),parameter :: ap_psn_limit=0.7_dp 107 real(dp),parameter :: ap_psn_1=0.9_dp*ap_psn_limit,ap_psn_2=1.1_dp*ap_psn_limit 108 real(dp),parameter :: fpi3=third*four_pi 109 real(dp),parameter :: psn_aa=69.7029_dp,psn_ba=-107.4927_dp,psn_bb=141.8458_dp 110 real(dp),parameter :: psn_ca=23.7182_dp,psn_cb=-33.6472_dp ,psn_cc=5.21152_dp 111 real(dp),parameter :: sk_a=-1.56_dp,sk_b=0.1324_dp,sk_c=-4.092_dp,sk_d=51.96_dp,sk_e=0.7207_dp 112 real(dp),parameter :: rsfac=0.6203504908994000_dp 113 real(dp) :: arse,brse,crse,darse,dbrse,dcrse,d2arse,d2brse,d2crse 114 real(dp) :: d2eps,deps,dexc,dexcdg,dexc_p,d2expgga,dexpgga,d2invrs,dinvrs,d2kf,dkf,d2nqtf2,dnqtf2 115 real(dp) :: drse,drsp,d2exc,d2exc_p,d2rse,d2rsp,d2sqr,dsqr 116 real(dp) :: eexp,eps,exc,exc_p,expgga,invf,dinvf,d2invf,invrhoe,invrhop,invrs,invsqr 117 real(dp) :: kf,logrs,nqtf2,opr2,ratio_ap,ratio_psn,rhoe,rhop,rse,rsp,sqr 118 character(len=500) :: msg 119 !arrays 120 real(dp),allocatable :: rsepts(:),rsppts(:) 121 122 ! ************************************************************************* 123 124 gga=(ngr==npt) 125 need_dvxce=present(dvxce) 126 need_dvxcp=present(dvxcp) 127 128 if (gga.and.ixcpositron==2) then 129 msg = 'xcpositron: GGA not yet implemented for ixcpositron=2 !' 130 MSG_ERROR(msg) 131 end if 132 if (posdensity0_limit.and.ixcpositron==2) then 133 msg = 'xcpositron: ixcpositron=2 cannot be treated in the zero positron density limit !' 134 MSG_ERROR(msg) 135 end if 136 if (abs(ixcpositron)/=1.and.ixcpositron/=11.and.ixcpositron/=2.and.ixcpositron/=3.and.ixcpositron/=31) then 137 msg = 'xcpositron: unknown electron-positron correlation !' 138 MSG_ERROR(msg) 139 end if 140 141 !Compute density radii for rhor_el, rhor_pos 142 ABI_ALLOCATE(rsepts,(npt)) 143 call invcb(rhoer(:),rsepts,npt) 144 rsepts(:)=rsfac*rsepts(:) 145 if (ixcpositron==2) then 146 ABI_ALLOCATE(rsppts,(npt)) 147 call invcb(rhopr(:),rsppts,npt) 148 rsppts(:)=rsfac*rsppts(:) 149 end if 150 151 !Loop over grid points 152 !---------------------------------------------------- 153 do ipt=1,npt 154 155 rhoe=rhoer(ipt) 156 rhop=rhopr(ipt) 157 exc=zero;dexc=zero;d2exc=zero;dexcdg=zero 158 159 rse=rsepts(ipt) 160 invrhoe=one/rhoe 161 drse=-third*rse*invrhoe 162 if (need_dvxce) d2rse= four/nine*rse*invrhoe**2 163 164 ! Arponen & Pajane parametrization for electron 165 if (ixcpositron/=11.and.ixcpositron/=31) then 166 if (rse<0.302_dp) then 167 invrs=one/rse;invsqr=sqrt(invrs);logrs=log(rse) 168 exc =ap_a1*invsqr+(ap_b1*logrs+ap_c1)*logrs+ap_d1 169 dexc=drse*invrs*(-half*ap_a1*invsqr+two*ap_b1*logrs+ap_c1) 170 if (need_dvxce) d2exc=(d2rse/drse-drse*invrs)*dexc+drse**2*invrs**2*(quarter*ap_a1*invsqr+two*ap_b1) 171 else if (rse>=0.302_dp.and.rse<=0.56_dp) then 172 invrs=one/rse 173 exc =ap_a2+ap_b2*invrs**2 174 dexc=-drse*ap_b2*two*invrs**3 175 if (need_dvxce) d2exc=d2rse/drse*dexc+six*drse**2*ap_b2*invrs**4 176 else if (rse>0.56_dp.and.rse<=8.0_dp) then 177 invrs=one/(rse+2.5_dp) 178 dinvrs=-drse*invrs**2 179 ! jmb : d2rse initialized only if need_dvxce = .True. 180 if (need_dvxce) d2invrs=-d2rse*invrs**2-two*invrs*drse**2 181 exc =ap_a3+ap_b3*invrs**2+ap_c3*invrs 182 dexc=two*ap_b3*invrs*dinvrs+ap_c3*dinvrs 183 if (need_dvxce) d2exc=two*ap_b3*dinvrs**2+(two*ap_b3*invrs+ap_c3)*d2invrs 184 else 185 exc =ap_a4*rhoe**2+ap_b4*rhoe+ap_c4 186 dexc =two*ap_a4*rhoe+ap_b4 187 if (need_dvxce) d2exc=two*ap_a4 188 end if 189 190 ! Sterne & Kaiser parametrization for electron 191 else 192 eexp=exp(-(rse+sk_c)**2/sk_d) 193 opr2=(one+rse**2) 194 arse=atan(rse) 195 exc = sk_a/sqrt(arse)+sk_b*eexp+sk_e 196 dexc= -(two*sk_b*eexp*(sk_c+rse)/sk_d + sk_a/(two*opr2*sqrt(arse)**3))*drse 197 if (need_dvxce) d2exc=-(two*sk_b*eexp*(sk_c+rse)/sk_d + sk_a/(two*opr2*arse**1.5_dp))*d2rse & 198 & +(two*sk_b*eexp*(two*sk_c**2-sk_d+four*sk_c*rse+two*rse**2)/sk_d**2 & 199 & +sk_a*(three+four*rse*arse)/(four*opr2**2*sqrt(arse)**5))*drse**2 200 end if 201 202 ! Puska, Seitsonen and Nieminen parametrization for positron 203 if (ixcpositron==2.and.rse>=ap_psn_1) then 204 rsp=rsppts(ipt) 205 invrhop=one/rhop 206 drsp=-third*rsp*invrhop 207 if (need_dvxcp) d2rsp= four/nine*rsp*invrhop**2 208 exc_p=zero;dexc_p=zero;d2exc_p=zero 209 if (rsp<0.302_dp) then 210 invrs=one/rsp;invsqr=sqrt(invrs);logrs=log(rsp) 211 exc_p =ap_a1*invsqr+(ap_b1*logrs+ap_c1)*logrs+ap_d1 212 dexc_p=drsp*invrs*(-half*ap_a1*invsqr+two*ap_b1*logrs+ap_c1) 213 if (need_dvxcp) d2exc_p=(d2rsp/drsp-drsp*invrs)*dexc_p+drsp**2*invrs**2*(quarter*ap_a1*invsqr+two*ap_b1) 214 else if (rsp>=0.302_dp.and.rsp<=0.56_dp) then 215 invrs=one/rsp 216 exc_p =ap_a2+ap_b2*invrs**2 217 dexc_p=-drsp*ap_b2*two*invrs**3 218 if (need_dvxcp) d2exc_p=d2rsp/drsp*dexc_p+six*drsp**2*ap_b2*invrs**4 219 else if (rsp>0.56_dp.and.rsp<=8.0_dp) then 220 invrs=one/(rsp+2.5_dp) 221 dinvrs=-drsp*invrs**2 222 ! jmb : d2rsp initialized only if need_dvxcp = .True.* 223 if (need_dvxcp) d2invrs=-d2rsp*invrs**2-two*invrs*drsp**2 224 exc_p =ap_a3+ap_b3*invrs**2+ap_c3*invrs 225 dexc_p=two*ap_b3*invrs*dinvrs+ap_c3*dinvrs 226 if (need_dvxcp) d2exc_p=two*ap_b3*dinvrs**2+(two*ap_b3*invrs+ap_c3)*d2invrs 227 else 228 exc_p =ap_a4*rhop**2+ap_b4*rhop+ap_c4 229 dexc_p =two*ap_a4*rhop+ap_b4 230 if (need_dvxcp) d2exc_p=two*ap_a4 231 end if 232 end if 233 234 ! GGA correction 235 if (gga) then 236 kf=(three*pi*pi*rhoe)**third 237 nqtf2=(rhoe*sqrt(four*kf/pi))**2 238 eps=grhoe2(ipt)/nqtf2 239 if (eps<zero) then 240 MSG_ERROR('xcpositron: problem, negative GGA espilon !') 241 end if 242 expgga=exp(-alpha_gga*eps*third) 243 244 dkf=pi*pi/(sqrt(three*pi*pi*rhoe)**third) 245 d2kf=-two*pi*pi*pi*pi*(three*pi*pi*rhoe)**(-5.0_dp/3.0_dp) 246 sqr=sqrt(four*kf/pi) 247 dsqr=(four*dkf/pi)/(two*sqr) 248 d2sqr=two/(pi*sqr*dkf)*(d2kf*sqr-dsqr*dkf) 249 nqtf2=(rhoe*sqr)**two 250 dnqtf2=two*(sqr+rhoe*dsqr)*rhoe*sqr 251 d2nqtf2=two*(rhoe*sqr*(two*dsqr+rhoe*d2sqr) & 252 & +sqr*(sqr+rhoe*dsqr) & 253 & +rhoe*(sqr+rhoe*dsqr) ) 254 deps=-grhoe2(ipt)*dnqtf2/(nqtf2**two) 255 d2eps=-grhoe2(ipt)/(nqtf2*nqtf2*dnqtf2)*(d2nqtf2*nqtf2*nqtf2-two*nqtf2*dnqtf2*dnqtf2) 256 dexpgga=-alpha_gga*third*deps*expgga 257 d2expgga=-alpha_gga*third*(d2eps*expgga+deps*dexpgga) 258 259 exc = exc *expgga 260 dexc=(dexc*expgga+exc*dexpgga) 261 if (need_dvxce) d2exc=d2exc*expgga+two*dexc*dexpgga+exc*d2expgga 262 if (abs(grhoe2(ipt))<1.e24_dp) dexcdg=-exc*alpha_gga*two_thirds/nqtf2 263 end if 264 265 ! Computation of XC energy, potentials and kernels 266 ! Zero positron density limit 267 if (ixcpositron/=2.or.rse<ap_psn_1) then 268 fnxc(ipt)=rhop*exc 269 vxce(ipt)=rhop*dexc 270 vxcp(ipt)=exc 271 if (need_dvxce) dvxce(ipt)=rhop*d2exc 272 if (need_dvxcp) dvxcp(ipt)=zero 273 if (gga) vxcegr(ipt)=rhop*dexcdg 274 else 275 ! Puska, Seitsonen and Nieminen functional 276 arse=psn_aa+psn_ba*rse+psn_ca*rse**2 277 brse=psn_ba+psn_bb*rse+psn_cb*rse**2 278 crse=psn_ca+psn_cb*rse+psn_cc*rse**2 279 darse=(psn_ba+two*psn_ca*rse)*drse 280 dbrse=(psn_bb+two*psn_cb*rse)*drse 281 dcrse=(psn_cb+two*psn_cc*rse)*drse 282 invf=arse+brse*rsp+crse*rsp**2+invrhop/exc+invrhoe/exc_p 283 fnxc(ipt)=one/invf 284 dinvf=darse+dbrse*rsp+dcrse*rsp**2-invrhop*dexc/exc**2-invrhoe**2/exc_p 285 vxce(ipt)=-dinvf/invf**2 286 if (need_dvxce) then 287 d2arse=darse*d2rse/drse+two*psn_ca*drse**2 288 d2brse=dbrse*d2rse/drse+two*psn_cb*drse**2 289 d2crse=dcrse*d2rse/drse+two*psn_cc*drse**2 290 d2invf=d2arse+d2brse*rsp+d2crse*rsp**2 & 291 & +invrhop*(two*dexc**2/exc-d2exc)/exc**2+two*invrhoe**3/exc_p 292 dvxce(ipt)=(two*dinvf**2/invf-d2invf)/invf**2 293 end if 294 dinvf=(brse+two*crse*rsp)*drsp-invrhop**2/exc-invrhoe*dexc_p/exc_p**2 295 vxcp(ipt)=-dinvf/invf**2 296 if (need_dvxcp) then 297 d2invf=two*crse*drsp+(brse+two*crse*rsp)*d2rsp & 298 & +two*invrhop**3/exc+invrhoe*(two*dexc_p**2/exc_p-d2exc_p)/exc_p**2 299 dvxcp(ipt)=(two*dinvf**2/invf-d2invf)/invf**2 300 end if 301 ! For small rse, use pure Arponen/Pajanne functional 302 ! Around the limit (rse=0.7, see PSN paper), switch smoothly from PSN to AP 303 if (rse>=ap_psn_1.and.rse<=ap_psn_2) then 304 ratio_psn=(rse-ap_psn_1)/(ap_psn_2-ap_psn_1);ratio_ap=one-ratio_psn 305 fnxc(ipt)=ratio_psn*fnxc(ipt)+ratio_ap*rhop*exc 306 vxce(ipt)=ratio_psn*vxce(ipt)+ratio_ap*rhop*dexc 307 vxcp(ipt)=ratio_psn*vxcp(ipt)+ratio_ap*exc 308 if (need_dvxce) dvxce(ipt)=ratio_psn*dvxce(ipt)+ratio_ap*rhop*d2exc 309 if (need_dvxcp) dvxcp(ipt)=ratio_psn*dvxcp(ipt) 310 end if 311 end if 312 313 ! Debug statements: use polynomial functionals 314 if (idebug>0) then 315 if (idebug==4) then ! order 4 316 fnxc(ipt)=tol3*((rhop**4+rhoe**4)/12._dp+(rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe**2) 317 vxce(ipt)=tol3*((rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe+rhop*rhoe**2) 318 vxcp(ipt)=tol3*((rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe+rhop*rhoe**2) 319 if (need_dvxce) dvxce(ipt)=tol3*(rhop**3/3._dp+rhop*rhoe**2+rhop**2+two*rhop*rhoe) 320 if (need_dvxcp) dvxcp(ipt)=tol3*(rhoe**3/3._dp+rhoe*rhop**2+rhoe**2+two*rhop*rhoe) 321 end if 322 if (idebug==3) then ! order 3 323 fnxc(ipt)=tol3*((rhop**3+rhoe**3)*third+rhop**2*rhoe+rhop*rhoe**2) 324 vxce(ipt)=tol3*(rhop+rhoe)**2 325 vxcp(ipt)=tol3*(rhop+rhoe)**2 326 if (need_dvxce) dvxce(ipt)=tol3*two*rhoe 327 if (need_dvxcp) dvxcp(ipt)=tol3*two*rhop 328 end if 329 if (idebug==2) then ! order 2 330 fnxc(ipt)=tol3*(rhop+rhoe)**2 331 vxce(ipt)=tol3*two*(rhop+rhoe) 332 vxcp(ipt)=tol3*two*(rhop+rhoe) 333 if (need_dvxce) dvxce(ipt)=tol3*two 334 if (need_dvxcp) dvxcp(ipt)=tol3*two 335 end if 336 if (idebug==1) then ! order 1 337 fnxc(ipt)=tol3*(rhop+rhoe) 338 vxce(ipt)=tol3 339 vxcp(ipt)=tol3 340 if (need_dvxce) dvxce(ipt)=zero 341 if (need_dvxcp) dvxcp(ipt)=zero 342 end if 343 end if 344 345 end do ! ipt 346 347 ABI_DEALLOCATE(rsepts) 348 if (ixcpositron==2) then 349 ABI_DEALLOCATE(rsppts) 350 end if 351 352 !Convert everything in Hartree units 353 fnxc(:)=half*fnxc(:) 354 vxce(:)=half*vxce(:) 355 vxcp(:)=half*vxcp(:) 356 if (need_dvxce) dvxce(:)=half*dvxce(:) 357 if (need_dvxcp) dvxcp(:)=half*dvxcp(:) 358 if (gga) vxcegr(:)=half*vxcegr(:) 359 360 end subroutine xcpositron