1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_vs98.F 7C> Implementation of the VS98 correlation functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief The Voorhis and Scuseria correlation functional 17C> 18C> The VS98 functional [1,2] is a meta-GGA. This routine implements 19C> the correlation component. 20C> 21C> Due to the form of the meta-GGAs we need to screen on the kinetic 22C> energy density to ensure that LDA will be obtained when the kinetic 23C> energy density goes to zero [3]. 24C> 25C> ### References ### 26C> 27C> [1] T van Voorhis, GE Scuseria, 28C> "A novel form for the exchange-correlation energy functional", 29C> J.Chem.Phys. <b>109</b>, 400-410 (1998), DOI: 30C> <a href="https://doi.org/10.1063/1.476577"> 31C> 10.1063/1.476577</a>. 32C> 33C> [2] T van Voorhis, GE Scuseria, 34C> Erratum: "A novel form for the exchange-correlation energy 35C> functional", 36C> J.Chem.Phys. <b>129</b>, 219901-219901 (2008), DOI: 37C> <a href="https://doi.org/10.1063/1.3005348"> 38C> 10.1063/1.3005348</a>. 39C> 40C> [3] J. Gräfenstein, D. Izotov, D. Cremer, 41C> "Avoiding singularity problems associated with meta-GGA exchange 42C> and correlation functionals containing the kinetic energy 43C> density", J. Chem. Phys. <b>127</b>, 214103 (2007), DOI: 44C> <a href="https://doi.org/10.1063/1.2800011"> 45C> 10.1063/1.2800011</a>. 46C> 47c VSXC correlation functional 48c META GGA 49C utilizes ingredients: 50c rho - density 51c delrho - gradient of density 52c tau (tauN)- K.S kinetic energy density 53c ijzy - 1 VS98 54c ijzy - 2 M06-L 55c ijzy - 3 M06-HF 56c ijzy - 4 M06 57c ijzy - 5 M06-2X 58c 59#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 60#if defined(NWAD_PRINT) 61 Subroutine nwxc_c_vs98_p(param, tol_rho, ipol, nq, wght, rho, 62 & rgamma, tau, func) 63#else 64 Subroutine nwxc_c_vs98(param, tol_rho, ipol, nq, wght, rho, 65 & rgamma, tau, func) 66#endif 67#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 68 Subroutine nwxc_c_vs98_d2(param, tol_rho, ipol, nq, wght, rho, 69 & rgamma, tau, func) 70#else 71 Subroutine nwxc_c_vs98_d3(param, tol_rho, ipol, nq, wght, rho, 72 & rgamma, tau, func) 73#endif 74c 75c$Id$ 76c 77c Reference 78c [a] T. V. Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998). 79c [b] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006). 80 81c 82#include "nwad.fh" 83c 84 implicit none 85c 86#include "intf_nwxc_c_lsda.fh" 87#include "intf_nwxc_vs98ss.fh" 88c 89#include "nwxc_param.fh" 90c 91c Input and other parameters 92c 93#if defined(NWAD_PRINT) 94#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 95 type(nwad_dble)::param(*) !< [Input] Parameters of functional (table 1) 96 type(nwad_dble)::r7, r8, r9, r10, r11, r12 97 type(nwad_dble)::r13, r14, r15, r16, r17, r18 98#else 99 double precision param(*) !< [Input] Parameters of functional (table 1) 100 double precision r7, r8, r9, r10, r11, r12 101 double precision r13, r14, r15, r16, r17, r18 102#endif 103#else 104 double precision param(*) !< [Input] Parameters of functional (table 1) 105 !< - param(1): \f$ a_{\sigma\sigma'} \f$ 106 !< - param(2): \f$ b_{\sigma\sigma'} \f$ 107 !< - param(3): \f$ c_{\sigma\sigma'} \f$ 108 !< - param(4): \f$ d_{\sigma\sigma'} \f$ 109 !< - param(5): \f$ e_{\sigma\sigma'} \f$ 110 !< - param(6): \f$ f_{\sigma\sigma'} \f$ 111 !< - param(7): \f$ a_{\sigma\sigma} \f$ 112 !< - param(8): \f$ b_{\sigma\sigma} \f$ 113 !< - param(9): \f$ c_{\sigma\sigma} \f$ 114 !< - param(10): \f$ d_{\sigma\sigma} \f$ 115 !< - param(11): \f$ e_{\sigma\sigma} \f$ 116 !< - param(12): \f$ f_{\sigma\sigma} \f$ 117 double precision r7, r8, r9, r10, r11, r12 118 double precision r13, r14, r15, r16, r17, r18 119#endif 120 double precision tol_rho !< [Input] The lower limit on the density 121 integer nq !< [Input] The number of points 122 integer ipol !< [Input] The number of spin channels 123 double precision wght !< [Input] The weight of the functional 124c 125c Charge Density 126c 127 type(nwad_dble)::rho(nq,*) !< [Input] The density 128c 129c Charge Density Gradient Norm 130c 131 type(nwad_dble)::rgamma(nq,*) !< [Input] The density gradient norm 132c 133c Kinetic Energy Density 134c 135 type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density 136c 137c Functional values 138c 139 type(nwad_dble)::func(*) !< [Output] The functional value 140c 141c Sampling Matrices for the XC Potential 142c 143c double precision Amat(nq,*) !< [Output] Derivative wrt density 144c double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma 145c double precision Mmat(nq,*) !< [Output] Derivative wrt tau 146c 147c Threshold parameters 148c 149 double precision DTol,F1, F2, F3, F4, gab, cf 150 Data F1/1.0d0/,F2/2.0d0/, 151 & F3/3.0d0/,F4/4.0d0/,gab/0.00304966d0/, 152 & cf/9.115599720d0/ 153c 154c Local 155c 156 integer n 157 double precision tauN 158 159c call to the vs98css subroutine 160 type(nwad_dble)::PA,GAA,TA,FA,EUA,ZA,ChiA 161 double precision FPA,FGA,FTA,EUEGA,EUPA,ChiAP,ChiAG, 162 & ZAP,ZAT 163 type(nwad_dble)::PB,GBB,TB,FB,EUB,ZB,ChiB 164 double precision FPB,FGB,FTB,EUEGB,EUPB,ChiBP,ChiBG, 165 & ZBP,ZBT 166c 167 type(nwad_dble)::RS,Zeta,gcab,P,PotLC,EUEG,ZAB 168 double precision Pi, F43, F13, Pi34, F6, 169 &RSP,dZdA,dZdB,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ 170 type(nwad_dble)::XAB, kab, xk, zk 171 double precision dgdx,dgdz,dgdPA,dgdGA,dgdTA,dgdPB,dgdGB,dgdTB 172 double precision EUEGPA,EUEGPB 173 174 175c 176c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 177c 178c DTol=1.0d-7 179 dtol=tol_rho 180C Parameters for VS98 / M06-L / M06-HF / M06 / M06-2X 181 r7= param(1) 182 r8= param(2) 183 r9= param(3) 184 r10= param(4) 185 r11= param(5) 186 r12= param(6) 187 r13= param(7) 188 r14= param(8) 189 r15= param(9) 190 r16= param(10) 191 r17= param(11) 192 r18= param(12) 193c 194 Pi = F4*ATan(F1) 195 F6=6.0d0 196 F43 = F4 / F3 197 Pi34 = F3 / (F4*Pi) 198 F13 = F1 / F3 199c 200 do 20 n = 1, nq 201 if (ipol.eq.1) then 202 if (rho(n,R_A).lt.DTol) goto 20 203 else 204 if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20 205 endif 206 if (ipol.eq.1) then 207c 208c get the density, gradient, and tau for the alpha spin 209c from the total 210c 211 PA = rho(n,R_A)/F2 212 PB = 0.0d0 213 EUA = 0.0d0 214 ChiA = 0.0d0 215 EUB = 0.0d0 216 ChiB = 0.0d0 217 GAA = rgamma(n,G_AA)/4.0d0 218c In the bc95css subroutine, we use 2*TA as the tau, so we do 219c not divide the tau by 2 here 220 221 TA = tau(n,T_A) 222 if(ta.lt.dtol) goto 20 223 224#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 225#if defined(NWAD_PRINT) 226 Call nwxc_vs98ss_p(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 227 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 228 & r13,r14,r15,r16,r17,r18) 229#else 230 Call nwxc_vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 231 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 232 & r13,r14,r15,r16,r17,r18) 233#endif 234#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 235 Call nwxc_vs98ss_d2(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 236 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 237 & r13,r14,r15,r16,r17,r18) 238#else 239 Call nwxc_vs98ss_d3(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 240 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 241 & r13,r14,r15,r16,r17,r18) 242#endif 243 PB = PA 244 GBB = GAA 245 TB = TA 246 FB = FA 247 FPB = FPA 248 FGB = FGA 249 FTB = FTA 250 EUB = EUA 251 ZB = ZA 252 ChiB = ChiA 253 EUPB = EUPA 254 ChiBP = ChiAP 255 ChiBG = ChiAG 256 ZBP = ZAP 257 ZBT = ZAT 258 259 func(n) = func(n) + 2.0d0*FA*wght 260c Amat(n,D1_RA) = Amat(n,D1_RA) + FPA*wght 261c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + FGA*wght 262c Mmat(n,D1_TA) = Mmat(n,D1_TA) + FTA*wght 263 264 265c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 266 else ! ipol=2 267c 268c ======> SPIN-UNRESTRICTED <====== 269c 270c 271c alpha 272c 273 274 PA = 0.0d0 275 GAA = 0.0d0 276 TA = 0.0d0 277 EUA = 0.0d0 278 ChiA = 0.0d0 279 ZA = 0.0d0 280 if (rho(n,R_A).le.0.5d0*DTol) go to 25 281 PA = rho(n,R_A) 282c 283c In the bc95css subroutine, we use 2*TA as the tau 284c 285 if (tau(n,T_A).le.0.5d0*DTol) go to 25 286 GAA = rgamma(n,G_AA) 287 TA = tau(n,T_A)*2.0d0 288 289#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 290#if defined(NWAD_PRINT) 291 Call nwxc_vs98ss_p(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 292 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 293 & r13,r14,r15,r16,r17,r18) 294#else 295 Call nwxc_vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 296 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 297 & r13,r14,r15,r16,r17,r18) 298#endif 299#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 300 Call nwxc_vs98ss_d2(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 301 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 302 & r13,r14,r15,r16,r17,r18) 303#else 304 Call nwxc_vs98ss_d3(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA, 305 & ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT, 306 & r13,r14,r15,r16,r17,r18) 307#endif 308 func(n) = func(n) + FA*wght 309c Amat(n,D1_RA) = Amat(n,D1_RA) + FPA*wght 310c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + FGA*wght 311c 2*0.5=1.0 for Mmat 312c Mmat(n,D1_TA) = Mmat(n,D1_TA) + FTA*wght 313 314c 315c In the vs98ss subroutine, we use 2*TA as the tau, 316c 317c 318c Beta 319c 320 25 continue 321 PB = 0.0d0 322 GBB = 0.0d0 323 TB = 0.0d0 324 EUB = 0.0d0 325 ChiB = 0.0d0 326 ZB = 0.0d0 327 328 if(rho(n,R_B).le.0.5d0*DTol) go to 30 329 330 PB = rho(n,R_B) 331 332 if(tau(n,T_B).le.0.5d0*DTol) go to 30 333 GBB = rgamma(n,G_BB) 334 TB = tau(n,T_B)*2.0d0 335 336#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 337#if defined(NWAD_PRINT) 338 Call nwxc_vs98ss_p(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB, 339 & ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT, 340 & r13,r14,r15,r16,r17,r18) 341#else 342 Call nwxc_vs98ss(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB, 343 & ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT, 344 & r13,r14,r15,r16,r17,r18) 345#endif 346#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 347 Call nwxc_vs98ss_d2(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB, 348 & ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT, 349 & r13,r14,r15,r16,r17,r18) 350#else 351 Call nwxc_vs98ss_d3(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB, 352 & ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT, 353 & r13,r14,r15,r16,r17,r18) 354#endif 355 func(n) = func(n) + FB*wght 356c Amat(n,D1_RB) = Amat(n,D1_RB) + FPB*wght 357c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + FGB*wght 358c Mmat(n,D1_TB) = Mmat(n,D1_TB) + FTB*wght 359 360#if 0 361 write (0,'(A,3F20.6)') "BAmat Cmat Mmat",FPB,FGB,FTB 362#endif 363 endif 364 30 continue 365 P = PA + PB 366 If(PA.gt.DTol.or.PB.gt.DTol) then 367 RS = (Pi34/P) ** F13 368c RSP = -RS/(F3*P) 369 Zeta = (PA-PB)/P 370c dZdA = (F1-Zeta)/P 371c dZdB = (-F1-Zeta)/P 372#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 373#if defined(NWAD_PRINT) 374 Call nwxc_c_lsda_p(tol_rho,RS,Zeta,PotLC, 375 & dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 376#else 377 Call nwxc_c_lsda(tol_rho,RS,Zeta,PotLC, 378 & dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 379#endif 380#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 381 Call nwxc_c_lsda_d2(tol_rho,RS,Zeta,PotLC, 382 & dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 383#else 384 Call nwxc_c_lsda_d3(tol_rho,RS,Zeta,PotLC, 385 & dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 386#endif 387 EUEG = P*PotLC - EUA - EUB 388 ZAB = ZA + ZB 389 XAB = ChiA+ChiB 390 kab = F1 + gab*(XAB+ZAB) 391 xk = XAB/kab 392 zk = ZAB/kab 393#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 394#if defined(NWAD_PRINT) 395 call nwxc_gvt4_p(gcab,dgdx,dgdz,xk,zk,kab,gab,gab, 396 & r7,r8,r9,r10,r11,r12) 397#else 398 call nwxc_gvt4(gcab,dgdx,dgdz,xk,zk,kab,gab,gab, 399 & r7,r8,r9,r10,r11,r12) 400#endif 401#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 402 call nwxc_gvt4_d2(gcab,dgdx,dgdz,xk,zk,kab,gab,gab, 403 & r7,r8,r9,r10,r11,r12) 404#else 405 call nwxc_gvt4_d3(gcab,dgdx,dgdz,xk,zk,kab,gab,gab, 406 & r7,r8,r9,r10,r11,r12) 407#endif 408 func(n) = func(n) + gcab*EUEG*wght 409c dgdPA = dgdx*ChiAP + dgdz*ZAP 410c dgdGA = dgdx*ChiAG 411c dgdTA = dgdz*ZAT 412c dgdPB = dgdx*ChiBP + dgdz*ZBP 413c dgdGB = dgdx*ChiBG 414c dgdTB = dgdz*ZBT 415c EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA 416c EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB 417c if (ipol.eq.1) then 418c Amat(n,D1_RA) = Amat(n,D1_RA) 419c & + (EUEGPA*gcab + EUEG*dgdPA)*wght 420c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + EUEG*dgdGA*wght 421c Mmat(n,D1_TA) = Mmat(n,D1_TA) + EUEG*dgdTA*wght 422c else 423c Amat(n,D1_RA) = Amat(n,D1_RA) 424c & + (EUEGPA*gcab + EUEG*dgdPA)*wght 425c Amat(n,D1_RB) = Amat(n,D1_RB) 426c & + (EUEGPB*gcab + EUEG*dgdPB)*wght 427c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + EUEG*dgdGA*wght 428c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + EUEG*dgdGB*wght 429c Mmat(n,D1_TA) = Mmat(n,D1_TA) + EUEG*dgdTA*wght 430c Mmat(n,D1_TB) = Mmat(n,D1_TB) + EUEG*dgdTB*wght 431c endif 432 endIf 433c write (*,*) "Amat(n,1),Cmat(n,1),Mmat(n,1)",Amat(n,1),Cmat(n,1) 434c & ,Mmat(n,1),ipol 435c stop 43620 continue 437 end 438C> 439C> \brief Compute the same-spin part of the VS98 correlation functional 440C> 441C> This routine evaluates the same-spin part of the VS98 functional for 442C> 1 grid point and 1 spin-case. 443C> 444#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 445#if defined(NWAD_PRINT) 446 Subroutine nwxc_vs98ss_p(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi, 447 + EUEGP,ChiP,ChiG,ZP,ZT, 448 + r13,r14,r15,r16,r17,r18) 449#else 450 Subroutine nwxc_vs98ss(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi, 451 + EUEGP,ChiP,ChiG,ZP,ZT, 452 + r13,r14,r15,r16,r17,r18) 453#endif 454#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 455 Subroutine nwxc_vs98ss_d2(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi, 456 + EUEGP,ChiP,ChiG,ZP,ZT, 457 + r13,r14,r15,r16,r17,r18) 458#else 459 Subroutine nwxc_vs98ss_d3(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi, 460 + EUEGP,ChiP,ChiG,ZP,ZT, 461 + r13,r14,r15,r16,r17,r18) 462#endif 463#include "nwad.fh" 464 Implicit none 465c 466#include "intf_nwxc_c_lsda.fh" 467#include "intf_nwxc_gvt4.fh" 468C 469C Compute the same-spin part of the vs98 correlation functional for one grid 470C point and one spin-case. 471C 472 473 double precision tol_rho 474#if defined(NWAD_PRINT) 475#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 476 type(nwad_dble)::r13, r14, r15, r16, r17, r18 477#else 478 double precision r13, r14, r15, r16, r17, r18 479#endif 480#else 481 double precision r13, r14, r15, r16, r17, r18 482#endif 483 type(nwad_dble)::PX, GX, TX, F, EUEG, Chi, RS, Z, P, E, PotLC 484 type(nwad_dble)::rhoo, rho53, rho83 485 type(nwad_dble)::kc,xk,zk,d,gc, Zeta 486 double precision FP, FG, FT, DTol, ZP, ZT 487 double precision EUEGP, ChiP, ChiG, cf, gcc 488 double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11 489 double precision Pi, Pi34, F13, F23, F43, F53, F83, F113 490 double precision RSP, DX, DZ, dgdP, dgdG, dgdT 491 double precision DP, DG, DT 492 double precision F4o3, dgdx, dgdz 493 double precision d2LdSS, d2LdSZ, d2LdZZ, dLdS, dLdZ 494 495 Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/, 496 $ F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/, 497 $ gcc/0.00515088d0/,cf/9.115599720d0/ 498 499 500 F4o3 = 4.0d0/3.0d0 501 dtol=tol_rho 502 If(PX.le.DTol) then 503 EUEG = Zero 504 Chi = Zero 505 EUEGP = Zero 506 ChiP = Zero 507 ChiG = Zero 508 PX = Zero 509 GX = Zero 510 TX = Zero 511 F = Zero 512 FP = Zero 513 FG = Zero 514 FT = Zero 515 Z = Zero 516 ZP = Zero 517 ZT = Zero 518 else 519 Pi = F4*ATan(F1) 520 Pi34 = F3 / (F4*Pi) 521 F13 = F1 / F3 522 F23 = F2 / F3 523 F43 = F2 * F23 524 F53 = F5 / F3 525 F83 = F8 / F3 526 F113 = F11 / F3 527 rhoo = PX 528c rrho = 1.0d0/rhoo 529c rho43 = rhoo**F4o3 530c rho13 = rho43*rrho 531 rho53 = rhoo**F53 532 rho83 = rho53*rhoo 533 534 RS = (Pi34/PX) ** F13 535 Zeta = F1 536#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 537#if defined(NWAD_PRINT) 538 Call nwxc_c_lsda_p(tol_rho, 539 A RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 540#else 541 Call nwxc_c_lsda(tol_rho, 542 A RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 543#endif 544#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 545 Call nwxc_c_lsda_d2(tol_rho, 546 A RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 547#else 548 Call nwxc_c_lsda_d3(tol_rho, 549 A RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 550#endif 551 EUEG = PX*PotLC 552 Chi = GX/rho83 553 Z = (TX/rho53) - cf 554 kc = F1 + gcc*(Chi + Z) 555 xk = Chi/kc 556 zk = Z/kc 557 D = F1 - Chi/(F4*(Z + cf)) 558#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 559#if defined(NWAD_PRINT) 560 call nwxc_gvt4_p(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc, 561 & r13,r14,r15,r16,r17,r18) 562#else 563 call nwxc_gvt4(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc, 564 & r13,r14,r15,r16,r17,r18) 565#endif 566#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 567 call nwxc_gvt4_d2(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc, 568 & r13,r14,r15,r16,r17,r18) 569#else 570 call nwxc_gvt4_d3(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc, 571 & r13,r14,r15,r16,r17,r18) 572#endif 573 E = D*EUEG*gc 574c write (*,*) "Chi, Z, gc", CHi, Z, gc 575 F = E 576c 577c RSP = -RS/(F3*Px) 578c ChiG = F1/PX**F83 579c ChiP = -F83*Chi/PX 580c ZP = -F53 * TX/rho83 581c ZT = F1/rho53 582c DZ = Chi/(F4*(Z + cf)*(Z + cf)) 583c DX = -F1/(F4*(Z + cf)) 584c DP = DZ*ZP + DX*ChiP 585c DG = DX*ChiG 586c DT = DZ*ZT 587c dgdP = dgdx*ChiP + dgdz*ZP 588c dgdG = dgdx*ChiG 589c dgdT = dgdz*ZT 590c EUEGP = PotLC + PX*dLdS*RSP 591c FP = DP*EUEG*gc + D*EUEGP*gc + D*EUEG*dgdP 592c FG = DG*EUEG*gc + D*EUEG*dgdG 593c FT = DT*EUEG*gc + D*EUEG*dgdT 594 Endif 595 Return 596 End 597#ifndef NWAD_PRINT 598#define NWAD_PRINT 599c 600c Compile source again for Maxima 601c 602#include "nwxc_c_vs98.F" 603#endif 604#ifndef SECOND_DERIV 605#define SECOND_DERIV 606c 607c Compile source again for the 2nd derivative case 608c 609#include "nwxc_c_vs98.F" 610#endif 611#ifndef THIRD_DERIV 612#define THIRD_DERIV 613c 614c Compile source again for the 3rd derivative case 615c 616#include "nwxc_c_vs98.F" 617#endif 618#undef NWAD_PRINT 619C> @} 620