1c
2c     Modified to handle second derivatives while reusing code
3c
4c     BGJ - 8/98
5c
6#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
7      Subroutine xc_camxlsd(tol_rho, fac, lfac, nlfac, rho, Amat, nq,
8     &                    ipol, Ex, qwght, ldew, func)
9#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
10c     For locations of 2nd derivatives of functionals in array
11#include "dft2drv.fh"
12      Subroutine xc_camxlsd_d2(tol_rho, fac, lfac, nlfac, rho, Amat,
13     &                       Amat2, nq, ipol, Ex, qwght, ldew, func)
14#else
15c     For locations of 3rd derivatives of functionals in array
16#include "dft3drv.fh"
17      Subroutine xc_camxlsd_d3(tol_rho, fac, lfac, nlfac, rho, Amat,
18     1                       Amat2, Amat3, nq, ipol, Ex, qwght, ldew,
19     2                       func)
20#endif
21c
22C$Id$
23c
24      Implicit none
25#include "errquit.fh"
26c
27#include "stdio.fh"
28c
29      integer nq, ipol
30      double precision fac, Ex
31      logical ldew, lfac, nlfac
32      double precision func(*)  ! value of the functional [output]
33c
34c     Charge Density
35c
36      double precision rho(nq,(ipol*(ipol+1))/2)
37c
38c     Quadrature Weights
39c
40      double precision qwght(nq)
41c
42c     Partial First Derivatives of the Exchange Energy Functional
43c
44      double precision Amat(nq,ipol)
45      double precision Etmp,Atmp,Ctmp,A2tmp,C2tmp,C3tmp
46c
47#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
48c
49c     Partial Second Derivatives of the Exchange Energy Functional
50c
51c      double precision Amat2(nq,*)
52      double precision Amat2(nq,NCOL_AMAT2)
53#endif
54#ifdef THIRD_DERIV
55c
56c     Partial Third Derivatives of the Exchange Energy Functional
57c
58      double precision Amat3(nq,NCOL_AMAT3)
59      double precision A3tmp, C4tmp, C5tmp, C6tmp
60      double precision rhom23
61#endif
62c
63c     Compute the partial derivatives of the exchange functional of Dirac.
64c
65      double precision P1, P2, P3, P4, tol_rho
66c
67c     P1 =       -(3/PI)**(1/3)
68c     P2 = -(3/4)*(3/PI)**(1/3)
69c     P3 =       -(6/PI)**(1/3)
70c     P4 = -(3/4)*(6/PI)**(1/3)
71c
72      Parameter (P1 = -0.9847450218426959D+00)
73      Parameter (P2 = -0.7385587663820219D+00)
74      Parameter (P3 = -0.1240700981798799D+01)
75      Parameter (P4 = -0.9305257363490993D+00)
76      double precision rho13, rho32, rho33, one_third
77      Parameter (one_third = 1.d0/3.d0)
78      double precision two_ninth
79      Parameter (two_ninth = 2.0d0/9.0d0)
80      integer n
81c
82      if (ipol.eq.1)then
83c
84c        ======> SPIN-RESTRICTED <======
85c
86         do 10 n = 1, nq
87            if (rho(n,1).gt.tol_rho)then
88             rho13=rho(n,1)**one_third
89             Etmp = rho(n,1)*rho13*P2*fac
90             if(ldew)func(n) = func(n) + rho(n,1)*rho13*fac*P2
91             Atmp = rho13*P1*fac
92             Ctmp = 0.0d0
93#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
94             A2tmp = (rho13/rho(n,1))*2.0d0*one_third*P1*fac
95             C2tmp = 0.0d0
96             C3tmp = 0.0d0
97#endif
98#ifdef THIRD_DERIV
99             rhom23 = rho13/rho(n,1)
100             A3tmp = (rhom23/rho(n,1))*(-4.0d0)*two_ninth*P1*fac
101             C4tmp = 0.0d0
102             C5tmp = 0.0d0
103             C6tmp = 0.0d0
104#endif
105#ifdef THIRD_DERIV
106             call xc_att_xc_d3(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
107     &           C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp)
108c
109             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
110c
111             Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp
112#elif defined(SECOND_DERIV)
113             call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
114     &           C2tmp,C3tmp)
115c
116             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
117#else
118             call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
119#endif
120             Ex = Ex + qwght(n)*Etmp
121             Amat(n,1) = Amat(n,1) + Atmp
122            endif
123   10    continue
124c
125      else
126c
127c        ======> SPIN-UNRESTRICTED <======
128c
129         do 20 n = 1,nq
130           if (rho(n,1).gt.tol_rho)then
131             rho32=0.0d0
132             rho33=0.0d0
133             if (rho(n,2).gt.tol_rho) rho32=rho(n,2)**one_third
134             if (rho(n,3).gt.tol_rho) rho33=rho(n,3)**one_third
135c ---- alpha ----
136             Etmp = rho32*rho(n,2)*P4*fac
137c
138             Atmp = P3*rho32*fac
139             Ctmp = 0.0d0
140#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
141             A2tmp = 0.0d0
142             C2tmp = 0.0d0
143             C3tmp = 0.0d0
144c
145             if (rho(n,2).gt.tol_rho) then
146               A2tmp = one_third*P3*rho32/rho(n,2)*fac
147             endif
148#endif
149#ifdef THIRD_DERIV
150             A3tmp = 0.0d0
151             C4tmp = 0.0d0
152             C5tmp = 0.0d0
153             C6tmp = 0.0d0
154c
155             if (rho(n,2).gt.tol_rho) then
156               A3tmp = -two_ninth*P3*rho32/(rho(n,2)**2)*fac
157             endif
158#endif
159#ifdef THIRD_DERIV
160             if (rho(n,2).gt.tol_rho) then
161               call xc_att_xc_d3(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
162     &           C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp)
163             endif
164c
165             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
166c
167             Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp
168#elif defined(SECOND_DERIV)
169             if (rho(n,2).gt.tol_rho) then
170               call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
171     &           C2tmp,C3tmp)
172             endif
173c
174             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
175#else
176             if (rho(n,2).gt.tol_rho) then
177               call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
178             endif
179#endif
180             Amat(n,1) = Amat(n,1) + Atmp
181             Ex = Ex + qwght(n)*Etmp
182c ---- beta ----
183             Etmp = rho33*rho(n,3)*P4*fac
184             Atmp = P3*rho33*fac
185             Ctmp = 0.0d0
186#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
187             A2tmp = 0.0d0
188             C2tmp = 0.0d0
189             C3tmp = 0.0d0
190c
191             if (rho(n,3).gt.tol_rho) then
192               A2tmp = one_third*P3*rho33/rho(n,3)*fac
193             end if
194#endif
195#ifdef THIRD_DERIV
196             A3tmp = 0.0d0
197             C4tmp = 0.0d0
198             C5tmp = 0.0d0
199             C6tmp = 0.0d0
200c
201             if (rho(n,3).gt.tol_rho) then
202               A3tmp = -two_ninth*P3*rho33/(rho(n,3)**2)*fac
203             endif
204#endif
205#ifdef THIRD_DERIV
206             if (rho(n,3).gt.tol_rho) then
207               call xc_att_xc_d3(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
208     &           C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp)
209             endif
210c
211             Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
212c
213             Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + A3tmp
214#elif defined(SECOND_DERIV)
215             if (rho(n,3).gt.tol_rho) then
216               call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
217     &           C2tmp,C3tmp)
218             end if
219c
220             Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
221#else
222             if (rho(n,3).gt.tol_rho) then
223               call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
224             endif
225#endif
226             Amat(n,2) = Amat(n,2) + Atmp
227             Ex = Ex + qwght(n)*Etmp
228c
229             if (ldew)func(n) = func(n) + ( rho32*rho(n,2) +
230     &                                      rho33*rho(n,3)   )*P4*fac
231            endif
232   20    continue
233c
234      endif
235c
236      return
237      end
238c
239#ifndef SECOND_DERIV
240#define SECOND_DERIV
241c     Compile source again for the 2nd derivative case
242#include "xc_camxlsd.F"
243#endif
244c
245#ifndef THIRD_DERIV
246#define THIRD_DERIV
247c     Compile source again for the 3rd derivative case
248#include "xc_camxlsd.F"
249#endif
250