1c 2c$Id$ 3c 4#include "dft2drv.fh" 5c Tao,Perdew, Staroverov, Scuseria exchange functional 6c META GGA 7C utilizes ingredients: 8c rho - density 9c delrho - gradient of density 10c tau - K.S kinetic energy density 11c tauW - von Weiszacker kinetic energy density 12c tauU - uniform-gas KE density 13c References: 14c [a] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria 15c PRL 91, 146401 (2003). 16c [b] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria 17c JCP 120, 6898 (2004). 18 Subroutine xc_xtpss03(tol_rho, fac, rho, delrho, 19 & Amat, Cmat, nq, ipol, Ex, 20 & qwght, ldew, func, tau,Mmat) 21 implicit none 22c 23 double precision fac, Ex 24 integer nq, ipol 25 logical ldew 26 double precision func(*) ! value of the functional [output] 27c 28c Charge Density & Its Cube Root 29c 30 double precision rho(nq,ipol*(ipol+1)/2) 31c 32c Charge Density Gradient 33c 34 double precision delrho(nq,3,ipol) 35c 36c Quadrature Weights 37c 38 double precision qwght(nq) 39c 40c Sampling Matrices for the XC Potential & Energy 41c 42 double precision Amat(nq,ipol), Cmat(nq,*) 43c 44c kinetic energy density or tau 45c 46 double precision tau(nq,ipol), Mmat(nq,*) 47 double precision tol_rho 48c 49 integer ispin,cmatpos 50c 51 if (ipol.eq.1 )then 52c 53c SPIN-RESTRICTED 54c Ex = Ex[n] 55c 56 call xc_xtpss03_cs(tol_rho, fac, rho, delrho, 57 & Amat, Cmat, nq, Ex, 1d0, 58 & qwght, ldew, func, tau,Mmat) 59 else 60c 61c SPIN-UNRESTRICTED 62c Ex = Ex[2n_up]/2 + Ex[2n_down]/2 63 64 do ispin=1,2 65 if (ispin.eq.1) cmatpos=D1_GAA 66 if (ispin.eq.2) cmatpos=D1_GBB 67 call xc_xtpss03_cs(tol_rho, fac, 68 R rho(1,ispin+1), delrho(1,1,ispin), 69 & Amat(1,ispin), Cmat(1,cmatpos), 70 & nq, Ex, 2d0, 71 & qwght, ldew, func, 72 T tau(1,ispin),Mmat(1,ispin)) 73 enddo 74 75 endif 76 return 77 end 78 Subroutine xc_xtpss03_cs(tol_rho, fac, rho, delrho, 79 & Amat, Cmat, nq, Ex, facttwo, 80 & qwght, ldew, func, tau,Mmat) 81 implicit none 82c 83 double precision fac, Ex 84 integer nq 85 logical ldew 86 double precision func(*) ! value of the functional [output] 87c 88c Charge Density & Its Cube Root 89c 90 double precision rho(*) 91c 92c Charge Density Gradient 93c 94 double precision delrho(nq,3) 95c 96c Quadrature Weights 97c 98 double precision qwght(nq) 99c 100c Sampling Matrices for the XC Potential & Energy 101c 102 double precision Amat(nq), Cmat(nq) 103c 104c kinetic energy density or tau 105c 106 double precision tau(nq,*), Mmat(nq) 107c 108 double precision facttwo ! 2 for o.s. 1 for c.s. 109c 110 double precision tol_rho, pi 111 integer n 112 double precision rrho, rho43, rho13, gamma 113 double precision tauN, tauW, tauU 114 115 double precision p, qtil, x, al, mt, z 116 double precision F83, F23, F53, F43, F13 117 double precision G920 118 double precision b,c,e,es 119 double precision C1, C2, C3 120 double precision kap, mu 121 double precision xb,xc,xd 122 double precision x1,x2,x3,x4,x5,x6,x7 123 double precision P32, Ax 124c functional derivatives below FFFFFFFFFFFF 125 double precision dzdn, dpdn, daldn, dqtdn 126 double precision til1, til2 127 double precision dtil2dn, dtil1dn 128 double precision ax1, bx1, dx1dn 129 double precision dx2dn 130 double precision dxddn, dxcdn, dx3dn 131 double precision dx4dn, dx5dn, dx6dn, dx7dn 132 double precision dxdn 133 double precision xmany, dxmanydn 134 double precision dmtdn, derivn 135 136 double precision dzdg, dpdg, daldg, dqtdg 137 double precision dtil2dg 138 double precision dx1dg, dx2dg 139 double precision dxcdg, dxddg,dx3dg 140 double precision dx4dg, dx5dg, dx6dg, dx7dg 141 double precision dxmanydg, dxdg, dmtdg, derivg 142 143 double precision dzdt, daldt, dqtdt 144 double precision dx1dt, dx2dt, dx3dt 145 double precision dx5dt 146 double precision dxmanydt, dxdt, dmtdt, derivt 147 double precision afact2 148 double precision rhoval 149 150c functional derivatives above FFFFFFFFFFFF 151 152 parameter(kap=0.8040d0, mu=0.21951d0) 153 parameter (F43=4.d0/3.d0, F13=1.d0/3.d0) 154 parameter (F83=8.d0/3.d0, F23=2.d0/3.d0, F53=5.d0/3.d0) 155 parameter (G920 =9.d0/20.d0 ) 156 157 parameter(b=0.40d0, c=1.59096d0, e=1.537d0) 158 parameter (C1 = 10.d0/81.d0, 159 & C2 = 146.d0/2025.d0, 160 & C3 = -73.d0/405.d0 ) 161c 162 pi=acos(-1d0) 163 Ax = (-0.75d0)*(3d0/pi)**F13 164 P32 = (3.d0*pi**2)**F23 165 es=dsqrt(e) 166 afact2=1d0/facttwo 167c 168 do n = 1, nq 169 if (rho(n).ge.tol_rho) then 170 171c rho43= n*e_x^unif=exchange energy per electron for uniform electron gas 172c = n* Ax*n^(1/3) or n*C*n^(1/3) 173 174 rhoval=rho(n)*facttwo 175 rho43 = Ax*rhoval**F43 ! Ax*n^4/3 176 rrho = 1d0/rhoval ! reciprocal of rho 177 rho13 = rho43*rrho 178 179C Below we just sum up the LDA contribution to the functional 180 Ex = Ex + rho43*qwght(n)*fac*afact2 181 if (ldew) func(n)= func(n) + rho43*fac*afact2 182 Amat(n) = Amat(n) + F43*rho13*fac 183 184c 185 gamma = delrho(n,1)*delrho(n,1) + 186 & delrho(n,2)*delrho(n,2) + 187 & delrho(n,3)*delrho(n,3) 188 gamma=gamma*facttwo*facttwo 189 tauN = tau(n,1)*facttwo 190 tauW=0.125d0*gamma*rrho 191 tauU=0.3d0*P32*rhoval**F53 192c 193c Evaluate the Fx, i.e. mt(x) = Fx - 1 (LDA bit already there) 194c 195 p=gamma/(rhoval**F83*P32*4.d0) 196 if(abs(p).gt.tol_rho) then 197 z=tauW/tauN 198 al=(tauN-tauW)/tauU 199c al=dabs(al) 200 if(al.lt.0d0) al=0d0 201 202 qtil=(G920*(al-1.d0)/((1.d0+b*al*(al-1.d0))**.5d0)) + 203 + F23*p 204 205 xb=(c*z**2)/( (1+z**2)**2 ) 206 x1=(C1 + xb)*p 207 x2=C2*qtil*qtil 208 xc=C3*qtil 209 xd=(0.5d0*(.6d0*z)**2 + .5d0*p*p)**.5d0 210 x3=xc*xd 211 x4=C1*C1*p*p/kap 212 x5=2.d0*es*C1*(.6d0*z)**2 213 x6= e*mu*p*p*p 214 x7 = (1.d0+es*p)**(-2.d0) 215 216 x=(x1+x2+x3 +x4 +x5+x6)*x7 217 218 if (abs(x).lt.tol_rho) write(0,*) ' x for fx ',x 219 220c functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF 221 222C Derivatives wrt n, density below 223 dzdn=-z*rrho 224 dpdn = -p*rrho*F83 225 daldn=F53*( -p*dzdn/z**2 +dpdn*(-1.d0+1.d0/z) ) 226 227 til1=al-1.d0 228 til2=(1.d0+b*al*(al-1.d0))**(-0.5d0) 229 dtil1dn=daldn 230 dtil2dn=b*daldn*(2.d0*al-1d0)* 231 & (-.5d0)*(til2**3) 232 dqtdn = G920*(til2*dtil1dn+til1*dtil2dn)+F23*dpdn 233 234 ax1=c*p*z*z 235 bx1=(1+z*z)**(-2.d0) 236 dx1dn=(x1/p)*dpdn + 2d0*c*p*z*dzdn/((1d0+z*z)**3)*(1d0-z*z) 237 dx2dn=2.d0*C2*qtil*dqtdn 238 239 dxddn=.5d0/xd*( (3d0/5d0)**2*z*dzdn + 240 + p*dpdn) 241 dxcdn=C3*dqtdn 242 dx3dn=xc*dxddn+xd*dxcdn 243 244 dx4dn=(2.d0*x4/p)*dpdn 245 dx5dn=(2.d0*x5/z)*dzdn 246 dx6dn=(3.d0*x6/p)*dpdn 247 dx7dn=-2.d0*es*dpdn/(1.d0+es*p)**3 248 249 xmany=x1+x2+x3 +x4 +x5+x6 250 dxmanydn= dx1dn+dx2dn+dx3dn+dx4dn+dx5dn+dx6dn 251 dxdn=x7*dxmanydn+xmany*dx7dn 252C Derivatives wrt n, density above 253 254C Derivatives wrt gamma, below 255 256 dpdg=p/gamma 257 dzdg=z/gamma 258 daldg=(al/p)*dpdg-F53*(p/(z*z))*dzdg 259 260 dtil2dg=-0.5d0*daldg*b*(2.d0*al-1d0)*til2**3.d0 261 dqtdg=G920*(til1*dtil2dg+til2*daldg)+F23*dpdg 262 dx1dg=(x1/p)*dpdg + 2d0*c*p*z*dzdg/((1d0+z*z)**3)*(1d0-z*z) 263 264 dx2dg=C2*2.d0*qtil*dqtdg 265 266 dxcdg=C3*dqtdg 267 dxddg=.5d0/xd*( (3d0/5d0)**2*z*dzdg + 268 + p*dpdg) 269 dx3dg=xc*dxddg+xd*dxcdg 270 271 dx4dg=(2.d0*x4/p)*dpdg 272 dx5dg=(2.d0*x5/z)*dzdg 273 dx6dg=(3.d0*x6/p)*dpdg 274 275 dx7dg=-2.d0*es*dpdg*(1.d0+p*es)**(-3.d0) 276 277 dxmanydg= dx1dg+dx2dg+dx3dg+dx4dg+dx5dg+dx6dg 278 dxdg=x7*dxmanydg+xmany*dx7dg 279 280C Derivatives wrt tau, below 281c ttttttttttttttttttttttttttttttttttttttttttttttttt 282 dzdt= -z/tauN 283 daldt=1.d0/tauU 284 285 dqtdt=g920*daldt*til2*(1d0- 286 - 0.5d0*b*til1*til2*til2*(2d0*al-1d0)) 287 288 dx1dt=c*p*dzdt*2d0*z*(1d0-z*z)/((1.d0+z*z)**3) 289 dx2dt=2*c2*qtil*dqtdt 290 dx3dt=x3*(dqtdt/qtil + 291 & 0.5d0*(3d0/5d0)**2*z*dzdt/(xd*xd)) 292 dx5dt=2d0*(x5/z)*dzdt 293 294 dxmanydt= dx1dt+dx2dt+dx3dt+dx5dt 295 dxdt=x7*dxmanydt 296c ttttttttttttttttttttttttttttttttttttttttttttttttttt 297 298 mt = kap - kap/(1.d0 + x/kap) 299 300 Ex = Ex + mt*rho43*qwght(n)*fac*afact2 301 if (ldew) func(n)= func(n) + mt*rho43*fac*afact2 302 303 dmtdn=dxdn/(1.d0+x/kap)**2 304 derivn=mt*F43*rho13+rho43*dmtdn 305 306 dmtdg=dxdg/(1.d0+x/kap)**2 307 derivg = rho43*dmtdg 308 309 dmtdt=dxdt/(1.d0+x/kap)**2 310 derivt = rho43*dmtdt 311 Amat(n) = Amat(n) + derivn*fac 312c 313c 4x factor comes from gamma_aa = gamma_total/4 314c 315 Cmat(n)= Cmat(n) + 2d0*derivg*fac 316 Mmat(n)= Mmat(n) +0.5d0*derivt*fac 317 endif 318 endif 319 enddo 320 return 321 end 322 323 Subroutine xc_xtpss03_d2() 324 call errquit(' xtpss03: d2 not coded ',0,0) 325 return 326 end 327 328 329