1#ifndef SECOND_DERIV 2 Subroutine xc_xpw6(tol_rho, fac, lfac, nlfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,ijzy) 4#else 5 Subroutine xc_xpw6_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 6 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 7 & qwght,ldew,func,ijzy) 8#endif 9c 10C$Id$ 11c 12 implicit none 13c 14#include "dft2drv.fh" 15c 16 double precision fac, Ex 17 integer nq, ipol,ijzy 18 logical lfac, nlfac,ldew 19 double precision func(*) ! value of the functional [output] 20c 21c Charge Density & Its Cube Root 22c 23 double precision rho(nq,ipol*(ipol+1)/2) 24c 25c Charge Density Gradient 26c 27 double precision delrho(nq,3,ipol) 28c 29c Quadrature Weights 30c 31 double precision qwght(nq) 32c 33c Sampling Matrices for the XC Potential & Energy 34c 35 double precision Amat(nq,ipol), Cmat(nq,*) 36c 37c 38c Compute the partial derivatives of the exchange functional of Perdew91. 39c 40c Becke & Perdew Parameters 41c 42 double precision DPOW 43 double precision BETA, tol_rho, CPW91 44 45#ifdef SECOND_DERIV 46c 47c Second Derivatives of the Exchange Energy Functional 48c 49 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 50#endif 51c 52c References: 53c 54c 55c*************************************************************************** 56 57 integer n 58 double precision gamma 59c 60 if (ijzy.eq.2) then 61c 62c parameters for PWB6K 63c 64 BETA = 0.00539d0 65 CPW91 = 1.7077d0 66 DPOW = 4.0876d0 67 else 68c 69c parameters for PW6B95 70c 71 BETA = 0.00538d0 72 CPW91 = 1.7382d0 73 DPOW = 3.8901d0 74 endif 75 if (ipol.eq.1 )then 76c 77c ======> SPIN-RESTRICTED <====== 78c 79 do 10 n = 1, nq 80 if (rho(n,1).lt.tol_rho) goto 10 81 gamma = delrho(n,1,1)*delrho(n,1,1) + 82 & delrho(n,2,1)*delrho(n,2,1) + 83 & delrho(n,3,1)*delrho(n,3,1) 84 gamma=0.25d0*gamma 85#ifndef SECOND_DERIV 86 call xc_xpw6core(DPOW,BETA,CPW91,n,1, 87 & rho(n,1),gamma,qwght(n),func(n), 88 & tol_rho, fac, lfac, nlfac, 89 & Amat, Cmat, nq, ipol, Ex, ldew) 90#else 91 call xc_xpw6core_d2(DPOW,BETA,CPW91,n,1, 92 & rho(n,1),gamma,qwght(n),func(n), 93 & tol_rho, fac, lfac, nlfac, 94 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 95 & ldew) 96#endif 97c write (*,*) "Ex",Ex,Amat(n,1),Cmat(n,1) 98c stop 99 10 continue 100c 101 else 102c 103c ======> SPIN-UNRESTRICTED <====== 104c 105 do 20 n = 1, nq 106 if (rho(n,1).lt.tol_rho) goto 20 107c 108c Spin alpha: 109c 110 if (rho(n,2).gt.tol_rho) then 111 gamma = delrho(n,1,1)*delrho(n,1,1) + 112 & delrho(n,2,1)*delrho(n,2,1) + 113 & delrho(n,3,1)*delrho(n,3,1) 114#ifndef SECOND_DERIV 115 call xc_xpw6core(DPOW,BETA,CPW91,n,1, 116 & rho(n,2),gamma,qwght(n),func(n), 117 & tol_rho, fac, lfac, nlfac, 118 & Amat, Cmat, nq, ipol, Ex, ldew) 119#else 120 call xc_xpw6core_d2(DPOW,BETA,CPW91,n,1, 121 & rho(n,2),gamma,qwght(n),func(n), 122 & tol_rho, fac, lfac, nlfac, 123 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 124 & ldew) 125#endif 126 endif 127c 128c Spin beta: 129c 130 if (rho(n,3).gt.tol_rho) then 131 132 gamma = delrho(n,1,2)*delrho(n,1,2) + 133 & delrho(n,2,2)*delrho(n,2,2) + 134 & delrho(n,3,2)*delrho(n,3,2) 135#ifndef SECOND_DERIV 136 call xc_xpw6core(DPOW,BETA,CPW91,n,2, 137 & rho(n,3),gamma,qwght(n),func(n), 138 & tol_rho, fac, lfac, nlfac, 139 & Amat, Cmat, nq, ipol, Ex, ldew) 140#else 141 call xc_xpw6core_d2(DPOW,BETA,CPW91,n,2, 142 & rho(n,3),gamma,qwght(n),func(n), 143 & tol_rho, fac, lfac, nlfac, 144 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 145 & ldew) 146#endif 147 endif 148c 149 20 continue 150c 151 endif 152c 153 return 154 end 155 156 157 158 159#ifndef SECOND_DERIV 160 Subroutine xc_xpw6core(DPOW,BETA,CPW91,n,ispin, 161 . rho,gamma,qwght,func, 162 & tol_rho, fac, lfac, nlfac, 163 & Amat, Cmat, nq, ipol, Ex, ldew) 164#else 165 Subroutine xc_xpw6core_d2(DPOW,BETA,CPW91,n,ispin, 166 . rho,gamma,qwght,func, 167 & tol_rho, fac, lfac, nlfac, 168 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 169 & ldew) 170#endif 171c 172C$Id$ 173c 174 implicit none 175c 176#include "dft2drv.fh" 177c 178 double precision fac, Ex 179 integer nq, ipol 180 logical lfac, nlfac,ldew 181 double precision func ! value of the functional [output] 182 double precision rho 183 double precision qwght 184 integer ispin ! alpha=1; beta=2 185c 186c Sampling Matrices for the XC Potential & Energy 187c 188 double precision Amat(nq,ipol), Cmat(nq,*) 189c 190c 191c Compute the partial derivatives of the exchange functional of Perdew91. 192c 193c Becke & Perdew Parameters 194c 195 double precision DPOW 196 double precision BETA, tol_rho, AX, pi, 197 & CPW91,BETAPW91,big 198 199 Parameter (big=1d4) 200 201#ifdef SECOND_DERIV 202c 203c Second Derivatives of the Exchange Energy Functional 204c 205 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 206 double precision rhom23, f2x, d2den,d2num 207#endif 208c 209c References: 210c 211c 212c*************************************************************************** 213 214 integer n 215 double precision x, sinhm1,dsinhm1 216 double precision rho13, rho43, Xa, Ha, denom, num, 217 & fprimex, d1num,d1den, 218 & gamma,fx,x2a,ten6xd,expo,bbx2 219 integer D0R,D1G,D2RR,D2RG,D2GG 220 221c 222 sinhm1(x)=log(x+dsqrt(1d0+x*x)) 223 dsinhm1(x)=1d0/dsqrt(1d0+x*x) 224 pi=acos(-1.d0) 225 BETAPW91=(pi*36.d0)**(-5.d0/3.d0)*5.d0 226 AX=-(0.75d0/pi)**(1.d0/3.d0)*1.5d0 227 if(ispin.eq.1) then 228 D0R=1 229 D1G=D1_GAA 230 D2RR=D2_RA_RA 231 D2RG=D2_RA_GAA 232 D2GG=D2_GAA_GAA 233 else 234 D0R=2 235 D1G=D1_GBB 236 D2RR=D2_RB_RB 237 D2RG=D2_RB_GBB 238 D2GG=D2_GBB_GBB 239 endif 240 241 rho13 = (rho*ipol/2d0)**(1.d0/3.d0) 242 rho43 = rho13**4 243 if (gamma.gt.tol_rho**2)then 244 xa = sqrt(gamma)/rho43 245 x2a=xa*xa 246 ten6xd=Xa**DPOW*1.d-6 247 expo=0d0 248 if(CPW91*x2a.lt.big) expo=exp(-CPW91*x2a) 249 bbx2=(BETA-BETAPW91)*x2a 250 Ha = sinhm1(Xa) 251 denom = 1.d0/(1.d0 + 6d0*(beta*xa)*ha-ten6xd/ax) 252 num = -BETA*x2a+bbx2*expo+ten6xd 253 fx=num*denom 254 d1num=-2.d0*xa*(beta-bbx2*expo*(1d0/x2a-CPW91))+ 255 + ten6xd/xa*dpow 256 d1den=6.d0*beta*(ha + xa*dsinhm1(xa)) - 257 - ten6xd/ax/xa*dpow 258 fprimex=(d1num - d1den*fx)*denom 259 else 260 gamma = 0.d0 261 Xa = 0.d0 262 fx=0.d0 263 fprimex=0.d0 264 denom=0d0 265 d1den=0d0 266 expo=0d0 267 ten6xd=0d0 268 x2a=0d0 269 bbx2=0d0 270 endif 271c 272 if (lfac)then 273 Ex = Ex + 2d0/ipol*rho43*AX*qwght*fac 274 if(ldew)func = func + 2.d0/ipol*rho43*AX*fac 275 Amat(n,D0R) = Amat(n,D0R) + (4.d0/3.d0)*rho13*AX*fac 276 endif 277c 278 if (nlfac)then 279 Ex = Ex + 2d0/ipol*rho43*fx*qwght*fac 280 if(ldew)func = func + 2.d0/ipol*rho43*fx*fac 281 Amat(n,D0R) = Amat(n,D0R) + 282 + (4.d0/3.d0)*rho13*(fx-xa*fprimex)*fac 283 if (xa.gt.tol_rho) then 284 Cmat(n,D1G)=Cmat(n,D1G)+ 285 + .5d0*fprimex/sqrt(gamma)*fac 286 endif 287 endif 288c 289#ifdef SECOND_DERIV 290 rhom23 = 2d0/ipol*rho13/rho 291 if (lfac)then 292 Amat2(n,D2RR) = Amat2(n,D2RR) + 293 & (4d0/9d0)*rhom23*Ax*fac 294 endif 295 if (nlfac)then 296 if(gamma.gt.tol_rho**2)then 297 d2num=-2d0*beta + 298 + 2d0*bbx2*expo* 299 * (1d0/x2a-5d0*CPW91+2d0*CPW91*CPW91*x2a)+ 300 + ten6xd/x2a*dpow*(dpow-1d0) 301 d2den=6.d0*beta*dsinhm1(xa)*(2d0-x2a/(1d0+x2a)) - 302 - ten6xd/ax/x2a*dpow*(dpow-1d0) 303 f2x = denom*(d2num -fx*d2den - 2d0*d1den*fprimex) 304 else 305 f2x=0d0 306 endif 307 Amat2(n,D2RR) = Amat2(n,D2RR) 308 & + (4d0/9d0)*rhom23*(fx-xa*fprimex+4d0*x2a*f2x)*fac 309 Cmat2(n,D2RG) = Cmat2(n,D2RG) 310 & - (4d0/ipol/3d0)*(rhom23**2/rho)*f2x*fac 311 if (xa.gt.tol_rho) then 312 Cmat2(n,D2GG) = Cmat2(n,D2GG) 313 & - 0.25d0*gamma**(-1.5d0)*(fprimex-xa*f2x)*fac 314 endif 315 endif 316#endif 317c 318c 319 return 320 end 321#ifndef SECOND_DERIV 322#define SECOND_DERIV 323c 324c Compile source again for the 2nd derivative case 325c 326#include "xc_xpw6.F" 327#endif 328