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