Branch data Line data Source code
1 : : !!****m* ABINIT/m_xcpositron
2 : : !! NAME
3 : : !! m_xcpositron
4 : : !!
5 : : !! FUNCTION
6 : : !! Compute electron-positron correlation potentials and energy density.
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (GJ,MT)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_xcpositron
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_abicore
27 : :
28 : : use m_numeric_tools, only : invcb
29 : :
30 : : implicit none
31 : :
32 : : private
33 : : !!***
34 : :
35 : : public :: xcpositron
36 : : !!***
37 : :
38 : : contains
39 : : !!***
40 : :
41 : : !!****f* ABINIT/xcpositron
42 : : !! NAME
43 : : !! xcpositron
44 : : !!
45 : : !! FUNCTION
46 : : !! Compute electron-positron correlation potentials and energy density.
47 : : !! Used electron-positron correlation functional is controlled by ixcpositron argument.
48 : : !! Returns Fxc, Vxc_pos, Vxc_el from input rhor_pos and rhor_el for positron and electrons.
49 : : !!
50 : : !! INPUTS
51 : : !! grhoe2(ngr)=square of the gradient of electronic density rhoe (needed for GGA)
52 : : !! ixcpositron=type of electron-positron correlation functional:
53 : : !! 1: LDA zero positron density limit parametrized by Arponen & Pajanne
54 : : !! and provided by Boronski & Nieminen [1,2]
55 : : !! 11: LDA zero positron density limit parametrized by Arponen & Pajanne
56 : : !! and fitted by Sterne & Kaiser [1,3]
57 : : !! 2: LDA electron-positron correlation
58 : : !! provided by Puska, Seitsonen, and Nieminen [1,4]
59 : : !! 3: GGA zero positron density limit parametrized by Arponen & Pajanne
60 : : !! and provided by Boronski & Nieminen [1,2,5]
61 : : !! 31: GGA zero positron density limit parametrized by Arponen & Pajanne
62 : : !! and fitted by Sterne & Kaiser [1,3,5]
63 : : !! See references below
64 : : !! ngr=size of grho2 array (0 if LDA, npt if GGA)
65 : : !! npt=number of real space points on which density is provided
66 : : !! posdensity0_limit=True if we are in the zero positron density limit
67 : : !! rhoer(npt)=electron density (bohr^-3)
68 : : !! rhopr(npt)=positron density (bohr^-3)
69 : : !!
70 : : !! OUTPUT
71 : : !! fnxc(npt)=correlation energy per unit volume fxc
72 : : !! vxce(npt)=correlation potential for electron dfxc/drhoe (hartree)
73 : : !! vxcp(npt)=correlation potential for positron dfxc/drhop (hartree)
74 : : !! vxcegr(ngr)= 1/|gradRhoe| dfxc/d|gradRhoe| (empty if LDA, i.e. ngr=0)
75 : : !! Optional outputs:
76 : : !! dvxce(npt)=partial second derivatives of the xc energy wr to the electronic density
77 : : !! dvxce(:)=dVxce/dRhoe
78 : : !! dvxcp(npt)=partial second derivatives of the xc energy wr to the positronic density
79 : : !! dvxcp(:)=dVxcp/drhop
80 : : !!
81 : : !! NOTES
82 : : !! References for electron-positron correlation functionals:
83 : : !! [1] J. Arponen and E. Pajanne, Ann. Phys. (N.Y.) 121, 343 (1979) [[cite:Arponen1979a]].
84 : : !! [2] E. Boronski and R.M. Nieminen, Phys. Rev. B 34, 3820 (1986) [[cite:Boronski1986]].
85 : : !! [3] P.A. Sterne and J.H. Kaiser, Phys. Rev. B 43, 13892 (1991) [[cite:Sterne1991]].
86 : : !! [4] M.J. Puska, A.P. Seitsonen and R.M. Nieminen, Phys. Rev. B 52, 10947 (1994) [[cite:Puska1994]].
87 : : !! [5] B. Barbiellini, M.J. Puska, T. Torsti and R.M.Nieminen, Phys. Rev. B 51, 7341 (1995) [[cite:Barbiellini1995]]
88 : : !!
89 : : !! SOURCE
90 : :
91 : 138520 : subroutine xcpositron(fnxc,grhoe2,ixcpositron,ngr,npt,posdensity0_limit,rhoer,rhopr,vxce,vxcegr,vxcp,&
92 : 7427 : & dvxce,dvxcp) ! optional arguments
93 : :
94 : : !Arguments ------------------------------------
95 : : !scalars
96 : : integer,intent(in) :: ixcpositron,ngr,npt
97 : : logical,intent(in) :: posdensity0_limit
98 : : !arrays
99 : : real(dp),intent(in) :: grhoe2(ngr),rhoer(npt),rhopr(npt)
100 : : real(dp),intent(out) :: fnxc(npt),vxce(npt),vxcegr(ngr),vxcp(npt)
101 : : real(dp),intent(out),optional :: dvxce(npt),dvxcp(npt)
102 : :
103 : : !Local variables-------------------------------
104 : : !scalars
105 : : integer,parameter :: idebug=0
106 : : integer :: ipt
107 : : logical :: gga,need_dvxce,need_dvxcp
108 : : real(dp),parameter :: alpha_gga=0.22_dp
109 : : real(dp),parameter :: ap_a1=-1.56_dp,ap_b1=0.051_dp,ap_c1=-0.081_dp,ap_d1=1.14_dp
110 : : real(dp),parameter :: ap_a2=-0.92305_dp,ap_b2=-0.05459_dp
111 : : real(dp),parameter :: ap_a3=-0.6298_dp,ap_b3=-13.15111_dp,ap_c3=2.8655_dp
112 : : real(dp),parameter :: ap_a4=-179856.2768_dp,ap_b4=186.4207_dp,ap_c4=-0.524_dp
113 : : real(dp),parameter :: ap_psn_limit=0.7_dp
114 : : real(dp),parameter :: ap_psn_1=0.9_dp*ap_psn_limit,ap_psn_2=1.1_dp*ap_psn_limit
115 : : real(dp),parameter :: fpi3=third*four_pi
116 : : real(dp),parameter :: psn_aa=69.7029_dp,psn_ba=-107.4927_dp,psn_bb=141.8458_dp
117 : : real(dp),parameter :: psn_ca=23.7182_dp,psn_cb=-33.6472_dp ,psn_cc=5.21152_dp
118 : : real(dp),parameter :: sk_a=-1.56_dp,sk_b=0.1324_dp,sk_c=-4.092_dp,sk_d=51.96_dp,sk_e=0.7207_dp
119 : : real(dp),parameter :: rsfac=0.6203504908994000_dp
120 : : real(dp) :: arse,brse,crse,darse,dbrse,dcrse,d2arse,d2brse,d2crse
121 : : real(dp) :: d2eps,deps,dexc,dexcdg,dexc_p,d2expgga,dexpgga,d2invrs,dinvrs,d2kf,dkf,d2nqtf2,dnqtf2
122 : : real(dp) :: drse,drsp,d2exc,d2exc_p,d2rse,d2rsp,d2sqr,dsqr
123 : : real(dp) :: eexp,eps,exc,exc_p,expgga,invf,dinvf,d2invf,invrhoe,invrhop,invrs,invsqr
124 : : real(dp) :: kf,logrs,nqtf2,opr2,ratio_ap,ratio_psn,rhoe,rhop,rse,rsp,sqr
125 : : character(len=500) :: msg
126 : : !arrays
127 : 69260 : real(dp),allocatable :: rsepts(:),rsppts(:)
128 : :
129 : : ! *************************************************************************
130 : :
131 : 69260 : gga=(ngr==npt)
132 : 69260 : need_dvxce=present(dvxce)
133 : 69260 : need_dvxcp=present(dvxcp)
134 : :
135 [ - + - - ]: 69260 : if (gga.and.ixcpositron==2) then
136 : 0 : msg = 'xcpositron: GGA not yet implemented for ixcpositron=2 !'
137 : 0 : ABI_ERROR(msg)
138 : : end if
139 [ + + - + ]: 69260 : if (posdensity0_limit.and.ixcpositron==2) then
140 : 0 : msg = 'xcpositron: ixcpositron=2 cannot be treated in the zero positron density limit !'
141 : 0 : ABI_ERROR(msg)
142 : : end if
143 [ + + ]: 69260 : if (abs(ixcpositron)/=1.and.ixcpositron/=11.and.ixcpositron/=2.and.ixcpositron/=3.and.ixcpositron/=31) then
144 : 0 : msg = 'xcpositron: unknown electron-positron correlation !'
145 : 0 : ABI_ERROR(msg)
146 : : end if
147 : :
148 : : !Compute density radii for rhor_el, rhor_pos
149 [ - + + - : 207780 : ABI_MALLOC(rsepts,(npt))
+ - + - -
+ - + ]
150 : 69260 : call invcb(rhoer(:),rsepts,npt)
151 [ + + ]: 147695946 : rsepts(:)=rsfac*rsepts(:)
152 [ + + ]: 69260 : if (ixcpositron==2) then
153 [ - + + - : 1232 : ABI_MALLOC(rsppts,(npt))
- + - + ]
154 : 616 : call invcb(rhopr(:),rsppts,npt)
155 [ + + ]: 716856 : rsppts(:)=rsfac*rsppts(:)
156 : : end if
157 : :
158 : : !Loop over grid points
159 : : !----------------------------------------------------
160 [ + + ]: 147695946 : do ipt=1,npt
161 : :
162 : 147626686 : rhoe=rhoer(ipt)
163 : 147626686 : rhop=rhopr(ipt)
164 : 147626686 : exc=zero;dexc=zero;d2exc=zero;dexcdg=zero
165 : :
166 : 147626686 : rse=rsepts(ipt)
167 : 147626686 : invrhoe=one/rhoe
168 : 147626686 : drse=-third*rse*invrhoe
169 [ + + ]: 147626686 : if (need_dvxce) d2rse= four/nine*rse*invrhoe**2
170 : :
171 : : ! Arponen & Pajane parametrization for electron
172 [ + + ]: 147626686 : if (ixcpositron/=11.and.ixcpositron/=31) then
173 [ + + ]: 143164406 : if (rse<0.302_dp) then
174 : 30553159 : invrs=one/rse;invsqr=sqrt(invrs);logrs=log(rse)
175 : 30553159 : exc =ap_a1*invsqr+(ap_b1*logrs+ap_c1)*logrs+ap_d1
176 : 30553159 : dexc=drse*invrs*(-half*ap_a1*invsqr+two*ap_b1*logrs+ap_c1)
177 [ + + ]: 30553159 : if (need_dvxce) d2exc=(d2rse/drse-drse*invrs)*dexc+drse**2*invrs**2*(quarter*ap_a1*invsqr+two*ap_b1)
178 [ + - + + ]: 112611247 : else if (rse>=0.302_dp.and.rse<=0.56_dp) then
179 : 2970406 : invrs=one/rse
180 : 2970406 : exc =ap_a2+ap_b2*invrs**2
181 : 2970406 : dexc=-drse*ap_b2*two*invrs**3
182 [ + + ]: 2970406 : if (need_dvxce) d2exc=d2rse/drse*dexc+six*drse**2*ap_b2*invrs**4
183 [ + - ]: 109640841 : else if (rse>0.56_dp.and.rse<=8.0_dp) then
184 : 109640841 : invrs=one/(rse+2.5_dp)
185 : 109640841 : dinvrs=-drse*invrs**2
186 : : ! jmb : d2rse initialized only if need_dvxce = .True.
187 [ + + ]: 109640841 : if (need_dvxce) d2invrs=-d2rse*invrs**2-two*invrs*drse**2
188 : 109640841 : exc =ap_a3+ap_b3*invrs**2+ap_c3*invrs
189 : 109640841 : dexc=two*ap_b3*invrs*dinvrs+ap_c3*dinvrs
190 [ + + ]: 109640841 : if (need_dvxce) d2exc=two*ap_b3*dinvrs**2+(two*ap_b3*invrs+ap_c3)*d2invrs
191 : : else
192 : 0 : exc =ap_a4*rhoe**2+ap_b4*rhoe+ap_c4
193 : 0 : dexc =two*ap_a4*rhoe+ap_b4
194 [ # # ]: 0 : if (need_dvxce) d2exc=two*ap_a4
195 : : end if
196 : :
197 : : ! Sterne & Kaiser parametrization for electron
198 : : else
199 : 4462280 : eexp=exp(-(rse+sk_c)**2/sk_d)
200 : 4462280 : opr2=(one+rse**2)
201 : 4462280 : arse=atan(rse)
202 : 4462280 : exc = sk_a/sqrt(arse)+sk_b*eexp+sk_e
203 : 4462280 : dexc= -(two*sk_b*eexp*(sk_c+rse)/sk_d + sk_a/(two*opr2*sqrt(arse)**3))*drse
204 [ + + ]: 4462280 : if (need_dvxce) d2exc=-(two*sk_b*eexp*(sk_c+rse)/sk_d + sk_a/(two*opr2*arse**1.5_dp))*d2rse &
205 : : & +(two*sk_b*eexp*(two*sk_c**2-sk_d+four*sk_c*rse+two*rse**2)/sk_d**2 &
206 : 648000 : & +sk_a*(three+four*rse*arse)/(four*opr2**2*sqrt(arse)**5))*drse**2
207 : : end if
208 : :
209 : : ! Puska, Seitsonen and Nieminen parametrization for positron
210 [ + + + + ]: 147626686 : if (ixcpositron==2.and.rse>=ap_psn_1) then
211 : 559720 : rsp=rsppts(ipt)
212 : 559720 : invrhop=one/rhop
213 : 559720 : drsp=-third*rsp*invrhop
214 [ - + ]: 559720 : if (need_dvxcp) d2rsp= four/nine*rsp*invrhop**2
215 : 559720 : exc_p=zero;dexc_p=zero;d2exc_p=zero
216 [ - + ]: 559720 : if (rsp<0.302_dp) then
217 : 0 : invrs=one/rsp;invsqr=sqrt(invrs);logrs=log(rsp)
218 : 0 : exc_p =ap_a1*invsqr+(ap_b1*logrs+ap_c1)*logrs+ap_d1
219 : 0 : dexc_p=drsp*invrs*(-half*ap_a1*invsqr+two*ap_b1*logrs+ap_c1)
220 [ # # ]: 0 : if (need_dvxcp) d2exc_p=(d2rsp/drsp-drsp*invrs)*dexc_p+drsp**2*invrs**2*(quarter*ap_a1*invsqr+two*ap_b1)
221 [ + - - + ]: 559720 : else if (rsp>=0.302_dp.and.rsp<=0.56_dp) then
222 : 0 : invrs=one/rsp
223 : 0 : exc_p =ap_a2+ap_b2*invrs**2
224 : 0 : dexc_p=-drsp*ap_b2*two*invrs**3
225 [ # # ]: 0 : if (need_dvxcp) d2exc_p=d2rsp/drsp*dexc_p+six*drsp**2*ap_b2*invrs**4
226 [ + + ]: 559720 : else if (rsp>0.56_dp.and.rsp<=8.0_dp) then
227 : 331912 : invrs=one/(rsp+2.5_dp)
228 : 331912 : dinvrs=-drsp*invrs**2
229 : : ! jmb : d2rsp initialized only if need_dvxcp = .True.*
230 [ - + ]: 331912 : if (need_dvxcp) d2invrs=-d2rsp*invrs**2-two*invrs*drsp**2
231 : 331912 : exc_p =ap_a3+ap_b3*invrs**2+ap_c3*invrs
232 : 331912 : dexc_p=two*ap_b3*invrs*dinvrs+ap_c3*dinvrs
233 [ - + ]: 331912 : if (need_dvxcp) d2exc_p=two*ap_b3*dinvrs**2+(two*ap_b3*invrs+ap_c3)*d2invrs
234 : : else
235 : 227808 : exc_p =ap_a4*rhop**2+ap_b4*rhop+ap_c4
236 : 227808 : dexc_p =two*ap_a4*rhop+ap_b4
237 [ - + ]: 227808 : if (need_dvxcp) d2exc_p=two*ap_a4
238 : : end if
239 : : end if
240 : :
241 : : ! GGA correction
242 [ - + ]: 147626686 : if (gga) then
243 : 0 : kf=(three*pi*pi*rhoe)**third
244 : 0 : nqtf2=(rhoe*sqrt(four*kf/pi))**2
245 : 0 : eps=grhoe2(ipt)/nqtf2
246 [ # # ]: 0 : if (eps<zero) then
247 : 0 : ABI_ERROR('xcpositron: problem, negative GGA espilon !')
248 : : end if
249 : 0 : expgga=exp(-alpha_gga*eps*third)
250 : :
251 : 0 : dkf=pi*pi/(sqrt(three*pi*pi*rhoe)**third)
252 : 0 : d2kf=-two*pi*pi*pi*pi*(three*pi*pi*rhoe)**(-5.0_dp/3.0_dp)
253 : 0 : sqr=sqrt(four*kf/pi)
254 : 0 : dsqr=(four*dkf/pi)/(two*sqr)
255 : 0 : d2sqr=two/(pi*sqr*dkf)*(d2kf*sqr-dsqr*dkf)
256 : 0 : nqtf2=(rhoe*sqr)**two
257 : 0 : dnqtf2=two*(sqr+rhoe*dsqr)*rhoe*sqr
258 : : d2nqtf2=two*(rhoe*sqr*(two*dsqr+rhoe*d2sqr) &
259 : : & +sqr*(sqr+rhoe*dsqr) &
260 : 0 : & +rhoe*(sqr+rhoe*dsqr) )
261 : 0 : deps=-grhoe2(ipt)*dnqtf2/(nqtf2**two)
262 : 0 : d2eps=-grhoe2(ipt)/(nqtf2*nqtf2*dnqtf2)*(d2nqtf2*nqtf2*nqtf2-two*nqtf2*dnqtf2*dnqtf2)
263 : 0 : dexpgga=-alpha_gga*third*deps*expgga
264 : 0 : d2expgga=-alpha_gga*third*(d2eps*expgga+deps*dexpgga)
265 : :
266 : 0 : exc = exc *expgga
267 : 0 : dexc=(dexc*expgga+exc*dexpgga)
268 [ # # ]: 0 : if (need_dvxce) d2exc=d2exc*expgga+two*dexc*dexpgga+exc*d2expgga
269 [ # # ]: 0 : if (abs(grhoe2(ipt))<1.e24_dp) dexcdg=-exc*alpha_gga*two_thirds/nqtf2
270 : : end if
271 : :
272 : : ! Computation of XC energy, potentials and kernels
273 : : ! Zero positron density limit
274 [ + + + + ]: 147626686 : if (ixcpositron/=2.or.rse<ap_psn_1) then
275 : 147066966 : fnxc(ipt)=rhop*exc
276 : 147066966 : vxce(ipt)=rhop*dexc
277 : 147066966 : vxcp(ipt)=exc
278 [ + + ]: 147066966 : if (need_dvxce) dvxce(ipt)=rhop*d2exc
279 [ - + ]: 147066966 : if (need_dvxcp) dvxcp(ipt)=zero
280 [ - + ]: 147066966 : if (gga) vxcegr(ipt)=rhop*dexcdg
281 : : else
282 : : ! Puska, Seitsonen and Nieminen functional
283 : 559720 : arse=psn_aa+psn_ba*rse+psn_ca*rse**2
284 : 559720 : brse=psn_ba+psn_bb*rse+psn_cb*rse**2
285 : 559720 : crse=psn_ca+psn_cb*rse+psn_cc*rse**2
286 : 559720 : darse=(psn_ba+two*psn_ca*rse)*drse
287 : 559720 : dbrse=(psn_bb+two*psn_cb*rse)*drse
288 : 559720 : dcrse=(psn_cb+two*psn_cc*rse)*drse
289 : 559720 : invf=arse+brse*rsp+crse*rsp**2+invrhop/exc+invrhoe/exc_p
290 : 559720 : fnxc(ipt)=one/invf
291 : 559720 : dinvf=darse+dbrse*rsp+dcrse*rsp**2-invrhop*dexc/exc**2-invrhoe**2/exc_p
292 : 559720 : vxce(ipt)=-dinvf/invf**2
293 [ - + ]: 559720 : if (need_dvxce) then
294 : 0 : d2arse=darse*d2rse/drse+two*psn_ca*drse**2
295 : 0 : d2brse=dbrse*d2rse/drse+two*psn_cb*drse**2
296 : 0 : d2crse=dcrse*d2rse/drse+two*psn_cc*drse**2
297 : : d2invf=d2arse+d2brse*rsp+d2crse*rsp**2 &
298 : 0 : & +invrhop*(two*dexc**2/exc-d2exc)/exc**2+two*invrhoe**3/exc_p
299 : 0 : dvxce(ipt)=(two*dinvf**2/invf-d2invf)/invf**2
300 : : end if
301 : 559720 : dinvf=(brse+two*crse*rsp)*drsp-invrhop**2/exc-invrhoe*dexc_p/exc_p**2
302 : 559720 : vxcp(ipt)=-dinvf/invf**2
303 [ - + ]: 559720 : if (need_dvxcp) then
304 : : d2invf=two*crse*drsp+(brse+two*crse*rsp)*d2rsp &
305 : 0 : & +two*invrhop**3/exc+invrhoe*(two*dexc_p**2/exc_p-d2exc_p)/exc_p**2
306 : 0 : dvxcp(ipt)=(two*dinvf**2/invf-d2invf)/invf**2
307 : : end if
308 : : ! For small rse, use pure Arponen/Pajanne functional
309 : : ! Around the limit (rse=0.7, see PSN paper), switch smoothly from PSN to AP
310 [ + - + + ]: 559720 : if (rse>=ap_psn_1.and.rse<=ap_psn_2) then
311 : 157052 : ratio_psn=(rse-ap_psn_1)/(ap_psn_2-ap_psn_1);ratio_ap=one-ratio_psn
312 : 157052 : fnxc(ipt)=ratio_psn*fnxc(ipt)+ratio_ap*rhop*exc
313 : 157052 : vxce(ipt)=ratio_psn*vxce(ipt)+ratio_ap*rhop*dexc
314 : 157052 : vxcp(ipt)=ratio_psn*vxcp(ipt)+ratio_ap*exc
315 [ - + ]: 157052 : if (need_dvxce) dvxce(ipt)=ratio_psn*dvxce(ipt)+ratio_ap*rhop*d2exc
316 [ - + ]: 157052 : if (need_dvxcp) dvxcp(ipt)=ratio_psn*dvxcp(ipt)
317 : : end if
318 : : end if
319 : :
320 : : ! Debug statements: use polynomial functionals
321 : 69260 : if (idebug>0) then
322 : : if (idebug==4) then ! order 4
323 : : fnxc(ipt)=tol3*((rhop**4+rhoe**4)/12._dp+(rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe**2)
324 : : vxce(ipt)=tol3*((rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe+rhop*rhoe**2)
325 : : vxcp(ipt)=tol3*((rhop**3*rhoe+rhop*rhoe**3)/3._dp+rhop**2*rhoe+rhop*rhoe**2)
326 : : if (need_dvxce) dvxce(ipt)=tol3*(rhop**3/3._dp+rhop*rhoe**2+rhop**2+two*rhop*rhoe)
327 : : if (need_dvxcp) dvxcp(ipt)=tol3*(rhoe**3/3._dp+rhoe*rhop**2+rhoe**2+two*rhop*rhoe)
328 : : end if
329 : : if (idebug==3) then ! order 3
330 : : fnxc(ipt)=tol3*((rhop**3+rhoe**3)*third+rhop**2*rhoe+rhop*rhoe**2)
331 : : vxce(ipt)=tol3*(rhop+rhoe)**2
332 : : vxcp(ipt)=tol3*(rhop+rhoe)**2
333 : : if (need_dvxce) dvxce(ipt)=tol3*two*rhoe
334 : : if (need_dvxcp) dvxcp(ipt)=tol3*two*rhop
335 : : end if
336 : : if (idebug==2) then ! order 2
337 : : fnxc(ipt)=tol3*(rhop+rhoe)**2
338 : : vxce(ipt)=tol3*two*(rhop+rhoe)
339 : : vxcp(ipt)=tol3*two*(rhop+rhoe)
340 : : if (need_dvxce) dvxce(ipt)=tol3*two
341 : : if (need_dvxcp) dvxcp(ipt)=tol3*two
342 : : end if
343 : : if (idebug==1) then ! order 1
344 : : fnxc(ipt)=tol3*(rhop+rhoe)
345 : : vxce(ipt)=tol3
346 : : vxcp(ipt)=tol3
347 : : if (need_dvxce) dvxce(ipt)=zero
348 : : if (need_dvxcp) dvxcp(ipt)=zero
349 : : end if
350 : : end if
351 : :
352 : : end do ! ipt
353 : :
354 [ - + ]: 69260 : ABI_FREE(rsepts)
355 [ + + ]: 69260 : if (ixcpositron==2) then
356 [ - + ]: 616 : ABI_FREE(rsppts)
357 : : end if
358 : :
359 : : !Convert everything in Hartree units
360 [ + + ]: 147695946 : fnxc(:)=half*fnxc(:)
361 [ + + ]: 147695946 : vxce(:)=half*vxce(:)
362 [ + + ]: 147695946 : vxcp(:)=half*vxcp(:)
363 [ + + + + ]: 45692915 : if (need_dvxce) dvxce(:)=half*dvxce(:)
364 [ - + - - ]: 69260 : if (need_dvxcp) dvxcp(:)=half*dvxcp(:)
365 [ - + - - ]: 69260 : if (gga) vxcegr(:)=half*vxcegr(:)
366 : :
367 [ + + - + : 145947 : end subroutine xcpositron
- + ]
368 : : !!***
369 : :
370 : : end module m_xcpositron
371 : : !!***
|
