1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3C> \ingroup nwxc
4C> @{
5C>
6C> \file nwxc_c_lyp.F
7C> The LYP correlation functional
8C>
9C> @}
10#endif
11#endif
12C>
13C> \ingroup nwxc_priv
14C> @{
15C>
16C> \brief Evaluate the LYP correlation functional
17C>
18C> Evaluate the  LYP correlation functional [1-4].
19C>
20C> ### References ###
21C>
22C> [1] C. Lee, W. Yang, R.G. Parr,
23C> "Development of the Colle-Salvetti correlation-energy formula into
24C> a functional of the electron density", Phys. Rev. B <b>37</b>,
25C> 785-789 (1988), DOI:
26C> <a href="https://doi.org/10.1103/PhysRevB.37.785">
27C> 10.1103/PhysRevB.37.785</a>.
28C>
29C> [2]  R. Colle, O. Salvetti,
30C> "Approximate calculation of the correlation energy for the closed
31C> shells", Theor. Chim. Acta <b>37</b>, 329-334 (1975), DOI:
32C> <a href="https://doi.org/10.1007/BF01028401">
33C> 10.1007/BF01028401</a>.
34C>
35C> [3] B. Miehlich, A. Savin, H. Stoll, H. Preuss,
36C> "Results obtained with the correlation energy density functionals of
37C> Becke, and Lee, Yang and Parr", Chem. Phys. Lett. <b>157</b>, 200-206
38C> (1989), DOI: <a href="https://doi.org/10.1016/0009-2614(89)87234-3">
39C> 10.1016/0009-2614(89)87234-3</a>.
40C>
41C> [4] B.G. Johnson, P.M.W. Gill, J.A. Pople,
42C> "The performance of a family of density functional methods",
43C> J. Chem. Phys. <b>98</b>, 5612-5626 (1993), DOI:
44C> <a href="https://doi.org/10.1063/1.464906">
45C> 10.1063/1.464906</a>.
46C>
47#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
48#if defined(NWAD_PRINT)
49      Subroutine nwxc_c_lyp_p(tol_rho, ipol, nq, wght, rho, rgamma,
50     &                      func)
51#else
52      Subroutine nwxc_c_lyp(tol_rho, ipol, nq, wght, rho, rgamma,
53     &                      func)
54#endif
55#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
56      Subroutine nwxc_c_lyp_d2(tol_rho, ipol, nq, wght, rho,
57     &                         rgamma, func)
58#else
59      Subroutine nwxc_c_lyp_d3(tol_rho, ipol, nq, wght, rho,
60     &                         rgamma, func)
61#endif
62c
63C$Id$
64c
65#include "nwad.fh"
66      implicit none
67c
68#include "nwxc_param.fh"
69c
70c     Input and other parameters
71c
72      double precision tol_rho !< [Input] The lower limit on the density
73      integer ipol             !< [Input] The number of spin channels
74      integer nq               !< [Input] The number of points
75      double precision wght    !< [Input] The weight of the functional
76c
77c     Charge Density
78c
79      type(nwad_dble)::rho(nq,*)    !< [Input] The density
80c
81c     Charge Density Gradient
82c
83      type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients
84c
85c     Sampling Matrices for the XC Potential
86c
87      type(nwad_dble)::func(nq)    !< [Output] The value of the functional
88c     double precision Amat(nq,*)   !< [Output] The derivative wrt rho
89c     double precision Cmat(nq,*)   !< [Output] The derivative wrt rgamma
90#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
91c
92c     Sampling Matrices for the XC Kernel
93c
94c     double precision Amat2(nq,*)  !< [Output] The 2nd derivative wrt rho
95c     double precision Cmat2(nq,*)  !< [Output] The 2nd derivative wrt rgamma
96                                    !< and possibly rho
97#endif
98#if defined(THIRD_DERIV)
99c
100c     Sampling Matrices for the XC Kernel
101c
102c     double precision Amat3(nq,*)  !< [Output] The 3rd derivative wrt rho
103c     double precision Cmat3(nq,*)  !< [Output] The 3rd derivative wrt rgamma
104                                    !< and possibly rho
105#endif
106      double precision F13, F43, F113, F83, F53, F19, F79, P1,
107     &                 A, B, C, D
108c
109      Parameter (F13 = 1.D0/3.D0, F43 = 4.D0*F13, F113 = 11.D0*F13,
110     &           F83 = 8.D0*F13, F53 = 5.D0*F13, F19 = 1.D0/9.D0,
111     &           F79 = 7.D0*F19)
112c
113c     P1 = 2**(11/3)*(3/10)*(3*PI**2)**(2/3)
114c
115      Parameter (P1 = 0.3646239897876487D+02)
116c
117c     Colle-Salvetti Empirical Parameters
118c
119      Parameter (A = 0.04918D0)
120      Parameter (B = 0.13200D0)
121      Parameter (C = 0.25330D0)
122      Parameter (D = 0.34900D0)
123c
124c     Compute the partial derivatives of the correlation functional of
125c     Lee, Yang and Parr.
126c
127c     References:
128c
129c     Colle & Salvetti, Theor. Chim. Acta 37, 329 (1975)
130c     Lee, Yang & Parr, Phys. Rev. B 37, 785 (1988)
131c     Miehlich, Savin, Stoll & Preuss, Chem. Phys. Lett. 157, 200 (1989)
132c     Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993)
133c
134      integer n
135      double precision c1, c2, ab, ratcd
136      type(nwad_dble)::rrho, rhoa, rhob, d1t(2)
137      type(nwad_dble)::rhoa83, rhob83
138      type(nwad_dble)::rrho2, rhoa2, rhob2, rhoab, rho2
139      type(nwad_dble)::h1, h2, h3, om, t, tm_in
140      type(nwad_dble)::de, de11, de47
141      type(nwad_dble)::zero
142      type(nwad_dble)::xrarb
143c
144      double precision d1xrarb(2)
145      double precision d1tm_in(2)
146c
147      type(nwad_dble)::gaa, gab, gbb
148      type(nwad_dble)::f1, f2, d1f1(2), d1f2(2), f, d1f(5),
149     &     d2fgaa(2), d2fgab(2), d2fgbb(2)
150c
151c     Coefficients of first two terms in LYP functional and other
152c     commonly occurring factors
153c
154      zero = 0.0d0
155      c1 = -4.0d0*a
156      c2 = -P1*a*b
157      ab = a*b
158      ratcd = c/d
159c
160      if (ipol.eq.1)then
161c
162c        ======> SPIN-RESTRICTED <======
163c
164         do 10 n = 1, nq
165            if (rho(n,R_T).lt.tol_rho)goto 10
166            rrho = 1.0d0/rho(n,R_T)
167            rhoa = 0.5d0*rho(n,R_T)
168            rrho2 = rrho*rrho
169            rho2 = 1.0d0/rrho2
170            rhoa2 = rhoa*rhoa
171            rhoab = rhoa2
172c           rhoa53 = rhoa**F53
173            rhoa83 = rhoa**F83
174c
175            h2 = d*rrho**F13
176c           d1h2 = -F13*h2*rrho
177c
178            h3 = ratcd*h2
179c           d1h3 = ratcd*d1h2
180c
181            h1 = 1d0/(1d0+h2)
182c           d1h1 = -h1*h1*d1h2
183c
184            om = exp(-h3)*h1*rrho**F113
185c           t = d1h3+h1*d1h2+F113*rrho
186c           d1om = -om*t
187c
188            de = h3+h1*h2
189c           d1de = d1h3 + d1h1*h2 + h1*d1h2
190c
191            f1 = h1*rhoab*rrho
192c           d1f1(1) = d1h1*rhoab*rrho - h1*rhoab*rrho2
193c           d1f1(1) = d1f1(1) + h1*rhoa*rrho
194c
195            xrarb = rhoab*(rhoa83+rhoa83)
196            f2 = om*xrarb
197c           d1xrarb(1) = rhoa*(F113*rhoa83+rhoa83)
198c           d1f2(1) = d1om*xrarb + om*d1xrarb(1)
199c
200c           gaa =(delrho(n,1,1)*delrho(n,1,1) +
201c    &            delrho(n,2,1)*delrho(n,2,1) +
202c    &            delrho(n,3,1)*delrho(n,3,1))*0.25d0
203            gaa = rgamma(n,G_TT)*0.25d0
204c
205            de11 = de - 11.0d0
206            de47 = 47.0d0 - 7.0d0*de
207c
208c Daniel (10-30-12): tm_in is what I call Qi (the inside term)
209            tm_in = 1.0d0 - 3.0d0*de - de11*rhoa*rrho
210c Daniel (10-23-12): "t" is what I call Q or S.
211            t = F19*rhoab*tm_in - rhoa2
212c Daniel (10-30-12): Derivatives of the inside term
213c           d1tm_in(1) = -(3.0d0+rhoa*rrho)*d1de-de11*rhoa*rrho2
214c           d1tm_in(2) = -(3.0d0+rhoa*rrho)*d1de+de11*rhoa*rrho2
215c Daniel (10-23-12): d1t(1) is the derivative with respect to rhoa,
216c and d1t(2) is the derivative with respect to rhob.
217c           d1t(1) = F19*( rhoa*tm_in + rhoab*d1tm_in(1) )
218c           d1t(2) = F19*( rhoa*tm_in + rhoab*d1tm_in(2) )
219c    &             - 2.0d0*rhoa
220c
221            d1f(3) = -ab*om*t
222c
223c           d2fgaa(1) = -ab*( d1om*t + om*d1t(1) )
224c           d2fgaa(2) = -ab*( d1om*t + om*d1t(2) )
225c
226c Daniel (10-23-12): "t" is what I call R.
227            t = F19*rhoab*de47-F43*rho2
228c           d1t(1) = F19*rhoa*de47 - F79*rhoab*d1de - F83*rho(n,R_T)
229            d1f(4) = -ab*om*t
230c           d2fgab(1) = -ab*(d1om*t+om*d1t(1))
231c
232c           d2fgbb(1) = d2fgaa(2)
233c
234            f = c1*f1 + c2*f2 + gaa*(2d0*d1f(3) + d1f(4))
235c           d1f(1) = c1*d1f1(1) + c2*d1f2(1)
236c    &             + gaa*(d2fgaa(1) + d2fgab(1) + d2fgbb(1))
237c
238            func(n) = func(n) + f*wght
239c           Amat(n,D1_RA) = Amat(n,D1_RA) + d1f(1)*wght
240c           Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*wght
241c           Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(4)*wght
242c
243   10    continue
244      else
245c
246c        ======> SPIN-UNRESTRICTED <======
247c
248         do 20 n = 1,nq
249            if (rho(n,R_A)+rho(n,R_B).lt.tol_rho)goto 20
250            rrho = 1.0d0/(rho(n,R_A)+rho(n,R_B))
251c           rhoa = max(zero,rho(n,R_A))
252c           rhob = max(zero,rho(n,R_B))
253            if (rho(n,R_A).ge.0.5d0*tol_rho) then
254              rhoa = rho(n,R_A)
255            else
256              rhoa = 0.0d0
257              rrho = 1.0d0/rho(n,R_B)
258            endif
259            if (rho(n,R_B).ge.0.5d0*tol_rho) then
260              rhob = rho(n,R_B)
261            else
262              rhob = 0.0d0
263              rrho = 1.0d0/rho(n,R_A)
264            endif
265            rrho2 = rrho*rrho
266            rho2 = 1d0/rrho2
267            rhoa2 = rhoa*rhoa
268            rhob2 = rhob*rhob
269            rhoab = rhoa*rhob
270c           rhoa53 = rhoa**F53
271c           rhob53 = rhob**F53
272            rhoa83 = rhoa**F83
273            rhob83 = rhob**F83
274c
275            h2 = d*rrho**F13
276c           d1h2 = -F13*h2*rrho
277c
278            h3 = ratcd*h2
279c           d1h3 = ratcd*d1h2
280c
281            h1 = 1d0/(1d0+h2)
282c           d1h1 = -h1*h1*d1h2
283c
284            om = exp(-h3)*h1*rrho**F113
285c           t = d1h3+h1*d1h2+F113*rrho
286c           d1om = -om*t
287c
288            de = h3+h1*h2
289c           d1de = d1h3 + d1h1*h2 + h1*d1h2
290c
291c Daniel (3-21-13): f1 is J/(-4*a) in my notes.
292            f1 = h1*rhoab*rrho
293c           d1f1(1) = d1h1*rhoab*rrho - h1*rhoab*rrho2
294c           d1f1(2) = d1f1(1)
295c           d1f1(1) = d1f1(1) + h1*rhob*rrho
296c           d1f1(2) = d1f1(2) + h1*rhoa*rrho
297c
298c Daniel (10-30-12): Define xrarb here
299            xrarb = rhoab*(rhoa83+rhob83)
300c
301            f2 =om*xrarb
302c
303c           d1xrarb(1) = rhob*(F113*rhoa83+rhob83)
304c           d1xrarb(2) = rhoa*(F113*rhob83+rhoa83)
305c
306c           d1f2(1) = d1om*xrarb
307c           d1f2(2) = d1f2(1)
308c           d1f2(1) = d1f2(1) + om*d1xrarb(1)
309c           d1f2(2) = d1f2(2) + om*d1xrarb(2)
310c
311            gaa = rgamma(n,G_AA)
312            gbb = rgamma(n,G_BB)
313            gab = rgamma(n,G_AB)
314            if (rho(n,R_A).lt.0.5d0*tol_rho) then
315              gaa = 0.0d0
316              gab = 0.0d0
317            endif
318            if (rho(n,R_B).lt.0.5d0*tol_rho) then
319              gbb = 0.0d0
320              gab = 0.0d0
321            endif
322c
323            de11 = de - 11d0
324            de47 = 47d0 - 7d0*de
325c Daniel (10-30-12): tm_in is what I call Qi (the inside term)
326            tm_in = 1.0d0 - 3.0d0*de - de11*rhoa*rrho
327            t = F19*rhoab*tm_in-rhob2
328c Daniel (10-30-12): Derivatives of the inside term
329c           d1tm_in(1) = -(3.0d0+rhoa*rrho)*d1de-de11*rhob*rrho2
330c           d1tm_in(2) = -(3.0d0+rhoa*rrho)*d1de+de11*rhoa*rrho2
331c           d1t(1) = F19*(rhob*tm_in + rhoab*d1tm_in(1))
332c           d1t(2) = F19*(rhoa*tm_in + rhoab*d1tm_in(2)) - 2.0d0*rhob
333c
334            d1f(3) = -ab*om*t
335c           d2fgaa(1) = -ab*(d1om*t+om*d1t(1))
336c           d2fgaa(2) = -ab*(d1om*t+om*d1t(2))
337c
338            t = F19*rhoab*de47-F43*rho2
339c           d1t(1) = F19*rhob*de47 - F79*rhoab*d1de
340c    &             - F83*(rho(n,R_A)+rho(n,R_B))
341c           d1t(2) = F19*rhoa*de47 - F79*rhoab*d1de
342c    &             - F83*(rho(n,R_A)+rho(n,R_B))
343            d1f(4) = -ab*om*t
344c           d2fgab(1) = -ab*(d1om*t+om*d1t(1))
345c           d2fgab(2) = -ab*(d1om*t+om*d1t(2))
346c Daniel (3-21-13): tm_in is what I call Si (the inside term)
347            tm_in = 1.0d0 - 3.0d0*de - de11*rhob*rrho
348            t = F19*rhoab*tm_in - rhoa2
349c Daniel (10-30-12): Derivatives of the inside term
350c           d1tm_in(1) = -(3.0d0+rhob*rrho)*d1de + de11*rhob*rrho2
351c           d1tm_in(2) = -(3.0d0+rhob*rrho)*d1de - de11*rhoa*rrho2
352c           d1t(1) = F19*( rhob*tm_in + rhoab*d1tm_in(1) )
353c    1             - 2.0d0*rhoa
354c           d1t(2) = F19*( rhoa*tm_in + rhoab*d1tm_in(2) )
355            d1f(5) = -ab*om*t
356c           d2fgbb(1) = -ab*( d1om*t + om*d1t(1) )
357c           d2fgbb(2) = -ab*( d1om*t + om*d1t(2) )
358c
359            f = c1*f1 + c2*f2 + gaa*d1f(3) + gab*d1f(4) + gbb*d1f(5)
360c           d1f(1) = c1*d1f1(1) + c2*d1f2(1)
361c    &             + gaa*d2fgaa(1) + gab*d2fgab(1) + gbb*d2fgbb(1)
362c           d1f(2) = c1*d1f1(2) + c2*d1f2(2)
363c    &             + gaa*d2fgaa(2) + gab*d2fgab(2) + gbb*d2fgbb(2)
364c
365            func(n) = func(n) + f*wght
366c           Amat(n,D1_RA) = Amat(n,D1_RA) + d1f(1)*wght
367c           Amat(n,D1_RB) = Amat(n,D1_RB) + d1f(2)*wght
368c           Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + d1f(3)*wght
369c           Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + d1f(4)*wght
370c           Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + d1f(5)*wght
371c
372   20    continue
373      endif
374      return
375      end
376#ifndef NWAD_PRINT
377#define NWAD_PRINT
378c
379c     Compile source again for Maxima
380c
381#include "nwxc_c_lyp.F"
382#endif
383#ifndef SECOND_DERIV
384#define SECOND_DERIV
385c
386c     Compile source again for the 2nd derivative case
387c
388#include "nwxc_c_lyp.F"
389#endif
390#ifndef THIRD_DERIV
391#define THIRD_DERIV
392c
393c     Compile source again for the 3rd derivative case
394c
395#include "nwxc_c_lyp.F"
396#endif
397#undef NWAD_PRINT
398C> @}
399