TABLE OF CONTENTS
ABINIT/xcspol [ Functions ]
NAME
xcspol
FUNCTION
Spin-polarized exchange and correlation, parameterized by Mike Teter of Corning Incorporated.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
nspden=number of spin-density components npts= number of points to be computed order=its absolute value gives the maximal derivative of Exc to be computed. rspts(npts)=Seitz electron radius (bohr) zeta(npts)=$(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$=degree of polarization (ignored if nspden=1, in which case zeta should be 0)
OUTPUT
if(abs(order)>1) dvxc(npts,1+nspden)= (Hartree*bohr^3)
if(nspden=1 .and. order==2): dvxc(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty
if(nspden=1 .and. order==-2): also compute dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$
if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\uparrow)$,
dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$
exc(npts)=exchange-correlation energy density (hartree)
vxc(npts,nspden)=xc potent. (d($\rho$*exc)/d($\rho\uparrow$)) and d/d($\rho\downarrow$) (ha)
(only overall potential d($\rho$*exc)/d($\rho$) returned in vxc(1) for nspden=1)
ndvxc= size of dvxc(npts,ndvxc)
Normalization: Exc=$\int(exc(r)*\rho(r) d^3 r)$ for $\rho$(r)=electron density.
TODO
To be added later d2vxc=derivative $d^2 (Vxc)/d(rho)^2$ (hartree*bohr^6)
NOTES
This form is based on Mike Teter s rational polynomial exc=-(a0+a1*rs+a2*rs**2+a3*rs**3)/(b1*rs+b2*rs**2+b3*rs**3+b4*rs**4) where the parameters are fit to reproduce (in this case) the Perdew-Wang parameterization of the correlation energy given in Phys. Rev. B 45, 13244-13249 (1992). Each parameter is interpolated between zeta=0 and 1 by a_i(zeta)=a_i(0)+(a_i(1)-a_i(0))*f_x(zeta) and f_x(zeta)=[(1+zeta)$^{4/3}$+(1-zeta)$^{4/3}$-2]/(2*(2$^{1/3}$-1)). Beware : in this expression, zeta is actually replaced by zeta*alpha_zeta, where alpha_zeta is very close to 1, but slightly lower. This is to remove the singularity in the derivatives when abs(zeta) is 1 Below, a_i(1)-a_i(0) is called "da" for delta a, same for b s. rs = $(3/(4\pi))^{1/3} * \rho(r)^{-1/3}$ zeta = $(\rho\uparrow-\rho\downarrow)/(\rho\uparrow+\rho\downarrow)$ b1 must be 1 and a0 must be $(3/4)(3/(2\pi))^{2/3}$.
PARENTS
drivexc
CHILDREN
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 78 subroutine xcspol(exc,npts,nspden,order,rspts,vxc,zeta,ndvxc,& !Mandatory arguments 79 & dvxc) !Optional arguments 80 81 use defs_basis 82 use m_profiling_abi 83 use m_errors 84 85 !This section has been created automatically by the script Abilint (TD). 86 !Do not modify the following lines by hand. 87 #undef ABI_FUNC 88 #define ABI_FUNC 'xcspol' 89 !End of the abilint section 90 91 implicit none 92 93 !Arguments ------------------------------------ 94 !scalars 95 integer,intent(in) :: ndvxc,npts,nspden,order 96 !arrays 97 real(dp),intent(in) :: rspts(npts),zeta(npts) 98 real(dp),intent(out) :: exc(npts),vxc(npts,nspden) 99 real(dp),intent(out),optional :: dvxc(npts,ndvxc) 100 101 !Local variables------------------------------- 102 !The generation of density from rs needs rsfac and rsfac^(-3) : 103 !rsfac=(3/(4 Pi))^(1/3) ; rsfacm3=4pi/3 104 !Mike Teter s parameters of 8 April 1993. 105 !New parameters which accomodate spin polarization (fit to P-W) 106 !Paramagnetic limit:a0p,...b4p 107 !(a0=(3/4)(3/(2 Pi))^(2/3)) (note that b1=1 is fixed) 108 !Differences, ferromagnetic - paramagnetic (delta params):da1,da2,da3,db1,db2,db3,db4 109 !scalars 110 integer :: ipts 111 real(dp),parameter :: a0p=.4581652932831429_dp,a1p=2.217058676663745_dp 112 real(dp),parameter :: a2p=0.7405551735357053_dp,a3p=0.01968227878617998_dp 113 real(dp),parameter :: alpha_zeta=one-1.0d-6,b1p=one,b2p=4.504130959426697_dp 114 real(dp),parameter :: b3p=1.110667363742916_dp,b4p=0.02359291751427506_dp 115 real(dp),parameter :: da0=.119086804055547_dp,da1=0.6157402568883345_dp 116 real(dp),parameter :: da2=0.1574201515892867_dp,da3=0.003532336663397157_dp 117 real(dp),parameter :: db1=zero,db2=0.2673612973836267_dp 118 real(dp),parameter :: db3=0.2052004607777787_dp,db4=0.004200005045691381_dp 119 real(dp),parameter :: ft=4._dp/3._dp,rsfac=0.6203504908994000_dp 120 real(dp),parameter :: rsfacm3=rsfac**(-3) 121 real(dp) :: a0,a1,a2,a3,b1,b2,b3,b4,d1,d1m1,d2d1drs2,d2d1drsdf,d2excdf2 122 real(dp) :: d2excdrs2,d2excdrsdf,d2excdz2,d2fxcdz2,d2n1drs2,d2n1drsdf,dd1df 123 real(dp) :: dd1drs,dexcdf,dexcdrs,dexcdz,dfxcdz,dn1df,dn1drs,dvxcdrs 124 real(dp) :: dvxcpdrho,dvxcpdz,excipt,fact,fxc,n1 125 real(dp) :: rhom1,rs,vxcp,zet,zetm,zetm_third 126 real(dp) :: zetp,zetp_third 127 character(len=500) :: message 128 !no_abirules 129 !Set a minimum rho below which terms are 0 130 real(dp),parameter :: rhotol=1.d-28 131 !real(dp) :: delta,rho,rho_dn,rho_dnm,rho_dnp,rho_up,rho_upm,rho_upp,zeta_mean 132 133 ! ************************************************************************* 134 135 !Checks the compatibility between the presence of dvxc and ndvxc 136 if(ndvxc /=0 .neqv. present(dvxc))then 137 message = 'If ndvxc/=0 there must be the optional argument dvxc' 138 MSG_BUG(message) 139 end if 140 141 !Checks the compatibility between the inputs and the presence of the optional arguments 142 if(abs(order) <= 1 .and. ndvxc /= 0)then 143 write(message, '(4a,i0)' )ch10,& 144 & 'The order chosen does not need the presence',ch10,& 145 & 'of the vector dvxc, that is needed only with |order|>1 , while we have',order 146 MSG_BUG(message) 147 end if 148 149 if(nspden == 1 .and. ndvxc /=0 .and. ndvxc /= 2)then 150 write(message,'(a,i0)')' Once nspden=1 we must have ndvxc=2, while we have',ndvxc 151 MSG_BUG(message) 152 end if 153 154 if(nspden == 2 .and. ndvxc /=0 .and. ndvxc /= 3)then 155 write(message, '(a,i0)' )' Once nspden=2 we must have ndvxc=3, while we have',ndvxc 156 MSG_BUG(message) 157 end if 158 159 160 !Although fact is parameter value, some compilers are not able to evaluate 161 !it at compile time. 162 fact=one/(two**(four*third)-two) 163 164 !DEBUG 165 !Finite-difference debugging, do not take away 166 !debug=1 167 !zeta_mean=0.1_dp 168 !delta=0.0001 169 !if(debug==1)then 170 !do ipts=1,npts,5 171 !rho=ipts*0.01_dp 172 !rho_up=rho*(one+zeta_mean)*half 173 !rho_dn=rho*(one-zeta_mean)*half 174 !rho_upp=rho_up+delta 175 !rho_upm=rho_up-delta 176 !rho_dnp=rho_dn+delta 177 !rho_dnm=rho_dn-delta 178 !First possibility : vary rho up , and then rho down 179 !zeta(ipts )=(rho_up -rho_dn )/(rho_up +rho_dn ) 180 !zeta(ipts+1)=(rho_upp-rho_dn )/(rho_upp+rho_dn ) 181 !zeta(ipts+2)=(rho_upm-rho_dn )/(rho_upm+rho_dn ) 182 !zeta(ipts+3)=(rho_up -rho_dnp)/(rho_up +rho_dnp) 183 !zeta(ipts+4)=(rho_up -rho_dnm)/(rho_up +rho_dnm) 184 !rspts(ipts )=rsfac*(rho_up +rho_dn )**(-third) 185 !rspts(ipts+1)=rsfac*(rho_upp+rho_dn )**(-third) 186 !rspts(ipts+2)=rsfac*(rho_upm+rho_dn )**(-third) 187 !rspts(ipts+3)=rsfac*(rho_up +rho_dnp)**(-third) 188 !rspts(ipts+4)=rsfac*(rho_up +rho_dnm)**(-third) 189 !DEBUGBUG : another possibility : vary rho and zeta 190 !zeta(ipts+1)=zeta(ipts ) 191 !zeta(ipts+2)=zeta(ipts ) 192 !zeta(ipts+3)=zeta(ipts )+delta 193 !zeta(ipts+4)=zeta(ipts )-delta 194 !rspts(ipts+1)=rsfac*(rho+delta)**(-third) 195 !rspts(ipts+2)=rsfac*(rho-delta )**(-third) 196 !rspts(ipts+3)=rspts(ipts ) 197 !rspts(ipts+4)=rspts(ipts ) 198 !ENDDEBUGBUG 199 !end do 200 !end if 201 !nspden=2 202 !order=2 203 !ENDDEBUG 204 205 if (nspden==1) then 206 ! separate cases with respect to order 207 if(order==-2) then 208 ! No spin-polarization so skip steps related to zeta not 0 209 do ipts=1,npts 210 211 rs=rspts(ipts) 212 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 213 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 214 d1m1=one/d1 215 216 ! Exchange-correlation energy 217 excipt=-n1*d1m1 218 exc(ipts)=excipt 219 220 ! Exchange-correlation potential 221 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 222 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 223 224 ! dexcdrs is d(exc)/d(rs) 225 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 226 vxc(ipts,1)=excipt-third*rs*dexcdrs 227 228 ! If the exchange-correlation kernel is needed 229 230 d2n1drs2=2._dp*a2p+rs*(6._dp*a3p) 231 d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p)) 232 ! d2excdrs2 is d2(exc)/d(rs)2 233 d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 234 dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2) 235 ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs) 236 dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs 237 238 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 239 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 240 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 241 dexcdf=-(dn1df+excipt*dd1df)*d1m1 242 ! d2(fxc)/d(zeta)2 243 d2fxcdz2=ft*third*(alpha_zeta**2)*2._dp*fact 244 ! d2(exc)/d(zeta)2 245 d2excdz2=d2fxcdz2*dexcdf 246 rhom1=rsfacm3*rs**3 247 dvxc(ipts,2)= dvxc(ipts,1) - d2excdz2*rhom1 248 end do 249 else if(order**2>1) then 250 ! No spin-polarization so skip steps related to zeta not 0 251 do ipts=1,npts 252 253 rs=rspts(ipts) 254 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 255 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 256 d1m1=one/d1 257 258 ! Exchange-correlation energy 259 excipt=-n1*d1m1 260 exc(ipts)=excipt 261 262 ! Exchange-correlation potential 263 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 264 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 265 266 ! dexcdrs is d(exc)/d(rs) 267 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 268 vxc(ipts,1)=excipt-third*rs*dexcdrs 269 270 ! If the exchange-correlation kernel is needed 271 d2n1drs2=2._dp*a2p+rs*(6._dp*a3p) 272 d2d1drs2=2._dp*b2p+rs*(6._dp*b3p+rs*(12._dp*b4p)) 273 ! d2excdrs2 is d2(exc)/d(rs)2 274 d2excdrs2=-(d2n1drs2+2._dp*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 275 dvxcdrs=third*(2.0_dp*dexcdrs-rs*d2excdrs2) 276 ! And d(vxc)/d(rho)=(-rs/(3*rho))*d(vxc)/d(rs) 277 dvxc(ipts,1)= -rs**4*rsfacm3*third*dvxcdrs 278 279 end do 280 else 281 ! No spin-polarization so skip steps related to zeta not 0 282 do ipts=1,npts 283 284 rs=rspts(ipts) 285 n1=a0p+rs*(a1p+rs*(a2p+rs*a3p)) 286 d1=rs*(b1p+rs*(b2p+rs*(b3p+rs*b4p))) 287 d1m1=one/d1 288 289 ! Exchange-correlation energy 290 excipt=-n1*d1m1 291 exc(ipts)=excipt 292 293 ! Exchange-correlation potential 294 dn1drs=a1p+rs*(2._dp*a2p+rs*(3._dp*a3p)) 295 dd1drs=b1p+rs*(2._dp*b2p+rs*(3._dp*b3p+rs*(4._dp*b4p))) 296 297 ! dexcdrs is d(exc)/d(rs) 298 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 299 vxc(ipts,1)=excipt-third*rs*dexcdrs 300 end do 301 302 end if 303 304 305 ! Allows for nspden==1, in the case of testing nspden=1 against nspden=2 306 else if (nspden<=2) then 307 308 309 ! DEBUG 310 ! do not take away : allows to compare nspden=1 and nspden=2 coding 311 ! if (nspden==1)then 312 ! zeta(:)=zero 313 ! end if 314 ! ENDDEBUG 315 ! separate cases with respect to order 316 if(abs(order)>1) then 317 ! Allow for spin polarization. This part could be optimized for speed. 318 do ipts=1,npts 319 320 rs=rspts(ipts) 321 zet=zeta(ipts) 322 zetp=one+zet*alpha_zeta 323 zetm=one-zet*alpha_zeta 324 zetp_third=zetp**third 325 zetm_third=zetm**third 326 ! Exchange energy spin interpolation function f(zeta) 327 fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact 328 329 a0=a0p+fxc*da0 330 a1=a1p+fxc*da1 331 a2=a2p+fxc*da2 332 a3=a3p+fxc*da3 333 b1=b1p+fxc*db1 334 b2=b2p+fxc*db2 335 b3=b3p+fxc*db3 336 b4=b4p+fxc*db4 337 338 n1= a0+rs*(a1+rs*(a2+rs*a3)) 339 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 340 d1m1=one/d1 341 342 ! Exchange-correlation energy 343 excipt=-n1*d1m1 344 exc(ipts)=excipt 345 346 ! Exchange-correlation potential 347 dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3)) 348 dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 349 ! dexcdrs is d(exc)/d(rs) 350 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 351 352 ! Only vxcp contributes when paramagnetic 353 vxcp=excipt-third*rs*dexcdrs 354 355 ! d(fxc)/d(zeta) (which is 0 at zeta=0) 356 dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact 357 358 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 359 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 360 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 361 362 ! dexcdz is d(exc)/d(zeta) 363 dexcdf=-(dn1df+excipt*dd1df)*d1m1 364 dexcdz=dfxcdz*dexcdf 365 366 ! Compute Vxc for both spin channels 367 368 vxc(ipts,1)=vxcp - (zet-one)*dexcdz 369 vxc(ipts,2)=vxcp - (zet+one)*dexcdz 370 371 ! DEBUG Allow to check the variation of rho and zeta 372 ! vxc(ipts,1)=vxcp 373 ! vxc(ipts,2)=dexcdz 374 ! ENDDEBUG 375 ! Compute second derivative with respect to rho 376 d2n1drs2=2._dp*a2+rs*(6._dp*a3) 377 d2d1drs2=2._dp*b2+rs*(6._dp*b3+rs*(12._dp*b4)) 378 ! d2excdrs2 is d2(exc)/d(rs)2 379 d2excdrs2=-(d2n1drs2+two*dexcdrs*dd1drs+excipt*d2d1drs2)*d1m1 380 dvxcdrs=third*(two*dexcdrs-rs*d2excdrs2) 381 ! And d(vxc)/d(rho) paramagnetic =(-rs/(3*rho))*d(vxcp)/d(rs) 382 ! remember : 1/rho=(4pi/3)*rs**3=rsfacm3*rs**3 383 rhom1=rsfacm3*rs**3 384 dvxcpdrho= -rs*rhom1*third * dvxcdrs 385 386 ! Compute mixed second derivative with respect to rho and zeta 387 d2n1drsdf=da1+rs*(2._dp*da2+rs*(3._dp*da3)) 388 d2d1drsdf=db1+rs*(2._dp*db2+rs*(3._dp*db3+rs*(4._dp*db4))) 389 ! d2excdrsdf is d2(exc)/d(rs)df 390 d2excdrsdf=-(d2n1drsdf+dexcdrs*dd1df+dexcdf*dd1drs+excipt*d2d1drsdf)*d1m1 391 ! d(vxc)/d(zeta) paramagnetic 392 dvxcpdz=dexcdz-third*rs*dfxcdz*d2excdrsdf 393 394 ! Compute second derivative with respect to zeta 395 ! the second derivative of n1 and d1 wrt f vanishes 396 d2excdf2=-(two*dexcdf*dd1df)*d1m1 397 ! d2(fxc)/d(zeta)2 398 d2fxcdz2=ft*third*(alpha_zeta**2)*(zetp_third**(-2)+zetm_third**(-2))*fact 399 ! d2(exc)/d(zeta)2 400 d2excdz2=d2fxcdz2*dexcdf+dfxcdz**2*d2excdf2 401 402 ! Compute now the three second derivatives of the Exc energy with respect 403 ! to : wrt twice spin-up ; wrt spin-up and spin-dn ; wrt twice spin-down 404 dvxc(ipts,1)= dvxcpdrho & 405 & +two*rhom1*( one-zet)*(dvxcpdz-dexcdz) & 406 & +d2excdz2*rhom1*(one-zet)**2 407 dvxc(ipts,2)= dvxcpdrho & 408 & +two*rhom1*( -zet)*(dvxcpdz-dexcdz) & 409 & +d2excdz2*rhom1*(one-zet)*(-one-zet) 410 ! if(nspden==2)then 411 dvxc(ipts,3)= dvxcpdrho & 412 & +two*rhom1*(-one-zet)*(dvxcpdz-dexcdz) & 413 & +d2excdz2*rhom1*(-one-zet)**2 414 ! else 415 ! ! For testing purposes, need the spin-averaged quantity 416 ! dvxc(ipts,1)= ( dvxc(ipts,1) + dvxc(ipts,2) ) * half 417 ! end if 418 419 ! DEBUG Allow to check the variation of rho and zeta 420 ! dvxc(ipts,1)=dvxcpdrho 421 ! dvxc(ipts,2)=d2excdz2 422 ! dvxc(ipts,3)=dvxcpdz 423 ! ENDDEBUG 424 end do 425 else 426 ! Allow for spin polarization. This part could be optimized for speed. 427 do ipts=1,npts 428 429 rs=rspts(ipts) 430 zet=zeta(ipts) 431 zetp=one+zet*alpha_zeta 432 zetm=one-zet*alpha_zeta 433 zetp_third=zetp**third 434 zetm_third=zetm**third 435 ! Exchange energy spin interpolation function f(zeta) 436 fxc=( zetp*zetp_third + zetm*zetm_third - two ) *fact 437 438 a0=a0p+fxc*da0 439 a1=a1p+fxc*da1 440 a2=a2p+fxc*da2 441 a3=a3p+fxc*da3 442 b1=b1p+fxc*db1 443 b2=b2p+fxc*db2 444 b3=b3p+fxc*db3 445 b4=b4p+fxc*db4 446 447 n1= a0+rs*(a1+rs*(a2+rs*a3)) 448 d1=rs*(b1+rs*(b2+rs*(b3+rs*b4))) 449 d1m1=one/d1 450 451 ! Exchange-correlation energy 452 excipt=-n1*d1m1 453 exc(ipts)=excipt 454 455 ! Exchange-correlation potential 456 dn1drs=a1+rs*(2._dp*a2+rs*(3._dp*a3)) 457 dd1drs=b1+rs*(2._dp*b2+rs*(3._dp*b3+rs*(4._dp*b4))) 458 ! dexcdrs is d(exc)/d(rs) 459 dexcdrs=-(dn1drs+excipt*dd1drs)*d1m1 460 461 ! Only vxcp contributes when paramagnetic 462 vxcp=excipt-third*rs*dexcdrs 463 464 ! d(fxc)/d(zeta) (which is 0 at zeta=0) 465 dfxcdz=ft*alpha_zeta*(zetp_third-zetm_third)*fact 466 467 ! dn1df=d(n1)/d(fxc) and dd1df=d(d1)/d(fxc) 468 dn1df=da0+rs*(da1+rs*(da2+rs*da3)) 469 dd1df=rs*(db1+rs*(db2+rs*(db3+rs*db4))) 470 471 ! dexcdz is d(exc)/d(zeta) 472 dexcdf=-(dn1df+excipt*dd1df)*d1m1 473 dexcdz=dfxcdz*dexcdf 474 475 ! Compute Vxc for both spin channels 476 477 vxc(ipts,1)=vxcp - (zet-one)*dexcdz 478 vxc(ipts,2)=vxcp - (zet+one)*dexcdz 479 480 ! DEBUG Allow to check the variation of rho and zeta 481 ! vxc(ipts,1)=vxcp 482 ! vxc(ipts,2)=dexcdz 483 ! ENDDEBUG 484 end do 485 end if 486 else 487 488 ! Disallowed value for nspden 489 write(message, '(3a,i0)' )& 490 & ' Argument nspden must be 1 or 2; ',ch10,& 491 & ' Value provided as argument was ',nspden 492 MSG_BUG(message) 493 end if 494 495 !DEBUG 496 !Finite-difference debugging, do not take away 497 !if(debug==1)then 498 !write(std_out,*)' delta =',delta 499 !do ipts=1,npts,5 500 !rho=(rspts(ipts)/rsfac)**(-3) 501 !write(std_out,'(a,i5,a,2es16.8)' ) ' Point number',ipts,' with rho,zeta=',rho,zeta(ipts) 502 !write(std_out,'(3es16.8)' )exc(ipts)*rho,vxc(ipts,1),vxc(ipts,2) 503 !write(std_out,'(3es16.8)' )dvxc(ipts,1),dvxc(ipts,3),dvxc(ipts,2) 504 !write(std_out,'(3es16.8)' )exc(ipts)*rho,& 505 !& ( exc(ipts+1)*(rho+delta) - exc(ipts+2)*(rho-delta) )/2._dp/delta,& 506 !& ( exc(ipts+3)*(rho+delta) - exc(ipts+4)*(rho-delta) )/2._dp/delta 507 !write(std_out,'(4es16.8)' )& 508 !& ( vxc(ipts+1,1) - vxc(ipts+2,1) )/2._dp/delta,& 509 !& ( vxc(ipts+3,2) - vxc(ipts+4,2) )/2._dp/delta,& 510 !& ( vxc(ipts+3,1) - vxc(ipts+4,1) )/2._dp/delta,& 511 !& ( vxc(ipts+1,2) - vxc(ipts+2,2) )/2._dp/delta 512 !end do 513 !stop 514 !end if 515 !ENDDEBUG 516 517 !DEBUG 518 !if(order==-2)then 519 !write(std_out,*)' xcspol : ipts,npts ',ipts,npts 520 !write(std_out,*)dvxcdrs,d2excdz2,d2fxcdz2,dexcdf 521 !write(std_out,*)rhom1 522 !write(std_out,*)dvxc(1000,1),dvxc(1000,2) 523 !stop 524 !end if 525 !ENDDEBUG 526 527 end subroutine xcspol