1c     SSB-D exchange functional part 1
2c         (the one that depends on s)
3c
4c     References:
5c     [a] J.P. Perdew, K. Burke, and M. Ernzerhof, PRL 77, 3865 (1996).
6c     [b] M. Swart, M. Sola, and F.M. Bickelhaupt, JCP 131, 094103 (2009).
7c
8#ifndef SECOND_DERIV
9      Subroutine xc_ssbD_1(tol_rho, fac, lfac, nlfac, rho, delrho,
10     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
11#else
12      Subroutine xc_ssbD_1_d2(tol_rho, fac, lfac, nlfac, rho, delrho,
13     &                        Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
14     &                        qwght,ldew,func)
15#endif
16c
17c$Id$
18c
19      implicit none
20c
21#include "dft2drv.fh"
22c
23      double precision fac, Ex
24      integer nq, ipol
25      logical lfac, nlfac,ldew
26      double precision func(*)  ! value of the functional [output]
27c
28c     Charge Density & Its Cube Root
29c
30      double precision rho(nq,ipol*(ipol+1)/2)
31c
32c     Charge Density Gradient
33c
34      double precision delrho(nq,3,ipol)
35c
36c     Quadrature Weights
37c
38      double precision qwght(nq)
39c
40c     Sampling Matrices for the XC Potential & Energy
41c
42      double precision amat(nq,ipol), cmat(nq,*)
43#ifdef SECOND_DERIV
44      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
45#endif
46c
47      double precision tol_rho, pi
48      double precision rA, rB, rC, rD, rE, rU
49      double precision C, Cs
50      double precision F43, F13
51#ifdef SECOND_DERIV
52      double precision F73
53#endif
54      parameter (rA=1.079966d0, rB=0.197465d0, rC=0.272729d0)
55      parameter (rE=5.873645d0, rU=-0.749940d0)
56      parameter (rD=rB*(1.0d0-rU))
57c
58      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
59#ifdef SECOND_DERIV
60      parameter (F73=7.d0/3.d0)
61#endif
62c
63      integer n
64      double precision rrho, rho43, rho13, gamma, gam12, s, d1s(2),
65     &      d, g, gp, d1g(2), facSSB
66#ifdef SECOND_DERIV
67      double precision rhom23, d2s(3), gpp, d2g(3), gssb2
68#endif
69      double precision gssb0,gssb1
70      gssb0(s)= rB*s*s/(1d0+rC*s*s)
71     +               - rD*s*s/(1d0+rE*s**4)
72      gssb1(s)= 2d0*rB*s/(1d0+rC*s*s)**2 +
73     +         (2d0*rD*rE*s**5 - 2d0*rD*s)/(1d0+rE*s**4)**2
74#ifdef SECOND_DERIV
75      gssb2(s)= 8d0*rB/(1d0+rC*s*s)**3 - 6d0*rB/(1d0+rC*s*s)**2
76     +       + 36d0*rD/(1d0+rE*s**4)**2 - 32d0*rD/(1d0+rE*s**4)**3
77     -        - 6d0*rD/(1d0+rE*s**4)
78#endif
79c
80      pi = acos(-1.d0)
81      C = -3d0/(4d0*pi)*(3d0*pi*pi)**F13
82      Cs = 0.5d0/(3d0*pi*pi)**F13
83      Cs = Cs * C               ! account for including C in rho43
84c
85      if (ipol.eq.1 )then
86c
87c        ======> SPIN-RESTRICTED <======
88c
89#ifdef IFCV81
90CDEC$ NOSWP
91#endif
92         do 10 n = 1, nq
93            if (rho(n,1).lt.tol_rho) goto 10
94            rho43 = C*rho(n,1)**F43
95            rrho = 1d0/rho(n,1)
96            rho13 = F43*rho43*rrho
97#ifdef SECOND_DERIV
98            rhom23 = F13*rho13*rrho
99#endif
100c
101            gamma = delrho(n,1,1)*delrho(n,1,1) +
102     &              delrho(n,2,1)*delrho(n,2,1) +
103     &              delrho(n,3,1)*delrho(n,3,1)
104            gam12 = dsqrt(gamma)
105            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 10
106c
107
108            s = Cs*gam12/rho43
109            d1s(1) = -F43*s*rrho
110            d1s(2) = 0.5d0*s/gamma
111c
112c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
113c
114            g=gssb0(s)
115            gp=gssb1(s)
116c
117            d1g(1) = gp*d1s(1)
118            d1g(2) = gp*d1s(2)
119            Ex = Ex + rho43*g*qwght(n)*fac
120            if(ldew)func(n) = func(n) + rho43*g*fac
121            Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac
122            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*rho43*d1g(2)*fac
123#ifdef SECOND_DERIV
124            d2s(1) = -F73*d1s(1)*rrho
125            d2s(2) = -F43*d1s(2)*rrho
126            d2s(3) = -0.5d0*d1s(2)/gamma
127            gpp=gssb2(s)
128            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
129            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
130            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
131            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
132     &           +(rhom23*g
133     &           + 2.d0*rho13*d1g(1)
134     &           + rho43*d2g(1))*fac*2d0
135            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
136     &           +(rho13*d1g(2)
137     &           + rho43*d2g(2))*fac*4d0
138            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
139     &           + rho43*d2g(3)*fac*8d0
140#endif
141 10      continue
142c
143      else
144c
145c        ======> SPIN-UNRESTRICTED <======
146c
147#ifdef IFCV81
148CDEC$ NOSWP
149#endif
150         do 20 n = 1, nq
151            if (rho(n,1).lt.tol_rho) goto 20
152c
153c     Alpha
154c
155            if (rho(n,2).lt.tol_rho) goto 25
156            rho43 = C*(2d0*rho(n,2))**F43
157            rrho = 0.5d0/rho(n,2)
158            rho13 = F43*rho43*rrho
159#ifdef SECOND_DERIV
160            rhom23 = F13*rho13*rrho
161#endif
162            gamma = delrho(n,1,1)*delrho(n,1,1) +
163     &              delrho(n,2,1)*delrho(n,2,1) +
164     &              delrho(n,3,1)*delrho(n,3,1)
165            gam12 = 2d0*dsqrt(gamma)
166            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 25
167c
168            s = Cs*gam12/rho43
169            d1s(1) = -F43*s*rrho
170            d1s(2) = 0.5d0*s/gamma
171c
172c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
173c
174
175            g=gssb0(s)
176            gp=gssb1(s)
177c
178            d1g(1) = gp*d1s(1)
179            d1g(2) = gp*d1s(2)
180            Ex = Ex + rho43*g*qwght(n)*fac*0.5d0
181            if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0
182            Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac
183            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 0.5d0*rho43*d1g(2)*fac
184#ifdef SECOND_DERIV
185            d2s(1) = -F73*d1s(1)*rrho
186            d2s(2) = -F43*d1s(2)*rrho
187            d2s(3) = -0.5d0*d1s(2)/gamma
188            gpp=gssb2(s)
189            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
190            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
191            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
192            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
193     &           +(rhom23*g
194     &           + 2.d0*rho13*d1g(1)
195     &           + rho43*d2g(1))*fac*2d0
196            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
197     &           +(rho13*d1g(2)
198     &           + rho43*d2g(2))*fac
199            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
200     &           + rho43*d2g(3)*fac*0.5d0
201#endif
202c
203c     Beta
204c
205 25         continue
206            if (rho(n,3).lt.tol_rho) goto 20
207            rho43 = C*(2d0*rho(n,3))**F43
208            rrho = 0.5d0/rho(n,3)
209            rho13 = F43*rho43*rrho
210#ifdef SECOND_DERIV
211            rhom23 = F13*rho13*rrho
212#endif
213            gamma = delrho(n,1,2)*delrho(n,1,2) +
214     &              delrho(n,2,2)*delrho(n,2,2) +
215     &              delrho(n,3,2)*delrho(n,3,2)
216            gam12 = 2d0*dsqrt(gamma)
217            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 20
218c
219            s = Cs*gam12/rho43
220            d1s(1) = -F43*s*rrho
221            d1s(2) = 0.5d0*s/gamma
222c
223c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
224c
225            g=gssb0(s)
226            gp=gssb1(s)
227c
228            d1g(1) = gp*d1s(1)
229            d1g(2) = gp*d1s(2)
230            Ex = Ex + rho43*g*qwght(n)*fac*0.5d0
231            if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0
232            Amat(n,2) = Amat(n,2) + (rho13*g+rho43*d1g(1))*fac
233            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 0.5d0*rho43*d1g(2)*fac
234#ifdef SECOND_DERIV
235            d2s(1) = -F73*d1s(1)*rrho
236            d2s(2) = -F43*d1s(2)*rrho
237            d2s(3) = -0.5d0*d1s(2)/gamma
238            gpp=gssb2(s)
239            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
240            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
241            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
242            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
243     &           +(rhom23*g
244     &           + 2.d0*rho13*d1g(1)
245     &           + rho43*d2g(1))*fac*2d0
246            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
247     &           +(rho13*d1g(2)
248     &           + rho43*d2g(2))*fac
249            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB)
250     &           + rho43*d2g(3)*fac*0.5d0
251#endif
252c
253 20      continue
254      endif
255c
256      return
257      end
258#ifndef SECOND_DERIV
259#define SECOND_DERIV
260c
261c     Compile source again for the 2nd derivative case
262c
263#include "xc_ssbD_1.F"
264#endif
265