1#ifndef SECOND_DERIV
2      Subroutine xc_op(tol_rho, whichf,
3     &     fac, lfac, nlfac, rho, delrho,
4     &                      Amat, Cmat, nq, ipol, Ec, qwght,ldew,func)
5#else
6      Subroutine xc_op_d2(tol_rho, whichf,
7     &     fac, lfac, nlfac, rho, delrho,
8     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ec,
9     &                         qwght,ldew,func)
10#endif
11c
12C$Id$
13c
14      implicit none
15c
16#include "dft2drv.fh"
17c
18      double precision tol_rho, fac, Ec
19      character*4 whichf
20      integer nq, ipol
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 QABOP,QAB88OP,QABPBOP
48      Parameter (QAB88OP=2.3670D0,QABPBOP=2.3789D0)
49c
50c References:
51c    Tsuneda, Suzumura, Hirao, JCP 110, 10664 (1999)
52c    Tsuneda, Suzumura, Hirao, JCP 111, 5656 (1999)
53c
54c***************************************************************************
55c
56      integer n
57      double precision rho13, rho43, gamma, x
58      double precision kalpha,kbeta, rho13a, rho13b,rhoa,rhob
59      double precision banb, hbab, hbabx
60      double precision dhdab,dhdabx,dkadra,dkbdrb,dkadxa,dkbdxb,
61     A     dbabdra,dbabdrb,dbabdga,dbabdgb,dkadga,dkbdgb,
62     A     dbabdka,dbabdkb
63c
64      hbabx(x) = (1.5214d0*x + 0.5764d0)/
65     /           (x**2*(x**2+1.1284d0*x+0.3183d0))
66      dhdabx(x) = -(4.5642d0*x**4+5.7391d0*x**3+
67     +     2.4355*x**2+0.3669d0*x)/
68     /           ((x**4+1.1284d0*x**3+0.3183d0*x**2)**2)
69c
70      if(whichf.eq.'be88') then
71         QABOP=QAB88OP
72      endif
73      if(whichf.eq.'pb96') then
74         QABOP=QABPBOP
75      endif
76      if (ipol.eq.1) then
77c
78c        ======> SPIN-RESTRICTED <======
79c
80         do 10 n = 1, nq
81            if (rho(n,1).lt.tol_rho) goto 10
82c
83c           Spin alpha:
84c
85            rhoa=rho(n,1)*0.5d0
86            rho13a = (rhoa)**(1.d0/3.d0)
87            rho43 = rho13a**4
88            gamma = delrho(n,1,1)*delrho(n,1,1) +
89     &              delrho(n,2,1)*delrho(n,2,1) +
90     &              delrho(n,3,1)*delrho(n,3,1)
91            gamma = 0.25d0 * gamma
92            if (dsqrt(gamma).gt.tol_rho)then
93               x = dsqrt(gamma) / rho43
94               call xc_kop(tol_rho,whichf,x,
95     &              kalpha, dkadxa)
96               dkadra = -(4d0/3d0)*x*dkadxa/rhoa
97               dkadga = (dkadxa/rho43)*0.5d0/dsqrt(gamma)
98            else
99               x=0d0
100               call xc_kop(tol_rho,whichf,x,
101     &              kalpha, dkadxa)
102               dkadra = 0d0
103               dkadga = 0d0
104            endif
105c
106c
107
108            banb = qabop * rho13a * kalpha *0.5d0
109
110            if(banb.ne.0) then
111               dbabdra = banb*0.5d0*
112     /              (1d0/(3d0*rhoa)+dkadra/kalpha)
113
114               dbabdga = banb/kalpha*dkadga*0.5d0
115
116               hbab = hbabx(banb)
117               dhdab = dhdabx(banb)
118            else
119               dbabdra =0d0
120               dbabdga =0d0
121               hbab = 0d0
122               dhdab = 0d0
123            endif
124
125            Ec = Ec - rhoa**2*hbab*qwght(n)*fac
126            if(ldew)func(n) = func(n) - rhoa**2*hbab*fac
127            Amat(n,1) = Amat(n,1) -
128     -           (rhoa*hbab + rhoa**2*dhdab*dbabdra)*fac
129
130c
131            if (x.gt.tol_rho) then
132                Cmat(n,D1_GAA) = Cmat(n,D1_GAA) -
133     -              rhoa**2*dhdab*dbabdga*fac
134             endif
135c
136 10      continue
137c
138      else
139c
140c        ======> SPIN-UNRESTRICTED <======
141c
142         do 20 n = 1, nq
143            if (rho(n,1).lt.tol_rho) goto 20
144            if (rho(n,2).ge.tol_rho*0.5d0)  then
145c
146c           Spin alpha:
147c
148               rhoa=rho(n,2)
149               rho13a = rhoa**(1.d0/3.d0)
150               rho43 = rho13a*rhoa
151               gamma = delrho(n,1,1)*delrho(n,1,1) +
152     &              delrho(n,2,1)*delrho(n,2,1) +
153     &              delrho(n,3,1)*delrho(n,3,1)
154               if (dsqrt(gamma).gt.tol_rho)then
155                  x = dsqrt(gamma) / rho43
156                  call xc_kop(tol_rho,whichf,x,
157     &                 kalpha, dkadxa)
158
159                  dkadra = -(4d0/3d0)*x*dkadxa/rhoa
160                  dkadga = dkadxa*0.5d0/(rho43*dsqrt(gamma))
161               else
162                  x = 0d0
163               endif
164            else
165               rhoa=0d0
166               rho13a=0d0
167               x = 0d0
168            endif
169            if(x.eq.0d0) then
170               call xc_kop(tol_rho,whichf,x,
171     &              kalpha, dkadxa)
172               dkadra = 0d0
173               dkadga = 0d0
174            endif
175c
176c           Spin beta:
177c
178            if (rho(n,3).ge.tol_rho*0.5d0) then
179c
180               rhob=rho(n,3)
181               rho13b = rhob**(1.d0/3.d0)
182               rho43 = rho13b*rhob
183               gamma = delrho(n,1,2)*delrho(n,1,2) +
184     &              delrho(n,2,2)*delrho(n,2,2) +
185     &              delrho(n,3,2)*delrho(n,3,2)
186               if (dsqrt(gamma).gt.tol_rho)then
187                  x = dsqrt(gamma) / rho43
188                  call xc_kop(tol_rho,whichf,x,
189     &                 kbeta, dkbdxb)
190
191                  dkbdrb = -(4d0/3d0)*x*dkbdxb/rhob
192                  dkbdgb = dkbdxb*0.5d0/(rho43*dsqrt(gamma))
193               else
194                  x = 0d0
195               endif
196            else
197               if(rho13a.eq.0) goto 20
198               rhob=0d0
199               rho13b=0d0
200               x=0d0
201            endif
202            if(x.eq.0d0) then
203               call xc_kop(tol_rho,whichf,x,
204     &              kbeta, dkbdxb)
205               dkbdrb = 0d0
206               dkbdgb=  0d0
207            endif
208
209            banb = qabop*(rho13a*kalpha*rho13b*kbeta)/
210     /           (rho13a*kalpha+rho13b*kbeta)
211
212            if(banb.ne.0) then
213               dbabdra = banb*kbeta*rho13b/
214     /              (rho13a*kalpha+rho13b*kbeta)*
215     /              (1d0/(3d0*rhoa)+dkadra/kalpha)
216               dbabdrb = banb*kalpha*rho13a/
217     /              (rho13a*kalpha+rho13b*kbeta)*
218     /              (1d0/(3d0*rhob)+dkbdrb/kbeta)
219
220               dbabdga = banb*rho13b*kbeta/
221     /              ((rho13a*kalpha+rho13b*kbeta)*kalpha)*
222     *              dkadga
223               dbabdgb = banb*rho13a*kalpha/
224     /              ((rho13a*kalpha+rho13b*kbeta)*kbeta)*
225     *              dkbdgb
226
227               hbab = hbabx(banb)
228               dhdab = dhdabx(banb)
229            else
230               dbabdra =0d0
231               dbabdrb =0d0
232               dbabdga =0d0
233               dbabdgb =0d0
234               hbab = 0d0
235               dhdab = 0d0
236            endif
237
238            Ec = Ec - rhoa*rhob*hbab*qwght(n)*fac
239            if (ldew) func(n) = func(n) - rhoa*rhob*hbab*fac
240            Amat(n,1) = Amat(n,1) -
241     -           (rhob*hbab + rhoa*rhob*dhdab*dbabdra)*fac
242            Amat(n,2) = Amat(n,2) -
243     -           (rhoa*hbab + rhoa*rhob*dhdab*dbabdrb)*fac
244c
245c
246            if (x.gt.tol_rho) then
247               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) -
248     -              rhoa*rhob*dhdab*dbabdga*fac
249               Cmat(n,D1_GBB) = Cmat(n,D1_GBB) -
250     -              rhoa*rhob*dhdab*dbabdgb*fac
251            endif
252
253c
254c
255 20      continue
256c
257      endif
258c
259      return
260      end
261#ifndef SECOND_DERIV
262#define SECOND_DERIV
263c
264c     Compile source again for the 2nd derivative case
265c
266#include "xc_op.F"
267#endif
268