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.

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 ]

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

[ Top ] [ Modules ]

NAME

  m_gammapositron

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2022 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