1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_x_sogga.F 7C> The SOGGA, SOGGA11 and SOGGA-X exchange functionals 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief Evaluate the SOGGA, SOGGA11 and SOGGA11-X exchange functionals 17C> 18C> The SOGGA, SOGGA11 and SOGGA11-X functionals are GGA functionals that 19C> are formulated and optimized such that the coefficient of the 20C> gradient correction term matches that of the appropriate expansion ofC> the electronic energy [1,2,3]. 21C> 22C> ### References ### 23C> 24C> [1] Y. Zhao, D.G. Truhlar, 25C> "Construction of a generalized gradient approximation by 26C> restoring the density-gradient expansion and enforcing a tight 27C> Lieb-Oxford bound", J. Chem. Phys. <b>128</b> (2008) 184109, 28C> DOI: 29C> <a href="https://doi.org/10.1063/1.2912068"> 30C> 10.1063/1.2912068</a>. 31C> 32C> [2] R. Peverati, Y. Zhao, D.G. Truhlar, 33C> "Generalized gradient approximation that recovers the 34C> second-order density-gradient expansion with optimized 35C> across-the-board performance", J. Phys. Chem. Lett. <b>2</b> 36C> (2011) 1991-1997, DOI: 37C> <a href="https://doi.org/10.1021/jz200616w"> 38C> 10.1021/jz200616w</a>. 39C> 40C> [3] R. Peverati, D.G. Truhlar, 41C> "Communication: A global hybrid generalized gradient 42C> approximation to the exchange-correlation functional that 43C> satisfies the second-order density-gradient constraint and has 44C> broad applicability in chemistry", J. Chem. Phys. <b>135</b> 45C> (2011) 191102, DOI: 46C> <a href="https://doi.org/10.1063/1.3663871"> 47C> 10.1063/1.3663871</a>. 48C> 49#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 50#if defined(NWAD_PRINT) 51 Subroutine nwxc_x_sogga_p(param, tol_rho, ipol, nq, wght, rho, 52 & rgamma, ffunc) 53#else 54 Subroutine nwxc_x_sogga(param, tol_rho, ipol, nq, wght, rho, 55 & rgamma, ffunc) 56#endif 57#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 58 Subroutine nwxc_x_sogga_d2(param, tol_rho, ipol, nq, wght, rho, 59 & rgamma, ffunc) 60#else 61 Subroutine nwxc_x_sogga_d3(param, tol_rho, ipol, nq, wght, rho, 62 & rgamma, ffunc) 63#endif 64c 65c$Id$ 66c 67c**********************************************************************c 68c c 69c SOGGA11X evaluates the exchange part of the SOGGA, SOGGA11 c 70c and SOGGA11-X functionals on the grid. c 71c c 72c a) Zhao and Truhlar, J.Chem.Phys., 128, 184109 (2008) c 73c b) Peverati, Zhao and Truhlar, J.Phys.Chem.Lett, 2, 1991 (2011) c 74c c) Peverati and Truhlar, J.Chem.Phys, 135, 191102 (2011) c 75c c 76c ijzy = 1 - SOGGA functional (a) - it requires PBE correlation c 77c ijzy = 2 - SOGGA11 functional (b) c 78c ijzy = 3 - SOGGA11-X functional (c) c 79c c 80c Coded by Roberto Peverati (12/11) c 81c c 82c**********************************************************************c 83c 84#include "nwad.fh" 85c 86 implicit none 87#include "nwxc_param.fh" 88c 89c Input and other parameters 90c 91#if defined(NWAD_PRINT) 92#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 93 type(nwad_dble)::param(*)!< [Input] Parameters of the functional, 94 type(nwad_dble)::CxA0,CxA1,CxA2,CxA3,CxA4,CxA5 95 type(nwad_dble)::CxB0,CxB1,CxB2,CxB3,CxB4,CxB5 96 type(nwad_dble)::FA0 97 type(nwad_dble)::FB0 98#else 99 double precision param(*)!< [Input] Parameters of the functional, 100 double precision CxA0,CxA1,CxA2,CxA3,CxA4,CxA5 101 double precision CxB0,CxB1,CxB2,CxB3,CxB4,CxB5 102 double precision FA0 103 double precision FB0 104#endif 105#else 106 double precision param(*)!< [Input] Parameters of the functional, 107 !< see Table 1 of [1] and Table 1 of [2]. 108 !< - param(1): \f$ a_0 \f$ 109 !< - param(2): \f$ a_1 \f$ 110 !< - param(3): \f$ a_2 \f$ 111 !< - param(4): \f$ a_3 \f$ 112 !< - param(5): \f$ a_4 \f$ 113 !< - param(6): \f$ a_5 \f$ 114 !< - param(7): \f$ b_0 \f$ 115 !< - param(8): \f$ b_1 \f$ 116 !< - param(9): \f$ b_2 \f$ 117 !< - param(10): \f$ b_3 \f$ 118 !< - param(11): \f$ b_4 \f$ 119 !< - param(12): \f$ b_5 \f$ 120 double precision CxA0,CxA1,CxA2,CxA3,CxA4,CxA5 121 double precision CxB0,CxB1,CxB2,CxB3,CxB4,CxB5 122 double precision FA0 123 double precision FB0 124#endif 125 double precision tol_rho !< [Input] The lower limit on the density 126 integer ipol !< [Input] The number of spin channels 127 integer nq !< [Input] The number of points 128 double precision wght !< [Input] The weight of the functional 129c 130c Charge Density 131c 132 type(nwad_dble)::rho(nq,*) !< [Input] The density 133c 134c Charge Density Gradient 135c 136 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 137c 138c Sampling Matrices for the XC Potential 139c 140 type(nwad_dble)::ffunc(nq) !< [Output] The value of the functional 141c double precision Amat(nq,*) !< [Output] The derivative wrt rho 142c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 143c 144 double precision pi 145c 146c Intermediate derivative results, etc. 147c 148 integer n 149c 150 type(nwad_dble)::rho43, rho13, rhoo 151c 152 double precision AS, ASO, AX, DELOCDR 153 double precision DFA1DG, DFA1DR, DFA1DY 154 double precision DFA2DG, DFA2DR, DFA2DY 155 double precision DFA3DG, DFA3DR, DFA3DY 156 double precision DFA4DG, DFA4DR, DFA4DY 157 double precision DFA5DG, DFA5DR, DFA5DY 158 double precision DFB1DG, DFB1DR, DFB1DY 159 double precision DFB2DG, DFB2DR, DFB2DY 160 double precision DFB3DG, DFB3DR, DFB3DY 161 double precision DFB4DG, DFB4DR, DFB4DY 162 double precision DFB5DG, DFB5DR, DFB5DY 163 double precision DFEXPDPON, DFFRACDPON, DFGGAXDG, DFGGAXDR 164 double precision DYDG, DYDR, DTOL 165 double precision MU 166 type(nwad_dble)::FA1, FA2, FA3, FA4, FA5 167 type(nwad_dble)::FB1, FB2, FB3, FB4, FB5 168c type(nwad_dble)::Gam12, Gam, X, FEXP, FFRAC, FGGAX, S, Y, ELOC 169 type(nwad_dble)::Gam, X, FEXP, FFRAC, FGGAX, Y, ELOC 170 type(nwad_dble)::PON 171c 172 double precision f1,f2,f3,f4,f5,f8 173 double precision F1o3,F4o3,F48 174 parameter( F1=1.0D+00, F2=2.0D+00, F3=3.0D+00, 175 $ F4=4.0D+00, F5=5.0D+00, F8=8.0D+00, 176 $ F48=48.0D+00) 177c 178 pi=acos(-1d0) 179c 180c if (ijzy.eq.1) then 181c SOGGA 182c CxA0 = 0.5d0 183c CxA1 = 0.276d0 184c CxA2 = 0d0 185c CxA3 = 0d0 186c CxA4 = 0d0 187c CxA5 = 0d0 188c CxB0 = 0.5d0 189c CxB1 = 0.276d0 190c CxB2 = 0d0 191c CxB3 = 0d0 192c CxB4 = 0d0 193c CxB5 = 0d0 194c elseif (ijzy.eq.2) then 195c SOGGA11 196c CxA0 = 0.50000d0 197c CxA1 = -2.95535d0 198c CxA2 = 15.7974d0 199c CxA3 = -91.1804d0 200c CxA4 = 96.2030d0 201c CxA5 = 0.18683d0 202c CxB0 = 0.50000d0 203c CxB1 = 3.50743d0 204c CxB2 = -12.9523d0 205c CxB3 = 49.7870d0 206c CxB4 = -33.2545d0 207c CxB5 = -11.1396d0 208c elseif (ijzy.eq.3) then 209c SOGGA11-X 210c CxA0 = 2.99250d-01 211c CxA1 = 3.21638d+00 212c CxA2 = -3.55605d+00 213c CxA3 = 7.65852d+00 214c CxA4 = -1.12830d+01 215c CxA5 = 5.25813d+00 216c CxB0 = 2.99250d-01 217c CxB1 = -2.88595d+00 218c CxB2 = 3.23617d+00 219c CxB3 = -2.45393d+00 220c CxB4 = -3.75495d+00 221c CxB5 = 3.96613d+00 222c endif 223 CxA0 = param(1) 224 CxA1 = param(2) 225 CxA2 = param(3) 226 CxA3 = param(4) 227 CxA4 = param(5) 228 CxA5 = param(6) 229 CxB0 = param(7) 230 CxB1 = param(8) 231 CxB2 = param(9) 232 CxB3 = param(10) 233 CxB4 = param(11) 234 CxB5 = param(12) 235c 236 DTol = tol_rho 237 F1o3 = F1/F3 238 F4o3 = F4/F3 239 Pi = ACos(-F1) 240 AsO = (F48*PI*PI)**F1o3 241 As = F1/AsO 242 Ax = -(F3/F2) * (F4o3*Pi)**(-F1o3) 243 mu = 0.2236536053d0 244c 245 if (ipol.eq.1 )then 246c 247c ======> SPIN-RESTRICTED <====== 248c or 249c SPIN-UNPOLARIZED 250c 251c 252 do 10 n = 1, nq 253 if (rho(n,R_T).lt.DTol) goto 10 254 rhoo = rho(n,R_T)/F2 255 rho43 = rhoo**F4o3 256 rho13 = rho43/rhoo 257 Gam = rgamma(n,G_TT)/F4 258c Gam =(delrho(n,1,1)*delrho(n,1,1) + 259c & delrho(n,2,1)*delrho(n,2,1) + 260c & delrho(n,3,1)*delrho(n,3,1))/F4 261c Gam12 = sqrt(Gam) 262c if(gam12.lt.dtol) goto 10 263c 264 Eloc = Ax*Rho43 265c x = Gam12/Rho43 266c s = As*x 267c y = s*s 268 y = As*As*Gam/(Rho43*Rho43) 269 PON = mu*y 270 Ffrac = F1-F1/(F1+PON) 271 Fexp = F1-exp(-PON) 272 fa0 = CxA0 273 fa1 = CxA1 *Ffrac 274 fa2 = CxA2 *Ffrac**F2 275 fa3 = CxA3 *Ffrac**F3 276 fa4 = CxA4 *Ffrac**F4 277 fa5 = CxA5 *Ffrac**F5 278 fb0 = CxB0 279 fb1 = CxB1 *Fexp 280 fb2 = CxB2 *Fexp**F2 281 fb3 = CxB3 *Fexp**F3 282 fb4 = CxB4 *Fexp**F4 283 fb5 = CxB5 *Fexp**F5 284c 285 Fggax = fa0+fa1+fa2+fa3+fa4+fa5 + 286 $ fb0+fb1+fb2+fb3+fb4+fb5 287C 288C 1st derivatives. 289C 290 291c dElocdR=Ax*F4o3*Rho13 292c dydR = -(F8/F3)*y/Rhoo 293c dydG = y/Gam 294c dFfracdPON = F1/((F1+PON)**F2) 295c dFexpdPON = exp(-PON) 296c dfa1dy = CxA1 *mu*dFfracdPON 297c dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON 298c dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON 299c dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON 300c dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON 301c dfa1dR = dfa1dy *dydR 302c dfa2dR = dfa2dy *dydR 303c dfa3dR = dfa3dy *dydR 304c dfa4dR = dfa4dy *dydR 305c dfa5dR = dfa5dy *dydR 306c dfa1dG = dfa1dy *dydG 307c dfa2dG = dfa2dy *dydG 308c dfa3dG = dfa3dy *dydG 309c dfa4dG = dfa4dy *dydG 310c dfa5dG = dfa5dy *dydG 311c dfb1dy = CxB1 *mu*dFexpdPON 312c dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON 313c dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON 314c dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON 315c dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON 316c dfb1dR = dfb1dy *dydR 317c dfb2dR = dfb2dy *dydR 318c dfb3dR = dfb3dy *dydR 319c dfb4dR = dfb4dy *dydR 320c dfb5dR = dfb5dy *dydR 321c dfb1dG = dfb1dy *dydG 322c dfb2dG = dfb2dy *dydG 323c dfb3dG = dfb3dy *dydG 324c dfb4dG = dfb4dy *dydG 325c dfb5dG = dfb5dy *dydG 326c 327c dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR + 328c $ dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR 329c 330c dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG + 331c $ dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG 332c 333c Ex = Ex +F2*(Eloc*Fggax)*qwght(n) 334 ffunc(n)=ffunc(n)+F2*(Eloc*Fggax)*wght 335c Amat(n,D1_RA)=Amat(n,D1_RA) 336c $ +(dElocdR*Fggax+Eloc*dFggaxdR)*wght 337c Cmat(n,D1_GAA)=Cmat(n,D1_GAA)+Eloc*dFggaxdG*wght 338c 33910 continue 340c 341c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 342 else ! ipol=2 343c 344c ======> SPIN-UNRESTRICTED <====== 345 346c 347c use spin density functional theory ie n-->2n 348c Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta] 349c 350c Alpha ALPHA ALPHA 351c 352 do 20 n = 1, nq 353 if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20 354 if (rho(n,R_A).lt.DTol) goto 25 355 rhoo = rho(n,R_A) 356 rho43 = rhoo**F4o3 357 rho13 = rho43/rhoo 358 Gam = rgamma(n,G_AA) 359c Gam =(delrho(n,1,1)*delrho(n,1,1) + 360c & delrho(n,2,1)*delrho(n,2,1) + 361c & delrho(n,3,1)*delrho(n,3,1)) 362c Gam12 = sqrt(Gam) 363c if(gam12.lt.dtol) goto 25 364c 365 Eloc = Ax*Rho43 366c x = Gam12/Rho43 367c s = As*x 368c y = s*s 369 y = As*As*Gam/(Rho43*Rho43) 370 PON = mu*y 371 Ffrac = F1-F1/(F1+PON) 372 Fexp = F1-exp(-PON) 373 fa0 = CxA0 374 fa1 = CxA1 *Ffrac 375 fa2 = CxA2 *Ffrac**F2 376 fa3 = CxA3 *Ffrac**F3 377 fa4 = CxA4 *Ffrac**F4 378 fa5 = CxA5 *Ffrac**F5 379 fb0 = CxB0 380 fb1 = CxB1 *Fexp 381 fb2 = CxB2 *Fexp**F2 382 fb3 = CxB3 *Fexp**F3 383 fb4 = CxB4 *Fexp**F4 384 fb5 = CxB5 *Fexp**F5 385c 386 Fggax = fa0+fa1+fa2+fa3+fa4+fa5 + 387 $ fb0+fb1+fb2+fb3+fb4+fb5 388C 389C 1st derivatives. 390C 391c dElocdR=Ax*F4o3*Rho13 392c dydR = -(F8/F3)*y/Rhoo 393c dydG = y/Gam 394c dFfracdPON = F1/((F1+PON)**F2) 395c dFexpdPON = exp(-PON) 396c dfa1dy = CxA1 *mu*dFfracdPON 397c dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON 398c dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON 399c dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON 400c dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON 401c dfa1dR = dfa1dy *dydR 402c dfa2dR = dfa2dy *dydR 403c dfa3dR = dfa3dy *dydR 404c dfa4dR = dfa4dy *dydR 405c dfa5dR = dfa5dy *dydR 406c dfa1dG = dfa1dy *dydG 407c dfa2dG = dfa2dy *dydG 408c dfa3dG = dfa3dy *dydG 409c dfa4dG = dfa4dy *dydG 410c dfa5dG = dfa5dy *dydG 411c dfb1dy = CxB1 *mu*dFexpdPON 412c dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON 413c dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON 414c dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON 415c dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON 416c dfb1dR = dfb1dy *dydR 417c dfb2dR = dfb2dy *dydR 418c dfb3dR = dfb3dy *dydR 419c dfb4dR = dfb4dy *dydR 420c dfb5dR = dfb5dy *dydR 421c dfb1dG = dfb1dy *dydG 422c dfb2dG = dfb2dy *dydG 423c dfb3dG = dfb3dy *dydG 424c dfb4dG = dfb4dy *dydG 425c dfb5dG = dfb5dy *dydG 426c 427c dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR + 428c $ dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR 429c 430c dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG + 431c $ dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG 432c 433c Ex = Ex + (Eloc*Fggax)*qwght(n) 434 ffunc(n)=ffunc(n)+ Eloc*Fggax*wght 435c Amat(n,D1_RA) = Amat(n,D1_RA) + dElocdR*Fggax*wght 436c $ + Eloc*dFggaxdR*wght 437c Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + Eloc*dFggaxdG*wght 438c 43925 continue 440c 441c Beta BETA BETA 442c 443 if (rho(n,R_B).lt.DTol) goto 20 444 rhoo = rho(n,R_B) 445 rho43 = rhoo**F4o3 446 rho13 = rho43/rhoo 447c 448 Gam = rgamma(n,G_BB) 449c Gam =(delrho(n,1,2)*delrho(n,1,2) + 450c & delrho(n,2,2)*delrho(n,2,2) + 451c & delrho(n,3,2)*delrho(n,3,2)) 452c Gam12 = sqrt(Gam) 453c if(gam12.lt.dtol) goto 20 454c 455 Eloc = Ax*Rho43 456c x = Gam12/Rho43 457c s = As*x 458c y = s*s 459 y = As*As*Gam/(Rho43*Rho43) 460 PON = mu*y 461 Ffrac = F1-F1/(F1+PON) 462 Fexp = F1-exp(-PON) 463 fa0 = CxA0 464 fa1 = CxA1 *Ffrac 465 fa2 = CxA2 *Ffrac**F2 466 fa3 = CxA3 *Ffrac**F3 467 fa4 = CxA4 *Ffrac**F4 468 fa5 = CxA5 *Ffrac**F5 469 fb0 = CxB0 470 fb1 = CxB1 *Fexp 471 fb2 = CxB2 *Fexp**F2 472 fb3 = CxB3 *Fexp**F3 473 fb4 = CxB4 *Fexp**F4 474 fb5 = CxB5 *Fexp**F5 475c 476 Fggax = fa0+fa1+fa2+fa3+fa4+fa5 + 477 $ fb0+fb1+fb2+fb3+fb4+fb5 478C 479C 1st derivatives. 480C 481 482c dElocdR=Ax*F4o3*Rho13 483c dydR = -(F8/F3)*y/Rhoo 484c dydG = y/Gam 485c dFfracdPON = F1/((F1+PON)**F2) 486c dFexpdPON = exp(-PON) 487c dfa1dy = CxA1 *mu*dFfracdPON 488c dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON 489c dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON 490c dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON 491c dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON 492c dfa1dR = dfa1dy *dydR 493c dfa2dR = dfa2dy *dydR 494c dfa3dR = dfa3dy *dydR 495c dfa4dR = dfa4dy *dydR 496c dfa5dR = dfa5dy *dydR 497c dfa1dG = dfa1dy *dydG 498c dfa2dG = dfa2dy *dydG 499c dfa3dG = dfa3dy *dydG 500c dfa4dG = dfa4dy *dydG 501c dfa5dG = dfa5dy *dydG 502c dfb1dy = CxB1 *mu*dFexpdPON 503c dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON 504c dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON 505c dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON 506c dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON 507c dfb1dR = dfb1dy *dydR 508c dfb2dR = dfb2dy *dydR 509c dfb3dR = dfb3dy *dydR 510c dfb4dR = dfb4dy *dydR 511c dfb5dR = dfb5dy *dydR 512c dfb1dG = dfb1dy *dydG 513c dfb2dG = dfb2dy *dydG 514c dfb3dG = dfb3dy *dydG 515c dfb4dG = dfb4dy *dydG 516c dfb5dG = dfb5dy *dydG 517c 518c dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR + 519c $ dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR 520c 521c dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG + 522c $ dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG 523c 524c Ex = Ex + (Eloc*Fggax)*qwght(n) 525 ffunc(n)=ffunc(n)+ Eloc*Fggax*wght 526 527c Amat(n,D1_RB) = Amat(n,D1_RB) + dElocdR*Fggax*wght 528c $ + Eloc*dFggaxdR*wght 529c 530c Cmat(n,D1_GBB)= Cmat(n,D1_GBB) + Eloc*dFggaxdG*wght 531c 53220 continue 533 endif 534 return 535 end 536#ifndef NWAD_PRINT 537#define NWAD_PRINT 538c 539c compile again for 2nd derivatives 540c 541#include "nwxc_x_sogga.F" 542#endif 543#ifndef SECOND_DERIV 544#define SECOND_DERIV 545c 546c compile again for 2nd derivatives 547c 548#include "nwxc_x_sogga.F" 549#endif 550#ifndef THIRD_DERIV 551#define THIRD_DERIV 552c 553c compile again for 3rd derivatives 554c 555#include "nwxc_x_sogga.F" 556#endif 557#undef NWAD_PRINT 558C> 559C> @} 560