1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_x_pw6.F 7C> The PW6 exchange functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief Evaluate the PW6 exchange functional 17C> 18C> The PW6 exchange functional [1] is a 6 parameter functional based 19C> on the Perdew-Wang-91 exchange functional and the Becke-95 20C> correlation functional. This subroutine implements just the exchange 21C> part but the parameters were co-optimized with the correlation 22C> functional. Also this subroutine just implements the gradient 23C> correction term (and not the LDA term). 24C> 25C> ### References ### 26C> 27C> [1] Y. Zhao, D. G. Truhlar, 28C> "Design of density functionals that are broadly accurate for 29C> thermochemistry, thermochemical kinetics, and nonbonded 30C> interactions", J. Phys. Chem. A <b>109</b> 5656-5667 (2005), 31C> DOI: 32C> <a href="https://doi.org/10.1021/jp050536c"> 33C> 10.1021/jp050536c</a>. 34C> 35#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 36#if defined(NWAD_PRINT) 37 Subroutine nwxc_x_pw6_p(param, tol_rho, ipol, nq, wght, rho, 38 & rgamma, func) 39#else 40 Subroutine nwxc_x_pw6(param, tol_rho, ipol, nq, wght, rho, rgamma, 41 & func) 42#endif 43#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 44 Subroutine nwxc_x_pw6_d2(param, tol_rho, ipol, nq, wght, 45 & rho, rgamma, func) 46#else 47 Subroutine nwxc_x_pw6_d3(param, tol_rho, ipol, nq, wght, 48 & rho, rgamma, func) 49#endif 50c 51C$Id$ 52c 53#include "nwad.fh" 54c 55 implicit none 56c 57#include "nwxc_param.fh" 58c 59c Input and other parameters 60c 61#if defined(NWAD_PRINT) 62#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 63 type(nwad_dble)::param(*)!< [Input] Parameters of functional 64 type(nwad_dble)::DPOW, BETA, CPW91 65#else 66 double precision param(*)!< [Input] Parameters of functional 67 double precision DPOW, BETA, CPW91 68#endif 69#else 70 double precision param(*)!< [Input] Parameters of functional 71 !< (see [1] Eq.(5) and Table 2): 72 !< - param(1): \f$ b \f$ 73 !< - param(2): \f$ c \f$ 74 !< - param(3): \f$ d \f$ 75 double precision DPOW, BETA, CPW91 76#endif 77 double precision tol_rho !< [Input] The lower limit on the density 78 integer nq !< [Input] The number of points 79 integer ipol !< [Input] The number of spin channels 80 double precision wght !< [Input] The weight of the functional 81c 82c Charge Density 83c 84 type(nwad_dble)::rho(nq,*) !< [Input] The density 85c 86c Charge Density Gradient 87c 88 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 89c 90c Sampling Matrices for the XC Potential & Energy 91c 92 type(nwad_dble)::func(nq) !< [Output] The value of the functional 93c double precision amat(nq,*) !< [Output] The derivative wrt rho 94c double precision cmat(nq,*) !< [Output] The derivative wrt rgamma 95#ifdef SECOND_DERIV 96c double precision Amat2(nq,*) !< [Output] The 2nd derivative wrt rho 97c double precision Cmat2(nq,*) !< [Output] The 2nd derivative wrt rgamma 98c !< and possibly rho 99#endif 100c 101c Compute the partial derivatives of the exchange functional of Perdew91. 102c 103c Becke & Perdew Parameters 104c 105 106c Data BETA/0.00538d0/, CPW91/1.7382d0/,DPOW/3.8901d0/! PW6B95 107 108c 109c*************************************************************************** 110 111 integer n 112 type(nwad_dble)::gamma 113c 114c parameters for PWB6K 115c 116c if (ijzy.eq.2) then 117c BETA = 0.00539d0 118c CPW91= 1.7077d0 119c DPOW= 4.0876d0 120c endif 121c 122 BETA = param(1) 123 CPW91 = param(2) 124 DPOW = param(3) 125c 126 if (ipol.eq.1 )then 127c 128c ======> SPIN-RESTRICTED <====== 129c 130 do 10 n = 1, nq 131 if (rho(n,R_T).lt.tol_rho) goto 10 132 gamma = rgamma(n,G_TT) 133c gamma = delrho(n,1,1)*delrho(n,1,1) + 134c & delrho(n,2,1)*delrho(n,2,1) + 135c & delrho(n,3,1)*delrho(n,3,1) 136 gamma=0.25d0*gamma 137#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 138#if defined(NWAD_PRINT) 139 call nwxc_x_pw6core_p(DPOW,BETA,CPW91,n,1, 140 & rho(n,R_T),gamma,func(n), 141 & tol_rho, wght, 142 & nq, ipol) 143#else 144 call nwxc_x_pw6core(DPOW,BETA,CPW91,n,1, 145 & rho(n,R_T),gamma,func(n), 146 & tol_rho, wght, 147 & nq, ipol) 148#endif 149#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 150 call nwxc_x_pw6core_d2(DPOW,BETA,CPW91,n,1, 151 & rho(n,R_T),gamma,func(n), 152 & tol_rho, wght, 153 & nq, ipol) 154#else 155 call nwxc_x_pw6core_d3(DPOW,BETA,CPW91,n,1, 156 & rho(n,R_T),gamma,func(n), 157 & tol_rho, wght, 158 & nq, ipol) 159#endif 160c write (*,*) "Ex",Ex,Amat(n,1),Cmat(n,1) 161c stop 162 10 continue 163c 164 else 165c 166c ======> SPIN-UNRESTRICTED <====== 167c 168 do 20 n = 1, nq 169 if (rho(n,R_A)+rho(n,R_B).lt.tol_rho) goto 20 170c 171c Spin alpha: 172c 173 if (rho(n,R_A).gt.tol_rho) then 174 gamma = rgamma(n,G_AA) 175c gamma = delrho(n,1,1)*delrho(n,1,1) + 176c & delrho(n,2,1)*delrho(n,2,1) + 177c & delrho(n,3,1)*delrho(n,3,1) 178#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 179#if defined(NWAD_PRINT) 180 call nwxc_x_pw6core_p(DPOW,BETA,CPW91,n,1, 181 & rho(n,R_A),gamma,func(n), 182 & tol_rho, wght, 183 & nq, ipol) 184#else 185 call nwxc_x_pw6core(DPOW,BETA,CPW91,n,1, 186 & rho(n,R_A),gamma,func(n), 187 & tol_rho, wght, 188 & nq, ipol) 189#endif 190#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 191 call nwxc_x_pw6core_d2(DPOW,BETA,CPW91,n,1, 192 & rho(n,R_A),gamma,func(n), 193 & tol_rho, wght, 194 & nq, ipol) 195#else 196 call nwxc_x_pw6core_d3(DPOW,BETA,CPW91,n,1, 197 & rho(n,R_A),gamma,func(n), 198 & tol_rho, wght, 199 & nq, ipol) 200#endif 201 endif 202c 203c Spin beta: 204c 205 if (rho(n,R_B).gt.tol_rho) then 206 gamma = rgamma(n,G_BB) 207c gamma = delrho(n,1,2)*delrho(n,1,2) + 208c & delrho(n,2,2)*delrho(n,2,2) + 209c & delrho(n,3,2)*delrho(n,3,2) 210#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 211#if defined(NWAD_PRINT) 212 call nwxc_x_pw6core_p(DPOW,BETA,CPW91,n,2, 213 & rho(n,R_B),gamma,func(n), 214 & tol_rho, wght, 215 & nq, ipol) 216#else 217 call nwxc_x_pw6core(DPOW,BETA,CPW91,n,2, 218 & rho(n,R_B),gamma,func(n), 219 & tol_rho, wght, 220 & nq, ipol) 221#endif 222#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 223 call nwxc_x_pw6core_d2(DPOW,BETA,CPW91,n,2, 224 & rho(n,R_B),gamma,func(n), 225 & tol_rho, wght, 226 & nq, ipol) 227#else 228 call nwxc_x_pw6core_d3(DPOW,BETA,CPW91,n,2, 229 & rho(n,R_B),gamma,func(n), 230 & tol_rho, wght, 231 & nq, ipol) 232#endif 233 endif 234c 235 20 continue 236c 237 endif 238c 239 return 240 end 241 242 243 244 245#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 246#if defined(NWAD_PRINT) 247 Subroutine nwxc_x_pw6core_p(DPOW,BETA,CPW91,n,ispin, 248 . rho,gamma,func, tol_rho, wght, nq, ipol) 249#else 250 Subroutine nwxc_x_pw6core(DPOW,BETA,CPW91,n,ispin, 251 . rho,gamma,func, tol_rho, wght, nq, ipol) 252#endif 253#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 254 Subroutine nwxc_x_pw6core_d2(DPOW,BETA,CPW91,n,ispin, 255 . rho,gamma,func, tol_rho, wght, nq, ipol) 256#else 257 Subroutine nwxc_x_pw6core_d3(DPOW,BETA,CPW91,n,ispin, 258 . rho,gamma,func, tol_rho, wght, nq, ipol) 259#endif 260c 261C$Id$ 262c 263#include "nwad.fh" 264c 265 implicit none 266c 267#include "nwxc_param.fh" 268c 269 double precision wght 270 integer nq, ipol 271 type(nwad_dble)::func ! value of the functional [output] 272 type(nwad_dble)::rho 273 integer ispin ! alpha=1; beta=2 274c 275c Sampling Matrices for the XC Potential & Energy 276c 277c double precision Amat(nq,ipol), Cmat(nq,*) 278c 279c 280c Compute the partial derivatives of the exchange functional of Perdew91. 281c 282c Becke & Perdew Parameters 283c 284#if defined(NWAD_PRINT) 285#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 286 type(nwad_dble)::DPOW, BETA, CPW91 287#else 288 double precision DPOW, BETA, CPW91 289#endif 290#else 291 double precision DPOW, BETA, CPW91 292#endif 293 double precision tol_rho, AX, pi, BETAPW91,big 294 295 Parameter (big=1d4) 296 297#ifdef SECOND_DERIV 298c 299c Second Derivatives of the Exchange Energy Functional 300c 301c double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 302c double precision rhom23, f2x, d2den,d2num 303#endif 304c 305c References: 306c 307c 308c*************************************************************************** 309 310 integer n 311 type(nwad_dble)::x, sinhm1,dsinhm1 312 type(nwad_dble)::rho13, rho43, Xa, Ha, denom, num, 313 & fprimex, d1num,d1den, 314 & gamma,fx,x2a,ten6xd,expo,bbx2 315 integer D0R,D1G,D2RR,D2RG,D2GG 316 317c 318c sinhm1(x)=log(x+sqrt(1d0+x*x)) 319c dsinhm1(x)=1d0/dsqrt(1d0+x*x) 320c 321 pi=acos(-1.d0) 322 BETAPW91=(pi*36.d0)**(-5.d0/3.d0)*5.d0 323 AX=-(0.75d0/pi)**(1.d0/3.d0)*1.5d0 324 if(ispin.eq.1) then 325 D0R=D1_RA 326 D1G=D1_GAA 327 D2RR=D2_RA_RA 328 D2RG=D2_RA_GAA 329 D2GG=D2_GAA_GAA 330 else 331 D0R=D1_RB 332 D1G=D1_GBB 333 D2RR=D2_RB_RB 334 D2RG=D2_RB_GBB 335 D2GG=D2_GBB_GBB 336 endif 337 338 rho13 = (rho*(ipol/2d0))**(1.d0/3.d0) 339 rho43 = rho13**4.0d0 340 if (gamma.gt.tol_rho**2)then 341 xa = sqrt(gamma)/rho43 342 x2a=xa*xa 343 ten6xd=Xa**DPOW*1.d-6 344 expo=0d0 345 if(CPW91*x2a.lt.big) expo=exp(-CPW91*x2a) 346 bbx2=(BETA-BETAPW91)*x2a 347 Ha = asinh(Xa) 348 denom = 1.d0/(1.d0 + 6d0*(beta*xa)*ha-ten6xd/ax) 349 num = -BETA*x2a+bbx2*expo+ten6xd 350 fx=num*denom 351c d1num=-2.d0*xa*(beta-bbx2*expo*(1d0/x2a-CPW91))+ 352c + ten6xd/xa*dpow 353c d1den=6.d0*beta*(ha + xa*dsinhm1(xa)) - 354c - ten6xd/ax/xa*dpow 355c fprimex=(d1num - d1den*fx)*denom 356 else 357 gamma = 0.d0 358 Xa = 0.d0 359 fx=0.d0 360 fprimex=0.d0 361 denom=0d0 362 d1den=0d0 363 expo=0d0 364 ten6xd=0d0 365 x2a=0d0 366 bbx2=0d0 367 endif 368c 369c if (lfac)then 370c Ex = Ex + 2d0/ipol*rho43*AX*qwght*fac 371c func = func + 2.d0/ipol*rho43*AX*fac 372c Amat(n,D0R) = Amat(n,D0R) + (4.d0/3.d0)*rho13*AX*fac 373c endif 374c 375c if (nlfac)then 376c Ex = Ex + 2d0/ipol*rho43*fx*qwght*fac 377 func = func + 2.d0/ipol*rho43*fx*wght 378c Amat(n,D0R) = Amat(n,D0R) + 379c + (4.d0/3.d0)*rho13*(fx-xa*fprimex)*wght 380c if (xa.gt.tol_rho) then 381c Cmat(n,D1G)=Cmat(n,D1G)+ 382c + .5d0*fprimex/sqrt(gamma)*wght 383c endif 384c endif 385c 386#ifdef SECOND_DERIV 387c rhom23 = 2d0/ipol*rho13/rho 388c if (lfac)then 389c Amat2(n,D2RR) = Amat2(n,D2RR) + 390c & (4d0/9d0)*rhom23*Ax*fac 391c endif 392c if (nlfac)then 393c if(gamma.gt.tol_rho**2)then 394c d2num=-2d0*beta + 395c + 2d0*bbx2*expo* 396c * (1d0/x2a-5d0*CPW91+2d0*CPW91*CPW91*x2a)+ 397c + ten6xd/x2a*dpow*(dpow-1d0) 398c d2den=6.d0*beta*dsinhm1(xa)*(2d0-x2a/(1d0+x2a)) - 399c - ten6xd/ax/x2a*dpow*(dpow-1d0) 400c f2x = denom*(d2num -fx*d2den - 2d0*d1den*fprimex) 401c else 402c f2x=0d0 403c endif 404c Amat2(n,D2RR) = Amat2(n,D2RR) 405c & + (4d0/9d0)*rhom23*(fx-xa*fprimex+4d0*x2a*f2x)*wght 406c Cmat2(n,D2RG) = Cmat2(n,D2RG) 407c & - (4d0/ipol/3d0)*(rhom23**2/rho)*f2x*wght 408c if (xa.gt.tol_rho) then 409c Cmat2(n,D2GG) = Cmat2(n,D2GG) 410c & - 0.25d0*gamma**(-1.5d0)*(fprimex-xa*f2x)*wght 411c endif 412c endif 413#endif 414c 415c 416 return 417 end 418#ifndef NWAD_PRINT 419#define NWAD_PRINT 420c 421c Compile source again for Maxima 422c 423#include "nwxc_x_pw6.F" 424#endif 425#ifndef SECOND_DERIV 426#define SECOND_DERIV 427c 428c Compile source again for the 2nd derivative case 429c 430#include "nwxc_x_pw6.F" 431#endif 432#ifndef THIRD_DERIV 433#define THIRD_DERIV 434c 435c Compile source again for the 3rd derivative case 436c 437#include "nwxc_x_pw6.F" 438#endif 439#undef NWAD_PRINT 440C> 441C> @} 442