1#ifndef SECOND_DERIV
2      Subroutine xc_xpw91(tol_rho, fac, lfac, nlfac, rho, delrho,
3     &                      Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
4#else
5      Subroutine xc_xpw91_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
46C
47C
48
49#ifdef SECOND_DERIV
50c
51c     Second Derivatives of the Exchange Energy Functional
52c
53      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
54#endif
55c
56c References:
57c
58c
59c***************************************************************************
60
61      integer n
62      double precision gamma
63c
64      if (ipol.eq.1 )then
65c
66c        ======> SPIN-RESTRICTED <======
67c
68         do 10 n = 1, nq
69            if (rho(n,1).lt.tol_rho) goto 10
70            gamma = delrho(n,1,1)*delrho(n,1,1) +
71     &           delrho(n,2,1)*delrho(n,2,1) +
72     &           delrho(n,3,1)*delrho(n,3,1)
73            gamma=0.25d0*gamma
74#ifndef SECOND_DERIV
75            call xc_xpw91core(DPOW,BETA,n,1,
76     &           rho(n,1),gamma,qwght(n),func(n),
77     &           tol_rho, fac, lfac, nlfac,
78     &           Amat, Cmat, nq, ipol, Ex, ldew)
79#else
80               call xc_xpw91core_d2(DPOW,BETA,n,1,
81     &           rho(n,1),gamma,qwght(n),func(n),
82     &           tol_rho, fac, lfac, nlfac,
83     &           Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
84     &           ldew)
85#endif
86
87   10    continue
88c
89      else
90c
91c        ======> SPIN-UNRESTRICTED <======
92c
93         do 20 n = 1, nq
94           if (rho(n,1).lt.tol_rho) goto 20
95c
96c           Spin alpha:
97c
98           if (rho(n,2).gt.tol_rho) then
99              gamma =    delrho(n,1,1)*delrho(n,1,1) +
100     &             delrho(n,2,1)*delrho(n,2,1) +
101     &             delrho(n,3,1)*delrho(n,3,1)
102#ifndef SECOND_DERIV
103               call xc_xpw91core(DPOW,BETA,n,1,
104     &           rho(n,2),gamma,qwght(n),func(n),
105     &              tol_rho, fac, lfac, nlfac,
106     &              Amat, Cmat, nq, ipol, Ex, ldew)
107#else
108               call xc_xpw91core_d2(DPOW,BETA,n,1,
109     &           rho(n,2),gamma,qwght(n),func(n),
110     &              tol_rho, fac, lfac, nlfac,
111     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
112     &                         ldew)
113#endif
114            endif
115c
116c           Spin beta:
117c
118            if (rho(n,3).gt.tol_rho) then
119
120            gamma =   delrho(n,1,2)*delrho(n,1,2) +
121     &           delrho(n,2,2)*delrho(n,2,2) +
122     &           delrho(n,3,2)*delrho(n,3,2)
123#ifndef SECOND_DERIV
124            call xc_xpw91core(DPOW,BETA,n,2,
125     &           rho(n,3),gamma,qwght(n),func(n),
126     &           tol_rho, fac, lfac, nlfac,
127     &           Amat, Cmat, nq, ipol, Ex, ldew)
128#else
129               call xc_xpw91core_d2(DPOW,BETA,n,2,
130     &           rho(n,3),gamma,qwght(n),func(n),
131     &           tol_rho, fac, lfac, nlfac,
132     &           Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
133     &           ldew)
134#endif
135            endif
136c
137   20    continue
138c
139      endif
140c
141      return
142      end
143#ifndef SECOND_DERIV
144#define SECOND_DERIV
145c
146c     Compile source again for the 2nd derivative case
147c
148#include "xc_xpw91.F"
149#endif
150