1#if defined(FUJITSU_VPP)
2!ocl scalar
3#endif
4#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
5      Subroutine xc_camb88(tol_rho, fac, lfac, nlfac, rho, delrho,
6     &                      Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
7#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
8      Subroutine xc_camb88_d2(tol_rho, fac, lfac, nlfac, rho, delrho,
9     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
10     &                         qwght,ldew,func)
11#else
12      Subroutine xc_camb88_d3(tol_rho, fac, lfac, nlfac, rho, delrho,
13     1                         Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3,
14     2                         nq, ipol, Ex, qwght,ldew,func)
15#endif
16c
17C$Id$
18c
19c     Coulomb attenuated Becke88 functional
20c
21      implicit none
22c
23#include "dft2drv.fh"
24#include "dft3drv.fh"
25c
26      double precision tol_rho, fac, Ex
27      integer nq, ipol
28      logical lfac, nlfac,ldew
29      double precision func(*)  ! value of the functional [output]
30c
31c     Charge Density
32c
33      double precision rho(nq,ipol*(ipol+1)/2)
34c
35c     Charge Density Gradient
36c
37      double precision delrho(nq,3,ipol)
38c
39c     Quadrature Weights
40c
41      double precision qwght(nq)
42c
43c     Sampling Matrices for the XC Potential
44c
45      double precision Amat(nq,ipol), Cmat(nq,*)
46c
47#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
48c
49c     Second Derivatives of the Exchange Energy Functional
50c
51      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
52      double precision A2tmp, C2tmp, C3tmp
53#endif
54#ifdef THIRD_DERIV
55c
56c     Third Derivatives of the Exchange Energy Functional
57c
58      double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3)
59      double precision A3tmp, C4tmp, C5tmp, C6tmp
60#endif
61c
62      double precision BETA
63      Parameter (BETA = 0.0042D0)
64c
65c References:
66c
67c    Becke, Phys. Rev. A 38, 3098 (1988)
68c    Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993)
69c
70c***************************************************************************
71c
72      integer n
73      double precision arcsinh, darcsinh, d2arcsinh
74      double precision C, rho13, rho43, gamma, x, g, gdenom, dg,
75     &     dgdenom, t, Etmp, Atmp, Ctmp
76      double precision gdenom2
77#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
78      double precision rhom23, d2g, d2gdenom
79      double precision gdenom3
80#endif
81c
82#ifdef THIRD_DERIV
83      double precision rhom53, d3g, d3gdenom
84      double precision gdenom4
85#endif
86c
87      arcsinh(x)=log(x+dsqrt(1d0+x*x))
88      darcsinh(x)=1d0/dsqrt(1d0+x*x)
89      d2arcsinh(x) = -x/dsqrt(1d0+x*x)**3
90c
91c     Uniform electron gas constant
92c
93      C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0)
94c
95      if (ipol.eq.1) then
96c
97c        ======> SPIN-RESTRICTED <======
98c
99         do 10 n = 1, nq
100            if (rho(n,1).lt.tol_rho) goto 10
101c
102c           Spin alpha:
103c
104            rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
105            rho43 = rho13**4
106            gamma = delrho(n,1,1)*delrho(n,1,1) +
107     &              delrho(n,2,1)*delrho(n,2,1) +
108     &              delrho(n,3,1)*delrho(n,3,1)
109            if (dsqrt(gamma).gt.tol_rho)then
110               gamma = 0.25d0 * gamma
111               x = dsqrt(gamma) / rho43
112            else
113               x = 0d0
114            endif
115c
116            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
117            gdenom2 = gdenom*gdenom
118            g = -BETA*x*x / gdenom
119            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
120c            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
121            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2
122c
123            Etmp = 0.d0
124            Atmp = 0.d0
125            Ctmp = 0.d0
126            if (lfac) then
127               Etmp = 2d0*rho43*C
128               Atmp = (4d0/3d0)*rho13*C
129            endif
130c
131            if (nlfac) then
132               Etmp = Etmp + 2d0*rho43*g
133               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)
134            endif
135c
136            if (x.gt.tol_rho) then
137               Ctmp = 0.5d0 * dg / sqrt(gamma)
138            endif
139c
140#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
141            A2tmp = 0d0
142            C2tmp = 0d0
143            C3tmp = 0d0
144            if(lfac) g = g + C           ! Add local contribution back to g
145            rhom23 = rho13 / (0.5d0*rho(n,1))
146            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
147            gdenom3 = gdenom2*gdenom
148            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2
149     &           + BETA*x*x*d2gdenom/gdenom2
150     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom3
151c
152c rr
153            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)
154c rg
155            C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g
156c gg
157            if (x.gt.tol_rho) then
158               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)
159            endif
160#endif
161#ifdef THIRD_DERIV
162            A3tmp = 0.0d0
163            C4tmp = 0.0d0
164            C5tmp = 0.0d0
165            C6tmp = 0.0d0
166c
167            rhom53 = rhom23 / (0.5d0*rho(n,1))
168c
169            d3gdenom = 6.0d0*BETA*
170     1          d2arcsinh(x)*( 3.0d0
171     2                       - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x))
172c
173            gdenom4 = gdenom3*gdenom
174c
175            d3g = 6.0d0*BETA*dgdenom/gdenom2
176     1          - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3
177     2          + 6.0d0*BETA*x*d2gdenom/gdenom2
178     3          + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4
179     4          - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3
180     5          + BETA*x*x*d3gdenom/gdenom2
181c
182c rrr
183            A3tmp = (8.0d0/27.0d0)*rhom53*(-g + x*dg
184     1                                    - 18.0d0*x*x*d2g
185     2                                    - 8.0d0*x*x*x*d3g)
186c
187c rrg
188            C4tmp = (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g
189     1                                           + 4.0d0*x*x*x*d3g)
190c
191c rgg
192            C5tmp = -(8.0d0/3.0d0)*(rhom23/rho(n,1)**3)/dsqrt(gamma)
193     1               *d3g
194c
195c ggg
196            if (x.gt.tol_rho) then
197              C6tmp = (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg
198     1                                              - 3.0d0*x*d2g
199     2                                              + x*x*d3g)
200            endif
201#endif
202c
203#ifdef THIRD_DERIV
204            call xc_att_xc_d3(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
205     &           C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp)
206c
207            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac
208            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac
209            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac
210c
211            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp*fac
212            Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + C4tmp*fac
213            Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) + C5tmp*
214     *           fac
215            Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + C6tmp*
216     *           fac
217#elif defined(SECOND_DERIV)
218            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
219     &           C2tmp,C3tmp)
220c
221            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac
222            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac
223            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac
224#else
225            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
226#endif
227            Ex = Ex + qwght(n)*Etmp*fac
228            if(ldew) func(n) = func(n) + Etmp*fac
229            Amat(n,1) = Amat(n,1) + Atmp*fac
230            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp*fac
231 10      continue
232c
233      else
234c
235c        ======> SPIN-UNRESTRICTED <======
236c
237         do 20 n = 1, nq
238            if (rho(n,1).lt.tol_rho) goto 20
239            if (rho(n,2).lt.tol_rho) goto 25
240c
241c           Spin alpha:
242c
243            rho13 = rho(n,2)**(1.d0/3.d0)
244            rho43 = rho13*rho(n,2)
245            gamma = delrho(n,1,1)*delrho(n,1,1) +
246     &              delrho(n,2,1)*delrho(n,2,1) +
247     &              delrho(n,3,1)*delrho(n,3,1)
248            if (dsqrt(gamma).gt.tol_rho)then
249               x = dsqrt(gamma) / rho43
250            else
251               x = 0d0
252            endif
253c
254            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
255            g = -BETA*x*x / gdenom
256            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
257            gdenom2 = gdenom*gdenom
258            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2
259c
260            Etmp = 0.d0
261            Atmp = 0.d0
262            Ctmp = 0.d0
263            if (lfac) then
264               Etmp = rho43*C
265               Atmp = (4d0/3d0)*rho13*C
266            endif
267c
268            if (nlfac) then
269               Etmp = Etmp + rho43*g
270               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)
271            endif
272c
273            if (x.gt.tol_rho) then
274               Ctmp = 0.5d0*dg / sqrt(gamma)
275            endif
276c
277#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
278            if (lfac) g = g + C           ! Add local contribution back to g
279            rhom23 = rho13 / rho(n,2)
280            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
281            gdenom3 = gdenom2*gdenom
282            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2
283     &           + BETA*x*x*d2gdenom/gdenom2
284     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom3
285c            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
286c     &           + BETA*x*x*d2gdenom/gdenom**2
287c     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
288c
289            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)
290c            C2tmp =  (2d0/3d0)*(rhom23**2/rho(n,2))*d2g
291            C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,2))*d2g
292            if (x.gt.tol_rho) then
293               C3tmp = -0.25d0*gamma**(-1.5d0)*(dg-x*d2g)
294            endif
295#endif
296#ifdef THIRD_DERIV
297            rhom53 = rhom23 / rho(n,2)
298c
299            d3gdenom = 6.0d0*BETA*
300     1          d2arcsinh(x)*( 3.0d0
301     2                       - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x))
302c
303            gdenom4 = gdenom3*gdenom
304c
305            d3g = 6.0d0*BETA*dgdenom/gdenom2
306     1          - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3
307     2          + 6.0d0*BETA*x*d2gdenom/gdenom2
308     3          + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4
309     4          - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3
310     5          + BETA*x*x*d3gdenom/gdenom2
311c
312            A3tmp = (8.0d0/27.0d0)*rhom53*(-g + x*dg
313     1                                    - 18.0d0*x*x*d2g
314     2                                    - 8.0d0*x*x*x*d3g)
315c
316            C4tmp = (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g
317     1                                           + 4.0d0*x*x*x*d3g)
318c
319            C5tmp = -(1.0d0/3.0d0)*(rhom23/rho(n,2)**3)/dsqrt(gamma)
320     1               *d3g
321c
322            if (x.gt.tol_rho) then
323              C6tmp = (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg
324     1                                              - 3.0d0*x*d2g
325     2                                              + x*x*d3g)
326            endif
327#endif
328c
329#ifdef THIRD_DERIV
330            call xc_att_xc_d3(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
331     &           C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp)
332c
333            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac
334            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac
335            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac
336c
337            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp*fac
338            Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + C4tmp*fac
339            Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) + C5tmp*
340     *           fac
341            Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + C6tmp*
342     *           fac
343#elif defined(SECOND_DERIV)
344            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp,
345     &           C3tmp)
346            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp*fac
347            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp*fac
348            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp*fac
349#else
350           call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
351#endif
352            Ex = Ex + qwght(n)*Etmp*fac
353            if(ldew) func(n) = func(n) + Etmp*fac
354            Amat(n,1) = Amat(n,1) + Atmp*fac
355            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp*fac
356c
357 25         continue
358c
359c           Spin beta:
360c
361            if (rho(n,3).lt.tol_rho) goto 20
362c
363            rho13 = rho(n,3)**(1.d0/3.d0)
364            rho43 = rho13*rho(n,3)
365            gamma = delrho(n,1,2)*delrho(n,1,2) +
366     &              delrho(n,2,2)*delrho(n,2,2) +
367     &              delrho(n,3,2)*delrho(n,3,2)
368            if (dsqrt(gamma).gt.tol_rho)then
369               x = dsqrt(gamma) / rho43
370            else
371               x = 0d0
372            endif
373c
374            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
375            g = -BETA*x*x / gdenom
376            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
377            gdenom2 = gdenom*gdenom
378            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2
379c
380            Etmp = 0.d0
381            Atmp = 0.d0
382            Ctmp = 0.d0
383            if (lfac) then
384               Etmp = rho43*C
385               Atmp = (4d0/3d0)*rho13*C
386            endif
387c
388            if (nlfac) then
389               Etmp = Etmp + rho43*g
390               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)
391            endif
392c
393            if (x.gt.tol_rho) then
394               Ctmp = 0.5d0*dg / sqrt(gamma)
395            endif
396c
397#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
398            A2tmp = 0d0
399            C2tmp = 0d0
400            C3tmp = 0d0
401            if(lfac) g = g + C           ! Add local contribution back to g
402            rhom23 = rho13 / rho(n,3)
403            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
404            gdenom3 = gdenom2*gdenom
405            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2
406     &           + BETA*x*x*d2gdenom/gdenom2
407     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom3
408            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)
409            C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g
410            if (x.gt.tol_rho) then
411               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)
412            endif
413#endif
414#ifdef THIRD_DERIV
415            rhom53 = rhom23 / rho(n,3)
416c
417            d3gdenom = 6.0d0*BETA*
418     1          d2arcsinh(x)*( 3.0d0
419     2                       - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x))
420c
421            gdenom4 = gdenom3*gdenom
422c
423            d3g = 6.0d0*BETA*dgdenom/gdenom2
424     1          - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3
425     2          + 6.0d0*BETA*x*d2gdenom/gdenom2
426     3          + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4
427     4          - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3
428     5          + BETA*x*x*d3gdenom/gdenom2
429c
430            A3tmp = (8.0d0/27.0d0)*rhom53*(-g + x*dg
431     1                                    - 18.0d0*x*x*d2g
432     2                                    - 8.0d0*x*x*x*d3g)
433c
434            C4tmp = (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g
435     1                                           + 4.0d0*x*x*x*d3g)
436c
437            C5tmp = -(1.0d0/3.0d0)*(rhom23/rho(n,3)**3)/dsqrt(gamma)
438     1               *d3g
439c
440            if (x.gt.tol_rho) then
441              C6tmp = (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg
442     1                                              - 3.0d0*x*d2g
443     2                                              + x*x*d3g)
444            endif
445#endif
446c
447#ifdef THIRD_DERIV
448            call xc_att_xc_d3(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
449     &           C2tmp,C3tmp,A3tmp,C4tmp,C5tmp,C6tmp)
450c
451            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp*fac
452            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp*fac
453            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp*fac
454c
455            Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + A3tmp*fac
456            Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) + C4tmp*fac
457            Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) + C5tmp*
458     *           fac
459            Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) + C6tmp*
460     *           fac
461#elif defined(SECOND_DERIV)
462            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp,
463     &           C3tmp)
464            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp*fac
465            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp*fac
466            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp*fac
467#else
468            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
469#endif
470            Ex = Ex + qwght(n)*Etmp*fac
471            if(ldew) func(n) = func(n) + Etmp*fac
472            Amat(n,2) = Amat(n,2) + Atmp*fac
473            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp*fac
474 20      continue
475c
476      endif
477c
478      return
479      end
480#ifndef SECOND_DERIV
481#define SECOND_DERIV
482c
483c     Compile source again for the 2nd derivative case
484c
485#include "xc_camb88.F"
486#endif
487c
488#ifndef THIRD_DERIV
489#define THIRD_DERIV
490c
491c     Compile source again for the 3rd derivative case
492c
493#include "xc_camb88.F"
494#endif
495