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