1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3C> \ingroup nwxc
4C> @{
5C>
6C> \file nwxc_x_camlsd.F
7C> The CAM-LSD exchange functional
8C>
9C> @}
10#endif
11#endif
12C> \ingroup nwxc_priv
13C> @{
14C>
15C> \brief Evaluate the CAM-LSD exchange functional
16C>
17C> Evaluate the CAM-LSD functional [1,2]. This routine is
18C> also used to implement CAM-B3LYP.
19C>
20C> ### References ###
21C>
22C> [1] T. Yanai, D.P. Tew, N.C. Handy,
23C> "A new hybrid exchange-correlation functional using the Coulomb-attenuating
24C> method (CAM-B3LYP)",
25C> Chem. Phys. Lett. <b>393</b>, 51-57 (2004), DOI:
26C> <a href="https://doi.org/10.1016/j.cplett.2004.06.011">
27C> 10.1016/j.cplett.2004.06.011</a>.
28C>
29C> [2] A.D. Becke,
30C> "Density-functional exchange-energy approximation with correct
31C> asymptotic behavior",
32C> Phys. Rev. A <b>38</b>, 3098-3100 (1998), DOI:
33C> <a href="https://doi.org/10.1103/PhysRevA.38.3098">
34C> 10.1103/PhysRevA.38.3098</a>.
35C>
36c
37c     Modified to handle second derivatives while reusing code
38c
39c     BGJ - 8/98
40c
41#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
42#if defined(NWAD_PRINT)
43      Subroutine nwxc_x_camlsd_p(param, tol_rho, ipol, nq, wght, rho,
44     +                           func)
45#else
46      Subroutine nwxc_x_camlsd(param, tol_rho, ipol, nq, wght, rho,
47     +                         func)
48#endif
49#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
50c     For locations of 2nd derivatives of functionals in array
51      Subroutine nwxc_x_camlsd_d2(param, tol_rho, ipol, nq, wght, rho,
52     +                            func)
53#else
54      Subroutine nwxc_x_camlsd_d3(param, tol_rho, ipol, nq, wght, rho,
55     +                            func)
56#endif
57c
58C$Id$
59c
60#include "nwad.fh"
61      Implicit none
62c
63#include "nwxc_param.fh"
64c
65#if defined(NWAD_PRINT)
66#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
67      type(nwad_dble)::param(*)!< [Input] Parameters of functional
68#else
69      double precision param(*)!< [Input] Parameters of functional
70#endif
71#else
72      double precision param(*)!< [Input] Parameters of functional
73                               !< - param(1): \f$ \alpha_{CAM} \f$
74                               !< - param(2): \f$ \beta_{CAM} \f$
75                               !< - param(3): \f$ \omega_{CAM} \f$
76#endif
77      double precision tol_rho  !< [Input] The lower limit on the density
78      integer nq                !< [Input] The number of points
79      integer ipol              !< [Input] The number of spin channels
80      double precision wght     !< [Input] The weight of the functional
81c
82c     Charge Density
83c
84      type(nwad_dble)::rho(nq,*) !< [Input] The density
85c
86c     The Exchange Energy Functional
87c
88      type(nwad_dble)::func(nq)  !< [Output] The value of the functional
89c
90c     Partial First Derivatives of the Exchange Energy Functional
91c
92c     double precision Amat(nq,*) !< [Output] 1st order partial derivatives
93c
94#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
95c
96c     Partial Second Derivatives of the Exchange Energy Functional
97c
98c     double precision Amat2(nq,*) !< [Output] 2nd order partial derivatives
99#endif
100#if defined(THIRD_DERIV)
101c
102c     Partial Third Order Derivatives of the Exchange Energy Functional
103c
104c     double precision Amat3(nq,*) !< [Output] 3rd order partial derivatives
105#endif
106c
107c     Compute the partial derivatives of the exchange functional of Dirac.
108c
109      double precision Atmp,Ctmp,A2tmp,C2tmp,C3tmp
110      type(nwad_dble):: Etmp
111      double precision A3tmp, C4tmp, C5tmp, C6tmp
112      double precision rhom23
113      double precision P1, P2, P3, P4
114c
115c     P1 =       -(3/PI)**(1/3)
116c     P2 = -(3/4)*(3/PI)**(1/3)
117c     P3 =       -(6/PI)**(1/3)
118c     P4 = -(3/4)*(6/PI)**(1/3)
119c
120      Parameter (P1 = -0.9847450218426959D+00)
121      Parameter (P2 = -0.7385587663820219D+00)
122      Parameter (P3 = -0.1240700981798799D+01)
123      Parameter (P4 = -0.9305257363490993D+00)
124      double precision  one_third,two_ninth
125      type(nwad_dble):: rho13, rho32, rho33
126      Parameter (one_third = 1.d0/3.d0)
127      Parameter (two_ninth = 2.d0/9.d0)
128      integer n
129c
130      if (ipol.eq.1)then
131c
132c        ======> SPIN-RESTRICTED <======
133c
134         do 10 n = 1, nq
135            if (rho(n,R_T).gt.tol_rho)then
136             rho13=rho(n,R_T)**one_third
137             Etmp = rho(n,R_T)*rho13*P2*wght
138c            Atmp = rho13*P1*wght
139c            Ctmp = 0.d0
140#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
141c            A2tmp = (rho13/rho(n,R_T))*2.0d0*one_third*P1*wght
142c            C2tmp = 0.d0
143c            C3tmp = 0.d0
144#endif
145#if defined(THIRD_DERIV)
146c            rhom23 = rho13/rho(n,R_T)
147c            A3tmp = (rhom23/rho(n,R_T))*(-4.0d0)*two_ninth*P1*wght
148c            C4tmp = 0.0d0
149c            C5tmp = 0.0d0
150c            C6tmp = 0.0d0
151#endif
152#if defined(THIRD_DERIV)
153             call nwxc_x_att_d3(param,tol_rho,rho(n,R_T),ipol,
154     &            Etmp)
155c
156c            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
157c
158c            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp
159#elif defined(SECOND_DERIV)
160
161             call nwxc_x_att_d2(param,tol_rho,rho(n,R_T),ipol,
162     &            Etmp)
163c            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
164#else
165#if defined(NWAD_PRINT)
166             call nwxc_x_att_p(param,tol_rho,rho(n,R_T),ipol,
167     &            Etmp)
168#else
169             call nwxc_x_att(param,tol_rho,rho(n,R_T),ipol,
170     &            Etmp)
171#endif
172#endif
173             func(n) = func(n) + Etmp
174c            Amat(n,D1_RA) = Amat(n,D1_RA) + Atmp
175            endif
176   10    continue
177c
178      else
179c
180c        ======> SPIN-UNRESTRICTED <======
181c
182         do 20 n = 1,nq
183             if (rho(n,R_A).gt.0.5d0*tol_rho) then
184               rho32=rho(n,R_A)**one_third
185             else
186               rho32=0.0d0
187             endif
188             if (rho(n,R_B).gt.0.5d0*tol_rho) then
189               rho33=rho(n,R_B)**one_third
190             else
191               rho33=0.0d0
192             endif
193c
194             if (rho(n,R_A).gt.0.5d0*tol_rho) then
195               Etmp = rho32*rho(n,R_A)*P4*wght
196             else
197               Etmp = 0.0d0
198             endif
199c            Atmp = P3*rho32*wght
200c            Ctmp = 0.d0
201#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
202c            A2tmp = 0.d0
203c            C2tmp = 0.d0
204c            C3tmp = 0.d0
205c            if (rho(n,R_A).gt.0.5d0*tol_rho) then
206c              A2tmp = one_third*P3*rho32/rho(n,R_A)*wght
207c            end if
208#endif
209#if defined(THIRD_DERIV)
210c            A3tmp = 0.0d0
211c            C4tmp = 0.0d0
212c            C5tmp = 0.0d0
213c            C6tmp = 0.0d0
214c
215c            if (rho(n,R_A).gt.0.5d0*tol_rho) then
216c              A3tmp = -two_ninth*P3*rho32/(rho(n,R_A)**2)*wght
217c            endif
218#endif
219#if defined(THIRD_DERIV)
220             if (rho(n,R_A).gt.0.5d0*tol_rho) then
221               call nwxc_x_att_d3(param,tol_rho,rho(n,R_A),ipol,
222     &              Etmp)
223             endif
224c
225c            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
226c
227c            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp
228#elif defined(SECOND_DERIV)
229             if (rho(n,R_A).gt.0.5d0*tol_rho) then
230               call nwxc_x_att_d2(param,tol_rho,rho(n,R_A),ipol,
231     &              Etmp)
232             end if
233c            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
234#else
235#if defined(NWAD_PRINT)
236             if (rho(n,R_A).gt.0.5d0*tol_rho) then
237               call nwxc_x_att_p(param,tol_rho,rho(n,R_A),ipol,
238     &              Etmp)
239             endif
240#else
241             if (rho(n,R_A).gt.0.5d0*tol_rho) then
242               call nwxc_x_att(param,tol_rho,rho(n,R_A),ipol,
243     &              Etmp)
244             endif
245#endif
246#endif
247             func(n) = func(n) + Etmp
248c            Amat(n,D1_RA) = Amat(n,D1_RA) + Atmp
249c
250c            Beta spin channel
251c
252             if (rho(n,R_B).gt.0.5d0*tol_rho) then
253               Etmp = rho33*rho(n,R_B)*P4*wght
254             else
255               Etmp = 0.0d0
256             endif
257c            Atmp = P3*rho33*wght
258c            Ctmp = 0.d0
259#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
260c            A2tmp = 0.d0
261c            C2tmp = 0.d0
262c            C3tmp = 0.d0
263c            if (rho(n,R_B).gt.0.5d0*tol_rho) then
264c              A2tmp = one_third*P3*rho33/rho(n,R_B)*wght
265c            end if
266#endif
267#if defined(THIRD_DERIV)
268c            A3tmp = 0.0d0
269c            C4tmp = 0.0d0
270c            C5tmp = 0.0d0
271c            C6tmp = 0.0d0
272c
273c            if (rho(n,R_B).gt.0.5d0*tol_rho) then
274c              A3tmp = -two_ninth*P3*rho33/(rho(n,R_B)**2)*wght
275c            endif
276#endif
277#if defined(THIRD_DERIV)
278             if (rho(n,R_B).gt.0.5d0*tol_rho) then
279               call nwxc_x_att_d3(param,tol_rho,rho(n,R_B),ipol,
280     &              Etmp)
281             endif
282c
283c            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
284c
285c            Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + A3tmp
286#elif defined(SECOND_DERIV)
287             if (rho(n,R_B).gt.0.5d0*tol_rho) then
288               call nwxc_x_att_d2(param,tol_rho,rho(n,R_B),ipol,
289     &              Etmp)
290             end if
291c            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
292#else
293#if defined(NWAD_PRINT)
294             if (rho(n,R_B).gt.0.5d0*tol_rho) then
295               call nwxc_x_att_p(param,tol_rho,rho(n,R_B),ipol,
296     &              Etmp)
297             endif
298#else
299             if (rho(n,R_B).gt.0.5d0*tol_rho) then
300               call nwxc_x_att(param,tol_rho,rho(n,R_B),ipol,
301     &              Etmp)
302             endif
303#endif
304#endif
305             func(n) = func(n) + Etmp
306c            Amat(n,D1_RB) = Amat(n,D1_RB) + Atmp
307c
308c            func(n) = func(n) + ( rho32*rho(n,R_A) +
309c    &                             rho33*rho(n,R_B)   )*P4*wght
310   20    continue
311c
312      endif
313c
314      return
315      end
316#ifndef NWAD_PRINT
317#define NWAD_PRINT
318c
319c     Compile source again for the 2nd derivative case
320c
321#include "nwxc_x_camlsd.F"
322#endif
323#ifndef SECOND_DERIV
324#define SECOND_DERIV
325c
326c     Compile source again for the 2nd derivative case
327c
328#include "nwxc_x_camlsd.F"
329#endif
330#ifndef THIRD_DERIV
331#define THIRD_DERIV
332c
333c     Compile source again for the 3rd derivative case
334c
335#include "nwxc_x_camlsd.F"
336#endif
337#undef NWAD_PRINT
338C> @}
339