1#ifndef SECOND_DERIV
2      Subroutine xc_gill96(tol_rho, fac, lfac, nlfac, rho, delrho,
3     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
4#else
5#include "dft2drv.fh"
6      Subroutine xc_gill96_d2(tol_rho, fac, lfac, nlfac, rho, delrho,
7     &                        Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
8     &                        qwght,ldew,func)
9#endif
10c
11C$Id$
12c
13      implicit none
14c
15c
16      double precision tol_rho, fac, Ex
17      integer nq, ipol
18      logical lfac, nlfac,ldew
19      double precision func(*)  ! value of the functional [output]
20c
21c     Charge Density
22c
23      double precision rho(nq,ipol*(ipol+1)/2)
24c
25c     Charge Density Gradient
26c
27      double precision delrho(nq,3,ipol)
28c
29c     Quadrature Weights
30c
31      double precision qwght(nq)
32c
33c     Sampling Matrices for the XC Potential
34c
35      double precision Amat(nq,ipol), Cmat(nq,*)
36#ifdef SECOND_DERIV
37      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
38#endif
39c
40      double precision BETA, ONE3, FOUR3
41#ifdef SECOND_DERIV
42      double precision SEVEN3
43#endif
44      Parameter (BETA = 1d0/137d0)
45      Parameter (ONE3 = 1d0/3d0, FOUR3 = 4d0/3d0)
46#ifdef SECOND_DERIV
47      Parameter (SEVEN3 = 7d0/3d0)
48#endif
49c
50c References:
51c
52c    Gill , Mol. Phys. 89, 433 (1996)
53c
54c***************************************************************************
55c
56      integer n
57      double precision C, rhoval, rho13, rho43, gamma, x, d1x(2),
58     &     g, d1g
59#ifdef SECOND_DERIV
60      double precision d2x(3), d2g
61#endif
62c
63c
64c     Uniform electron gas constant
65c
66      C = -(1.5d0)*(0.75d0/acos(-1d0))**(ONE3)
67c
68      if (ipol.eq.1) then
69c
70c        ======> SPIN-RESTRICTED <======
71c
72         do 10 n = 1, nq
73            if (rho(n,1).lt.tol_rho) goto 10
74            rhoval = 0.5d0*rho(n,1)
75c
76c           Spin alpha:
77c
78            rho13 = rhoval**ONE3
79            rho43 = rho13*rhoval
80c     Include factor of 4/3 in rho13 since it always appears with it
81            rho13 = FOUR3*rho13
82c
83            if (lfac) then
84               Ex = Ex + 2d0*rho43*C*qwght(n)*fac
85               if(ldew)func(n) = func(n) + 2.d0*rho43*C*fac
86               Amat(n,1) = Amat(n,1) + rho13*C*fac
87#ifdef SECOND_DERIV
88               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
89     &              + ONE3*(rho13/rhoval)*C*fac
90#endif
91            endif
92c
93            gamma = delrho(n,1,1)*delrho(n,1,1) +
94     &              delrho(n,2,1)*delrho(n,2,1) +
95     &              delrho(n,3,1)*delrho(n,3,1)
96            if (nlfac.and.dsqrt(gamma).gt.tol_rho)then
97               gamma = 0.25d0*gamma
98               x = dsqrt(gamma)/rho43
99               d1x(1) = -FOUR3*x/rhoval
100               d1x(2) = 0.5d0*x/gamma
101               g = -BETA*x*sqrt(x)
102               d1g = -1.5d0*BETA*sqrt(x)
103c
104               Ex = Ex + 2d0*rho43*g*qwght(n)*fac
105               if(ldew)func(n) = func(n) + 2.d0*rho43*g*fac
106               Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g*d1x(1))*fac
107               Cmat(n,1) = Cmat(n,1) + rho43*d1g*d1x(2)*fac
108#ifdef SECOND_DERIV
109               d2g = 0.5d0*d1g/x
110               d2x(1) = -SEVEN3*d1x(1)/rhoval
111               d2x(2) = -FOUR3*d1x(2)/rhoval
112               d2x(3) = -0.5d0*d1x(2)/gamma
113               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
114     &              + ONE3*(rho13/rhoval)*g*fac
115     &              + 2.d0*rho13*d1g*d1x(1)*fac
116     &              + rho43*(d2g*d1x(1)*d1x(1)+d1g*d2x(1))*fac
117               Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
118     &              + rho13*d1g*d1x(2)*fac
119     &              + rho43*(d2g*d1x(1)*d1x(2)+d1g*d2x(2))*fac
120               Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
121     &              + rho43*(d2g*d1x(2)*d1x(2)+d1g*d2x(3))*fac
122#endif
123            endif
124c
125 10      continue
126c
127      else
128c
129c        ======> SPIN-UNRESTRICTED <======
130c
131         do 20 n = 1, nq
132            if (rho(n,1).lt.tol_rho) goto 20
133            if (rho(n,2).lt.tol_rho) goto 25
134c
135c           Spin alpha:
136c
137            rhoval = rho(n,2)
138            rho13 = rhoval**ONE3
139            rho43 = rho13*rhoval
140c     Include factor of 4/3 in rho13 since it always appears with it
141            rho13 = FOUR3*rho13
142c
143            if (lfac) then
144               Ex = Ex + rho43*C*qwght(n)*fac
145               if(ldew)func(n) = func(n) + rho43*C*fac
146               Amat(n,1) = Amat(n,1) + rho13*C*fac
147#ifdef SECOND_DERIV
148               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
149     &              + ONE3*(rho13/rhoval)*C*fac
150#endif
151            endif
152c
153            gamma = delrho(n,1,1)*delrho(n,1,1) +
154     &              delrho(n,2,1)*delrho(n,2,1) +
155     &              delrho(n,3,1)*delrho(n,3,1)
156            if (nlfac.and.dsqrt(gamma).gt.tol_rho)then
157               x = dsqrt(gamma)/rho43
158               d1x(1) = -FOUR3*x/rhoval
159               d1x(2) = 0.5d0*x/gamma
160               g = -BETA*x*sqrt(x)
161               d1g = -1.5d0*BETA*sqrt(x)
162c
163               Ex = Ex + rho43*g*qwght(n)*fac
164               if(ldew)func(n) = func(n) + rho43*g*fac
165               Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g*d1x(1))*fac
166               Cmat(n,1) = Cmat(n,1) + rho43*d1g*d1x(2)*fac
167#ifdef SECOND_DERIV
168               d2g = 0.5d0*d1g/x
169               d2x(1) = -SEVEN3*d1x(1)/rhoval
170               d2x(2) = -FOUR3*d1x(2)/rhoval
171               d2x(3) = -0.5d0*d1x(2)/gamma
172               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
173     &              + ONE3*(rho13/rhoval)*g*fac
174     &              + 2.d0*rho13*d1g*d1x(1)*fac
175     &              + rho43*(d2g*d1x(1)*d1x(1)+d1g*d2x(1))*fac
176               Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
177     &              + rho13*d1g*d1x(2)*fac
178     &              + rho43*(d2g*d1x(1)*d1x(2)+d1g*d2x(2))*fac
179               Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
180     &              + rho43*(d2g*d1x(2)*d1x(2)+d1g*d2x(3))*fac
181#endif
182            endif
183c
184c
185 25         continue
186c
187c           Spin beta:
188c
189            if (rho(n,3).lt.tol_rho) goto 20
190c
191            rhoval = rho(n,3)
192            rho13 = rhoval**ONE3
193            rho43 = rho13*rhoval
194c     Include factor of 4/3 in rho13 since it always appears with it
195            rho13 = FOUR3*rho13
196c
197            if (lfac) then
198               Ex = Ex + rho43*C*qwght(n)*fac
199               if(ldew)func(n) = func(n) + rho43*C*fac
200               Amat(n,2) = Amat(n,2) + rho13*C*fac
201#ifdef SECOND_DERIV
202               Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
203     &              + ONE3*(rho13/rhoval)*C*fac
204#endif
205            endif
206c
207            gamma = delrho(n,1,2)*delrho(n,1,2) +
208     &              delrho(n,2,2)*delrho(n,2,2) +
209     &              delrho(n,3,2)*delrho(n,3,2)
210            if (nlfac.and.dsqrt(gamma).gt.tol_rho)then
211               x = dsqrt(gamma)/rho43
212               d1x(1) = -FOUR3*x/rhoval
213               d1x(2) = 0.5d0*x/gamma
214               g = -BETA*x*sqrt(x)
215               d1g = -1.5d0*BETA*sqrt(x)
216c
217               Ex = Ex + rho43*g*qwght(n)*fac
218               if(ldew)func(n) = func(n) + rho43*g*fac
219               Amat(n,2) = Amat(n,2) + (rho13*g+rho43*d1g*d1x(1))*fac
220               Cmat(n,3) = Cmat(n,3) + rho43*d1g*d1x(2)*fac
221#ifdef SECOND_DERIV
222               d2g = 0.5d0*d1g/x
223               d2x(1) = -SEVEN3*d1x(1)/rhoval
224               d2x(2) = -FOUR3*d1x(2)/rhoval
225               d2x(3) = -0.5d0*d1x(2)/gamma
226               Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
227     &              + ONE3*(rho13/rhoval)*g*fac
228     &              + 2.d0*rho13*d1g*d1x(1)*fac
229     &              + rho43*(d2g*d1x(1)*d1x(1)+d1g*d2x(1))*fac
230               Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
231     &              + rho13*d1g*d1x(2)*fac
232     &              + rho43*(d2g*d1x(1)*d1x(2)+d1g*d2x(2))*fac
233               Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB)
234     &              + rho43*(d2g*d1x(2)*d1x(2)+d1g*d2x(3))*fac
235#endif
236            endif
237c
238c
239 20      continue
240c
241      endif
242c
243      return
244      end
245#ifndef SECOND_DERIV
246#define SECOND_DERIV
247c
248c     Compile source again for the 2nd derivative case
249c
250#include "xc_gill96.F"
251#endif
252