1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_sogga.F 7C> The SOGGA11 and SOGGA-X correlation functionals 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief Evaluate the SOGGA11 and SOGGA11-X correlation functionals 17C> 18C> The SOGGA11 and SOGGA11-X functionals are GGA functionals that are 19C> formulated and optimized such that the coefficient of the gradient 20C> correction term matches that of the appropriate expansion of the 21C> electronic energy [1,2]. 22C> 23C> ### References ### 24C> 25C> [1] R. Peverati, Y. Zhao, D.G. Truhlar, 26C> "Generalized gradient approximation that recovers the 27C> second-order density-gradient expansion with optimized 28C> across-the-board performance", J. Phys. Chem. Lett. <b>2</b> 29C> (2011) 1991-1997, DOI: 30C> <a href="https://doi.org/10.1021/jz200616w"> 31C> 10.1021/jz200616w</a>. 32C> 33C> [2] R. Peverati, D.G. Truhlar, 34C> "Communication: A global hybrid generalized gradient 35C> approximation to the exchange-correlation functional that 36C> satisfies the second-order density-gradient constraint and has 37C> broad applicability in chemistry", J. Chem. Phys. <b>135</b> 38C> (2011) 191102, DOI: 39C> <a href="https://doi.org/10.1063/1.3663871"> 40C> 10.1063/1.3663871</a>. 41C> 42#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 43#if defined(NWAD_PRINT) 44 subroutine nwxc_c_sogga_p(param, tol_rho, ipol, nq, wght, rho, 45 & rgamma, ffunc) 46#else 47 subroutine nwxc_c_sogga(param, tol_rho, ipol, nq, wght, rho, 48 & rgamma, ffunc) 49#endif 50#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 51 subroutine nwxc_c_sogga_d2(param, tol_rho, ipol, nq, wght, rho, 52 & rgamma, ffunc) 53#else 54 subroutine nwxc_c_sogga_d3(param, tol_rho, ipol, nq, wght, rho, 55 & rgamma, ffunc) 56#endif 57c 58c$Id$ 59c 60c 61c**********************************************************************c 62c c 63c xc_csogga evaluates the corelation part of the SOGGA11 and c 64c SOGGA11-X functionals on the grid. c 65c c 66c a) Peverati, Zhao and Truhlar, J.Phys.Chem.Lett, 2, 1991 (2011) c 67c b) Peverati and Truhlar, J.Chem.Phys, 135, 191102 (2011) c 68c c 69c ijzy = 1 - SOGGA11 functional (a) c 70c ijzy = 2 - SOGGA11-X functional (b) c 71c c 72c Coded by Roberto Peverati (12/11) c 73c c 74c**********************************************************************c 75c 76#include "nwad.fh" 77 implicit none 78c 79#include "intf_nwxc_c_lsda.fh" 80#include "intf_nwxc_GZeta.fh" 81c 82#include "nwxc_param.fh" 83c 84c Input and other parameters 85c 86#if defined(NWAD_PRINT) 87#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 88 type(nwad_dble)::param(*) !< [Input] Parameters of the functional, 89 type(nwad_dble)::CcA(0:5),CcB(0:5) 90 type(nwad_dble)::FA0, FB0 91#else 92 double precision param(*) !< [Input] Parameters of the functional, 93 double precision CcA(0:5),CcB(0:5) 94 double precision FA0, FB0 95#endif 96#else 97 double precision param(*)!< [Input] Parameters of the functional, 98 !< see Table 1 of [1] and Table 1 of [2]. 99 !< - param(1): \f$ a_0 \f$ 100 !< - param(2): \f$ a_1 \f$ 101 !< - param(3): \f$ a_2 \f$ 102 !< - param(4): \f$ a_3 \f$ 103 !< - param(5): \f$ a_4 \f$ 104 !< - param(6): \f$ a_5 \f$ 105 !< - param(7): \f$ b_0 \f$ 106 !< - param(8): \f$ b_1 \f$ 107 !< - param(9): \f$ b_2 \f$ 108 !< - param(10): \f$ b_3 \f$ 109 !< - param(11): \f$ b_4 \f$ 110 !< - param(12): \f$ b_5 \f$ 111 double precision CcA(0:5),CcB(0:5) 112 double precision FA0, FB0 113#endif 114 double precision tol_rho !< [Input] The lower limit on the density 115 integer ipol !< [Input] The number of spin channels 116 integer nq !< [Input] The number of points 117 double precision wght !< [Input] The weight of the functional 118c 119c Charge Density 120c 121 type(nwad_dble)::rho(nq,*) !< [Input] The density 122c 123c Charge Density Gradient 124c 125 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 126c 127c Sampling Matrices for the XC Potential 128c 129 type(nwad_dble)::ffunc(nq) !< [Output] The value of the functional 130c double precision Amat(nq,*) !< [Output] The derivative wrt rho 131c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 132c 133c Intermediate derivative results, etc. 134c 135 integer n 136c 137 type(nwad_dble)::RHOA,RHOB,rhoval,gammaval,RS,Zeta 138 double precision BETA, CC0 139 double precision DELOCDA, DELOCDB, DFA1DG, DFA1DR, DFA1DZ 140 double precision DFA2DG, DFA2DR, DFA2DZ, DFA3DG, DFA3DR, DFA3DZ 141 double precision DFA4DG, DFA4DR, DFA4DZ, DFA5DG, DFA5DR, DFA5DZ 142 double precision DFB1DG, DFB1DR, DFB1DZ, DFB2DG, DFB2DR, DFB2DZ 143 double precision DFB3DG, DFB3DR, DFB3DZ, DFB4DG, DFB4DR, DFB4DZ 144 double precision DFB5DG, DFB5DR, DFB5DZ 145 double precision DFEXPDPON, DFFRACDPON, DFGGACDA, DFGGACDB 146 double precision DFGGACDG, DFGGACDGX, DFGGACDR, DFGGACDZ 147 double precision DGDGX, DGDZ, DLDA, DLDB, DLDS, DLDZ 148 double precision DPONDG, DPONDR, DPONDZ, DSDR 149 double precision DT2DG, DT2DR, DT2DZ, DZDA, DZDB 150 double precision F1O3, F1O6, F7O6 151 type(nwad_dble)::FA1, FA2, FA3, FA4, FA5 152 type(nwad_dble)::FB1, FB2, FB3, FB4, FB5 153 type(nwad_dble)::FEXP, FFRAC 154 type(nwad_dble)::FGGAC 155 double precision FKFAC 156 type(nwad_dble)::G, G2, G3 157 double precision PI, PI34 158c type(nwad_dble)::PON, T, T2 159 type(nwad_dble)::PON, T2 160 type(nwad_dble)::RHO16, RHO76, GRHO, POTLC, ELOC 161 double precision SCFAC, SKFAC, XNU 162 double precision dummy 163 164 165 double precision F1, F2, F3, F4, F5, F6, F7 166 Save F1, F2, F3, F4, F5, F6, F7 167 DATA F1/1.0D+00/, F2/2.0D+00/, F3/3.0D+00/, 168 $ F4/4.0D+00/, F5/5.0D+00/, F6/6.0D+00/, 169 $ F7/7.0D+00/ 170c 171c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 172c 173c if (ijzy.eq.1) then 174c CcA(0) = 0.50000D+00 175c CcA(1) = -4.62334D+00 176c CcA(2) = 8.00410D+00 177c CcA(3) = -130.226D+00 178c CcA(4) = 38.2685D+00 179c CcA(5) = 69.5599D+00 180c CcB(0) = 0.50000D+00 181c CcB(1) = 3.62334D+00 182c CcB(2) = 9.36393D+00 183c CcB(3) = 34.5114D+00 184c CcB(4) = -18.5684D+00 185c CcB(5) = -0.16519D+00 186c elseif (ijzy.eq.2) then 187c CcA(0) = 5.00000d-01 188c CcA(1) = 7.82439d+01 189c CcA(2) = 2.57211d+01 190c CcA(3) = -1.38830d+01 191c CcA(4) = -9.87375d+00 192c CcA(5) = -1.41357d+01 193c CcB(0) = 5.00000d-01 194c CcB(1) = -7.92439d+01 195c CcB(2) = 1.63725d+01 196c CcB(3) = 2.08129d+00 197c CcB(4) = 7.50769d+00 198c CcB(5) = -1.01861d+01 199c endif 200 CcA(0) = param(1) 201 CcA(1) = param(2) 202 CcA(2) = param(3) 203 CcA(3) = param(4) 204 CcA(4) = param(5) 205 CcA(5) = param(6) 206 CcB(0) = param(7) 207 CcB(1) = param(8) 208 CcB(2) = param(9) 209 CcB(3) = param(10) 210 CcB(4) = param(11) 211 CcB(5) = param(12) 212 213 Pi = ACos(-F1) 214 Pi34 = F3/(F4*Pi) 215 F1o3 = F1/F3 216 F1o6 = F1/F6 217 F7o6 = F7/F6 218 XNu = 15.75592D0 219 CC0 = 0.004235D0 220 beta= XNu*CC0 221 SCfac=F1 222 223 224c 225c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 226c 227 do 20 n = 1, nq 228 229 if (ipol.eq.1) then 230 RHOA = rho(n,R_T)/F2 231 RHOB = RHOA 232 gammaval = rgamma(n,G_TT) 233c gammaval =(delrho(n,1,1)*delrho(n,1,1) + 234c & delrho(n,2,1)*delrho(n,2,1) + 235c & delrho(n,3,1)*delrho(n,3,1)) 236 else 237 RhoA = 0.0d0 238 RhoB = 0.0d0 239 gammaval = 0.0d0 240 if (rho(n,R_A).gt.0.5d0*tol_rho) then 241 RhoA = rho(n,R_A) 242 gammaval = gammaval + rgamma(n,G_AA) 243 endif 244 if (rho(n,R_B).gt.0.5d0*tol_rho) then 245 RhoB = rho(n,R_B) 246 gammaval = gammaval + rgamma(n,G_BB) 247 if (rho(n,R_A).gt.0.5d0*tol_rho) then 248 gammaval = gammaval + 2.0d0*rgamma(n,G_AB) 249 endif 250 endif 251c gammaval = delrho(n,1,1)*delrho(n,1,1) + 252c & delrho(n,1,2)*delrho(n,1,2) + 253c & delrho(n,2,1)*delrho(n,2,1) + 254c & delrho(n,2,2)*delrho(n,2,2) + 255c & delrho(n,3,1)*delrho(n,3,1) + 256c & delrho(n,3,2)*delrho(n,3,2) + 257c & 2.d0*(delrho(n,1,1)*delrho(n,1,2) + 258c & delrho(n,2,1)*delrho(n,2,2) + 259c & delrho(n,3,1)*delrho(n,3,2)) 260 endif 261 Rhoval = RhoA + RhoB 262c GRho = max(dsqrt(gammaval),tol_rho) 263c GRho = sqrt(gammaval) 264 GRho = gammaval 265 if (rhoval.le.tol_rho) goto 20 266 rS = (Pi34/Rhoval)**F1o3 267 Zeta = (RhoA-RhoB)/Rhoval 268#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 269#if defined(NWAD_PRINT) 270 Call nwxc_c_lsda_p(tol_rho, 271 $ RS,Zeta,PotLC,dLdS,dLdZ,dummy,dummy,dummy) 272 Call nwxc_GZeta_p(Zeta,G,dGdZ,dummy,dummy) 273#else 274 Call nwxc_c_lsda(tol_rho, 275 $ RS,Zeta,PotLC,dLdS,dLdZ,dummy,dummy,dummy) 276 Call nwxc_GZeta(Zeta,G,dGdZ,dummy,dummy) 277#endif 278#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 279 Call nwxc_c_lsda_d2(tol_rho, 280 $ RS,Zeta,PotLC,dLdS,dLdZ,dummy,dummy,dummy) 281 Call nwxc_GZeta_d2(Zeta,G,dGdZ,dummy,dummy) 282#else 283 Call nwxc_c_lsda_d3(tol_rho, 284 $ RS,Zeta,PotLC,dLdS,dLdZ,dummy,dummy,dummy) 285 Call nwxc_GZeta_d3(Zeta,G,dGdZ,dummy,dummy) 286#endif 287c 288 Eloc = SCfac*Rhoval*PotLC 289 FKFac = (F3*Pi*Pi)**F1o3 290 SKFac = Sqrt(FKFac/Pi)*F2 291 Rho16 = Rhoval**F1o6 292 Rho76 = Rhoval*Rho16 293c T = GRho/(F2*SKFac*Rho76*G) 294c T2 = T*T 295 T2 = GRho/((F2*SKFac*Rho76*G)**2.0d0) 296 G2 = G*G 297 G3 = G*G2 298 PON = (G3*beta*T2)/PotLC 299 Ffrac = F1-F1/(F1-PON) 300 Fexp = F1-exp(PON) 301 fa0 = CcA(0) 302 fa1 = CcA(1) *Ffrac 303 fa2 = CcA(2) *Ffrac**F2 304 fa3 = CcA(3) *Ffrac**F3 305 fa4 = CcA(4) *Ffrac**F4 306 fa5 = CcA(5) *Ffrac**F5 307 fb0 = CcB(0) 308 fb1 = CcB(1) *Fexp 309 fb2 = CcB(2) *Fexp**F2 310 fb3 = CcB(3) *Fexp**F3 311 fb4 = CcB(4) *Fexp**F4 312 fb5 = CcB(5) *Fexp**F5 313 FggaC = fa0+fa1+fa2+fa3+fa4+fa5+ 314 $ fb0+fb1+fb2+fb3+fb4+fb5 315C 316C 1st derivatives. 317C 318c dSdR = -(F1o3*rS/Rhoval) 319c dZdA = (F1-Zeta)/Rhoval 320c dZdB = (-F1-Zeta)/Rhoval 321c dLdA = dLdS*dSdR + dLdz*dZdA 322c dLdB = dLdS*dSdR + dLdz*dZdB 323c dElocdA = SCfac*PotLC + SCfac*Rhoval*dLdA 324c dElocdB = SCfac*PotLC + SCfac*Rhoval*dLdB 325c dT2dR = -F2*F7o6*T2/Rhoval 326c dT2dG = F2*T2/GRho 327c dT2dZ = -F2*dGdZ*T2/G 328c dPONdR = -(dLdS*dSdR*PON/PotLC)+dT2dR*PON/T2 329c dPONdG = dT2dG*PON/T2 330c dPONdZ = F3*dGdZ*PON/G-dLdZ*PON/PotLC+dT2dZ*PON/T2 331c dFfracdPON = -F1/((F1-PON)**F2) 332c dFexpdPON = -exp(PON) 333c dfa1dR = CcA(1) *dFfracdPON*dPONdR 334c dfa2dR = CcA(2) *(F2 *Ffrac) *dFfracdPON*dPONdR 335c dfa3dR = CcA(3) *(F3 *Ffrac**F2)*dFfracdPON*dPONdR 336c dfa4dR = CcA(4) *(F4 *Ffrac**F3)*dFfracdPON*dPONdR 337c dfa5dR = CcA(5) *(F5 *Ffrac**F4)*dFfracdPON*dPONdR 338c dfa1dG = CcA(1) *dFfracdPON*dPONdG 339c dfa2dG = CcA(2) *(F2 *Ffrac) *dFfracdPON*dPONdG 340c dfa3dG = CcA(3) *(F3 *Ffrac**F2)*dFfracdPON*dPONdG 341c dfa4dG = CcA(4) *(F4 *Ffrac**F3)*dFfracdPON*dPONdG 342c dfa5dG = CcA(5) *(F5 *Ffrac**F4)*dFfracdPON*dPONdG 343c dfa1dZ = CcA(1) *dFfracdPON*dPONdZ 344c dfa2dZ = CcA(2) *(F2 *Ffrac) *dFfracdPON*dPONdZ 345c dfa3dZ = CcA(3) *(F3 *Ffrac**F2)*dFfracdPON*dPONdZ 346c dfa4dZ = CcA(4) *(F4 *Ffrac**F3)*dFfracdPON*dPONdZ 347c dfa5dZ = CcA(5) *(F5 *Ffrac**F4)*dFfracdPON*dPONdZ 348c dfb1dR = CcB(1) *dFexpdPON*dPONdR 349c dfb2dR = CcB(2) *(F2 *Fexp) *dFexpdPON*dPONdR 350c dfb3dR = CcB(3) *(F3 *Fexp**F2) *dFexpdPON*dPONdR 351c dfb4dR = CcB(4) *(F4 *Fexp**F3) *dFexpdPON*dPONdR 352c dfb5dR = CcB(5) *(F5 *Fexp**F4) *dFexpdPON*dPONdR 353c dfb1dG = CcB(1) *dFexpdPON*dPONdG 354c dfb2dG = CcB(2) *(F2 *Fexp) *dFexpdPON*dPONdG 355c dfb3dG = CcB(3) *(F3 *Fexp**F2) *dFexpdPON*dPONdG 356c dfb4dG = CcB(4) *(F4 *Fexp**F3) *dFexpdPON*dPONdG 357c dfb5dG = CcB(5) *(F5 *Fexp**F4) *dFexpdPON*dPONdG 358c dfb1dZ = CcB(1) *dFexpdPON*dPONdZ 359c dfb2dZ = CcB(2) *(F2 *Fexp) *dFexpdPON*dPONdZ 360c dfb3dZ = CcB(3) *(F3 *Fexp**F2) *dFexpdPON*dPONdZ 361c dfb4dZ = CcB(4) *(F4 *Fexp**F3) *dFexpdPON*dPONdZ 362c dfb5dZ = CcB(5) *(F5 *Fexp**F4) *dFexpdPON*dPONdZ 363c 364c dFggaCdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR + 365c $ dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR 366c dFggaCdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG + 367c $ dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG 368c dFggaCdZ = dfa1dZ+dfa2dZ+dfa3dZ+dfa4dZ+dfa5dZ + 369c $ dfb1dZ+dfb2dZ+dfb3dZ+dfb4dZ+dfb5dZ 370c 371c dFggaCdA = dFggaCdR + dFggaCdZ*dZdA 372c dFggaCdB = dFggaCdR + dFggaCdZ*dZdB 373c dGdGx = F1 / (F2*GRho) 374c dFggaCdGx = dFggaCdG*dGdGx 375c 376c Ec = Ec+ (Eloc*FggaC)*qwght(n) 377 ffunc(n) = ffunc(n)+ (Eloc*FggaC)*wght 378c 379c Amat(n,D1_RA) = Amat(n,D1_RA) + dElocdA*FggaC*wght + 380c $ Eloc*dFggaCdA*wght 381c if (ipol.eq.2) then 382c Amat(n,D1_RB) = Amat(n,D1_RB) + dElocdB*FggaC*wght + 383c $ Eloc*dFggaCdB*wght 384c endif 385c if (ipol.eq.1) then 386c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Eloc*dFggaCdGx*wght 387c Cmat(n,D1_GAB) = Cmat(n,D1_GAB) +F2*Eloc*dFggaCdGx*wght 388c else 389c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Eloc*dFggaCdGx*wght 390c Cmat(n,D1_GAB) = Cmat(n,D1_GAB) +F2*Eloc*dFggaCdGx*wght 391c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Eloc*dFggaCdGx*wght 392c endif 393 20 continue 394 end 395c 396#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 397#if defined(NWAD_PRINT) 398 Subroutine nwxc_GZeta_p(Zeta,GZet,dGZdz,d2GZdz,d3GZdz) 399#else 400 Subroutine nwxc_GZeta(Zeta,GZet,dGZdz,d2GZdz,d3GZdz) 401#endif 402#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 403 Subroutine nwxc_GZeta_d2(Zeta,GZet,dGZdz,d2GZdz,d3GZdz) 404#else 405 Subroutine nwxc_GZeta_d3(Zeta,GZet,dGZdz,d2GZdz,d3GZdz) 406#endif 407#include "nwad.fh" 408 Implicit Real*8(A-H,O-Z) 409C 410C Evaluate G(Zeta) and its derivatives for DFT. 411C 412 type(nwad_dble)::Zeta,GZet 413 type(nwad_dble)::OMZ,OMZ2,OMZ3,OPZ,OPZ2,OPZ3 414 Save F0,F1,F2,F3,F4,F9,F27 415 Data F0/0.0d0/,F1/1.0d0/,F2/2.0d0/,F3/3.0d0/,F4/4.0d0/, 416 $ F9/9.0d0/,F27/27.0D0/ 417 418 F1o3 = F1/F3 419 F1o9 = F1/F9 420 F4o27= F4/F27 421c 422 OMZ3 = F0 423 OPZ3 = F0 424 GZet = F0 425c dGZdz = F0 426c d2GZdz = F0 427c d3GZdz = F0 428c 429 OMZ = F1-Zeta 430 OPZ = F1+Zeta 431 OMZ2 = OMZ**2.0d0 432 OPZ2 = OPZ**2.0d0 433 if (OMZ.gt.F0) then 434 OMZ3 = OMZ**(-F1o3) 435c d2GZdz = d2GZdz-OMZ3/OMZ 436c d3GZdz = d3GZdz-OMZ3/OMZ2 437 endif 438 GZet = GZet+OMZ*OMZ3 439c dGZdz = dGZdz-OMZ3 440 if (OPZ.gt.F0) then 441 OPZ3 = OPZ**(-F1o3) 442c d2GZdz = d2GZdz-OPZ3/OPZ 443c d3GZdz = d3GZdz+OPZ3/OPZ2 444 endif 445 GZet = GZet+OPZ*OPZ3 446c dGZdz = dGZdz+OPZ3 447c 448 GZet = GZet/F2 449c dGZdz = dGZdz*F1o3 450c d2GZdz = d2GZdz*F1o9 451c d3GZdz = d3GZdz*F4o27 452 Return 453 End 454#ifndef NWAD_PRINT 455#define NWAD_PRINT 456c 457c Compile source again for the 2nd derivative case 458c 459#include "nwxc_c_sogga.F" 460#endif 461#ifndef SECOND_DERIV 462#define SECOND_DERIV 463c 464c Compile source again for the 2nd derivative case 465c 466#include "nwxc_c_sogga.F" 467#endif 468#ifndef THIRD_DERIV 469#define THIRD_DERIV 470c 471c Compile source again for the 3rd derivative case 472c 473#include "nwxc_c_sogga.F" 474#endif 475#undef NWAD_PRINT 476C> 477C> @} 478