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