Line data Source code
1 : !!****m* ABINIT/m_xchcth
2 : !! NAME
3 : !! m_xchcth
4 : !!
5 : !! FUNCTION
6 : !! Treat XC GGA functional of HCTH type.
7 : !! See Hamprecht, Cohen, Tozer and Handy, J. Chem. Phys. 109, 6264 (1998) [[cite:Hamprecht1998]] for HCTH-93
8 : !! Boese, Doltsinis, Handy and Sprik, J. Chem. Phys. 112, 1670 (2000) [[cite:Boese2000]] for HCTH-120 and HCTH-147.
9 : !! Boese and Handy , J. Chem. Phys. 114, 5497 (2001) [[cite:Boese2001]] for HCTH-407.
10 : !!
11 : !! COPYRIGHT
12 : !! Copyright (C) 2002-2026 ABINIT group (XG,LG)
13 : !! This file is distributed under the terms of the
14 : !! GNU General Public License, see ~abinit/COPYING
15 : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : module m_xchcth
26 :
27 : use defs_basis
28 : use m_errors
29 : use m_abicore
30 :
31 : use m_numeric_tools, only : invcb
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: xchcth
39 : !!***
40 :
41 : contains
42 : !!***
43 :
44 : !!****f* ABINIT/xchcth
45 : !! NAME
46 : !! xchcth
47 : !!
48 : !! FUNCTION
49 : !! Treat XC GGA functional of HCTH type.
50 : !! See Hamprecht, Cohen, Tozer and Handy, J. Chem. Phys. 109, 6264 (1998) [[cite:Hamprecht1998]] for HCTH-93.
51 : !! Boese, Doltsinis, Handy and Sprik, J. Chem. Phys. 112, 1670 (2000) [[cite:Boese2000]] for HCTH-120 and HCTH-147.
52 : !! Boese and Handy , J. Chem. Phys. 114, 5497 (2001) [[cite:Boese2001]] for HCTH-407.
53 : !!
54 : !! For a series of values of the density and the square of the
55 : !! gradient of the density, return the associated Exc energy,
56 : !! potential, and, in case of response-function, functions needed
57 : !! to build the XC kernel.
58 : !!
59 : !! INPUTS
60 : !! ixc=number of the XC functional : 16 for HCTH-93, 17 for HCTH-120, 26 for HCTH-147 and 27 for HCTH-407.
61 : !! npts= number of points to be computed
62 : !! nspden=1 for unpolarized, 2 for spin-polarized
63 : !! grho2_updn(npts,2*nspden-1)=square of the gradient of the spin-up,
64 : !! and, if nspden==2, spin-down, and total density (Hartree/Bohr**2)
65 : !! order=maximal derivative of Exc to be computed
66 : !! (1 => energy and potential, or 2 => also XC kernel )
67 : !! Warning : order=2 not yet available
68 : !! rho_updn(npts,nspden)=spin-up and spin-down density (Hartree/bohr**3)
69 : !!
70 : !! OUTPUT
71 : !!
72 : !! dvxcdgr(npts,3)=partial derivative of the exchange-correlation
73 : !! energy (exci*$\rho$) with respect to the spin-up (dvxcdgr(:,1)),
74 : !! spin-down (dvxcdgr(:,2)) square of gradients of the density
75 : !!
76 : !! exci(npts)=exchange-correlation energy density (hartree)
77 : !! vxci(npts,nspden)=partial derivative of the exchange-correlation energy (exci*$\rho$)
78 : !! with respect to the spin-down (vxci(:,1)) and spin-up (vxci(:,2) densities
79 : !! Normalization: Exc=$\int (exc(r)*\rho (r) d^3 r)$ for $\rho$(r)=electron density.
80 : !!
81 : !! TODO
82 : !! Response function not coded yet, but part of it are already present
83 : !!
84 : !! SOURCE
85 :
86 1265 : subroutine xchcth(dvxcdgr,exci,grho2_updn,ixc,npts,nspden,order,rho_updn,vxci)
87 :
88 : !Arguments ------------------------------------
89 : !scalars
90 : integer,intent(in) :: ixc,npts,nspden,order
91 : !arrays
92 : real(dp),intent(in) :: grho2_updn(npts,2*nspden-1),rho_updn(npts,nspden)
93 : real(dp),intent(out) :: dvxcdgr(npts,2),exci(npts),vxci(npts,nspden)
94 :
95 : !Local variables-------------------------------
96 : !scalars
97 : integer,save :: initialized=0
98 : integer :: ipts,ispden
99 : real(dp),parameter :: alpha_zeta2=one-1e-6_dp,alpha_zeta=one-1e-6_dp
100 : real(dp),parameter :: fsec_inv=one/1.709921_dp,gammacab=0.006_dp
101 : real(dp),parameter :: gammacsig=0.2_dp,gammax=0.004_dp
102 : real(dp),parameter :: rsfac=0.6203504908994000_dp
103 : real(dp),save :: factf_zeta,factfp_zeta,sixpi2_1_3,sixpi2m1_3,sq_rsfac
104 : real(dp),save :: sq_rsfac_inv,threefourth_divpi,twom1_3
105 : real(dp) :: ccab0,ccab1,ccab2,ccab3,ccab4,ccsig0,ccsig1,ccsig2,ccsig3,ccsig4
106 : real(dp) :: coeffss,cxsig0,cxsig1,cxsig2,cxsig3,cxsig4,d2ecrs0_drs2
107 : real(dp) :: d2ecrs1_drs2,d2ecrs_drs2,d2ecrs_drsdzeta,d2ecrs_dzeta2
108 : real(dp) :: d2fzeta4_dzeta2,d2gcrs_drs2,d2macrs_drs2,decrs0_drs,decrs1_drho
109 : real(dp) :: decrs1_drs,decrs_drs,decrs_dzeta,dfzeta4_dzeta,dgcabdss
110 : real(dp) :: dgcrs_drs,dgcsigdss,dgxsigdss,divcab,divcsig,divx,dmacrs_drs
111 : real(dp) :: drhoecab_drhodn,drhoecab_drhoup,drhoecrs1_drhodn,drhoecrs1_drhoup
112 : real(dp) :: drsdrho,dssdg,dssdndg,dssdndrho,dssdrho,dssupdg,dssupdrho,ducabdss
113 : real(dp) :: ducsigdss,duxsigdss,ec0_a1,ec0_aa,ec0_b1,ec0_b2,ec0_b3,ec0_b4
114 : real(dp) :: ec0_den,ec0_log,ec0_q0,ec0_q1,ec0_q1p,ec0_q1pp,ec1_a1,ec1_aa
115 : real(dp) :: ec1_b1,ec1_b2,ec1_b3,ec1_b4,ec1_den,ec1_log,ec1_q0,ec1_q1,ec1_q1p
116 : real(dp) :: ec1_q1pp,ecrs,ecrs0,ecrs1,ex_lsd,exc,f_zeta,factfpp_zeta
117 : real(dp) :: fp_zeta,fpp_zeta,gcab,gcrs,gcsig,gxsig,mac_a1,mac_aa,mac_b1
118 : real(dp) :: mac_b2,mac_b3,mac_b4,mac_den,mac_log,mac_q0,mac_q1,mac_q1p
119 : real(dp) :: mac_q1pp,macrs,rho,rho_inv
120 : real(dp) :: rhoecab,rhoecrs1_dn,rhoecrs1_up,rhomo6,rhomot,rhoo6
121 : real(dp) :: rhotmo6,rhotmot,rhoto6,rhotot,rhotot_inv,rs,rsm1_2,sqr_rs,ss,ss_dn
122 : real(dp) :: ss_up,ssavg,ucab,ucsig,uxsig,vxcadd,zeta,zeta4,zetm_1_3
123 : real(dp) :: zetp_1_3
124 : character(len=500) :: message
125 : !arrays
126 1265 : real(dp),allocatable :: rho_updnm1_3(:,:),rhoarr(:),rhom1_3(:),zetm(:)
127 1265 : real(dp),allocatable :: zetmm1_3(:),zetp(:),zetpm1_3(:)
128 : !no_abirules
129 : !real(dp) :: delta,factor,grr,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean
130 :
131 : ! *************************************************************************
132 :
133 : !DEBUG
134 : !write(std_out,*)' xchcth : enter'
135 : !write(std_out,*)' nspden=',nspden
136 : !ENDDEBUG
137 :
138 1265 : if (order/=1) then
139 0 : write(message, '(a,i0)' )' Order must be 1 ; argument was ',order
140 0 : ABI_BUG(message)
141 : end if
142 :
143 1265 : if(initialized==0)then
144 4 : twom1_3=two**(-third)
145 4 : sixpi2_1_3=(six*pi**2)**third
146 4 : sixpi2m1_3=one/sixpi2_1_3
147 4 : threefourth_divpi=three_quarters*piinv
148 4 : factf_zeta= one / ( two**(four/three)-two )
149 4 : factfp_zeta= four_thirds * factf_zeta * alpha_zeta2
150 4 : sq_rsfac=sqrt(rsfac)
151 4 : sq_rsfac_inv=one/sq_rsfac
152 4 : initialized=1
153 : end if
154 :
155 1265 : if(ixc==16)then
156 : ! HCTH-93 parameters from table II of JCP 109, 6264 (1998). Note that there is
157 : ! an error in tables III of JCP 112, 1670 (2000) and JCP 114, 5497 (2001) for the coefficient ccab1
158 : cxsig0=1.09320_dp ; cxsig1=-0.744056_dp ; cxsig2=5.59920_dp ; cxsig3=-6.78549_dp ; cxsig4=4.49357_dp
159 : ccsig0=0.222601_dp; ccsig1=-0.0338622_dp; ccsig2=-0.0125170_dp; ccsig3=-0.802496_dp; ccsig4=1.55396_dp
160 : ccab0=0.729974_dp ; ccab1=3.35287_dp ; ccab2=-11.543_dp ; ccab3=8.08564_dp ; ccab4=-4.47857_dp
161 : else if(ixc==17)then
162 : ! HCTH-120 parameters from table III of JCP 112, 1670 (2000)
163 : ! Note the correction of the sign of cxsig1 and ccsig1, as there is a misprint in the HCTH paper !
164 : ! see the exchange of mail with MMarques 23 Dec 2008, and 6 Jan 2009.
165 : ! Now, the functional agrees with the libXC results
166 : cxsig0=1.09163_dp; cxsig1=-0.7472_dp; cxsig2=5.0783_dp; cxsig3=-4.1075_dp; cxsig4=1.1717_dp
167 : ccsig0=0.48951_dp; ccsig1=-0.2607_dp; ccsig2=0.4329_dp; ccsig3=-1.9925_dp; ccsig4=2.4853_dp
168 : ccab0=0.51473_dp ; ccab1=6.9298_dp ; ccab2=-24.707_dp; ccab3=23.110_dp ; ccab4=-11.323_dp
169 : ! Exactly same values than in lib_xc
170 : cxsig0=1.09163_dp; cxsig1=-0.747215_dp; cxsig2=5.07833_dp; cxsig3=-4.10746_dp; cxsig4=1.17173_dp
171 : ccsig0=0.489508_dp; ccsig1=-0.260699_dp; ccsig2=0.432917_dp; ccsig3=-1.99247_dp; ccsig4=2.48531_dp
172 : ccab0=0.51473_dp ; ccab1=6.92982_dp ; ccab2=-24.7073_dp; ccab3=23.1098_dp ; ccab4=-11.3234_dp
173 : else if(ixc==26)then
174 : ! HCTH-147 parameters from table III of JCP 112, 1670 (2000)
175 : cxsig0=1.09025_dp; cxsig1=-0.7992_dp; cxsig2=5.5721_dp ; cxsig3=-5.8676_dp; cxsig4=3.0454_dp
176 : ccsig0=0.56258_dp; ccsig1=-0.0171_dp; ccsig2=-1.3064_dp; ccsig3= 1.0575_dp; ccsig4=0.8854_dp
177 : ccab0=0.54235_dp ; ccab1=7.0146_dp ; ccab2=-28.382_dp ; ccab3=35.033_dp ; ccab4=-20.428_dp
178 : ! Exactly same values than in lib_xc
179 : cxsig0=1.09025_dp; cxsig1=-0.799194_dp; cxsig2=5.57212_dp ; cxsig3=-5.86760_dp; cxsig4=3.04544_dp
180 : ccsig0=0.562576_dp; ccsig1= 0.0171436_dp; ccsig2=-1.30636_dp; ccsig3= 1.05747_dp; ccsig4=0.885429_dp
181 : ccab0=0.542352_dp ; ccab1=7.01464_dp ; ccab2=-28.3822_dp ; ccab3=35.0329_dp ; ccab4=-20.4284_dp
182 : else if(ixc==27)then
183 : ! HCTH-407 parameters from table IV of JCP 114, 5497 (2001)
184 : cxsig0=1.08184_dp; cxsig1=-0.5183_dp; cxsig2=3.4256_dp; cxsig3=-2.6290_dp; cxsig4=2.2886_dp
185 : ccsig0=1.18777_dp; ccsig1=-2.4029_dp; ccsig2=5.6174_dp; ccsig3=-9.1792_dp; ccsig4=6.2480_dp
186 : ccab0=0.58908_dp ; ccab1=4.4237_dp ; ccab2=-19.222_dp; ccab3=42.572_dp ; ccab4=-42.005_dp
187 : ! Exactly same values than in lib_xc
188 : cxsig0=1.08184_dp; cxsig1=-0.518339_dp; cxsig2=3.42562_dp; cxsig3=-2.62901_dp; cxsig4=2.28855_dp
189 : ccsig0=1.18777_dp; ccsig1=-2.40292_dp; ccsig2=5.61741_dp; ccsig3=-9.17923_dp; ccsig4=6.24798_dp
190 : ccab0=0.589076_dp ; ccab1=4.42374_dp ; ccab2=-19.2218_dp; ccab3=42.5721_dp ; ccab4=-42.0052_dp
191 : else
192 0 : write(message, '(a,i0)' )' xchcth : ixc must be 16, 17, 26, or 27 ; argument was ',ixc
193 0 : ABI_BUG(message)
194 : end if
195 :
196 : !Parameters for the Perdew-Wang 92 LSD as well as LSD-RPA,
197 : !see Table I of Phys.Rev.B 45,13244 (1992) [[cite:Perdew1992a]]
198 1265 : ec0_aa=0.031091_dp ; ec1_aa=0.015545_dp ; mac_aa=0.016887_dp
199 1265 : ec0_a1=0.21370_dp ; ec1_a1=0.20548_dp ; mac_a1=0.11125_dp
200 1265 : ec0_b1=7.5957_dp ; ec1_b1=14.1189_dp ; mac_b1=10.357_dp
201 1265 : ec0_b2=3.5876_dp ; ec1_b2=6.1977_dp ; mac_b2=3.6231_dp
202 1265 : ec0_b3=1.6382_dp ; ec1_b3=3.3662_dp ; mac_b3=0.88026_dp
203 1265 : ec0_b4=0.49294_dp ; ec1_b4=0.62517_dp ; mac_b4=0.49671_dp
204 :
205 : !DEBUG
206 : !Finite-difference debugging, do not take away
207 : !Note : here work with collinear gradients. Might be generalized ...
208 : !debug=2 ! Choose 1 (rho grads) or 2 (grho grads)
209 : !factor=one
210 : !zeta_mean=0.98_dp
211 : !zeta_mean=zero
212 : !delta=0.000025*factor
213 : !delta=0.0000125*factor
214 : !if(debug/=0)then
215 : !do ipts=1,npts,5
216 : !rho=ipts*0.zero*factor
217 : !rho_up=rho*(one+zeta_mean)*half
218 : !rho_dn=rho*(one-zeta_mean)*half
219 : !rho_upp=rho_up+delta
220 : !rho_upm=rho_up-delta
221 : !rho_dnp=rho_dn+delta
222 : !rho_dnm=rho_dn-delta
223 : !! Here, vary rho
224 : !if(debug==1)then
225 : !rho_updn(ipts ,1)=rho_up ; rho_updn(ipts ,2)=rho_dn
226 : !rho_updn(ipts+1,1)=rho_upp; rho_updn(ipts+1,2)=rho_dn
227 : !rho_updn(ipts+2,1)=rho_upm; rho_updn(ipts+2,2)=rho_dn
228 : !rho_updn(ipts+3,1)=rho_up ; rho_updn(ipts+3,2)=rho_dnp
229 : !rho_updn(ipts+4,1)=rho_up ; rho_updn(ipts+4,2)=rho_dnm
230 : !grho2_updn(ipts:ipts+4,1)=(0.2_dp*factor)**2 ! grad2 of spin up density
231 : !grho2_updn(ipts:ipts+4,2)=(0.2_dp*factor)**2 ! grad2 of spin down density
232 : !grho2_updn(ipts:ipts+4,3)=(0.3_dp*factor)**2 ! grad2 of total density
233 : !else
234 : !! Here, vary grho (interchange rho and grho)
235 : !grho2_updn(ipts ,1)=rho_up**2 ; grho2_updn(ipts ,2)=rho_dn**2
236 : !grho2_updn(ipts+1,1)=rho_upp**2; grho2_updn(ipts+1,2)=rho_dn**2
237 : !grho2_updn(ipts+2,1)=rho_upm**2; grho2_updn(ipts+2,2)=rho_dn**2
238 : !grho2_updn(ipts+3,1)=rho_up**2 ; grho2_updn(ipts+3,2)=rho_dnp**2
239 : !grho2_updn(ipts+4,1)=rho_up**2 ; grho2_updn(ipts+4,2)=rho_dnm**2
240 : !grho2_updn(ipts ,3)=(ipts*0.zero*factor)**2
241 : !grho2_updn(ipts+1,3)=(ipts*0.zero*factor+delta)**2
242 : !grho2_updn(ipts+2,3)=(ipts*0.zero*factor-delta)**2
243 : !grho2_updn(ipts+3,3)=(ipts*0.zero*factor+delta)**2 ! identical to ipts+1
244 : !grho2_updn(ipts+4,3)=(ipts*0.zero*factor-delta)**2 ! identical to ipts+2
245 : !rho_updn(ipts:ipts+4,1)=0.2_dp*factor*(one+zeta_mean)*half ! spin up density
246 : !rho_updn(ipts:ipts+4,2)=0.2_dp*factor*(one-zeta_mean)*half ! spin down density
247 : !end if
248 : !end do
249 : !end if
250 : !Usual option :
251 : !nspden=2 ; order=2
252 : !GGA
253 : !nspden=2 ; order=1
254 : !Might take also, although finite difference later is meaningless
255 : !nspden=1 ; order=-2
256 : !ENDDEBUG
257 :
258 1265 : if(order**2 >1)then
259 1265 : factfpp_zeta= third * factfp_zeta * alpha_zeta2
260 : end if
261 :
262 3795 : ABI_MALLOC(rhoarr,(npts))
263 2530 : ABI_MALLOC(rhom1_3,(npts))
264 3795 : ABI_MALLOC(rho_updnm1_3,(npts,2))
265 2808 : do ispden=1,nspden
266 2808 : call invcb(rho_updn(:,ispden),rho_updnm1_3(:,ispden),npts)
267 : end do
268 1265 : if(nspden==1)then
269 3770715 : rhoarr(:)=two*rho_updn(:,1)
270 3770715 : rhom1_3(:)=twom1_3*rho_updnm1_3(:,1)
271 3770715 : rho_updnm1_3(:,2)=rho_updnm1_3(:,1)
272 : else
273 1032470 : rhoarr(:)=rho_updn(:,1)+rho_updn(:,2)
274 278 : call invcb(rhoarr,rhom1_3,npts)
275 556 : ABI_MALLOC(zetm,(npts))
276 556 : ABI_MALLOC(zetmm1_3,(npts))
277 556 : ABI_MALLOC(zetp,(npts))
278 556 : ABI_MALLOC(zetpm1_3,(npts))
279 1032470 : do ipts=1,npts
280 1032192 : rhotmot=rhom1_3(ipts)
281 1032192 : rhotot_inv=rhotmot*rhotmot*rhotmot
282 1032192 : zeta=(rho_updn(ipts,1)-rho_updn(ipts,2))*rhotot_inv
283 1032192 : zetp(ipts)=one+zeta*alpha_zeta
284 1032470 : zetm(ipts)=one-zeta*alpha_zeta
285 : end do
286 278 : call invcb(zetp,zetpm1_3,npts)
287 278 : call invcb(zetm,zetmm1_3,npts)
288 : end if
289 :
290 1265 : if (nspden==1) then
291 :
292 987 : if(order==-2) then
293 :
294 0 : do ipts=1,npts
295 :
296 : ! -----------------------------------------------------------------------
297 : ! First take care of the spin-split part of the functional
298 0 : exc=zero
299 0 : ispden=1
300 0 : rho =rho_updn(ipts,ispden)
301 0 : rhomot=rho_updnm1_3(ipts,ispden)
302 0 : rho_inv=rhomot*rhomot*rhomot
303 :
304 : ! Exchange part
305 0 : ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho
306 : ! Note that this definition differs from the PBE one
307 0 : coeffss=rho_inv*rho_inv*rhomot*rhomot
308 0 : ss=grho2_updn(ipts,ispden)*coeffss
309 0 : dssdrho=-eight*third*ss*rho_inv
310 0 : dssdg=two*coeffss
311 :
312 0 : divx=one/(one+gammax*ss)
313 0 : uxsig=gammax*ss*divx
314 0 : duxsigdss=gammax*divx*(one-ss*gammax*divx)
315 :
316 0 : gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4)))
317 : dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))&
318 0 : & *duxsigdss
319 :
320 0 : exc=exc+ex_lsd*rho*gxsig
321 0 : vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho)
322 0 : dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg
323 :
324 : ! Spin parallel correlation part
325 : ! Note that this definition is for fully spin-polarized quantities
326 0 : rs=rsfac*rhomot
327 0 : rhomo6=sqrt(rhomot)
328 0 : sqr_rs=sq_rsfac*rhomo6
329 0 : rhoo6=rho*rhomot*rhomot*rhomo6
330 0 : rsm1_2=sq_rsfac_inv*rhoo6
331 0 : drsdrho=-third*rs*rho_inv
332 :
333 0 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
334 0 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
335 0 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
336 0 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
337 : ! ec1_log=log( one + one / ec1_q1 )
338 0 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
339 0 : ecrs1=ec1_q0*ec1_log
340 0 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
341 0 : decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho
342 :
343 : ! Store the LSDA corr energy and density derivative
344 0 : rhoecrs1_up=rho*ecrs1
345 0 : drhoecrs1_drhoup=decrs1_drho
346 0 : ss_up=ss
347 0 : dssupdrho=dssdrho
348 0 : dssupdg=dssdg
349 :
350 0 : divcsig=one/(one+gammacsig*ss)
351 0 : ucsig=gammacsig*ss*divcsig
352 0 : ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig)
353 :
354 0 : gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4)))
355 : dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))&
356 0 : & *ducsigdss
357 :
358 : ! DEBUG
359 : ! gcsig=zero
360 : ! dgcsigdss=zero
361 : ! ENDDEBUG
362 0 : exc=exc+ecrs1*rho*gcsig
363 : vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+&
364 0 : & ecrs1*rho*dgcsigdss*dssdrho
365 0 : dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg
366 :
367 0 : rhoecrs1_dn=rhoecrs1_up
368 0 : drhoecrs1_drhodn=drhoecrs1_drhoup
369 0 : ss_dn=ss_up
370 0 : dssdndrho=dssupdrho
371 0 : dssdndg=dssupdg
372 0 : exc=exc*2
373 :
374 : ! -----------------------------------------------------------------------------
375 : ! Then takes care of the LSD correlation part of the functional
376 :
377 0 : rhotot=rhoarr(ipts)
378 0 : rhotmot=rhom1_3(ipts)
379 0 : rhotot_inv=rhotmot*rhotmot*rhotmot
380 0 : rhotmo6=sqrt(rhotmot)
381 0 : rhoto6=rhotot*rhotmot*rhotmot*rhotmo6
382 :
383 : ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f
384 0 : rs=rsfac*rhotmot
385 0 : sqr_rs=sq_rsfac*rhotmo6
386 0 : rsm1_2=sq_rsfac_inv*rhoto6
387 :
388 : ! Formulas A6-A8 of PW92LSD
389 0 : ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs)
390 0 : ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs)
391 0 : ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs)
392 0 : ec0_den=one/(ec0_q1*ec0_q1+ec0_q1)
393 : ! ec0_log=log( one + one / ec0_q1 )
394 0 : ec0_log=-log( ec0_q1*ec0_q1*ec0_den )
395 0 : ecrs0=ec0_q0*ec0_log
396 0 : decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den
397 :
398 0 : ec0_q1pp=half*ec0_aa*(-ec0_b1*rsm1_2**3+3.d0*ec0_b3*rsm1_2+8.d0*ec0_b4)
399 : d2ecrs0_drs2= 4.0d0*ec0_aa*ec0_a1*ec0_q1p*ec0_den &
400 : & -ec0_q0*ec0_q1pp*ec0_den &
401 0 : & +ec0_q0*ec0_q1p**2*ec0_den**2*(2.d0*ec0_q1+one)
402 :
403 :
404 :
405 0 : mac_q0=-2.0d0*mac_aa*(one+mac_a1*rs)
406 0 : mac_q1=2.0d0*mac_aa*(mac_b1*sqr_rs+mac_b2*rs+mac_b3*rs*sqr_rs+mac_b4*rs*rs)
407 0 : mac_q1p=mac_aa*(mac_b1*rsm1_2+2.d0*mac_b2+3.d0*mac_b3*sqr_rs+4.d0*mac_b4*rs)
408 0 : mac_den=one/(mac_q1*mac_q1+mac_q1)
409 0 : mac_log=-log( mac_q1*mac_q1*mac_den )
410 0 : macrs=mac_q0*mac_log
411 0 : dmacrs_drs= -2.0d0*mac_aa*mac_a1*mac_log - mac_q0*mac_q1p*mac_den
412 :
413 :
414 0 : ecrs=ecrs0
415 0 : decrs_drs=decrs0_drs
416 0 : decrs_dzeta=zero
417 :
418 0 : d2ecrs_drs2=d2ecrs0_drs2
419 :
420 0 : d2ecrs_dzeta2=alpha_zeta**2*(-macrs)
421 :
422 0 : zeta=zero
423 :
424 : ! At this point, the coding of the PW92 functional finishes.
425 :
426 : ! The correlation between different spin from HCTH is now computed
427 : ! First, the part without gradient correction factor
428 0 : rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up
429 0 : vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta
430 0 : drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup
431 0 : drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn
432 :
433 : ! Now, the gradient correction factor
434 0 : ssavg=half*(ss_up+ss_dn)
435 0 : divcab=one/(one+gammacab*ssavg)
436 0 : ucab=gammacab*ssavg*divcab
437 0 : ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab)
438 :
439 0 : gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4)))
440 : dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))&
441 0 : & *ducabdss
442 :
443 0 : exc=exc+rhoecab*gcab
444 :
445 0 : vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho
446 0 : dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg
447 : ! If non spin-polarized, treat spin down contribution now, similar to spin up
448 : ! vxci(ipts,2)=vxci(ipts,1)
449 0 : dvxcdgr(ipts,2)=dvxcdgr(ipts,1)
450 :
451 : ! Final division by the total density, to give the energy density
452 0 : exci(ipts)=exc*rhotot_inv
453 :
454 : end do ! ipts=1,npts
455 :
456 987 : else if(order**2>1) then
457 :
458 0 : do ipts=1,npts
459 :
460 : ! -----------------------------------------------------------------------
461 : ! First take care of the spin-split part of the functional
462 0 : exc=zero
463 0 : ispden=1
464 0 : rho =rho_updn(ipts,ispden)
465 0 : rhomot=rho_updnm1_3(ipts,ispden)
466 0 : rho_inv=rhomot*rhomot*rhomot
467 :
468 : ! Exchange part
469 0 : ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho
470 : ! Note that this definition differs from the PBE one
471 0 : coeffss=rho_inv*rho_inv*rhomot*rhomot
472 0 : ss=grho2_updn(ipts,ispden)*coeffss
473 0 : dssdrho=-eight*third*ss*rho_inv
474 0 : dssdg=two*coeffss
475 :
476 0 : divx=one/(one+gammax*ss)
477 0 : uxsig=gammax*ss*divx
478 0 : duxsigdss=gammax*divx*(one-ss*gammax*divx)
479 :
480 0 : gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4)))
481 : dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))&
482 0 : & *duxsigdss
483 :
484 0 : exc=exc+ex_lsd*rho*gxsig
485 0 : vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho)
486 0 : dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg
487 :
488 : ! Spin parallel correlation part
489 : ! Note that this definition is for fully spin-polarized quantities
490 0 : rs=rsfac*rhomot
491 0 : rhomo6=sqrt(rhomot)
492 0 : sqr_rs=sq_rsfac*rhomo6
493 0 : rhoo6=rho*rhomot*rhomot*rhomo6
494 0 : rsm1_2=sq_rsfac_inv*rhoo6
495 0 : drsdrho=-third*rs*rho_inv
496 :
497 0 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
498 0 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
499 0 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
500 0 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
501 : ! ec1_log=log( one + one / ec1_q1 )
502 0 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
503 0 : ecrs1=ec1_q0*ec1_log
504 0 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
505 0 : decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho
506 :
507 : ! Store the LSDA corr energy and density derivative
508 0 : rhoecrs1_up=rho*ecrs1
509 0 : drhoecrs1_drhoup=decrs1_drho
510 0 : ss_up=ss
511 0 : dssupdrho=dssdrho
512 0 : dssupdg=dssdg
513 :
514 0 : divcsig=one/(one+gammacsig*ss)
515 0 : ucsig=gammacsig*ss*divcsig
516 0 : ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig)
517 :
518 0 : gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4)))
519 : dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))&
520 0 : & *ducsigdss
521 :
522 : ! DEBUG
523 : ! gcsig=zero
524 : ! dgcsigdss=zero
525 : ! ENDDEBUG
526 0 : exc=exc+ecrs1*rho*gcsig
527 : vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+&
528 0 : & ecrs1*rho*dgcsigdss*dssdrho
529 0 : dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg
530 :
531 0 : rhoecrs1_dn=rhoecrs1_up
532 0 : drhoecrs1_drhodn=drhoecrs1_drhoup
533 0 : ss_dn=ss_up
534 0 : dssdndrho=dssupdrho
535 0 : dssdndg=dssupdg
536 0 : exc=exc*2
537 :
538 : ! -----------------------------------------------------------------------------
539 : ! Then takes care of the LSD correlation part of the functional
540 :
541 0 : rhotot=rhoarr(ipts)
542 0 : rhotmot=rhom1_3(ipts)
543 0 : rhotot_inv=rhotmot*rhotmot*rhotmot
544 0 : rhotmo6=sqrt(rhotmot)
545 0 : rhoto6=rhotot*rhotmot*rhotmot*rhotmo6
546 :
547 : ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f
548 0 : rs=rsfac*rhotmot
549 0 : sqr_rs=sq_rsfac*rhotmo6
550 0 : rsm1_2=sq_rsfac_inv*rhoto6
551 :
552 : ! Formulas A6-A8 of PW92LSD
553 0 : ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs)
554 0 : ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs)
555 0 : ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs)
556 0 : ec0_den=one/(ec0_q1*ec0_q1+ec0_q1)
557 : ! ec0_log=log( one + one / ec0_q1 )
558 0 : ec0_log=-log( ec0_q1*ec0_q1*ec0_den )
559 0 : ecrs0=ec0_q0*ec0_log
560 0 : decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den
561 0 : ec0_q1pp=half*ec0_aa*(-ec0_b1*rsm1_2**3+3.d0*ec0_b3*rsm1_2+8.d0*ec0_b4)
562 : d2ecrs0_drs2= 4.0d0*ec0_aa*ec0_a1*ec0_q1p*ec0_den &
563 : & -ec0_q0*ec0_q1pp*ec0_den &
564 0 : & +ec0_q0*ec0_q1p**2*ec0_den**2*(2.d0*ec0_q1+one)
565 :
566 :
567 0 : ecrs=ecrs0
568 0 : decrs_drs=decrs0_drs
569 0 : decrs_dzeta=zero
570 0 : d2ecrs_drs2=d2ecrs0_drs2
571 0 : zeta=zero
572 :
573 : ! At this point, the coding of the PW92 functional finishes.
574 :
575 : ! The correlation between different spin from HCTH is now computed
576 : ! First, the part without gradient correction factor
577 0 : rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up
578 0 : vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta
579 0 : drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup
580 0 : drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn
581 :
582 : ! Now, the gradient correction factor
583 0 : ssavg=half*(ss_up+ss_dn)
584 0 : divcab=one/(one+gammacab*ssavg)
585 0 : ucab=gammacab*ssavg*divcab
586 0 : ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab)
587 :
588 0 : gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4)))
589 : dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))&
590 0 : & *ducabdss
591 :
592 0 : exc=exc+rhoecab*gcab
593 :
594 0 : vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho
595 0 : dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg
596 : ! If non spin-polarized, treat spin down contribution now, similar to spin up
597 : ! vxci(ipts,2)=vxci(ipts,1)
598 0 : dvxcdgr(ipts,2)=dvxcdgr(ipts,1)
599 :
600 : ! Final division by the total density, to give the energy density
601 0 : exci(ipts)=exc*rhotot_inv
602 :
603 : end do ! ipts=1,npts
604 : else
605 :
606 3770715 : do ipts=1,npts
607 :
608 : ! -----------------------------------------------------------------------
609 : ! First take care of the spin-split part of the functional
610 3769728 : exc=zero
611 3769728 : ispden=1
612 3769728 : rho =rho_updn(ipts,ispden)
613 3769728 : rhomot=rho_updnm1_3(ipts,ispden)
614 3769728 : rho_inv=rhomot*rhomot*rhomot
615 :
616 : ! Exchange part
617 3769728 : ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho
618 : ! Note that this definition differs from the PBE one
619 3769728 : coeffss=rho_inv*rho_inv*rhomot*rhomot
620 3769728 : ss=grho2_updn(ipts,ispden)*coeffss
621 3769728 : dssdrho=-eight*third*ss*rho_inv
622 3769728 : dssdg=two*coeffss
623 :
624 3769728 : divx=one/(one+gammax*ss)
625 3769728 : uxsig=gammax*ss*divx
626 3769728 : duxsigdss=gammax*divx*(one-ss*gammax*divx)
627 :
628 3769728 : gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4)))
629 : dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))&
630 3769728 : & *duxsigdss
631 :
632 3769728 : exc=exc+ex_lsd*rho*gxsig
633 3769728 : vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho)
634 3769728 : dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg
635 :
636 : ! Spin parallel correlation part
637 : ! Note that this definition is for fully spin-polarized quantities
638 3769728 : rs=rsfac*rhomot
639 3769728 : rhomo6=sqrt(rhomot)
640 3769728 : sqr_rs=sq_rsfac*rhomo6
641 3769728 : rhoo6=rho*rhomot*rhomot*rhomo6
642 3769728 : rsm1_2=sq_rsfac_inv*rhoo6
643 3769728 : drsdrho=-third*rs*rho_inv
644 :
645 3769728 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
646 3769728 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
647 3769728 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
648 3769728 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
649 : ! ec1_log=log( one + one / ec1_q1 )
650 3769728 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
651 3769728 : ecrs1=ec1_q0*ec1_log
652 3769728 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
653 3769728 : decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho
654 :
655 : ! Store the LSDA corr energy and density derivative
656 3769728 : rhoecrs1_up=rho*ecrs1
657 3769728 : drhoecrs1_drhoup=decrs1_drho
658 3769728 : ss_up=ss
659 3769728 : dssupdrho=dssdrho
660 3769728 : dssupdg=dssdg
661 :
662 3769728 : divcsig=one/(one+gammacsig*ss)
663 3769728 : ucsig=gammacsig*ss*divcsig
664 3769728 : ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig)
665 :
666 3769728 : gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4)))
667 : dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))&
668 3769728 : & *ducsigdss
669 :
670 : ! DEBUG
671 : ! gcsig=zero
672 : ! dgcsigdss=zero
673 : ! ENDDEBUG
674 3769728 : exc=exc+ecrs1*rho*gcsig
675 : vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+&
676 3769728 : & ecrs1*rho*dgcsigdss*dssdrho
677 3769728 : dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg
678 :
679 3769728 : rhoecrs1_dn=rhoecrs1_up
680 3769728 : drhoecrs1_drhodn=drhoecrs1_drhoup
681 3769728 : ss_dn=ss_up
682 3769728 : dssdndrho=dssupdrho
683 3769728 : dssdndg=dssupdg
684 3769728 : exc=exc*2
685 :
686 : ! -----------------------------------------------------------------------------
687 : ! Then takes care of the LSD correlation part of the functional
688 :
689 3769728 : rhotot=rhoarr(ipts)
690 3769728 : rhotmot=rhom1_3(ipts)
691 3769728 : rhotot_inv=rhotmot*rhotmot*rhotmot
692 3769728 : rhotmo6=sqrt(rhotmot)
693 3769728 : rhoto6=rhotot*rhotmot*rhotmot*rhotmo6
694 :
695 : ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f
696 3769728 : rs=rsfac*rhotmot
697 3769728 : sqr_rs=sq_rsfac*rhotmo6
698 3769728 : rsm1_2=sq_rsfac_inv*rhoto6
699 :
700 : ! Formulas A6-A8 of PW92LSD
701 3769728 : ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs)
702 3769728 : ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs)
703 3769728 : ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs)
704 3769728 : ec0_den=one/(ec0_q1*ec0_q1+ec0_q1)
705 : ! ec0_log=log( one + one / ec0_q1 )
706 3769728 : ec0_log=-log( ec0_q1*ec0_q1*ec0_den )
707 3769728 : ecrs0=ec0_q0*ec0_log
708 3769728 : decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den
709 :
710 3769728 : ecrs=ecrs0
711 3769728 : decrs_drs=decrs0_drs
712 3769728 : decrs_dzeta=zero
713 3769728 : zeta=zero
714 :
715 : ! At this point, the coding of the PW92 functional finishes.
716 :
717 : ! The correlation between different spin from HCTH is now computed
718 : ! First, the part without gradient correction factor
719 3769728 : rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up
720 3769728 : vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta
721 3769728 : drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup
722 3769728 : drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn
723 :
724 : ! Now, the gradient correction factor
725 3769728 : ssavg=half*(ss_up+ss_dn)
726 3769728 : divcab=one/(one+gammacab*ssavg)
727 3769728 : ucab=gammacab*ssavg*divcab
728 3769728 : ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab)
729 :
730 3769728 : gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4)))
731 : dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))&
732 3769728 : & *ducabdss
733 :
734 3769728 : exc=exc+rhoecab*gcab
735 :
736 3769728 : vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho
737 3769728 : dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg
738 : ! If non spin-polarized, treat spin down contribution now, similar to spin up
739 : ! vxci(ipts,2)=vxci(ipts,1)
740 3769728 : dvxcdgr(ipts,2)=dvxcdgr(ipts,1)
741 :
742 : ! Final division by the total density, to give the energy density
743 3770715 : exci(ipts)=exc*rhotot_inv
744 :
745 : end do ! ipts=1,npts
746 : end if
747 :
748 278 : else if(nspden==2) then
749 :
750 278 : if(order**2>1) then
751 :
752 0 : do ipts=1,npts
753 :
754 : ! -----------------------------------------------------------------------
755 : ! First take care of the spin-split part of the functional
756 : exc=zero
757 0 : do ispden=1,nspden
758 0 : rho =rho_updn(ipts,ispden)
759 0 : rhomot=rho_updnm1_3(ipts,ispden)
760 0 : rho_inv=rhomot*rhomot*rhomot
761 :
762 : ! Exchange part
763 0 : ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho
764 : ! Note that this definition differs from the PBE one
765 0 : coeffss=rho_inv*rho_inv*rhomot*rhomot
766 0 : ss=grho2_updn(ipts,ispden)*coeffss
767 0 : dssdrho=-eight*third*ss*rho_inv
768 0 : dssdg=two*coeffss
769 :
770 0 : divx=one/(one+gammax*ss)
771 0 : uxsig=gammax*ss*divx
772 0 : duxsigdss=gammax*divx*(one-ss*gammax*divx)
773 :
774 0 : gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4)))
775 : dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))&
776 0 : & *duxsigdss
777 :
778 0 : exc=exc+ex_lsd*rho*gxsig
779 0 : vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho)
780 0 : dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg
781 :
782 : ! Spin parallel correlation part
783 : ! Note that this definition is for fully spin-polarized quantities
784 0 : rs=rsfac*rhomot
785 0 : rhomo6=sqrt(rhomot)
786 0 : sqr_rs=sq_rsfac*rhomo6
787 0 : rhoo6=rho*rhomot*rhomot*rhomo6
788 0 : rsm1_2=sq_rsfac_inv*rhoo6
789 0 : drsdrho=-third*rs*rho_inv
790 :
791 0 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
792 0 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
793 0 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
794 0 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
795 : ! ec1_log=log( one + one / ec1_q1 )
796 0 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
797 0 : ecrs1=ec1_q0*ec1_log
798 0 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
799 0 : decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho
800 :
801 : ! Store the LSDA corr energy and density derivative
802 0 : if(ispden==1)then
803 0 : rhoecrs1_up=rho*ecrs1
804 0 : drhoecrs1_drhoup=decrs1_drho
805 0 : ss_up=ss
806 0 : dssupdrho=dssdrho
807 0 : dssupdg=dssdg
808 : else
809 0 : rhoecrs1_dn=rho*ecrs1
810 0 : drhoecrs1_drhodn=decrs1_drho
811 0 : ss_dn=ss
812 0 : dssdndrho=dssdrho
813 0 : dssdndg=dssdg
814 : end if
815 :
816 0 : divcsig=one/(one+gammacsig*ss)
817 0 : ucsig=gammacsig*ss*divcsig
818 0 : ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig)
819 :
820 0 : gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4)))
821 : dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))&
822 0 : & *ducsigdss
823 :
824 : ! DEBUG
825 : ! gcsig=zero
826 : ! dgcsigdss=zero
827 : ! ENDDEBUG
828 0 : exc=exc+ecrs1*rho*gcsig
829 : vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+&
830 0 : & ecrs1*rho*dgcsigdss*dssdrho
831 0 : dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg
832 :
833 : end do
834 :
835 : ! -----------------------------------------------------------------------------
836 : ! Then takes care of the LSD correlation part of the functional
837 :
838 0 : rhotot=rhoarr(ipts)
839 0 : rhotmot=rhom1_3(ipts)
840 0 : rhotot_inv=rhotmot*rhotmot*rhotmot
841 0 : rhotmo6=sqrt(rhotmot)
842 0 : rhoto6=rhotot*rhotmot*rhotmot*rhotmo6
843 :
844 : ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f
845 0 : rs=rsfac*rhotmot
846 0 : sqr_rs=sq_rsfac*rhotmo6
847 0 : rsm1_2=sq_rsfac_inv*rhoto6
848 :
849 : ! Formulas A6-A8 of PW92LSD
850 0 : ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs)
851 0 : ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs)
852 0 : ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs)
853 0 : ec0_den=one/(ec0_q1*ec0_q1+ec0_q1)
854 : ! ec0_log=log( one + one / ec0_q1 )
855 0 : ec0_log=-log( ec0_q1*ec0_q1*ec0_den )
856 0 : ecrs0=ec0_q0*ec0_log
857 0 : decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den
858 0 : ec0_q1pp=half*ec0_aa*(-ec0_b1*rsm1_2**3+3.d0*ec0_b3*rsm1_2+8.d0*ec0_b4)
859 : d2ecrs0_drs2= 4.0d0*ec0_aa*ec0_a1*ec0_q1p*ec0_den &
860 : & -ec0_q0*ec0_q1pp*ec0_den &
861 0 : & +ec0_q0*ec0_q1p**2*ec0_den**2*(2.d0*ec0_q1+one)
862 :
863 0 : mac_q0=-2.0d0*mac_aa*(one+mac_a1*rs)
864 0 : mac_q1=2.0d0*mac_aa*(mac_b1*sqr_rs+mac_b2*rs+mac_b3*rs*sqr_rs+mac_b4*rs*rs)
865 0 : mac_q1p=mac_aa*(mac_b1*rsm1_2+2.d0*mac_b2+3.d0*mac_b3*sqr_rs+4.d0*mac_b4*rs)
866 0 : mac_den=one/(mac_q1*mac_q1+mac_q1)
867 0 : mac_log=-log( mac_q1*mac_q1*mac_den )
868 0 : macrs=mac_q0*mac_log
869 0 : dmacrs_drs= -2.0d0*mac_aa*mac_a1*mac_log - mac_q0*mac_q1p*mac_den
870 :
871 0 : zeta=(rho_updn(ipts,1)-rho_updn(ipts,2))*rhotot_inv
872 0 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
873 0 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
874 0 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
875 0 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
876 : ! ec1_log=log( one + one / ec1_q1 )
877 0 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
878 0 : ecrs1=ec1_q0*ec1_log
879 0 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
880 :
881 : ! alpha_zeta is introduced in order to remove singularities for fully
882 : ! polarized systems.
883 0 : zetp_1_3=(one+zeta*alpha_zeta)*zetpm1_3(ipts)**2
884 0 : zetm_1_3=(one-zeta*alpha_zeta)*zetmm1_3(ipts)**2
885 :
886 : f_zeta=( (one+zeta*alpha_zeta2)*zetp_1_3 + &
887 0 : & (one-zeta*alpha_zeta2)*zetm_1_3 - 2.0d0 ) * factf_zeta
888 0 : fp_zeta=( zetp_1_3 - zetm_1_3 ) * factfp_zeta
889 0 : zeta4=zeta**4
890 0 : gcrs=ecrs1-ecrs0+macrs*fsec_inv
891 : ! ecrs=ecrs0+f_zeta*(-macrs*(one-zeta4)*fsec_inv+(ecrs1-ecrs0)*zeta4)
892 0 : ecrs=ecrs0+f_zeta*(zeta4*gcrs-macrs*fsec_inv)
893 :
894 0 : dgcrs_drs=decrs1_drs-decrs0_drs+dmacrs_drs*fsec_inv
895 : ! decrs_drs=decrs0_drs+f_zeta*&
896 : ! & (-dmacrs_drs*(one-zeta4)*fsec_inv+(decrs1_drs-decrs0_drs)*zeta4)
897 0 : decrs_drs=decrs0_drs+f_zeta*(zeta4*dgcrs_drs-dmacrs_drs*fsec_inv)
898 0 : dfzeta4_dzeta=4.0d0*zeta**3*f_zeta+fp_zeta*zeta4
899 0 : decrs_dzeta=dfzeta4_dzeta*gcrs-fp_zeta*macrs*fsec_inv
900 :
901 0 : ec1_q1pp=half*ec1_aa*(-ec1_b1*rsm1_2**3+3.d0*ec1_b3*rsm1_2+8.d0*ec1_b4)
902 : d2ecrs1_drs2= 4.0d0*ec1_aa*ec1_a1*ec1_q1p*ec1_den &
903 : & -ec1_q0*ec1_q1pp*ec1_den &
904 0 : & +ec1_q0*ec1_q1p**2*ec1_den**2*(2.d0*ec1_q1+one)
905 :
906 0 : mac_q1pp=half*mac_aa*(-mac_b1*rsm1_2**3+3.d0*mac_b3*rsm1_2+8.d0*mac_b4)
907 : d2macrs_drs2= 4.0d0*mac_aa*mac_a1*mac_q1p*mac_den &
908 : & -mac_q0*mac_q1pp*mac_den &
909 0 : & +mac_q0*mac_q1p**2*mac_den**2*(2.d0*mac_q1+one)
910 :
911 0 : d2gcrs_drs2=d2ecrs1_drs2-d2ecrs0_drs2+d2macrs_drs2*fsec_inv
912 0 : fpp_zeta=(zetpm1_3(ipts)**2+zetmm1_3(ipts)**2) * factfpp_zeta
913 : d2fzeta4_dzeta2=12.0d0*zeta**2*f_zeta &
914 : & + 8.0d0*zeta**3*fp_zeta &
915 0 : & + zeta4 *fpp_zeta
916 :
917 : d2ecrs_drs2=d2ecrs0_drs2+&
918 0 : & f_zeta*(zeta4*d2gcrs_drs2-d2macrs_drs2*fsec_inv)
919 0 : d2ecrs_drsdzeta=dfzeta4_dzeta*dgcrs_drs-fp_zeta*dmacrs_drs*fsec_inv
920 0 : d2ecrs_dzeta2=d2fzeta4_dzeta2*gcrs-fpp_zeta*macrs*fsec_inv
921 :
922 :
923 : ! At this point, the coding of the PW92 functional finishes.
924 :
925 : ! The correlation between different spin from HCTH is now computed
926 : ! First, the part without gradient correction factor
927 0 : rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up
928 0 : vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta
929 0 : drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup
930 0 : drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn
931 :
932 : ! Now, the gradient correction factor
933 0 : ssavg=half*(ss_up+ss_dn)
934 0 : divcab=one/(one+gammacab*ssavg)
935 0 : ucab=gammacab*ssavg*divcab
936 0 : ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab)
937 :
938 0 : gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4)))
939 : dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))&
940 0 : & *ducabdss
941 :
942 0 : exc=exc+rhoecab*gcab
943 :
944 0 : vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho
945 0 : dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg
946 0 : vxci(ipts,2)=vxci(ipts,2)+drhoecab_drhodn*gcab+rhoecab*dgcabdss*half*dssdndrho
947 0 : dvxcdgr(ipts,2)=dvxcdgr(ipts,2)+rhoecab*dgcabdss*half*dssdndg
948 :
949 : ! Final division by the total density, to give the energy density
950 0 : exci(ipts)=exc*rhotot_inv
951 :
952 : end do ! ipts=1,npts
953 :
954 : else
955 :
956 1032470 : do ipts=1,npts
957 :
958 : ! -----------------------------------------------------------------------
959 : ! First take care of the spin-split part of the functional
960 : exc=zero
961 3096576 : do ispden=1,nspden
962 2064384 : rho =rho_updn(ipts,ispden)
963 2064384 : rhomot=rho_updnm1_3(ipts,ispden)
964 2064384 : rho_inv=rhomot*rhomot*rhomot
965 :
966 : ! Exchange part
967 2064384 : ex_lsd= - threefourth_divpi * sixpi2_1_3*rhomot*rhomot*rho
968 : ! Note that this definition differs from the PBE one
969 2064384 : coeffss=rho_inv*rho_inv*rhomot*rhomot
970 2064384 : ss=grho2_updn(ipts,ispden)*coeffss
971 2064384 : dssdrho=-eight*third*ss*rho_inv
972 2064384 : dssdg=two*coeffss
973 :
974 2064384 : divx=one/(one+gammax*ss)
975 2064384 : uxsig=gammax*ss*divx
976 2064384 : duxsigdss=gammax*divx*(one-ss*gammax*divx)
977 :
978 2064384 : gxsig=cxsig0+uxsig*(cxsig1+uxsig*(cxsig2+uxsig*(cxsig3+uxsig*cxsig4)))
979 : dgxsigdss=(cxsig1+uxsig*(two*cxsig2+uxsig*(three*cxsig3+uxsig*four*cxsig4)))&
980 2064384 : & *duxsigdss
981 :
982 2064384 : exc=exc+ex_lsd*rho*gxsig
983 2064384 : vxci(ipts,ispden)=ex_lsd*(four_thirds*gxsig+rho*dgxsigdss*dssdrho)
984 2064384 : dvxcdgr(ipts,ispden)=ex_lsd*rho*dgxsigdss*dssdg
985 :
986 : ! Spin parallel correlation part
987 : ! Note that this definition is for fully spin-polarized quantities
988 2064384 : rs=rsfac*rhomot
989 2064384 : rhomo6=sqrt(rhomot)
990 2064384 : sqr_rs=sq_rsfac*rhomo6
991 2064384 : rhoo6=rho*rhomot*rhomot*rhomo6
992 2064384 : rsm1_2=sq_rsfac_inv*rhoo6
993 2064384 : drsdrho=-third*rs*rho_inv
994 :
995 2064384 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
996 2064384 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
997 2064384 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
998 2064384 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
999 : ! ec1_log=log( one + one / ec1_q1 )
1000 2064384 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
1001 2064384 : ecrs1=ec1_q0*ec1_log
1002 2064384 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
1003 2064384 : decrs1_drho=ecrs1+decrs1_drs*drsdrho*rho
1004 :
1005 : ! Store the LSDA corr energy and density derivative
1006 2064384 : if(ispden==1)then
1007 1032192 : rhoecrs1_up=rho*ecrs1
1008 1032192 : drhoecrs1_drhoup=decrs1_drho
1009 1032192 : ss_up=ss
1010 1032192 : dssupdrho=dssdrho
1011 1032192 : dssupdg=dssdg
1012 : else
1013 1032192 : rhoecrs1_dn=rho*ecrs1
1014 1032192 : drhoecrs1_drhodn=decrs1_drho
1015 1032192 : ss_dn=ss
1016 1032192 : dssdndrho=dssdrho
1017 1032192 : dssdndg=dssdg
1018 : end if
1019 :
1020 2064384 : divcsig=one/(one+gammacsig*ss)
1021 2064384 : ucsig=gammacsig*ss*divcsig
1022 2064384 : ducsigdss=gammacsig*divcsig*(one-ss*gammacsig*divcsig)
1023 :
1024 2064384 : gcsig=ccsig0+ucsig*(ccsig1+ucsig*(ccsig2+ucsig*(ccsig3+ucsig*ccsig4)))
1025 : dgcsigdss=(ccsig1+ucsig*(two*ccsig2+ucsig*(three*ccsig3+ucsig*four*ccsig4)))&
1026 2064384 : & *ducsigdss
1027 :
1028 : ! DEBUG
1029 : ! gcsig=zero
1030 : ! dgcsigdss=zero
1031 : ! ENDDEBUG
1032 2064384 : exc=exc+ecrs1*rho*gcsig
1033 : vxci(ipts,ispden)=vxci(ipts,ispden)+decrs1_drho*gcsig+&
1034 2064384 : & ecrs1*rho*dgcsigdss*dssdrho
1035 3096576 : dvxcdgr(ipts,ispden)=dvxcdgr(ipts,ispden)+ecrs1*rho*dgcsigdss*dssdg
1036 :
1037 : end do
1038 :
1039 : ! -----------------------------------------------------------------------------
1040 : ! Then takes care of the LSD correlation part of the functional
1041 :
1042 1032192 : rhotot=rhoarr(ipts)
1043 1032192 : rhotmot=rhom1_3(ipts)
1044 1032192 : rhotot_inv=rhotmot*rhotmot*rhotmot
1045 1032192 : rhotmo6=sqrt(rhotmot)
1046 1032192 : rhoto6=rhotot*rhotmot*rhotmot*rhotmo6
1047 :
1048 : ! From now, the coding of the PW92 functional start. It is identical in xcpbe.f
1049 1032192 : rs=rsfac*rhotmot
1050 1032192 : sqr_rs=sq_rsfac*rhotmo6
1051 1032192 : rsm1_2=sq_rsfac_inv*rhoto6
1052 :
1053 : ! Formulas A6-A8 of PW92LSD
1054 1032192 : ec0_q0=-2.0d0*ec0_aa*(one+ec0_a1*rs)
1055 1032192 : ec0_q1=2.0d0*ec0_aa*(ec0_b1*sqr_rs+ec0_b2*rs+ec0_b3*rs*sqr_rs+ec0_b4*rs*rs)
1056 1032192 : ec0_q1p=ec0_aa*(ec0_b1*rsm1_2+2.d0*ec0_b2+3.d0*ec0_b3*sqr_rs+4.d0*ec0_b4*rs)
1057 1032192 : ec0_den=one/(ec0_q1*ec0_q1+ec0_q1)
1058 : ! ec0_log=log( one + one / ec0_q1 )
1059 1032192 : ec0_log=-log( ec0_q1*ec0_q1*ec0_den )
1060 1032192 : ecrs0=ec0_q0*ec0_log
1061 1032192 : decrs0_drs= -2.0d0*ec0_aa*ec0_a1*ec0_log - ec0_q0*ec0_q1p *ec0_den
1062 :
1063 1032192 : mac_q0=-2.0d0*mac_aa*(one+mac_a1*rs)
1064 1032192 : mac_q1=2.0d0*mac_aa*(mac_b1*sqr_rs+mac_b2*rs+mac_b3*rs*sqr_rs+mac_b4*rs*rs)
1065 1032192 : mac_q1p=mac_aa*(mac_b1*rsm1_2+2.d0*mac_b2+3.d0*mac_b3*sqr_rs+4.d0*mac_b4*rs)
1066 1032192 : mac_den=one/(mac_q1*mac_q1+mac_q1)
1067 1032192 : mac_log=-log( mac_q1*mac_q1*mac_den )
1068 1032192 : macrs=mac_q0*mac_log
1069 1032192 : dmacrs_drs= -2.0d0*mac_aa*mac_a1*mac_log - mac_q0*mac_q1p*mac_den
1070 :
1071 1032192 : zeta=(rho_updn(ipts,1)-rho_updn(ipts,2))*rhotot_inv
1072 1032192 : ec1_q0=-two*ec1_aa*(one+ec1_a1*rs)
1073 1032192 : ec1_q1=two*ec1_aa*(ec1_b1*sqr_rs+ec1_b2*rs+ec1_b3*rs*sqr_rs+ec1_b4*rs*rs)
1074 1032192 : ec1_q1p=ec1_aa*(ec1_b1*rsm1_2+two*ec1_b2+three*ec1_b3*sqr_rs+four*ec1_b4*rs)
1075 1032192 : ec1_den=one/(ec1_q1*ec1_q1+ec1_q1)
1076 : ! ec1_log=log( one + one / ec1_q1 )
1077 1032192 : ec1_log=-log( ec1_q1*ec1_q1*ec1_den )
1078 1032192 : ecrs1=ec1_q0*ec1_log
1079 1032192 : decrs1_drs= -two*ec1_aa*ec1_a1*ec1_log - ec1_q0*ec1_q1p *ec1_den
1080 :
1081 : ! alpha_zeta is introduced in order to remove singularities for fully
1082 : ! polarized systems.
1083 1032192 : zetp_1_3=(one+zeta*alpha_zeta)*zetpm1_3(ipts)**2
1084 1032192 : zetm_1_3=(one-zeta*alpha_zeta)*zetmm1_3(ipts)**2
1085 :
1086 : f_zeta=( (one+zeta*alpha_zeta2)*zetp_1_3 + &
1087 1032192 : & (one-zeta*alpha_zeta2)*zetm_1_3 - 2.0d0 ) * factf_zeta
1088 1032192 : fp_zeta=( zetp_1_3 - zetm_1_3 ) * factfp_zeta
1089 1032192 : zeta4=zeta**4
1090 1032192 : gcrs=ecrs1-ecrs0+macrs*fsec_inv
1091 : ! ecrs=ecrs0+f_zeta*(-macrs*(one-zeta4)*fsec_inv+(ecrs1-ecrs0)*zeta4)
1092 1032192 : ecrs=ecrs0+f_zeta*(zeta4*gcrs-macrs*fsec_inv)
1093 :
1094 1032192 : dgcrs_drs=decrs1_drs-decrs0_drs+dmacrs_drs*fsec_inv
1095 : ! decrs_drs=decrs0_drs+f_zeta*&
1096 : ! & (-dmacrs_drs*(one-zeta4)*fsec_inv+(decrs1_drs-decrs0_drs)*zeta4)
1097 1032192 : decrs_drs=decrs0_drs+f_zeta*(zeta4*dgcrs_drs-dmacrs_drs*fsec_inv)
1098 1032192 : dfzeta4_dzeta=4.0d0*zeta**3*f_zeta+fp_zeta*zeta4
1099 1032192 : decrs_dzeta=dfzeta4_dzeta*gcrs-fp_zeta*macrs*fsec_inv
1100 :
1101 : ! At this point, the coding of the PW92 functional finishes.
1102 :
1103 : ! The correlation between different spin from HCTH is now computed
1104 : ! First, the part without gradient correction factor
1105 1032192 : rhoecab=ecrs*rhotot-rhoecrs1_dn-rhoecrs1_up
1106 1032192 : vxcadd=ecrs-rs*third*decrs_drs-zeta*decrs_dzeta
1107 1032192 : drhoecab_drhoup=vxcadd+decrs_dzeta-drhoecrs1_drhoup
1108 1032192 : drhoecab_drhodn=vxcadd-decrs_dzeta-drhoecrs1_drhodn
1109 :
1110 : ! Now, the gradient correction factor
1111 1032192 : ssavg=half*(ss_up+ss_dn)
1112 1032192 : divcab=one/(one+gammacab*ssavg)
1113 1032192 : ucab=gammacab*ssavg*divcab
1114 1032192 : ducabdss=gammacab*divcab*(one-ssavg*gammacab*divcab)
1115 :
1116 1032192 : gcab=ccab0+ucab*(ccab1+ucab*(ccab2+ucab*(ccab3+ucab*ccab4)))
1117 : dgcabdss=(ccab1+ucab*(two*ccab2+ucab*(three*ccab3+ucab*four*ccab4)))&
1118 1032192 : & *ducabdss
1119 :
1120 1032192 : exc=exc+rhoecab*gcab
1121 :
1122 1032192 : vxci(ipts,1)=vxci(ipts,1)+drhoecab_drhoup*gcab+rhoecab*dgcabdss*half*dssupdrho
1123 1032192 : dvxcdgr(ipts,1)=dvxcdgr(ipts,1)+rhoecab*dgcabdss*half*dssupdg
1124 1032192 : vxci(ipts,2)=vxci(ipts,2)+drhoecab_drhodn*gcab+rhoecab*dgcabdss*half*dssdndrho
1125 1032192 : dvxcdgr(ipts,2)=dvxcdgr(ipts,2)+rhoecab*dgcabdss*half*dssdndg
1126 :
1127 : ! Final division by the total density, to give the energy density
1128 1032470 : exci(ipts)=exc*rhotot_inv
1129 :
1130 : end do ! ipts=1,npts
1131 :
1132 : end if
1133 :
1134 : else
1135 : ! Disallowed value for nspden
1136 : write(message, '(3a,i0)' )&
1137 0 : & ' Argument nspden must be 1 or 2; ',ch10,&
1138 0 : & ' Value provided as argument was ',nspden
1139 0 : ABI_BUG(message)
1140 : end if
1141 :
1142 : !DEBUG
1143 : !Finite-difference debugging, do not take away
1144 : !Beware that dvxcdgr(:,3) no longer exists
1145 : !if(debug/=0)then
1146 : !do ipts=1,npts,5
1147 :
1148 : !rho=rho_updn(ipts,1)+rho_updn(ipts,2)
1149 : !write(std_out,'(a,i5,a,es16.8)' ) ' Point number',ipts,' with rho=',rho
1150 :
1151 : !! For rho
1152 : !if(debug==1)then
1153 : !write(std_out,'(3es16.8)' )exci(ipts)*rho,vxci(ipts,1),vxci(ipts,2)
1154 : !else
1155 : !! For grho2
1156 : !write(std_out,'(4es16.8)' )exci(ipts)*rho,dvxcdgr(ipts,1),&
1157 : !& dvxcdgr(ipts,2),dvxcdgr(ipts,3)
1158 : !end if
1159 :
1160 : !if(debug==1)then
1161 : !! For rho
1162 : !write(std_out,'(3es16.8)' )exci(ipts)*rho,&
1163 : !& ( exci(ipts+1)*(rho+delta) - exci(ipts+2)*(rho-delta) )/2.d0/delta,&
1164 : !& ( exci(ipts+3)*(rho+delta) - exci(ipts+4)*(rho-delta) )/2.d0/delta
1165 : !write(std_out,'(3es16.8)' )&
1166 : !& ( vxci(ipts+1,1) - vxci(ipts+2,1) )/2.d0/delta,&
1167 : !& ( vxci(ipts+3,1) - vxci(ipts+4,1) )/2.d0/delta,&
1168 : !& ( vxci(ipts+3,2) - vxci(ipts+4,2) )/2.d0/delta
1169 : !write(std_out,'(4es16.8)' )&
1170 : !& ( dvxcdgr(ipts+1,1) - dvxcdgr(ipts+2,1) )/2.d0/delta,&
1171 : !& ( dvxcdgr(ipts+3,2) - dvxcdgr(ipts+4,2) )/2.d0/delta,&
1172 : !& ( dvxcdgr(ipts+1,3) - dvxcdgr(ipts+2,3) )/2.d0/delta,&
1173 : !& ( dvxcdgr(ipts+3,3) - dvxcdgr(ipts+4,3) )/2.d0/delta
1174 : !else
1175 : !! For grho2 (should distinguish exchange and correlation ...)
1176 : !grr=sqrt(grho2_updn(ipts,1)) ! Analysis of exchange
1177 : !grr=sqrt(grho2_updn(ipts,3)) ! Analysis of correlation
1178 : !write(std_out,'(3es16.8)' )exci(ipts)*rho,&
1179 : !& ( exci(ipts+1)*rho - exci(ipts+2)*rho )/2.d0/delta/grr,&
1180 : !& ( exci(ipts+3)*rho - exci(ipts+4)*rho )/2.d0/delta/grr
1181 : !write(std_out,'(3es16.8)' )&
1182 : !& ( vxci(ipts+1,1) - vxci(ipts+2,1) )/2.d0/delta/grr,&
1183 : !& ( vxci(ipts+3,1) - vxci(ipts+4,1) )/2.d0/delta/grr,&
1184 : !& ( vxci(ipts+3,2) - vxci(ipts+4,2) )/2.d0/delta/grr
1185 : !write(std_out,'(4es16.8)' )&
1186 : !& ( dvxcdgr(ipts+1,1) - dvxcdgr(ipts+2,1) )/2.d0/delta/grr,&
1187 : !& ( dvxcdgr(ipts+3,2) - dvxcdgr(ipts+4,2) )/2.d0/delta/grr,&
1188 : !& ( dvxcdgr(ipts+1,3) - dvxcdgr(ipts+2,3) )/2.d0/delta/grr,&
1189 : !& ( dvxcdgr(ipts+3,3) - dvxcdgr(ipts+4,3) )/2.d0/delta/grr
1190 : !end if
1191 : !end do
1192 : !stop
1193 : !end if
1194 : !ENDDEBUG
1195 :
1196 1265 : ABI_FREE(rhoarr)
1197 1265 : ABI_FREE(rhom1_3)
1198 1265 : ABI_FREE(rho_updnm1_3)
1199 1265 : if(nspden==2) then
1200 278 : ABI_FREE(zetm)
1201 278 : ABI_FREE(zetmm1_3)
1202 278 : ABI_FREE(zetp)
1203 278 : ABI_FREE(zetpm1_3)
1204 : end if
1205 :
1206 : !DEBUG
1207 : !write(std_out,*)' xchcth : exit'
1208 : !write(std_out,*)' nspden=',nspden
1209 : !if(order==2)stop
1210 : !ENDDEBUG
1211 :
1212 1265 : end subroutine xchcth
1213 : !!***
1214 :
1215 : end module m_xchcth
1216 : !!***
|
