1#if defined(FUJITSU_VPP)
2!ocl scalar
3#endif
4#ifndef SECOND_DERIV
5      Subroutine xc_cams12x(tol_rho, fac, lfac, nlfac, rho, delrho,
6     &              Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,is12x)
7#else
8      Subroutine xc_cams12x_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,*)
39      double precision Atmp, Ctmp, Etmp
40      double precision A2tmp, C2tmp, C3tmp
41c
42#ifdef SECOND_DERIV
43c
44c     Second Derivatives of the Exchange Energy Functional
45c
46      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
47#endif
48c
49      double precision rB, rC, rD, rA, rK, rE, rH, rG2, rH2
50      double precision ums, vms, dudx, dvdx, d2udx2, d2vdx2
51c
52c References:
53c
54c    Swart, Chem. Phys. Lett. (2013) DOI:10.1016/j.cplett.2013.06.045.
55c
56c***************************************************************************
57c
58      integer n
59      double precision C, rho13, rho43, gamma, x, g, gdenom, dg
60      double precision dgdenom, t, x2, x3, x4, g1, g2
61      double precision g1h1, g2h1, g1h2, g2h2, g1h3, g3h1
62      double precision hdenom, dhdenom, d2hdenom, PI, rM
63      double precision gc4, dgc4, d2gc4
64      parameter (rM=60.770665d0)
65      parameter (PI = 3.1415926535897932385d0)
66#ifdef SECOND_DERIV
67      double precision rhom23, d2g, d2gdenom, d2g1, d2g2, d2g3
68#endif
69c
70      if (is12x.eq.1) then
71c
72cswar1      1.03323556     0.00417251     0.00115216     0.75700000     0.00000000
73cswar2      0.00706184     1.20250451     0.86124355     0.00000000     0.34485046
74cswar3      1.00000000     0.00000000     0.00000000     0.00000000     1.00000000     1.52420731
75c
76        rA = 1.03323556d0
77        rK = 0.757d0
78        rC = 0.00417251d0
79        rD = 0.00115216d0
80        rE = 0.00706184d0
81      elseif (is12x.eq.2) then
82c
83cswar1      1.02149642     0.00825905     0.00235804     0.75700000     0.25000000
84cswar2      0.00654977     1.08034183     0.37999939     0.00000000     0.35897845
85cswar3      1.00000000     0.00000000     0.00000000     0.00000000     1.00000000     0.48516891
86c
87        rA = 1.02149642d0
88        rK = 0.757d0
89        rC = 0.00825905d0
90        rD = 0.00235804d0
91        rE = 0.00654977d0
92      else
93        stop 'error in xc_cams12x.F'
94      endif
95      rB = 1d0 + rK - rA
96c
97c     Uniform electron gas constant
98c
99      C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0)
100c
101      if (ipol.eq.1) then
102c
103c        ======> SPIN-RESTRICTED <======
104c
105         do 10 n = 1, nq
106            if (rho(n,1).lt.tol_rho) goto 10
107c
108c           Spin alpha:
109c
110            rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
111            rho43 = rho13**4
112c
113            Etmp = 0.d0
114            Atmp = 0.d0
115            Ctmp = 0.d0
116            if (lfac) then
117               Etmp = rA * 2d0*rho43*C*fac
118               Atmp = rA * (4d0/3d0)*rho13*C*fac
119            endif
120c
121            gamma = delrho(n,1,1)*delrho(n,1,1) +
122     &              delrho(n,2,1)*delrho(n,2,1) +
123     &              delrho(n,3,1)*delrho(n,3,1)
124            if (dsqrt(gamma).gt.tol_rho)then
125               gamma = 0.25d0 * gamma
126               x = dsqrt(gamma) / rho43
127               x2 = x*x
128            else
129               x = 0d0
130               x2 = 0d0
131            endif
132c
133            gdenom = 1d0 + rC*x2 + rD*x2*x2
134            hdenom = 1d0 + rE*x2
135            ums = 1d0 - 1d0 / gdenom
136            vms = 1d0 - 1d0 / hdenom
137            g = C*rB*ums*vms
138c
139            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
140            dvdx = 2d0*rE*x/(hdenom**2)
141            dg = C*rB*(dudx*vms + ums*dvdx)
142c
143            if (nlfac) then
144               Etmp = Etmp + 2d0*rho43*g*fac
145               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
146            endif
147c
148            if (x.gt.tol_rho) then
149               Ctmp = 0.5d0 * dg / sqrt(gamma) * fac
150            endif
151c
152#ifdef SECOND_DERIV
153c
154c           Add local contribution back to g
155c
156            if(lfac) g = g + rA * C
157c
158            rhom23 = rho13 / (0.5d0*rho(n,1))
159
160            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
161     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
162            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
163            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
164c
165            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
166            C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac
167            if (x.gt.tol_rho) then
168               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
169            else
170               C3tmp = 0d0
171            endif
172c
173            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
174     &           C2tmp,C3tmp)
175            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
176            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
177            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
178#else
179            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
180#endif
181            Ex = Ex + qwght(n)*Etmp
182            if (ldew) func(n) = func(n) + Etmp
183            Amat(n,1) = Amat(n,1) + Atmp
184            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
185c
186 10      continue
187c
188      else
189c
190c        ======> SPIN-UNRESTRICTED <======
191c
192         do 20 n = 1, nq
193            if (rho(n,1).lt.tol_rho) goto 20
194            if (rho(n,2).lt.tol_rho) goto 25
195c
196c           Spin alpha:
197c
198            rho13 = rho(n,2)**(1.d0/3.d0)
199            rho43 = rho13*rho(n,2)
200c
201            Etmp = 0.d0
202            Atmp = 0.d0
203            Ctmp = 0.d0
204            if (lfac) then
205               Etmp = rA * rho43*C*fac
206               Atmp = rA * (4d0/3d0)*rho13*C*fac
207            endif
208c
209            gamma = delrho(n,1,1)*delrho(n,1,1) +
210     &              delrho(n,2,1)*delrho(n,2,1) +
211     &              delrho(n,3,1)*delrho(n,3,1)
212            if (dsqrt(gamma).gt.tol_rho)then
213               x = dsqrt(gamma) / rho43
214               x2 = x*x
215            else
216               x = 0d0
217               x2 = 0d0
218            endif
219c
220            gdenom = 1d0 + rC*x2 + rD*x2*x2
221            hdenom = 1d0 + rE*x2
222            ums = 1d0 - 1d0 / gdenom
223            vms = 1d0 - 1d0 / hdenom
224            g = C*rB*ums*vms
225c
226            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
227            dvdx = 2d0*rE*x/(hdenom**2)
228            dg = C*rB*(dudx*vms + ums*dvdx)
229c
230            if (nlfac) then
231               Etmp = Etmp + rho43*g*fac
232               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
233            endif
234c
235            if (x.gt.tol_rho) then
236               t = dg / sqrt(gamma) * fac
237               Ctmp = t * 0.5d0
238            endif
239c
240#ifdef SECOND_DERIV
241c
242c           Add local contribution back to g
243c
244            if (lfac) g = g + rA * C
245c
246            rhom23 = rho13 / rho(n,2)
247c
248            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
249     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
250            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
251            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
252c
253            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
254            C2tmp = - (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac
255            if (x.gt.tol_rho) then
256               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
257            else
258               C3tmp = 0d0
259            endif
260
261            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
262     &           C2tmp,C3tmp)
263            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
264            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
265            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
266#else
267            call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
268#endif
269            Ex = Ex + qwght(n)*Etmp
270            if (ldew) func(n) = func(n) + Etmp
271            Amat(n,1) = Amat(n,1) + Atmp
272            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
273c
274 25         continue
275c
276c           Spin beta:
277c
278            if (rho(n,3).lt.tol_rho) goto 20
279c
280            rho13 = rho(n,3)**(1.d0/3.d0)
281            rho43 = rho13*rho(n,3)
282c
283            Etmp = 0.d0
284            Atmp = 0.d0
285            Ctmp = 0.d0
286            if (lfac) then
287               Etmp = rA * rho43*C*fac
288               Atmp = rA * (4d0/3d0)*rho13*C*fac
289            endif
290c
291            gamma = delrho(n,1,2)*delrho(n,1,2) +
292     &              delrho(n,2,2)*delrho(n,2,2) +
293     &              delrho(n,3,2)*delrho(n,3,2)
294            if (dsqrt(gamma).gt.tol_rho)then
295               x = dsqrt(gamma) / rho43
296               x2 = x*x
297            else
298               x = 0d0
299               x2 = 0d0
300            endif
301c
302            gdenom = 1d0 + rC*x2 + rD*x2*x2
303            hdenom = 1d0 + rE*x2
304            ums = 1d0 - 1d0 / gdenom
305            vms = 1d0 - 1d0 / hdenom
306            g = C*rB*ums*vms
307c
308            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
309            dvdx = 2d0*rE*x/(hdenom**2)
310            dg = C*rB*(dudx*vms + ums*dvdx)
311c
312            if (nlfac) then
313               Etmp = Etmp + rho43*g*fac
314               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
315            endif
316c
317            if (x.gt.tol_rho) then
318               t = dg / sqrt(gamma) * fac
319               Ctmp = t * 0.5d0
320            endif
321c
322#ifdef SECOND_DERIV
323c
324c           Add local contribution back to g
325c
326            if (lfac) g = g + rA * C
327c
328            rhom23 = rho13 / rho(n,3)
329c
330            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
331     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
332            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
333            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
334
335c
336            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
337            C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac
338            if (x.gt.tol_rho) then
339               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
340            else
341               C3tmp = 0d0
342            endif
343c
344            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
345     &           C2tmp,C3tmp)
346            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
347            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp
348            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp
349#else
350            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
351#endif
352            Ex = Ex + qwght(n)*Etmp
353            if (ldew) func(n) = func(n) + Etmp
354            Amat(n,2) = Amat(n,2) + Atmp
355            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp
356c
357 20      continue
358c
359      endif
360c
361      return
362      end
363#ifndef SECOND_DERIV
364#define SECOND_DERIV
365c
366c     Compile source again for the 2nd derivative case
367c
368#include "xc_cams12x.F"
369#endif
370