1#if !defined SECOND_DERIV
2      subroutine xc_lb94_e(tol_rho, xfac, rho, delrho, Amat, nq,
3     I     ipol,  Ex, qwght,ldew,func,delrhoq)
4#elif defined(SECOND_DERIV)
5c     For locations of 2nd derivatives of functionals in array
6      subroutine xc_lb94_e_d2(tol_rho, xfac, rho, delrho, Amat,
7     I     Amat2, nq, ipol,  Ex, qwght,ldew,func,delrhoq)
8#include "dft2drv.fh"
9      implicit none
10#endif
11c
12      double precision tol_rho, xfac
13      double precision Ex
14      integer nq, ipol
15      logical ldew
16      double precision func(*),qwght(*)
17      double precision rho(nq,*)
18      double precision delrho(nq,3,*)
19      double precision Amat(nq,*)
20#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
21c
22c     Partial Second Derivatives of the Exchange Energy Functional
23c
24      double precision Amat2(nq,*)
25#endif
26      double precision delrhoq(nq,*) ! [in] 3*rho+ddot(delrho,qxyz)
27c
28      integer n,ii
29      double precision ex_lda   ! ignored
30c
31
32      ex_lda=0d0
33#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
34      call xc_dirac(tol_rho, 1d0, .true.,.false.,rho,
35     &     amat, nq, ipol,ex_lda,qwght,.false.,func)
36
37      call xc_lb94(tol_rho, 0d0, rho, delrho,
38     &     amat, nq, ipol)
39#endif
40#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
41      call xc_dirac_d2(tol_rho, 1d0, .true., .false., rho, Amat,
42     &                       Amat2, nq, ipol, ex_lda, qwght, .false.,
43     &     func)
44
45      call xc_lb94_d2(tol_rho, 0d0, rho, delrho,
46     &     amat, amat2, nq, ipol)
47#endif
48
49c
50c     eq. 12 of DOI:10.1080/00268979600100011
51c     or
52c     eq. 35 of DOI:10.1103/PhysRevA.51.170
53c     delrhoq is computed in xc_delrhodotr
54c
55#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
56      do ii=1,ipol
57         do n=1,nq
58            ex=ex+delrhoq(n,ii)*amat(n,ii)*qwght(n)
59         enddo
60      enddo
61#endif
62      return
63      end
64#if !defined SECOND_DERIV
65      Subroutine xc_lb94(tol_rho, fac, rho, delrho,
66     &                      Amat, nq, ipol)
67#elif defined(SECOND_DERIV)
68      Subroutine xc_lb94_d2(tol_rho, fac, rho, delrho,
69     &     Amat, Amat2,nq, ipol)
70#include "dft2drv.fh"
71#endif
72c
73C$Id$
74c
75      implicit none
76c
77c
78      double precision tol_rho, fac
79      integer nq, ipol
80c
81c     Charge Density
82c
83      double precision rho(nq,*)
84c
85c     Charge Density Gradient
86c
87      double precision delrho(nq,3,*)
88c
89c     Sampling Matrices for the XC Potential
90c
91      double precision Amat(nq,*)
92c
93#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
94c
95c     Partial Second Derivatives of the Exchange Energy Functional
96c
97      double precision Amat2(nq,*)
98#endif
99c
100      double precision BETA
101      Parameter (BETA = 0.05D0)
102c
103c References:
104c
105c    R. van Leeuwen & E. J. Baerends, Phys. Rev. A 49, 2421 (1994).
106c
107c***************************************************************************
108c
109      integer n
110      double precision arcsinh,darcsinh
111      double precision hrho
112      double precision rho13, rho43, gamma, x, g, gdenom
113c      double precisiobn rhom23
114#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
115      double precision dg,dgdenom
116#endif
117c
118      arcsinh(x)=log(x+dsqrt(1d0+x*x))
119      darcsinh(x)=1d0/dsqrt(1d0+x*x)
120c
121      if (ipol.eq.1) then
122c
123c        ======> SPIN-RESTRICTED <======
124c
125         do 10 n = 1, nq
126            if (rho(n,1).lt.tol_rho) goto 10
127c
128c           Spin alpha:
129c
130            hrho  = 0.5d0*rho(n,1)
131            rho13 = hrho**(1.d0/3.d0)
132            rho43 = rho13*hrho
133c            rhom23=rho13/hrho
134            gamma = delrho(n,1,1)*delrho(n,1,1) +
135     &              delrho(n,2,1)*delrho(n,2,1) +
136     &              delrho(n,3,1)*delrho(n,3,1)
137            if (dsqrt(gamma).gt.tol_rho)then
138               gamma = 0.25d0 * gamma
139               x = dsqrt(gamma) / rho43
140            else
141               x = 0d0
142            endif
143c
144            gdenom = 1d0 + 3d0*BETA*x*arcsinh(x)
145            g = -BETA*x*x / gdenom
146c
147            Amat(n,1) = Amat(n,1) + rho13*g*(1d0-fac)
148#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
149            dgdenom = 3d0*BETA*(arcsinh(x) + x*darcsinh(x))
150            dg = BETA*x*(x*dgdenom - 2d0*gdenom) /
151     G           (gdenom*gdenom)
152            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) +
153     R           BETA/(3d0*rho13*rho13)*(4d0*x*dg - g)
154#endif
155c
156 10      continue
157c
158      else
159c
160c        ======> SPIN-UNRESTRICTED <======
161c
162         do 20 n = 1, nq
163            if (rho(n,1).lt.tol_rho) goto 20
164            if (rho(n,2).lt.tol_rho) goto 25
165c
166c           Spin alpha:
167c
168            rho13 = rho(n,2)**(1.d0/3.d0)
169            rho43 = rho13*rho(n,2)
170            gamma = delrho(n,1,1)*delrho(n,1,1) +
171     &              delrho(n,2,1)*delrho(n,2,1) +
172     &              delrho(n,3,1)*delrho(n,3,1)
173            if (dsqrt(gamma).gt.tol_rho)then
174               x = dsqrt(gamma) / rho43
175            else
176               x = 0d0
177            endif
178c
179            gdenom = 1d0 + 3d0*BETA*x*arcsinh(x)
180            g = -BETA*x*x / gdenom
181c
182            Amat(n,1) = Amat(n,1) + rho13*g*(1d0-fac)
183c
184#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
185            dgdenom = 3d0*BETA*(arcsinh(x) + x*darcsinh(x))
186            dg = BETA*x*(x*dgdenom - 2d0*gdenom) /
187     G           (gdenom*gdenom)
188               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) +
189     R              BETA/(3d0*rho13*rho13)*(4d0*x*dg - g)
190#endif
191 25         continue
192c
193c           Spin beta:
194c
195            if (rho(n,3).lt.tol_rho) goto 20
196c
197            rho13 = rho(n,3)**(1.d0/3.d0)
198            rho43 = rho13*rho(n,3)
199            gamma = delrho(n,1,2)*delrho(n,1,2) +
200     &              delrho(n,2,2)*delrho(n,2,2) +
201     &              delrho(n,3,2)*delrho(n,3,2)
202            if (dsqrt(gamma).gt.tol_rho)then
203               x = dsqrt(gamma) / rho43
204            else
205               x = 0d0
206            endif
207c
208            gdenom = 1d0 + 3d0*BETA*x*arcsinh(x)
209            g = -BETA*x*x / gdenom
210c
211            Amat(n,2) = Amat(n,2) + rho13*g*(1d0-fac)
212#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
213            dgdenom = 3d0*BETA*(arcsinh(x) + x*darcsinh(x))
214            dg = BETA*x*(x*dgdenom - 2d0*gdenom) /
215     G           (gdenom*gdenom)
216            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) +
217     R              BETA/(3d0*rho13*rho13)*(4d0*x*dg - g)
218#endif
219c
220 20      continue
221c
222      endif
223c
224      return
225      end
226#ifndef SECOND_DERIV
227#define SECOND_DERIV
228c
229c     Compile source again for the 2nd derivative case
230c
231#include "xc_lb94.F"
232#endif
233