1#if defined(FUJITSU_VPP)
2!ocl scalar
3#endif
4#ifndef SECOND_DERIV
5      Subroutine xc_s12x(tol_rho, fac, lfac, nlfac, rho, delrho,
6     &              Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,is12x)
7#else
8      Subroutine xc_s12x_d2(tol_rho, fac, lfac, nlfac, rho, delrho,
9     &               Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
10     &                         qwght,ldew,func,is12x)
11#endif
12c
13C$Id$
14c
15      implicit none
16c
17#include "dft2drv.fh"
18c
19      double precision tol_rho, fac, Ex
20      integer nq, ipol, is12x
21      logical lfac, nlfac,ldew
22      double precision func(*)  ! value of the functional [output]
23c
24c     Charge Density
25c
26      double precision rho(nq,ipol*(ipol+1)/2)
27c
28c     Charge Density Gradient
29c
30      double precision delrho(nq,3,ipol)
31c
32c     Quadrature Weights
33c
34      double precision qwght(nq)
35c
36c     Sampling Matrices for the XC Potential
37c
38      double precision Amat(nq,ipol), Cmat(nq,*)
39c
40#ifdef SECOND_DERIV
41c
42c     Second Derivatives of the Exchange Energy Functional
43c
44      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
45#endif
46c
47      double precision rB, rC, rD, rA, rK, rE, rH, rG2, rH2
48      double precision ums, vms, dudx, dvdx, d2udx2, d2vdx2
49c
50c References:
51c
52c    Swart, Chem. Phys. Lett. (2013), DOI:10.1016/j.cplett.2013.06.045.
53c
54c***************************************************************************
55c
56      integer n
57      double precision C, rho13, rho43, gamma, x, g, gdenom, dg
58      double precision dgdenom, t, x2, x3, x4, g1, g2
59      double precision g1h1, g2h1, g1h2, g2h2, g1h3, g3h1
60      double precision hdenom, dhdenom, d2hdenom, PI, rM
61      double precision gc4, dgc4, d2gc4
62      parameter (rM=60.770665d0)
63      parameter (PI = 3.1415926535897932385d0)
64#ifdef SECOND_DERIV
65      double precision rhom23, d2g, d2gdenom, d2g1, d2g2, d2g3
66#endif
67c
68
69      if (is12x.eq.1) then
70c
71cswar1      1.03842032     0.00403198     0.00104596     0.75700000     0.00000000
72cswar2      0.00594635     1.17755954     0.84432515     0.00000000     0.00000000
73cswar3      1.00000000     0.00000000     0.00000000     0.00000000     1.00000000     0.00000000
74c
75        rA = 1.03842032d0
76        rK = 0.757d0
77        rC = 0.00403198d0
78        rD = 0.00104596d0
79        rE = 0.00594635d0
80      elseif (is12x.eq.2) then
81c
82cswar1      1.02543951     0.00761554     0.00211063     0.75700000     0.25000000
83cswar2      0.00604672     1.07735222     0.37705816     0.00000000     0.00000000
84cswar3      1.00000000     0.00000000     0.00000000     0.00000000     1.00000000     0.00000000
85c
86        rA = 1.02543951d0
87        rK = 0.757d0
88        rC = 0.00761554d0
89        rD = 0.00211063d0
90        rE = 0.00604672d0
91      else
92        stop 'error in xc_s12x.F'
93      endif
94      rB = 1d0 + rK - rA
95c
96c     Uniform electron gas constant
97c
98      C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0)
99c
100      if (ipol.eq.1) then
101c
102c        ======> SPIN-RESTRICTED <======
103c
104         do 10 n = 1, nq
105            if (rho(n,1).lt.tol_rho) goto 10
106c
107c           Spin alpha:
108c
109            rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
110            rho43 = rho13**4
111            gamma = delrho(n,1,1)*delrho(n,1,1) +
112     &              delrho(n,2,1)*delrho(n,2,1) +
113     &              delrho(n,3,1)*delrho(n,3,1)
114            if (dsqrt(gamma).gt.tol_rho)then
115               gamma = 0.25d0 * gamma
116               x = dsqrt(gamma) / rho43
117               x2 = x*x
118            else
119               x = 0d0
120               x2 = 0d0
121            endif
122c
123            gdenom = 1d0 + rC*x2 + rD*x2*x2
124            hdenom = 1d0 + rE*x2
125            ums = 1d0 - 1d0 / gdenom
126            vms = 1d0 - 1d0 / hdenom
127            g = C*rB*ums*vms
128c
129            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
130            dvdx = 2d0*rE*x/(hdenom**2)
131            dg = C*rB*(dudx*vms + ums*dvdx)
132c
133cmswart            if (lfac) then
134cmswart               Ex = Ex + rA*2d0*rho43*C*qwght(n)*fac
135cmswart               if(ldew)func(n) = func(n) + rA*2.d0*rho43*C*fac
136cmswart               Amat(n,1) = Amat(n,1) + rA*(4d0/3d0)*rho13*C*fac
137cmswart            endif
138c
139            if (nlfac) then
140               Ex = Ex + 2d0*rho43*g*qwght(n)*fac
141               if(ldew)func(n) = func(n) + 2.d0*rho43*g*fac
142               Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac
143            endif
144c
145            if (x.gt.tol_rho) then
146               t = 0.5d0 * dg / sqrt(gamma) * fac
147               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t
148            endif
149c
150#ifdef SECOND_DERIV
151cmswart            if(lfac) g = g + rA*C           ! Add local contribution back to g
152            rhom23 = rho13 / (0.5d0*rho(n,1))
153
154            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
155     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
156            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
157            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
158c
159            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
160     &           + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
161            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
162     &           - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac
163            if (x.gt.tol_rho) then
164               Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
165     &              - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
166            endif
167#endif
168c
169 10      continue
170c
171      else
172c
173c        ======> SPIN-UNRESTRICTED <======
174c
175         do 20 n = 1, nq
176            if (rho(n,1).lt.tol_rho) goto 20
177            if (rho(n,2).lt.tol_rho) goto 25
178c
179c           Spin alpha:
180c
181            rho13 = rho(n,2)**(1.d0/3.d0)
182            rho43 = rho13*rho(n,2)
183            gamma = delrho(n,1,1)*delrho(n,1,1) +
184     &              delrho(n,2,1)*delrho(n,2,1) +
185     &              delrho(n,3,1)*delrho(n,3,1)
186            if (dsqrt(gamma).gt.tol_rho)then
187               x = dsqrt(gamma) / rho43
188               x2 = x*x
189            else
190               x = 0d0
191               x2 = 0d0
192            endif
193c
194            gdenom = 1d0 + rC*x2 + rD*x2*x2
195            hdenom = 1d0 + rE*x2
196            ums = 1d0 - 1d0 / gdenom
197            vms = 1d0 - 1d0 / hdenom
198            g = C*rB*ums*vms
199c
200            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
201            dvdx = 2d0*rE*x/(hdenom**2)
202            dg = C*rB*(dudx*vms + ums*dvdx)
203c
204cmswart            if (lfac) then
205cmswart               Ex = Ex + rA*rho43*C*qwght(n)*fac
206cmswart               if (ldew)func(n) = func(n) + rA*rho43*C*fac
207cmswart               Amat(n,1) = Amat(n,1) + rA*(4d0/3d0)*rho13*C*fac
208cmswart            endif
209c
210            if (nlfac) then
211               Ex = Ex + rho43*g*qwght(n)*fac
212               if (ldew)func(n) = func(n) + rho43*g*fac
213               Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac
214            endif
215c
216            if (x.gt.tol_rho) then
217               t = dg / sqrt(gamma) * fac
218               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t * 0.5d0
219            endif
220c
221#ifdef SECOND_DERIV
222cmswart            if (lfac) g = g + rA*C           ! Add local contribution back to g
223            rhom23 = rho13 / rho(n,2)
224
225            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
226     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
227            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
228            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
229c
230            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
231     &           + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
232            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
233     &           - (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac
234            if (x.gt.tol_rho) then
235               Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
236     &              - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
237            endif
238#endif
239c
240 25         continue
241c
242c           Spin beta:
243c
244            if (rho(n,3).lt.tol_rho) goto 20
245c
246            rho13 = rho(n,3)**(1.d0/3.d0)
247            rho43 = rho13*rho(n,3)
248            gamma = delrho(n,1,2)*delrho(n,1,2) +
249     &              delrho(n,2,2)*delrho(n,2,2) +
250     &              delrho(n,3,2)*delrho(n,3,2)
251            if (dsqrt(gamma).gt.tol_rho)then
252               x = dsqrt(gamma) / rho43
253               x2 = x*x
254            else
255               x = 0d0
256               x2 = 0d0
257            endif
258c
259            gdenom = 1d0 + rC*x2 + rD*x2*x2
260            hdenom = 1d0 + rE*x2
261            ums = 1d0 - 1d0 / gdenom
262            vms = 1d0 - 1d0 / hdenom
263            g = C*rB*ums*vms
264c
265            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
266            dvdx = 2d0*rE*x/(hdenom**2)
267            dg = C*rB*(dudx*vms + ums*dvdx)
268c
269cmswart            if (lfac) then
270cmswart               Ex = Ex + rA*rho43*C*qwght(n)*fac
271cmswart               if (ldew)func(n) = func(n) + rA*rho43*C*fac
272cmswart               Amat(n,2) = Amat(n,2) + rA*(4d0/3d0)*rho13*C*fac
273cmswart            endif
274c
275            if (nlfac) then
276               Ex = Ex + rho43*g*qwght(n)*fac
277               if (ldew)func(n) = func(n) +rho43*g*fac
278               Amat(n,2) = Amat(n,2) + (4d0/3d0)*rho13*(g-x*dg)*fac
279            endif
280c
281            if (x.gt.tol_rho) then
282               t = dg / sqrt(gamma) * fac
283               Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + t * 0.5d0
284            endif
285c
286#ifdef SECOND_DERIV
287cmswart            if(lfac) g = g + rA*C           ! Add local contribution back to g
288            rhom23 = rho13 / rho(n,3)
289
290            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
291     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
292            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
293            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
294
295c
296            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
297     &           + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
298            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
299     &           - (2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac
300            if (x.gt.tol_rho) then
301               Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB)
302     &              - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
303            endif
304#endif
305c
306 20      continue
307c
308      endif
309c
310      return
311      end
312#ifndef SECOND_DERIV
313#define SECOND_DERIV
314c
315c     Compile source again for the 2nd derivative case
316c
317#include "xc_s12x.F"
318#endif
319