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