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

PARENTS

      gammapositron_fft,poslifetime,posratecore

CHILDREN

      invcb

SOURCE

 94 subroutine gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,npt,rhocore,rhoer,rhopr,usecore)
 95 
 96 
 97 !This section has been created automatically by the script Abilint (TD).
 98 !Do not modify the following lines by hand.
 99 #undef ABI_FUNC
100 #define ABI_FUNC 'gammapositron'
101 !End of the abilint section
102 
103  implicit none
104 
105 !Arguments ------------------------------------
106 !scalars
107  integer,intent(in) :: igamma,ngr,npt,usecore
108 !arrays
109  real(dp),intent(in) :: grhocore2(ngr*usecore),grhoe2(ngr),rhocore(npt*usecore),rhoer(npt),rhopr(npt)
110  real(dp),intent(out) :: gamma(npt,2)
111 
112 !Local variables-------------------------------
113 !scalars
114  integer :: iloop,ipt
115  logical :: gga
116  real(dp),parameter :: alpha_gga=0.22d0,rsfac=0.6203504908994000_dp
117  real(dp) :: aa,bb,cc,dg1
118  real(dp) :: drs,eps,expgga,g0,g1,g2,gg
119  real(dp) :: kf,kk,nqtf2,ratio1,ratio2,ratio3,rho1,rho2,rhoe,rhop,sqrs,rs,rse,rsp
120 !arrays
121  real(dp),allocatable :: grho2(:),rhor(:),rsepts(:),rsppts(:)
122 
123 ! *************************************************************************
124 
125  gga=(ngr==npt.and.igamma/=0)
126 
127  if (usecore/=0.and.usecore/=1) then
128    MSG_ERROR('Wrong value for usecore !')
129  end if
130  if (igamma/=0.and.igamma/=1.and.igamma/=2.and.igamma/=3.and.igamma/=4) then
131    MSG_ERROR('Unknown electron-positron correlation !')
132  end if
133 
134  ABI_ALLOCATE(rhor,(npt))
135  ABI_ALLOCATE(rsepts,(npt))
136  if (gga)  then
137    ABI_ALLOCATE(grho2,(npt))
138  end if
139 
140 !Eventually compute positronic density radii
141  if (igamma==1.or.igamma==4) then
142    ABI_ALLOCATE(rsppts,(npt))
143    call invcb(rhopr(:),rsppts,npt)
144    rsppts(:)=rsfac*rsppts(:)
145  end if
146 
147 !Loop: iloop=1: compute enhancement factor using total   electronic density
148 !iloop=2: compute enhancement factor using valence electronic density
149 !===================================================================================
150  do iloop=1,1+usecore
151 
152 !  Compute electronic density radii
153    if (iloop==1.and.usecore==1) then
154      rhor(1:npt)=rhoer(1:npt)+rhocore(1:npt)
155    else
156      rhor(1:npt)=rhoer(1:npt)
157    end if
158    call invcb(rhor(:),rsepts,npt)
159    rsepts(:)=rsfac*rsepts(:)
160 
161 !  Gradients for GGA
162    if (gga) then
163      if (iloop==1.and.usecore==1) then
164        grho2(1:npt)=grhoe2(1:npt)+grhocore2(1:npt)
165      else
166        grho2(1:npt)=grhoe2(1:npt)
167      end if
168    end if
169 
170 !  Loop over grid points
171    do ipt=1,npt
172 
173      rhoe=rhor(ipt)
174      rhop=rhopr(ipt)
175      rse =rsepts(ipt)
176      gg=zero
177 
178 !    Testing feature: gamma=1
179 !    -----------------------------------------------------------------------------------
180      if (igamma==0) then
181 
182        gg=one
183 
184 !    Boronski and Nieminen
185 !    -----------------------------------------------------------------------------------
186      else if (igamma==1) then
187 
188        rsp =rsppts(ipt)
189        if (rhoe>rhop) then
190          rho1=rhoe;rho2=rhop;rs=rse
191        else
192          rho1=rhop;rho2=rhoe;rs=rsp
193        end if
194        drs=-third*rs/rho1;sqrs=sqrt(rs)
195        ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1
196        g0=one+1.23_dp*rs+0.8295_dp*sqrs**3-1.26_dp*rs**2+0.3286_dp*sqrs**5+sixth   *rs**3
197        g1=one+0.51_dp*rs                  +0.65_dp*rs**2-0.51_dp  *sqrs**5+0.176_dp*rs**3
198        g2=one+0.60_dp*rs                  +0.63_dp*rs**2-0.48_dp  *sqrs**5+0.167_dp*rs**3
199        dg1=drs*(0.51_dp+two*0.65_dp*rs-2.5_dp*0.51_dp*sqrs**3+three*0.176_dp*rs**2)
200        kk=half*rho1*dg1
201        aa= two  *kk-six    *g1+eight  *g2-two *g0
202        bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0
203        cc=       kk-four   *g1+eight  *g2-four*g0
204        gg=g0+ratio3*aa+ratio2*bb+ratio1*cc
205 
206 !      Boronski and Nieminen RPA limit
207 !      -----------------------------------------------------------------------------------
208      else if (igamma==2) then
209 
210        rs=rse;sqrs=sqrt(rs)
211        gg=one !This is experimental to avoid divergences
212        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
213 
214 !      Sterne and Kaiser
215 !      -----------------------------------------------------------------------------------
216      else if (igamma==3) then
217 
218        rs=rse;sqrs=sqrt(rs)
219        gg=one !This is experimental to avoid divergences
220        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
221 
222 !      Puska, Seitsonen and Nieminen
223 !      -----------------------------------------------------------------------------------
224      else if (igamma==4) then
225 
226        rsp =rsppts(ipt)
227        if (rhoe>rhop) then
228          rho1=rhoe;rho2=rhop;rs=rse
229        else
230          rho1=rhop;rho2=rhoe;rs=rsp
231        end if
232        drs=-third*rs/rho1;sqrs=sqrt(rs)
233        ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1
234        g0=one+1.2300_dp*rs+0.9889_dp*sqrs**3-1.4820_dp*rs**2+0.3956_dp*sqrs**5+sixth*rs**3
235        g1=one+2.0286_dp*rs-3.3892_dp*sqrs**3+3.0547_dp*rs**2-1.0540_dp*sqrs**5+sixth*rs**3
236        g2=one+0.2499_dp*rs+0.2949_dp*sqrs**3+0.6944_dp*rs**2-0.5339_dp*sqrs**5+sixth*rs**3
237        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)
238        kk=half*rho1*dg1
239        aa= two  *kk-six    *g1+eight  *g2-two *g0
240        bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0
241        cc=       kk-four   *g1+eight  *g2-four*g0
242        gg=g0+ratio3*aa+ratio2*bb+ratio1*cc
243 
244      end if ! igamma
245 
246      if (gga) then
247        kf=(three*pi*pi*rhoe)**third
248        nqtf2=(rhoe*sqrt(four*kf/pi))**2
249        eps=grho2(ipt)/nqtf2
250        if (eps<zero) then
251          MSG_ERROR('  problem, negative GGA espilon !')
252        end if
253        expgga=exp(-alpha_gga*eps*third)
254        gg=one+(gg-one)*expgga
255      end if
256 
257 !    Store enhancement factor
258      gamma(ipt,iloop)=gg
259 
260    end do ! ipt
261  end do ! iloop
262 
263  ABI_DEALLOCATE(rhor)
264  ABI_DEALLOCATE(rsepts)
265  if (igamma==1.or.igamma==4)  then
266    ABI_DEALLOCATE(rsppts)
267  end if
268  if (gga)  then
269    ABI_DEALLOCATE(grho2)
270  end if
271 
272 !Case usecore=0 (no core density)
273  if (usecore==0) gamma(:,2)=gamma(:,1)
274 
275 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)

PARENTS

      posdoppler,poslifetime

CHILDREN

      gammapositron,xcden

SOURCE

319 subroutine gammapositron_fft(electronpositron,gamma,gprimd,igamma,mpi_enreg,&
320 &                            n3xccc,nfft,ngfft,rhor_e,rhor_p,xccc3d)
321 
322 
323 !This section has been created automatically by the script Abilint (TD).
324 !Do not modify the following lines by hand.
325 #undef ABI_FUNC
326 #define ABI_FUNC 'gammapositron_fft'
327 !End of the abilint section
328 
329  implicit none
330 
331 !Arguments ------------------------------------
332 !scalars
333  integer,intent(in) :: igamma,n3xccc,nfft
334  type(electronpositron_type),pointer :: electronpositron
335  type(MPI_type),intent(in) :: mpi_enreg
336 !arrays
337  integer,intent(in) :: ngfft(18)
338  real(dp),intent(in) :: gprimd(3,3),rhor_e(nfft),rhor_p(nfft),xccc3d(n3xccc)
339  real(dp),intent(out) :: gamma(nfft,2)
340 
341 !Local variables-------------------------------
342 !scalars
343  integer :: cplex,ishift,ngr,ngrad,nspden_ep,usecore
344 !arrays
345  real(dp),parameter :: qphon(3)=(/zero,zero,zero/)
346  real(dp),allocatable :: grhocore2(:),grhoe2(:),rhoc(:,:,:),rhoe(:,:,:)
347 
348 ! *************************************************************************
349 
350 !Several useful constants
351  usecore=n3xccc/nfft
352  cplex=1;ishift=0;ngrad=1;nspden_ep=1
353  if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
354  ngr=0;if (ngrad==2) ngr=nfft
355 
356 !Allocate several arrays
357  ABI_ALLOCATE(rhoe,(nfft,nspden_ep,ngrad**2))
358  ABI_ALLOCATE(grhoe2,(ngr))
359  ABI_ALLOCATE(grhocore2,(ngr*usecore))
360 
361 !Store electronic density and its gradients
362  call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,&
363 & mpi_enreg%paral_kgb,qphon,rhor_e,rhoe)
364 
365 !Compute squared gradient of the electronic density
366  if (ngrad==2) then
367    grhoe2(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2
368    if (usecore>0) then
369      ABI_ALLOCATE(rhoc,(nfft,1,ngrad**2))
370      call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,&
371 &     mpi_enreg%paral_kgb,qphon,xccc3d,rhoc)
372      grhocore2(:)=rhoc(:,1,2)**2+rhoc(:,1,3)**2+rhoc(:,1,4)**2
373      ABI_DEALLOCATE(rhoc)
374    end if
375  end if
376 
377 !Compute enhancement factor on FFT grid
378  call gammapositron(gamma,grhocore2,grhoe2,igamma,ngr,nfft,xccc3d,&
379 & rhoe(:,1,1),rhor_p,usecore)
380 
381 !Release temporary memory
382  ABI_DEALLOCATE(rhoe)
383  ABI_DEALLOCATE(grhoe2)
384  ABI_DEALLOCATE(grhocore2)
385 
386 end subroutine gammapositron_fft

ABINIT/m_gammapositron [ Modules ]

[ Top ] [ Modules ]

NAME

  m_gammapositron

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_gammapositron
27 
28  use defs_basis
29  use defs_abitypes
30  use m_abicore
31  use m_errors
32  use m_electronpositron
33 
34  use m_numeric_tools,   only : invcb
35  use m_xctk,            only : xcden
36 
37  implicit none
38 
39  private