Branch data Line data Source code
1 : : !!****m* ABINIT/m_gammapositron
2 : : !! NAME
3 : : !! m_gammapositron
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (MT,GJ)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
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
35 : : !!***
36 : :
37 : : public :: gammapositron
38 : : public :: gammapositron_fft
39 : : !!***
40 : :
41 : : contains
42 : : !!***
43 : :
44 : : !!****f* ABINIT/gammapositron
45 : : !! NAME
46 : : !! gammapositron
47 : : !!
48 : : !! FUNCTION
49 : : !! Compute positron electron-positron enhancement factor (contact density) used to compute positron lifetime.
50 : : !! Input is positronic rhop(r) and electronic rhoe(r) at a given set of points.
51 : : !!
52 : : !! INPUTS
53 : : !! grhocore2(ngr)=square of the gradient of core electronic density rhocore (needed for GGA)
54 : : !! grhoe2(ngr)=square of the gradient of valence electronic density rhoer (needed for GGA)
55 : : !! igamma=type of enhancement factor:
56 : : !! 1: Boronski and Nieminen [2]
57 : : !! 2: Boronski and Nieminen, RPA limit [2]
58 : : !! 3: Sterne and Kaiser [3]
59 : : !! 4: Puska, Seitsonen and Nieminen [4]
60 : : !! See references below
61 : : !! ngr=size of grho2 array (0 if LDA, npt if GGA)
62 : : !! npt=number of real space points on which density is provided
63 : : !! rhocore(npt*usecore)=core electron density (bohr^-3)
64 : : !! rhoer(npt) =valence electron density (bohr^-3)
65 : : !! rhopr(npt) =positron density (bohr^-3)
66 : : !! usecore =1 if core density is not zero
67 : : !!
68 : : !! OUTPUT
69 : : !! gamma(npt,2) =electron-positron enhancement factor,
70 : : !! gamma(:,1): using total electronic density
71 : : !! gamma(:,2): using valence electronic density
72 : : !!
73 : : !! NOTES
74 : : !! References for electron-positron correlation functionals:
75 : : !! [1] J. Arponen and E. Pajanne, Ann. Phys. (N.Y.) 121, 343 (1979) [[cite:Arponen1979a]].
76 : : !! [2] E. Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986) [[cite:Boronski1986]].
77 : : !! [3] P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991) [[cite:Sterne1991]].
78 : : !! [4] M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994) [[cite:Puska1994]].
79 : : !! [5] B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1995) [[cite:Barbiellini1995]]
80 : : !!
81 : : !! SOURCE
82 : :
83 : 8864 : 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 : 8864 : real(dp),allocatable :: grho2(:),rhor(:),rsepts(:),rsppts(:)
102 : :
103 : : ! *************************************************************************
104 : :
105 [ - + - - ]: 8864 : gga=(ngr==npt.and.igamma/=0)
106 : :
107 [ - + ]: 8864 : if (usecore/=0.and.usecore/=1) then
108 : 0 : ABI_ERROR('Wrong value for usecore !')
109 : : end if
110 [ - + ]: 8864 : if (igamma/=0.and.igamma/=1.and.igamma/=2.and.igamma/=3.and.igamma/=4) then
111 : 0 : ABI_ERROR('Unknown electron-positron correlation !')
112 : : end if
113 : :
114 [ - + + - : 26592 : ABI_MALLOC(rhor,(npt))
+ - + - -
+ - + ]
115 [ - + + - : 17728 : ABI_MALLOC(rsepts,(npt))
- + - + ]
116 [ - + ]: 8864 : if (gga) then
117 [ # # # # : 0 : ABI_MALLOC(grho2,(npt))
# # # # ]
118 : : end if
119 : :
120 : : !Eventually compute positronic density radii
121 [ + + ]: 8864 : if (igamma==1.or.igamma==4) then
122 [ - + + - : 2994 : ABI_MALLOC(rsppts,(npt))
- + - + ]
123 : 1497 : call invcb(rhopr(:),rsppts,npt)
124 [ + + ]: 3071800 : 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 [ + + ]: 25092 : do iloop=1,1+usecore
131 : :
132 : : ! Compute electronic density radii
133 [ + + + + ]: 16228 : if (iloop==1.and.usecore==1) then
134 [ + + ]: 11924050 : rhor(1:npt)=rhoer(1:npt)+rhocore(1:npt)
135 : : else
136 [ + + ]: 15799250 : rhor(1:npt)=rhoer(1:npt)
137 : : end if
138 : 16228 : call invcb(rhor(:),rsepts,npt)
139 [ + + ]: 27723300 : rsepts(:)=rsfac*rsepts(:)
140 : :
141 : : ! Gradients for GGA
142 [ - + ]: 16228 : if (gga) then
143 [ # # # # ]: 0 : if (iloop==1.and.usecore==1) then
144 [ # # ]: 0 : grho2(1:npt)=grhoe2(1:npt)+grhocore2(1:npt)
145 : : else
146 [ # # ]: 0 : grho2(1:npt)=grhoe2(1:npt)
147 : : end if
148 : : end if
149 : :
150 : : ! Loop over grid points
151 [ + + ]: 27732164 : do ipt=1,npt
152 : :
153 : 27707072 : rhoe=rhor(ipt)
154 : 27707072 : rhop=rhopr(ipt)
155 : 27707072 : rse =rsepts(ipt)
156 [ + + + + : 27707072 : 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 : 5819836 : rsp =rsppts(ipt)
169 [ + + ]: 5819836 : if (rhoe>rhop) then
170 : : rho1=rhoe;rho2=rhop;rs=rse
171 : : else
172 : 20818 : rho1=rhop;rho2=rhoe;rs=rsp
173 : : end if
174 : 5819836 : drs=-third*rs/rho1;sqrs=sqrt(rs)
175 : 5819836 : ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1
176 : 5819836 : 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 : 5819836 : g1=one+0.51_dp*rs +0.65_dp*rs**2-0.51_dp *sqrs**5+0.176_dp*rs**3
178 : 5819836 : g2=one+0.60_dp*rs +0.63_dp*rs**2-0.48_dp *sqrs**5+0.167_dp*rs**3
179 : 5819836 : 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 : 5819836 : kk=half*rho1*dg1
181 : 5819836 : aa= two *kk-six *g1+eight *g2-two *g0
182 : 5819836 : bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0
183 : 5819836 : cc= kk-four *g1+eight *g2-four*g0
184 : 5819836 : gg=g0+ratio3*aa+ratio2*bb+ratio1*cc
185 : :
186 : : ! Boronski and Nieminen RPA limit
187 : : ! -----------------------------------------------------------------------------------
188 : : else if (igamma==2) then
189 : :
190 : 17180256 : rs=rse;sqrs=sqrt(rs)
191 : 17180256 : gg=one !This is experimental to avoid divergences
192 [ + + ]: 17180256 : 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 : 4061856 : rs=rse;sqrs=sqrt(rs)
199 : 4061856 : gg=one !This is experimental to avoid divergences
200 [ + - ]: 4061856 : 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 : 243296 : rsp =rsppts(ipt)
207 [ + + ]: 243296 : if (rhoe>rhop) then
208 : : rho1=rhoe;rho2=rhop;rs=rse
209 : : else
210 : 5332 : rho1=rhop;rho2=rhoe;rs=rsp
211 : : end if
212 : 243296 : drs=-third*rs/rho1;sqrs=sqrt(rs)
213 : 243296 : ratio1=rho2/rho1;ratio2=ratio1*ratio1;ratio3=ratio2*ratio1
214 : 243296 : 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 : 243296 : 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 : 243296 : 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 : 243296 : 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 : 243296 : kk=half*rho1*dg1
219 : 243296 : aa= two *kk-six *g1+eight *g2-two *g0
220 : 243296 : bb=-three*kk+11.0_dp*g1-16.0_dp*g2+five*g0
221 : 243296 : cc= kk-four *g1+eight *g2-four*g0
222 : 243296 : gg=g0+ratio3*aa+ratio2*bb+ratio1*cc
223 : :
224 : : end if ! igamma
225 : :
226 [ - + ]: 27707072 : if (gga) then
227 : 0 : kf=(three*pi*pi*rhoe)**third
228 : 0 : nqtf2=(rhoe*sqrt(four*kf/pi))**2
229 : 0 : eps=grho2(ipt)/nqtf2
230 [ # # ]: 0 : if (eps<zero) then
231 : 0 : ABI_ERROR(' problem, negative GGA espilon !')
232 : : end if
233 : 0 : expgga=exp(-alpha_gga*eps*third)
234 : 0 : gg=one+(gg-one)*expgga
235 : : end if
236 : :
237 : : ! Store enhancement factor
238 : 27723300 : gamma(ipt,iloop)=gg
239 : :
240 : : end do ! ipt
241 : : end do ! iloop
242 : :
243 [ - + ]: 8864 : ABI_FREE(rhor)
244 [ - + ]: 8864 : ABI_FREE(rsepts)
245 [ + + ]: 8864 : if (igamma==1.or.igamma==4) then
246 [ - + ]: 1497 : ABI_FREE(rsppts)
247 : : end if
248 [ - + ]: 8864 : if (gga) then
249 [ # # ]: 0 : ABI_FREE(grho2)
250 : : end if
251 : :
252 : : !Case usecore=0 (no core density)
253 [ + + + + ]: 3882564 : if (usecore==0) gamma(:,2)=gamma(:,1)
254 : :
255 [ - + - + ]: 8864 : end subroutine gammapositron
256 : : !!***
257 : :
258 : : !!****f* ABINIT/gammapositron_fft
259 : : !! NAME
260 : : !! gammapositron_fft
261 : : !!
262 : : !! FUNCTION
263 : : !! Compute positron electron-positron enhancement factor on a real space FFT grid.
264 : : !!
265 : : !! INPUTS
266 : : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
267 : : !! gprimd(3,3)=dimensional primitive translations in reciprocal space (bohr^-1)
268 : : !! igamma=type of enhancement factor:
269 : : !! -1: gamma=one (test)
270 : : !! 1: Boronski and Nieminen [1]
271 : : !! 2: Boronski and Nieminen, RPA limit [1]
272 : : !! 3: Sterne and Kaiser [2]
273 : : !! 4: Puska, Seitsonen and Nieminen [3]
274 : : !! mpi_enreg=information about MPI parallelization
275 : : !! n3xccc=dimension of the xccc3d array (0 or nfft).
276 : : !! nfft=number of FFT grid points
277 : : !! ngfft(18)=contain all needed information about 3D FFT
278 : : !! rhor_e(nfft)=real space density electronic density (total density)
279 : : !! rhor_p(nfft)=real space density positronic density
280 : : !! xccc3d(n3xccc)=3D core electron density for XC core correction
281 : : !!
282 : : !! OUTPUT
283 : : !! gamma(nfft,2)=electron-positron enhancement factor,
284 : : !! gamma(:,1): using total electronic density
285 : : !! gamma(:,2): using valence electronic density
286 : : !!
287 : : !! NOTES
288 : : !! The input densities (rhor_e and rhor_p) should be positive everywhere
289 : : !! (call mkdenpos routine before entering this one)
290 : : !!
291 : : !! SOURCE
292 : :
293 : 184 : subroutine gammapositron_fft(electronpositron,gamma,gprimd,igamma,mpi_enreg,&
294 : 184 : & 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 : 184 : real(dp),allocatable :: grhocore2(:),grhoe2(:),rhoc(:,:,:),rhoe(:,:,:)
312 : :
313 : : ! *************************************************************************
314 : :
315 : : !Several useful constants
316 : 184 : usecore=n3xccc/nfft
317 : 184 : cplex=1;ishift=0;ngrad=1;nspden_ep=1
318 [ - + ]: 184 : if (electronpositron%ixcpositron==3.or.electronpositron%ixcpositron==31) ngrad=2
319 [ - + ]: 184 : ngr=0;if (ngrad==2) ngr=nfft
320 : :
321 : : !Allocate several arrays
322 [ - + + - : 736 : ABI_MALLOC(rhoe,(nfft,nspden_ep,ngrad**2))
- + + - +
- + - - +
- + ]
323 [ + - - - : 552 : ABI_MALLOC(grhoe2,(ngr))
+ - - + -
+ - + ]
324 [ + - - - : 552 : ABI_MALLOC(grhocore2,(ngr*usecore))
+ - - + -
+ - + ]
325 : :
326 : : !Store electronic density and its gradients
327 : 184 : 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 [ - + ]: 184 : if (ngrad==2) then
331 [ # # ]: 0 : grhoe2(:)=rhoe(:,1,2)**2+rhoe(:,1,3)**2+rhoe(:,1,4)**2
332 [ # # ]: 0 : if (usecore>0) then
333 [ # # # # : 0 : ABI_MALLOC(rhoc,(nfft,1,ngrad**2))
# # # # #
# # # ]
334 : 0 : call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden_ep,qphon,xccc3d,rhoc)
335 [ # # ]: 0 : grhocore2(:)=rhoc(:,1,2)**2+rhoc(:,1,3)**2+rhoc(:,1,4)**2
336 [ # # ]: 0 : 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 : 184 : & rhoe(:,1,1),rhor_p,usecore)
343 : :
344 : : !Release temporary memory
345 [ - + ]: 184 : ABI_FREE(rhoe)
346 [ - + ]: 184 : ABI_FREE(grhoe2)
347 [ - + ]: 184 : ABI_FREE(grhocore2)
348 : :
349 : 184 : end subroutine gammapositron_fft
350 : : !!***
351 : :
352 : : end module m_gammapositron
353 : : !!***
|
