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