1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_x_tpss03.F 7C> The TPSS exchange functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief Evaluate the TPSS exchange functional 17C> 18C> Evaluate the TPSS meta-GGA exchange functional [1,2]. 19C> 20C> Due to the form of the meta-GGAs we need to screen on the kinetic 21C> energy density to ensure that LDA will be obtained when the kinetic 22C> energy density goes to zero [3]. 23C> 24C> ### References ### 25C> 26C> [1] J. Tao, J.P. Perdew, V.N. Staveroverov, G.E. Scuseria, 27C> "Climbing the density functional ladder: Nonemperical meta- 28C> generalized gradient approximation designed for molecules 29C> and solids", 30C> Phys. Rev. Lett. <b>91</b>, 146401-146404 (2003), DOI: 31C> <a href="https://doi.org/10.1103/PhysRevLett.91.146401"> 32C> 10.1103/PhysRevLett.91.146401</a>. 33C> 34C> [2] J.P. Perdew, J. Tao, V.N. Staveroverov, G.E. Scuseria, 35C> "Meta-generalized gradient approximation: Explanation of a 36C> realistic nonempirical density functional", 37C> J. Chem. Phys. <b>120</b>, 6898-6911 (2004), DOI: 38C> <a href="https://doi.org/10.1063/1.1665298"> 39C> 10.1103/1.1665298</a>. 40C> 41C> [3] J. Gräfenstein, D. Izotov, D. Cremer, 42C> "Avoiding singularity problems associated with meta-GGA exchange 43C> and correlation functionals containing the kinetic energy 44C> density", J. Chem. Phys. <b>127</b>, 214103 (2007), DOI: 45C> <a href="https://doi.org/10.1063/1.2800011"> 46C> 10.1063/1.2800011</a>. 47C> 48c 49c$Id$ 50c 51c Tao,Perdew, Staroverov, Scuseria exchange functional 52c META GGA 53C utilizes ingredients: 54c rho - density 55c delrho - gradient of density 56c tau - K.S kinetic energy density 57c tauW - von Weiszacker kinetic energy density 58c tauU - uniform-gas KE density 59c References: 60c [a] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria 61c PRL 91, 146401 (2003). 62c [b] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria 63c JCP 120, 6898 (2004). 64#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 65#if defined(NWAD_PRINT) 66 Subroutine nwxc_x_tpss03_p(tol_rho, ipol, nq, wght, 67 & rho, rgamma, tau, func) 68#else 69 Subroutine nwxc_x_tpss03(tol_rho, ipol, nq, wght, 70 & rho, rgamma, tau, func) 71#endif 72#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 73 Subroutine nwxc_x_tpss03_d2(tol_rho, ipol, nq, wght, 74 & rho, rgamma, tau, func) 75#else 76 Subroutine nwxc_x_tpss03_d3(tol_rho, ipol, nq, wght, 77 & rho, rgamma, tau, func) 78#endif 79c 80#include "nwad.fh" 81c 82 implicit none 83c 84#include "nwxc_param.fh" 85c 86c Input and other parameters 87c 88 double precision tol_rho !< [Input] The lower limit on the density 89 integer nq !< [Input] The number of points 90 integer ipol !< [Input] The number of spin channels 91 double precision wght !< [Input] The weight of the functional 92c 93c Charge Density 94c 95 type(nwad_dble)::rho(nq,*) !< [Input] The density 96c 97c Charge Density Gradient 98c 99 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 100c 101c Kinetic Energy Density 102c 103 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 104c 105c The functional 106c 107 type(nwad_dble)::func(*) !< [Output] The value of the functional 108c 109c Sampling Matrices for the XC Potential & Energy 110c 111c double precision Amat(nq,*) !< [Output] The derivative wrt rho 112c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 113c double precision Mmat(nq,*) !< [Output] The derivative wrt tau 114c 115 integer ispin,cmatpos 116c 117 if (ipol.eq.1 )then 118c 119c SPIN-RESTRICTED 120c Ex = Ex[n] 121c 122#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 123#if defined(NWAD_PRINT) 124 call nwxc_x_tpss03_cs_p(1.0d0, tol_rho, ipol, nq, wght, 125 & rho, rgamma, tau, func) 126#else 127 call nwxc_x_tpss03_cs(1.0d0, tol_rho, ipol, nq, wght, 128 & rho, rgamma, tau, func) 129#endif 130#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 131 call nwxc_x_tpss03_cs_d2(1.0d0, tol_rho, ipol, nq, wght, 132 & rho, rgamma, tau, func) 133#else 134 call nwxc_x_tpss03_cs_d3(1.0d0, tol_rho, ipol, nq, wght, 135 & rho, rgamma, tau, func) 136#endif 137 else 138c 139c SPIN-UNRESTRICTED 140c Ex = Ex[2n_up]/2 + Ex[2n_down]/2 141 142#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 143#if defined(NWAD_PRINT) 144 call nwxc_x_tpss03_cs_p(2.0d0, tol_rho, ipol, nq, wght, 145 & rho(1,R_A), rgamma(1,G_AA), tau(1,T_A), 146 & func) 147 call nwxc_x_tpss03_cs_p(2.0d0, tol_rho, ipol, nq, wght, 148 & rho(1,R_B), rgamma(1,G_BB), tau(1,T_B), 149 & func) 150#else 151 call nwxc_x_tpss03_cs(2.0d0, tol_rho, ipol, nq, wght, 152 & rho(1,R_A), rgamma(1,G_AA), tau(1,T_A), 153 & func) 154 call nwxc_x_tpss03_cs(2.0d0, tol_rho, ipol, nq, wght, 155 & rho(1,R_B), rgamma(1,G_BB), tau(1,T_B), 156 & func) 157#endif 158#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 159 call nwxc_x_tpss03_cs_d2(2.0d0, tol_rho, ipol, nq, wght, 160 & rho(1,R_A), rgamma(1,G_AA), tau(1,T_A), 161 & func) 162 call nwxc_x_tpss03_cs_d2(2.0d0, tol_rho, ipol, nq, wght, 163 & rho(1,R_B), rgamma(1,G_BB), tau(1,T_B), 164 & func) 165#else 166 call nwxc_x_tpss03_cs_d3(2.0d0, tol_rho, ipol, nq, wght, 167 & rho(1,R_A), rgamma(1,G_AA), tau(1,T_A), 168 & func) 169 call nwxc_x_tpss03_cs_d3(2.0d0, tol_rho, ipol, nq, wght, 170 & rho(1,R_B), rgamma(1,G_BB), tau(1,T_B), 171 & func) 172#endif 173 endif 174 return 175 end 176#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 177#if defined(NWAD_PRINT) 178 Subroutine nwxc_x_tpss03_cs_p(facttwo, tol_rho, ipol, nq, wght, 179 & rho, rgamma, tau, func) 180#else 181 Subroutine nwxc_x_tpss03_cs(facttwo, tol_rho, ipol, nq, wght, 182 & rho, rgamma, tau, func) 183#endif 184#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 185 Subroutine nwxc_x_tpss03_cs_d2(facttwo, tol_rho, ipol, nq, wght, 186 & rho, rgamma, tau, func) 187#else 188 Subroutine nwxc_x_tpss03_cs_d3(facttwo, tol_rho, ipol, nq, wght, 189 & rho, rgamma, tau, func) 190#endif 191#include "nwad.fh" 192 implicit none 193c 194c Input and other parameters 195c 196 double precision facttwo !< [Input] Scale factor 197 !< - 1 for closed shell calculations 198 !< - 2 for open shell calculations 199 double precision tol_rho !< [Input] The lower limit on the density 200 integer nq !< [Input] The number of points 201 integer ipol !< [Input] The number of spin channels 202 double precision wght !< [Input] The weight of the functional 203c 204c Charge Density 205c 206 type(nwad_dble)::rho(nq) !< [Input] The density 207c 208c Charge Density Gradient 209c 210 type(nwad_dble)::rgamma(nq) !< [Input] The norm of the density gradients 211c 212c Kinetic Energy Density 213c 214 type(nwad_dble)::tau(nq) !< [Input] The kinetic energy density 215c 216c The functional 217c 218 type(nwad_dble)::func(*) !< [Output] The value of the functional 219c 220c Sampling Matrices for the XC Potential & Energy 221c 222c double precision Amat(nq) !< [Output] The derivative wrt rho 223c double precision Cmat(nq) !< [Output] The derivative wrt rgamma 224c double precision Mmat(nq) !< [Output] The derivative wrt tau 225c 226 double precision pi 227 integer n 228 type(nwad_dble)::rrho, rho43, gamma 229 type(nwad_dble)::tauN, tauW, tauU 230 231 type(nwad_dble):: p, qtil, x, al, mt, z 232 double precision F83, F23, F53, F43, F13 233 double precision G920 234 double precision b,c,e,es 235 double precision C1, C2, C3 236 double precision kap, mu 237 type(nwad_dble)::xb,xc,xd 238 type(nwad_dble)::x1,x2,x3,x4,x5,x6,x7 239 double precision P32, Ax 240c functional derivatives below FFFFFFFFFFFF 241 double precision dzdn, dpdn, daldn, dqtdn 242 double precision til1, til2 243 double precision dtil2dn, dtil1dn 244 double precision ax1, bx1, dx1dn 245 double precision dx2dn 246 double precision dxddn, dxcdn, dx3dn 247 double precision dx4dn, dx5dn, dx6dn, dx7dn 248 double precision dxdn 249 double precision xmany, dxmanydn 250 double precision dmtdn, derivn 251 252 double precision dzdg, dpdg, daldg, dqtdg 253 double precision dtil2dg 254 double precision dx1dg, dx2dg 255 double precision dxcdg, dxddg,dx3dg 256 double precision dx4dg, dx5dg, dx6dg, dx7dg 257 double precision dxmanydg, dxdg, dmtdg, derivg 258 259 double precision dzdt, daldt, dqtdt 260 double precision dx1dt, dx2dt, dx3dt 261 double precision dx5dt 262 double precision dxmanydt, dxdt, dmtdt, derivt 263 double precision afact2 264 type(nwad_dble)::rhoval 265 266c functional derivatives above FFFFFFFFFFFF 267 268 parameter(kap=0.8040d0, mu=0.21951d0) 269 parameter (F43=4.d0/3.d0, F13=1.d0/3.d0) 270 parameter (F83=8.d0/3.d0, F23=2.d0/3.d0, F53=5.d0/3.d0) 271 parameter (G920 =9.d0/20.d0 ) 272 273 parameter(b=0.40d0, c=1.59096d0, e=1.537d0) 274 parameter (C1 = 10.d0/81.d0, 275 & C2 = 146.d0/2025.d0, 276 & C3 = -73.d0/405.d0 ) 277c 278 pi=acos(-1d0) 279 Ax = (-0.75d0)*(3d0/pi)**F13 280 P32 = (3.d0*pi**2)**F23 281 es=dsqrt(e) 282 afact2=1d0/facttwo 283c 284 do n = 1, nq 285 rhoval=rho(n)*facttwo 286 if (rhoval.ge.tol_rho) then 287 288c rho43= n*e_x^unif=exchange energy per electron for uniform electron gas 289c = n* Ax*n^(1/3) or n*C*n^(1/3) 290 291 rho43 = Ax*rhoval**F43 ! Ax*n^4/3 292 rrho = 1d0/rhoval ! reciprocal of rho 293c rho13 = rho43*rrho 294 295C Below we just sum up the LDA contribution to the functional 296 func(n)= func(n) + rho43*wght*afact2 297c Amat(n) = Amat(n) + F43*rho13*wght 298 299c 300c gamma = delrho(n,1)*delrho(n,1) + 301c & delrho(n,2)*delrho(n,2) + 302c & delrho(n,3)*delrho(n,3) 303 gamma=rgamma(n) 304 gamma=gamma*facttwo*facttwo 305 tauN = tau(n)*facttwo 306 tauW=0.125d0*gamma*rrho 307 tauU=0.3d0*P32*rhoval**F53 308c 309c Evaluate the Fx, i.e. mt(x) = Fx - 1 (LDA bit already there) 310c 311 p=gamma/(rhoval**F83*P32*4.d0) 312 if (tauN.ge.tol_rho) then 313 z=tauW/tauN 314 else 315 z=0.0d0 316 endif 317 al=(tauN-tauW)/tauU 318c al=dabs(al) 319 if(al.lt.0d0) al=0d0 320 321 qtil=(G920*(al-1.d0)/((1.d0+b*al*(al-1.d0))**.5d0)) + 322 + F23*p 323 324 xb=(c*z**2)/( (1.0d0+z**2)**2) 325 x1=(C1 + xb)*p 326 x2=C2*qtil*qtil 327 xc=C3*qtil 328 if (gamma.lt.tol_rho) then 329 xd = 0.0d0 330 else 331 xd=(0.5d0*(.6d0*z)**2 + .5d0*p**2)**.5d0 332 endif 333 x3=xc*xd 334 x4=C1*C1*p**2/kap 335 x5=2.d0*es*C1*(.6d0*z)**2 336 x6= e*mu*p**3 337 x7 = (1.d0+es*p)**(-2) 338 339 x=(x1+x2+x3 +x4 +x5+x6)*x7 340 341 if (abs(x).lt.tol_rho) write(0,*) ' x for fx ',x 342 343c functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF 344 345C Derivatives wrt n, density below 346c dzdn=-z*rrho 347c dpdn = -p*rrho*F83 348c daldn=F53*( -p*dzdn/z**2 +dpdn*(-1.d0+1.d0/z) ) 349c 350c til1=al-1.d0 351c til2=(1.d0+b*al*(al-1.d0))**(-0.5d0) 352c dtil1dn=daldn 353c dtil2dn=b*daldn*(2.d0*al-1d0)* 354c & (-.5d0)*(til2**3) 355c dqtdn = G920*(til2*dtil1dn+til1*dtil2dn)+F23*dpdn 356c 357c ax1=c*p*z*z 358c bx1=(1+z*z)**(-2.d0) 359c dx1dn=(x1/p)*dpdn + 2d0*c*p*z*dzdn/((1d0+z*z)**3)*(1d0-z*z) 360c dx2dn=2.d0*C2*qtil*dqtdn 361c 362c dxddn=.5d0/xd*( (3d0/5d0)**2*z*dzdn + 363c + p*dpdn) 364c dxcdn=C3*dqtdn 365c dx3dn=xc*dxddn+xd*dxcdn 366c 367c dx4dn=(2.d0*x4/p)*dpdn 368c dx5dn=(2.d0*x5/z)*dzdn 369c dx6dn=(3.d0*x6/p)*dpdn 370c dx7dn=-2.d0*es*dpdn/(1.d0+es*p)**3 371c 372c xmany=x1+x2+x3 +x4 +x5+x6 373c dxmanydn= dx1dn+dx2dn+dx3dn+dx4dn+dx5dn+dx6dn 374c dxdn=x7*dxmanydn+xmany*dx7dn 375C Derivatives wrt n, density above 376 377C Derivatives wrt gamma, below 378 379c dpdg=p/gamma 380c dzdg=z/gamma 381c daldg=(al/p)*dpdg-F53*(p/(z*z))*dzdg 382c 383c dtil2dg=-0.5d0*daldg*b*(2.d0*al-1d0)*til2**3.d0 384c dqtdg=G920*(til1*dtil2dg+til2*daldg)+F23*dpdg 385c dx1dg=(x1/p)*dpdg + 2d0*c*p*z*dzdg/((1d0+z*z)**3)*(1d0-z*z) 386c 387c dx2dg=C2*2.d0*qtil*dqtdg 388c 389c dxcdg=C3*dqtdg 390c dxddg=.5d0/xd*( (3d0/5d0)**2*z*dzdg + 391c + p*dpdg) 392c dx3dg=xc*dxddg+xd*dxcdg 393c 394c dx4dg=(2.d0*x4/p)*dpdg 395c dx5dg=(2.d0*x5/z)*dzdg 396c dx6dg=(3.d0*x6/p)*dpdg 397c 398c dx7dg=-2.d0*es*dpdg*(1.d0+p*es)**(-3.d0) 399c 400c dxmanydg= dx1dg+dx2dg+dx3dg+dx4dg+dx5dg+dx6dg 401c dxdg=x7*dxmanydg+xmany*dx7dg 402 403C Derivatives wrt tau, below 404c ttttttttttttttttttttttttttttttttttttttttttttttttt 405c dzdt= -z/tauN 406c daldt=1.d0/tauU 407c 408c dqtdt=g920*daldt*til2*(1d0- 409c - 0.5d0*b*til1*til2*til2*(2d0*al-1d0)) 410c 411c dx1dt=c*p*dzdt*2d0*z*(1d0-z*z)/((1.d0+z*z)**3) 412c dx2dt=2*c2*qtil*dqtdt 413c dx3dt=x3*(dqtdt/qtil + 414c & 0.5d0*(3d0/5d0)**2*z*dzdt/(xd*xd)) 415c dx5dt=2d0*(x5/z)*dzdt 416c 417c dxmanydt= dx1dt+dx2dt+dx3dt+dx5dt 418c dxdt=x7*dxmanydt 419c ttttttttttttttttttttttttttttttttttttttttttttttttttt 420 421 mt = kap - kap/(1.d0 + x/kap) 422 423 func(n)= func(n) + mt*rho43*wght*afact2 424 425c dmtdn=dxdn/(1.d0+x/kap)**2 426c derivn=mt*F43*rho13+rho43*dmtdn 427 428c dmtdg=dxdg/(1.d0+x/kap)**2 429c derivg = rho43*dmtdg 430c 431c dmtdt=dxdt/(1.d0+x/kap)**2 432c derivt = rho43*dmtdt 433c Amat(n) = Amat(n) + derivn*wght 434c 435c 4x factor comes from gamma_aa = gamma_total/4 436c 437c Cmat(n)= Cmat(n) + 2d0*derivg*wght 438c Mmat(n)= Mmat(n) +0.5d0*derivt*wght 439 endif 440 enddo 441 return 442 end 443#ifndef NWAD_PRINT 444#define NWAD_PRINT 445c 446c Compile source again for the 2nd derivative case 447c 448#include "nwxc_x_tpss03.F" 449#endif 450#ifndef SECOND_DERIV 451#define SECOND_DERIV 452c 453c Compile source again for the 2nd derivative case 454c 455#include "nwxc_x_tpss03.F" 456#endif 457#ifndef THIRD_DERIV 458#define THIRD_DERIV 459c 460c Compile source again for the 3rd derivative case 461c 462#include "nwxc_x_tpss03.F" 463#endif 464#undef NWAD_PRINT 465C> 466C> @} 467