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.
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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 [1] 2: Boronski and Nieminen, RPA limit [1] 3: Sterne and Kaiser [2] 4: Puska, Seitsonen and Nieminen [3] 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] Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986). [2] P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991). [3] M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994). [4] B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1994)
PARENTS
gammapositron_fft,poslifetime,posratecore
CHILDREN
invcb
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,npt,rhocore,rhoer,rhopr,usecore) 60 61 use defs_basis 62 use m_profiling_abi 63 use m_errors 64 65 !This section has been created automatically by the script Abilint (TD). 66 !Do not modify the following lines by hand. 67 #undef ABI_FUNC 68 #define ABI_FUNC 'gammapositron' 69 use interfaces_41_xc_lowlevel 70 !End of the abilint section 71 72 implicit none 73 74 !Arguments ------------------------------------ 75 !scalars 76 integer,intent(in) :: igamma,ngr,npt,usecore 77 !arrays 78 real(dp),intent(in) :: grhocore2(ngr*usecore),grhoe2(ngr),rhocore(npt*usecore),rhoer(npt),rhopr(npt) 79 real(dp),intent(out) :: gamma(npt,2) 80 81 !Local variables------------------------------- 82 !scalars 83 integer :: iloop,ipt 84 logical :: gga 85 real(dp),parameter :: alpha_gga=0.22d0,rsfac=0.6203504908994000_dp 86 real(dp) :: aa,bb,cc,dg1 87 real(dp) :: drs,eps,expgga,g0,g1,g2,gg 88 real(dp) :: kf,kk,nqtf2,ratio1,ratio2,ratio3,rho1,rho2,rhoe,rhop,sqrs,rs,rse,rsp 89 !arrays 90 real(dp),allocatable :: grho2(:),rhor(:),rsepts(:),rsppts(:) 91 92 ! ************************************************************************* 93 94 gga=(ngr==npt.and.igamma/=0) 95 96 if (usecore/=0.and.usecore/=1) then 97 MSG_ERROR('Wrong value for usecore !') 98 end if 99 if (igamma/=0.and.igamma/=1.and.igamma/=2.and.igamma/=3.and.igamma/=4) then 100 MSG_ERROR('Unknown electron-positron correlation !') 101 end if 102 103 ABI_ALLOCATE(rhor,(npt)) 104 ABI_ALLOCATE(rsepts,(npt)) 105 if (gga) then 106 ABI_ALLOCATE(grho2,(npt)) 107 end if 108 109 !Eventually compute positronic density radii 110 if (igamma==1.or.igamma==4) then 111 ABI_ALLOCATE(rsppts,(npt)) 112 call invcb(rhopr(:),rsppts,npt) 113 rsppts(:)=rsfac*rsppts(:) 114 end if 115 116 !Loop: iloop=1: compute enhancement factor using total electronic density 117 !iloop=2: compute enhancement factor using valence electronic density 118 !=================================================================================== 119 do iloop=1,1+usecore 120 121 ! Compute electronic density radii 122 if (iloop==1.and.usecore==1) then 123 rhor(1:npt)=rhoer(1:npt)+rhocore(1:npt) 124 else 125 rhor(1:npt)=rhoer(1:npt) 126 end if 127 call invcb(rhor(:),rsepts,npt) 128 rsepts(:)=rsfac*rsepts(:) 129 130 ! Gradients for GGA 131 if (gga) then 132 if (iloop==1.and.usecore==1) then 133 grho2(1:npt)=grhoe2(1:npt)+grhocore2(1:npt) 134 else 135 grho2(1:npt)=grhoe2(1:npt) 136 end if 137 end if 138 139 ! Loop over grid points 140 do ipt=1,npt 141 142 rhoe=rhor(ipt) 143 rhop=rhopr(ipt) 144 rse =rsepts(ipt) 145 gg=zero 146 147 ! Testing feature: gamma=1 148 ! ----------------------------------------------------------------------------------- 149 if (igamma==0) then 150 151 gg=one 152 153 ! Boronski and Nieminen 154 ! ----------------------------------------------------------------------------------- 155 else if (igamma==1) then 156 157 rsp =rsppts(ipt) 158 if (rhoe>rhop) then 159 rho1=rhoe;rho2=rhop;rs=rse 160 else 161 rho1=rhop;rho2=rhoe;rs=rsp 162 end if 163 drs=-third*rs/rho1;sqrs=sqrt(rs) 164 ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1 165 g0=one+1.23_dp*rs+0.8295_dp*sqrs**3-1.26_dp*rs**2+0.3286_dp*sqrs**5+sixth *rs**3 166 g1=one+0.51_dp*rs +0.65_dp*rs**2-0.51_dp *sqrs**5+0.176_dp*rs**3 167 g2=one+0.60_dp*rs +0.63_dp*rs**2-0.48_dp *sqrs**5+0.167_dp*rs**3 168 dg1=drs*(0.51_dp+two*0.65_dp*rs-2.5_dp*0.51_dp*sqrs**3+three*0.176_dp*rs**2) 169 kk=half*rho1*dg1 170 aa= two *kk-six *g1+eight *g2-two *g0 171 bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0 172 cc= kk-four *g1+eight *g2-four*g0 173 gg=g0+ratio3*aa+ratio2*bb+ratio1*cc 174 175 ! Boronski and Nieminen RPA limit 176 ! ----------------------------------------------------------------------------------- 177 else if (igamma==2) then 178 179 rs=rse;sqrs=sqrt(rs) 180 gg=one !This is experimental to avoid divergences 181 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 182 183 ! Sterne and Kaiser 184 ! ----------------------------------------------------------------------------------- 185 else if (igamma==3) then 186 187 rs=rse;sqrs=sqrt(rs) 188 gg=one !This is experimental to avoid divergences 189 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 190 191 ! Puska, Seitsonen and Nieminen 192 ! ----------------------------------------------------------------------------------- 193 else if (igamma==4) then 194 195 rsp =rsppts(ipt) 196 if (rhoe>rhop) then 197 rho1=rhoe;rho2=rhop;rs=rse 198 else 199 rho1=rhop;rho2=rhoe;rs=rsp 200 end if 201 drs=-third*rs/rho1;sqrs=sqrt(rs) 202 ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1 203 g0=one+1.2300_dp*rs+0.9889_dp*sqrs**3-1.4820_dp*rs**2+0.3956_dp*sqrs**5+sixth*rs**3 204 g1=one+2.0286_dp*rs-3.3892_dp*sqrs**3+3.0547_dp*rs**2-1.0540_dp*sqrs**5+sixth*rs**3 205 g2=one+0.2499_dp*rs+0.2949_dp*sqrs**3+0.6944_dp*rs**2-0.5339_dp*sqrs**5+sixth*rs**3 206 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) 207 kk=half*rho1*dg1 208 aa= two *kk-six *g1+eight *g2-two *g0 209 bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0 210 cc= kk-four *g1+eight *g2-four*g0 211 gg=g0+ratio3*aa+ratio2*bb+ratio1*cc 212 213 end if ! igamma 214 215 if (gga) then 216 kf=(three*pi*pi*rhoe)**third 217 nqtf2=(rhoe*sqrt(four*kf/pi))**2 218 eps=grho2(ipt)/nqtf2 219 if (eps<zero) then 220 MSG_ERROR(' problem, negative GGA espilon !') 221 end if 222 expgga=exp(-alpha_gga*eps*third) 223 gg=one+(gg-one)*expgga 224 end if 225 226 ! Store enhancement factor 227 gamma(ipt,iloop)=gg 228 229 end do ! ipt 230 end do ! iloop 231 232 ABI_DEALLOCATE(rhor) 233 ABI_DEALLOCATE(rsepts) 234 if (igamma==1.or.igamma==4) then 235 ABI_DEALLOCATE(rsppts) 236 end if 237 if (gga) then 238 ABI_DEALLOCATE(grho2) 239 end if 240 241 !Case usecore=0 (no core density) 242 if (usecore==0) gamma(:,2)=gamma(:,1) 243 244 end subroutine gammapositron