1#ifndef SECOND_DERIV
2      Subroutine xc_xpw6(tol_rho, fac, lfac, nlfac, rho, delrho,
3     &                  Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,ijzy)
4#else
5      Subroutine xc_xpw6_d2(tol_rho, fac, lfac, nlfac, rho, delrho,
6     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
7     &                         qwght,ldew,func,ijzy)
8#endif
9c
10C$Id$
11c
12      implicit none
13c
14#include "dft2drv.fh"
15c
16      double precision fac, Ex
17      integer nq, ipol,ijzy
18      logical lfac, nlfac,ldew
19      double precision func(*)  ! value of the functional [output]
20c
21c     Charge Density & Its Cube Root
22c
23      double precision rho(nq,ipol*(ipol+1)/2)
24c
25c     Charge Density Gradient
26c
27      double precision delrho(nq,3,ipol)
28c
29c     Quadrature Weights
30c
31      double precision qwght(nq)
32c
33c     Sampling Matrices for the XC Potential & Energy
34c
35      double precision Amat(nq,ipol), Cmat(nq,*)
36c
37c
38c     Compute the partial derivatives of the exchange functional of Perdew91.
39c
40c     Becke & Perdew  Parameters
41c
42      double precision DPOW
43      double precision BETA,  tol_rho, CPW91
44
45#ifdef SECOND_DERIV
46c
47c     Second Derivatives of the Exchange Energy Functional
48c
49      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
50#endif
51c
52c References:
53c
54c
55c***************************************************************************
56
57      integer n
58      double precision gamma
59c
60      if (ijzy.eq.2) then
61c
62c       parameters for PWB6K
63c
64        BETA  = 0.00539d0
65        CPW91 = 1.7077d0
66        DPOW  = 4.0876d0
67      else
68c
69c       parameters for PW6B95
70c
71        BETA  = 0.00538d0
72        CPW91 = 1.7382d0
73        DPOW  = 3.8901d0
74      endif
75      if (ipol.eq.1 )then
76c
77c        ======> SPIN-RESTRICTED <======
78c
79         do 10 n = 1, nq
80            if (rho(n,1).lt.tol_rho) goto 10
81            gamma = delrho(n,1,1)*delrho(n,1,1) +
82     &           delrho(n,2,1)*delrho(n,2,1) +
83     &           delrho(n,3,1)*delrho(n,3,1)
84            gamma=0.25d0*gamma
85#ifndef SECOND_DERIV
86            call xc_xpw6core(DPOW,BETA,CPW91,n,1,
87     &           rho(n,1),gamma,qwght(n),func(n),
88     &           tol_rho, fac, lfac, nlfac,
89     &           Amat, Cmat, nq, ipol, Ex, ldew)
90#else
91               call xc_xpw6core_d2(DPOW,BETA,CPW91,n,1,
92     &           rho(n,1),gamma,qwght(n),func(n),
93     &           tol_rho, fac, lfac, nlfac,
94     &           Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
95     &           ldew)
96#endif
97c         write (*,*) "Ex",Ex,Amat(n,1),Cmat(n,1)
98c         stop
99   10    continue
100c
101      else
102c
103c        ======> SPIN-UNRESTRICTED <======
104c
105         do 20 n = 1, nq
106           if (rho(n,1).lt.tol_rho) goto 20
107c
108c           Spin alpha:
109c
110           if (rho(n,2).gt.tol_rho) then
111              gamma =    delrho(n,1,1)*delrho(n,1,1) +
112     &             delrho(n,2,1)*delrho(n,2,1) +
113     &             delrho(n,3,1)*delrho(n,3,1)
114#ifndef SECOND_DERIV
115               call xc_xpw6core(DPOW,BETA,CPW91,n,1,
116     &           rho(n,2),gamma,qwght(n),func(n),
117     &              tol_rho, fac, lfac, nlfac,
118     &              Amat, Cmat, nq, ipol, Ex, ldew)
119#else
120               call xc_xpw6core_d2(DPOW,BETA,CPW91,n,1,
121     &           rho(n,2),gamma,qwght(n),func(n),
122     &              tol_rho, fac, lfac, nlfac,
123     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
124     &                         ldew)
125#endif
126            endif
127c
128c           Spin beta:
129c
130            if (rho(n,3).gt.tol_rho) then
131
132            gamma =   delrho(n,1,2)*delrho(n,1,2) +
133     &           delrho(n,2,2)*delrho(n,2,2) +
134     &           delrho(n,3,2)*delrho(n,3,2)
135#ifndef SECOND_DERIV
136            call xc_xpw6core(DPOW,BETA,CPW91,n,2,
137     &           rho(n,3),gamma,qwght(n),func(n),
138     &           tol_rho, fac, lfac, nlfac,
139     &           Amat, Cmat, nq, ipol, Ex, ldew)
140#else
141               call xc_xpw6core_d2(DPOW,BETA,CPW91,n,2,
142     &           rho(n,3),gamma,qwght(n),func(n),
143     &           tol_rho, fac, lfac, nlfac,
144     &           Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
145     &           ldew)
146#endif
147            endif
148c
149   20    continue
150c
151      endif
152c
153      return
154      end
155
156
157
158
159#ifndef SECOND_DERIV
160      Subroutine xc_xpw6core(DPOW,BETA,CPW91,n,ispin,
161     .     rho,gamma,qwght,func,
162     &     tol_rho, fac, lfac, nlfac,
163     &                      Amat, Cmat, nq, ipol, Ex, ldew)
164#else
165      Subroutine xc_xpw6core_d2(DPOW,BETA,CPW91,n,ispin,
166     .     rho,gamma,qwght,func,
167     &     tol_rho, fac, lfac, nlfac,
168     &     Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
169     &     ldew)
170#endif
171c
172C$Id$
173c
174      implicit none
175c
176#include "dft2drv.fh"
177c
178      double precision fac, Ex
179      integer nq, ipol
180      logical lfac, nlfac,ldew
181      double precision func  ! value of the functional [output]
182      double precision rho
183      double precision qwght
184      integer ispin ! alpha=1; beta=2
185c
186c     Sampling Matrices for the XC Potential & Energy
187c
188      double precision Amat(nq,ipol), Cmat(nq,*)
189c
190c
191c     Compute the partial derivatives of the exchange functional of Perdew91.
192c
193c     Becke & Perdew  Parameters
194c
195      double precision DPOW
196      double precision BETA, tol_rho, AX, pi,
197     &     CPW91,BETAPW91,big
198
199      Parameter (big=1d4)
200
201#ifdef SECOND_DERIV
202c
203c     Second Derivatives of the Exchange Energy Functional
204c
205      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
206      double precision rhom23, f2x, d2den,d2num
207#endif
208c
209c References:
210c
211c
212c***************************************************************************
213
214      integer n
215      double precision x, sinhm1,dsinhm1
216      double precision rho13, rho43,  Xa,  Ha, denom, num,
217     &                  fprimex, d1num,d1den,
218     &        gamma,fx,x2a,ten6xd,expo,bbx2
219      integer D0R,D1G,D2RR,D2RG,D2GG
220
221c
222      sinhm1(x)=log(x+dsqrt(1d0+x*x))
223      dsinhm1(x)=1d0/dsqrt(1d0+x*x)
224      pi=acos(-1.d0)
225      BETAPW91=(pi*36.d0)**(-5.d0/3.d0)*5.d0
226      AX=-(0.75d0/pi)**(1.d0/3.d0)*1.5d0
227      if(ispin.eq.1) then
228         D0R=1
229         D1G=D1_GAA
230         D2RR=D2_RA_RA
231         D2RG=D2_RA_GAA
232         D2GG=D2_GAA_GAA
233      else
234         D0R=2
235         D1G=D1_GBB
236         D2RR=D2_RB_RB
237         D2RG=D2_RB_GBB
238         D2GG=D2_GBB_GBB
239      endif
240
241      rho13 = (rho*ipol/2d0)**(1.d0/3.d0)
242      rho43 = rho13**4
243      if (gamma.gt.tol_rho**2)then
244         xa = sqrt(gamma)/rho43
245         x2a=xa*xa
246         ten6xd=Xa**DPOW*1.d-6
247         expo=0d0
248         if(CPW91*x2a.lt.big) expo=exp(-CPW91*x2a)
249         bbx2=(BETA-BETAPW91)*x2a
250         Ha = sinhm1(Xa)
251         denom = 1.d0/(1.d0 + 6d0*(beta*xa)*ha-ten6xd/ax)
252         num = -BETA*x2a+bbx2*expo+ten6xd
253         fx=num*denom
254         d1num=-2.d0*xa*(beta-bbx2*expo*(1d0/x2a-CPW91))+
255     +        ten6xd/xa*dpow
256         d1den=6.d0*beta*(ha + xa*dsinhm1(xa)) -
257     -        ten6xd/ax/xa*dpow
258         fprimex=(d1num - d1den*fx)*denom
259      else
260         gamma = 0.d0
261         Xa = 0.d0
262         fx=0.d0
263         fprimex=0.d0
264         denom=0d0
265         d1den=0d0
266         expo=0d0
267         ten6xd=0d0
268         x2a=0d0
269         bbx2=0d0
270      endif
271c
272      if (lfac)then
273         Ex = Ex + 2d0/ipol*rho43*AX*qwght*fac
274         if(ldew)func = func + 2.d0/ipol*rho43*AX*fac
275         Amat(n,D0R) = Amat(n,D0R) + (4.d0/3.d0)*rho13*AX*fac
276      endif
277c
278      if (nlfac)then
279         Ex = Ex + 2d0/ipol*rho43*fx*qwght*fac
280         if(ldew)func = func + 2.d0/ipol*rho43*fx*fac
281         Amat(n,D0R) = Amat(n,D0R) +
282     +        (4.d0/3.d0)*rho13*(fx-xa*fprimex)*fac
283         if (xa.gt.tol_rho)  then
284            Cmat(n,D1G)=Cmat(n,D1G)+
285     +           .5d0*fprimex/sqrt(gamma)*fac
286         endif
287      endif
288c
289#ifdef SECOND_DERIV
290      rhom23 = 2d0/ipol*rho13/rho
291      if (lfac)then
292         Amat2(n,D2RR) = Amat2(n,D2RR) +
293     &        (4d0/9d0)*rhom23*Ax*fac
294      endif
295      if (nlfac)then
296      if(gamma.gt.tol_rho**2)then
297         d2num=-2d0*beta +
298     +        2d0*bbx2*expo*
299     *        (1d0/x2a-5d0*CPW91+2d0*CPW91*CPW91*x2a)+
300     +        ten6xd/x2a*dpow*(dpow-1d0)
301         d2den=6.d0*beta*dsinhm1(xa)*(2d0-x2a/(1d0+x2a)) -
302     -        ten6xd/ax/x2a*dpow*(dpow-1d0)
303         f2x = denom*(d2num -fx*d2den - 2d0*d1den*fprimex)
304         else
305            f2x=0d0
306         endif
307         Amat2(n,D2RR) = Amat2(n,D2RR)
308     &        + (4d0/9d0)*rhom23*(fx-xa*fprimex+4d0*x2a*f2x)*fac
309         Cmat2(n,D2RG) = Cmat2(n,D2RG)
310     &        - (4d0/ipol/3d0)*(rhom23**2/rho)*f2x*fac
311         if (xa.gt.tol_rho) then
312            Cmat2(n,D2GG) = Cmat2(n,D2GG)
313     &           - 0.25d0*gamma**(-1.5d0)*(fprimex-xa*f2x)*fac
314         endif
315      endif
316#endif
317c
318c
319      return
320      end
321#ifndef SECOND_DERIV
322#define SECOND_DERIV
323c
324c     Compile source again for the 2nd derivative case
325c
326#include "xc_xpw6.F"
327#endif
328