TABLE OF CONTENTS


ABINIT/gammapositron [ Functions ]

[ Top ] [ 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