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 : !!***
|
