TABLE OF CONTENTS
ABINIT/gammapositron [ Functions ]
NAME
gammapositron
FUNCTION
Compute positron electron-positron enhancement factor (contact density) used to compute positron lifetime. Input is positronic rhop(r) and electronic rhoe(r) at a given set of points.
INPUTS
grhocore2(ngr)=square of the gradient of core electronic density rhocore (needed for GGA) grhoe2(ngr)=square of the gradient of valence electronic density rhoer (needed for GGA) igamma=type of enhancement factor: 1: Boronski and Nieminen [2] 2: Boronski and Nieminen, RPA limit [2] 3: Sterne and Kaiser [3] 4: Puska, Seitsonen and Nieminen [4] 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 rhocore(npt*usecore)=core electron density (bohr^-3) rhoer(npt) =valence electron density (bohr^-3) rhopr(npt) =positron density (bohr^-3) usecore =1 if core density is not zero
OUTPUT
gamma(npt,2) =electron-positron enhancement factor, gamma(:,1): using total electronic density gamma(:,2): using valence electronic density
NOTES
References for electron-positron correlation functionals: [1] J. Arponen and E. Pajanne, Ann. Phys. (N.Y.) 121, 343 (1979) [[cite:Arponen1979a]]. [2] E. Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986) [[cite:Boronski1986]]. [3] P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991) [[cite:Sterne1991]]. [4] M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994) [[cite:Puska1994]]. [5] B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1995) [[cite:Barbiellini1995]]
SOURCE
83 subroutine gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,npt,rhocore,rhoer,rhopr,usecore) 84 85 !Arguments ------------------------------------ 86 !scalars 87 integer,intent(in) :: igamma,ngr,npt,usecore 88 !arrays 89 real(dp),intent(in) :: grhocore2(ngr*usecore),grhoe2(ngr),rhocore(npt*usecore),rhoer(npt),rhopr(npt) 90 real(dp),intent(out) :: gamma(npt,2) 91 92 !Local variables------------------------------- 93 !scalars 94 integer :: iloop,ipt 95 logical :: gga 96 real(dp),parameter :: alpha_gga=0.22d0,rsfac=0.6203504908994000_dp 97 real(dp) :: aa,bb,cc,dg1 98 real(dp) :: drs,eps,expgga,g0,g1,g2,gg 99 real(dp) :: kf,kk,nqtf2,ratio1,ratio2,ratio3,rho1,rho2,rhoe,rhop,sqrs,rs,rse,rsp 100 !arrays 101 real(dp),allocatable :: grho2(:),rhor(:),rsepts(:),rsppts(:) 102 103 ! ************************************************************************* 104 105 gga=(ngr==npt.and.igamma/=0) 106 107 if (usecore/=0.and.usecore/=1) then 108 ABI_ERROR('Wrong value for usecore !') 109 end if 110 if (igamma/=0.and.igamma/=1.and.igamma/=2.and.igamma/=3.and.igamma/=4) then 111 ABI_ERROR('Unknown electron-positron correlation !') 112 end if 113 114 ABI_MALLOC(rhor,(npt)) 115 ABI_MALLOC(rsepts,(npt)) 116 if (gga) then 117 ABI_MALLOC(grho2,(npt)) 118 end if 119 120 !Eventually compute positronic density radii 121 if (igamma==1.or.igamma==4) then 122 ABI_MALLOC(rsppts,(npt)) 123 call invcb(rhopr(:),rsppts,npt) 124 rsppts(:)=rsfac*rsppts(:) 125 end if 126 127 !Loop: iloop=1: compute enhancement factor using total electronic density 128 !iloop=2: compute enhancement factor using valence electronic density 129 !=================================================================================== 130 do iloop=1,1+usecore 131 132 ! Compute electronic density radii 133 if (iloop==1.and.usecore==1) then 134 rhor(1:npt)=rhoer(1:npt)+rhocore(1:npt) 135 else 136 rhor(1:npt)=rhoer(1:npt) 137 end if 138 call invcb(rhor(:),rsepts,npt) 139 rsepts(:)=rsfac*rsepts(:) 140 141 ! Gradients for GGA 142 if (gga) then 143 if (iloop==1.and.usecore==1) then 144 grho2(1:npt)=grhoe2(1:npt)+grhocore2(1:npt) 145 else 146 grho2(1:npt)=grhoe2(1:npt) 147 end if 148 end if 149 150 ! Loop over grid points 151 do ipt=1,npt 152 153 rhoe=rhor(ipt) 154 rhop=rhopr(ipt) 155 rse =rsepts(ipt) 156 gg=zero 157 158 ! Testing feature: gamma=1 159 ! ----------------------------------------------------------------------------------- 160 if (igamma==0) then 161 162 gg=one 163 164 ! Boronski and Nieminen 165 ! ----------------------------------------------------------------------------------- 166 else if (igamma==1) then 167 168 rsp =rsppts(ipt) 169 if (rhoe>rhop) then 170 rho1=rhoe;rho2=rhop;rs=rse 171 else 172 rho1=rhop;rho2=rhoe;rs=rsp 173 end if 174 drs=-third*rs/rho1;sqrs=sqrt(rs) 175 ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1 176 g0=one+1.23_dp*rs+0.8295_dp*sqrs**3-1.26_dp*rs**2+0.3286_dp*sqrs**5+sixth *rs**3 177 g1=one+0.51_dp*rs +0.65_dp*rs**2-0.51_dp *sqrs**5+0.176_dp*rs**3 178 g2=one+0.60_dp*rs +0.63_dp*rs**2-0.48_dp *sqrs**5+0.167_dp*rs**3 179 dg1=drs*(0.51_dp+two*0.65_dp*rs-2.5_dp*0.51_dp*sqrs**3+three*0.176_dp*rs**2) 180 kk=half*rho1*dg1 181 aa= two *kk-six *g1+eight *g2-two *g0 182 bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0 183 cc= kk-four *g1+eight *g2-four*g0 184 gg=g0+ratio3*aa+ratio2*bb+ratio1*cc 185 186 ! Boronski and Nieminen RPA limit 187 ! ----------------------------------------------------------------------------------- 188 else if (igamma==2) then 189 190 rs=rse;sqrs=sqrt(rs) 191 gg=one !This is experimental to avoid divergences 192 if (rs<=20._dp) gg=gg+1.23_dp*rs+0.8295_dp*sqrs**3-1.26_dp*rs**2+0.3286_dp*sqrs**5+sixth*rs**3 193 194 ! Sterne and Kaiser 195 ! ----------------------------------------------------------------------------------- 196 else if (igamma==3) then 197 198 rs=rse;sqrs=sqrt(rs) 199 gg=one !This is experimental to avoid divergences 200 if (rs<=20._dp) gg=gg+0.1512_dp*rs+2.414_dp*sqrs**3-2.01_dp*rs**2+0.4466_dp*sqrs**5+0.1667_dp*rs**3 201 202 ! Puska, Seitsonen and Nieminen 203 ! ----------------------------------------------------------------------------------- 204 else if (igamma==4) then 205 206 rsp =rsppts(ipt) 207 if (rhoe>rhop) then 208 rho1=rhoe;rho2=rhop;rs=rse 209 else 210 rho1=rhop;rho2=rhoe;rs=rsp 211 end if 212 drs=-third*rs/rho1;sqrs=sqrt(rs) 213 ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1 214 g0=one+1.2300_dp*rs+0.9889_dp*sqrs**3-1.4820_dp*rs**2+0.3956_dp*sqrs**5+sixth*rs**3 215 g1=one+2.0286_dp*rs-3.3892_dp*sqrs**3+3.0547_dp*rs**2-1.0540_dp*sqrs**5+sixth*rs**3 216 g2=one+0.2499_dp*rs+0.2949_dp*sqrs**3+0.6944_dp*rs**2-0.5339_dp*sqrs**5+sixth*rs**3 217 dg1=drs*(2.0286_dp-1.5_dp*3.3892_dp*sqrs+two*3.0547_dp*rs-2.5_dp*1.0540_dp*sqrs**3+three*sixth*rs**2) 218 kk=half*rho1*dg1 219 aa= two *kk-six *g1+eight *g2-two *g0 220 bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0 221 cc= kk-four *g1+eight *g2-four*g0 222 gg=g0+ratio3*aa+ratio2*bb+ratio1*cc 223 224 end if ! igamma 225 226 if (gga) then 227 kf=(three*pi*pi*rhoe)**third 228 nqtf2=(rhoe*sqrt(four*kf/pi))**2 229 eps=grho2(ipt)/nqtf2 230 if (eps<zero) then 231 ABI_ERROR(' problem, negative GGA espilon !') 232 end if 233 expgga=exp(-alpha_gga*eps*third) 234 gg=one+(gg-one)*expgga 235 end if 236 237 ! Store enhancement factor 238 gamma(ipt,iloop)=gg 239 240 end do ! ipt 241 end do ! iloop 242 243 ABI_FREE(rhor) 244 ABI_FREE(rsepts) 245 if (igamma==1.or.igamma==4) then 246 ABI_FREE(rsppts) 247 end if 248 if (gga) then 249 ABI_FREE(grho2) 250 end if 251 252 !Case usecore=0 (no core density) 253 if (usecore==0) gamma(:,2)=gamma(:,1) 254 255 end subroutine gammapositron
ABINIT/gammapositron_fft [ Functions ]
NAME
gammapositron_fft
FUNCTION
Compute positron electron-positron enhancement factor on a real space FFT grid.
INPUTS
electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1) igamma=type of enhancement factor: -1: gamma=one (test) 1: Boronski and Nieminen [1] 2: Boronski and Nieminen, RPA limit [1] 3: Sterne and Kaiser [2] 4: Puska, Seitsonen and Nieminen [3] mpi_enreg=information about MPI parallelization n3xccc=dimension of the xccc3d array (0 or nfft). nfft=number of FFT grid points ngfft(18)=contain all needed information about 3D FFT rhor_e(nfft)=real space density electronic density (total density) rhor_p(nfft)=real space density positronic density xccc3d(n3xccc)=3D core electron density for XC core correction
OUTPUT
gamma(nfft,2)=electron-positron enhancement factor, gamma(:,1): using total electronic density gamma(:,2): using valence electronic density
NOTES
The input densities (rhor_e and rhor_p) should be positive everywhere (call mkdenpos routine before entering this one)
SOURCE
293 subroutine gammapositron_fft(electronpositron,gamma,gprimd,igamma,mpi_enreg,& 294 & n3xccc,nfft,ngfft,rhor_e,rhor_p,xccc3d) 295 296 !Arguments ------------------------------------ 297 !scalars 298 integer,intent(in) :: igamma,n3xccc,nfft 299 type(electronpositron_type),pointer :: electronpositron 300 type(MPI_type),intent(in) :: mpi_enreg 301 !arrays 302 integer,intent(in) :: ngfft(18) 303 real(dp),intent(in) :: gprimd(3,3),rhor_e(nfft),rhor_p(nfft),xccc3d(n3xccc) 304 real(dp),intent(out) :: gamma(nfft,2) 305 306 !Local variables------------------------------- 307 !scalars 308 integer :: cplex,ishift,ngr,ngrad,nspden_ep,usecore 309 !arrays 310 real(dp),parameter :: qphon(3)=(/zero,zero,zero/) 311 real(dp),allocatable :: grhocore2(:),grhoe2(:),rhoc(:,:,:),rhoe(:,:,:) 312 313 ! ************************************************************************* 314 315 !Several useful constants 316 usecore=n3xccc/nfft 317 cplex=1;ishift=0;ngrad=1;nspden_ep=1 318 if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2 319 ngr=0;if (ngrad==2) ngr=nfft 320 321 !Allocate several arrays 322 ABI_MALLOC(rhoe,(nfft,nspden_ep,ngrad**2)) 323 ABI_MALLOC(grhoe2,(ngr)) 324 ABI_MALLOC(grhocore2,(ngr*usecore)) 325 326 !Store electronic density and its gradients 327 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,qphon,rhor_e,rhoe) 328 329 !Compute squared gradient of the electronic density 330 if (ngrad==2) then 331 grhoe2(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2 332 if (usecore>0) then 333 ABI_MALLOC(rhoc,(nfft,1,ngrad**2)) 334 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,qphon,xccc3d,rhoc) 335 grhocore2(:)=rhoc(:,1,2)**2+rhoc(:,1,3)**2+rhoc(:,1,4)**2 336 ABI_FREE(rhoc) 337 end if 338 end if 339 340 !Compute enhancement factor on FFT grid 341 call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,nfft,xccc3d,& 342 & rhoe(:,1,1),rhor_p,usecore) 343 344 !Release temporary memory 345 ABI_FREE(rhoe) 346 ABI_FREE(grhoe2) 347 ABI_FREE(grhocore2) 348 349 end subroutine gammapositron_fft
ABINIT/m_gammapositron [ Modules ]
NAME
m_gammapositron
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (MT,GJ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_gammapositron 22 23 use defs_basis 24 use m_abicore 25 use m_errors 26 use m_electronpositron 27 28 use defs_abitypes, only : MPI_type 29 use m_numeric_tools, only : invcb 30 use m_xctk, only : xcden 31 32 implicit none 33 34 private