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