1c     Perdew-Wang 86 exchange functional
2c
3c     References:
4c     [a] J. P. Perdew and W. Yue (actually, Y. Wang)
5c         Phys. Rev. B 33 (1986) 8800
6c
7#ifndef SECOND_DERIV
8      Subroutine xc_xpw86(tol_rho, fac, lfac, nlfac, rho, delrho,
9     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
10#else
11      Subroutine xc_xpw86_d2(tol_rho, fac, lfac, nlfac, rho, delrho,
12     &                        Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
13     &                        qwght,ldew,func)
14#endif
15c
16c$Id$
17c
18      implicit none
19c
20#include "dft2drv.fh"
21c
22      double precision tol_rho ! zero-density crit. (input)
23      double precision fac ! energy multiplier (input)
24      logical lfac, nlfac ! local and non-local labels (input)
25      logical ldew        ! compute funct(n) (input)
26      double precision Ex ! exchange energy (output)
27      integer nq  ! number of grid points (input)
28      integer ipol ! spinpolarization flag (input)
29      double precision func(*)  ! value of the functional [output]
30      double precision rho(nq,ipol*(ipol+1)/2) ! charge density (in)
31      double precision delrho(nq,3,ipol) ! gradient of rho (in)
32      double precision qwght(nq) ! quadrature weights (in)
33      double precision amat(nq,ipol) ! dex/drho matrix
34      double precision cmat(nq,*)  ! 2*dex/(dgrho)
35#ifdef SECOND_DERIV
36c     second derivatives? dkdc for now.
37      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
38#endif
39
40c      constants for PW86
41      double precision a, b, c, m1, Ax, F43, s_pfac
42      parameter(a = 1.296d0, b = 14d0, c = 0.2d0)
43      parameter(m1 = 1d0/15d0)
44      parameter(s_pfac = 6.18733545256027d0)
45      parameter(Ax = -0.738558766382022d0)
46      parameter(F43 = 4d0/3d0)
47
48      double precision pi, axr43, rh
49      double precision s, s2, s3, s4, s5, s6
50      double precision fs, sx, drho, drho2, dfds,oneminus
51      integer n
52c
53
54
55      pi = acos(-1.d0)
56
57c     spin-restricted
58      if (ipol.eq.1) then
59#ifdef IFCV81
60CDEC$ NOSWP
61#endif
62         do 10 n = 1, nq
63            if (rho(n,1).lt.tol_rho) goto 10
64
65            rh = rho(n,1)
66            drho2 = (delrho(n,1,1)*delrho(n,1,1) +
67     &              delrho(n,2,1)*delrho(n,2,1) +
68     &              delrho(n,3,1)*delrho(n,3,1))
69            drho = dsqrt(drho2)
70            if (drho.gt.tol_rho) then
71            s = drho / (s_pfac*rh**F43)
72            s2 = s * s
73            s3 = s2 * s
74            s4 = s3 * s
75            s5 = s4 * s
76            s6 = s5 * s
77            fs = (1 + a*s2 + b*s4 + c*s6)**m1
78            dfds = (2*a*s + 4*b*s3 + 6*c*s5)/(15.d0*fs**14)
79            oneminus=1d0-dfds/fs*s
80            else
81               s=0d0
82               fs=1d0
83               dfds=0d0
84               oneminus=1d0
85            endif
86
87            axr43 = Ax * rh**F43
88            sx = axr43 * fs
89            Ex = Ex + sx * qwght(n) * fac
90            if (ldew) func(n) = func(n) + sx * fac
91
92c     GC: dex/drho
93            Amat(n,1) = Amat(n,1) + F43*sx/rh * oneminus
94            if(drho.gt.tol_rho) then
95c     GC: 2*dex/d(grho2)
96               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + sx/fs*dfds*s/drho2
97            endif
98
99c     xxxx missing
100#ifdef SECOND_DERIV
101            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 0
102            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + 0
103            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + 0
104#endif
105 10      continue
106c
107      else
108c     spin unrestricted
109#ifdef IFCV81
110CDEC$ NOSWP
111#endif
112         do 100 n = 1, nq
113            if (rho(n,1).lt.tol_rho) goto 100
114
115c           alpha
116            if (rho(n,2).lt.tol_rho) goto 150
117            rh = rho(n,2) * 2
118            drho2 = (delrho(n,1,1)*delrho(n,1,1) +
119     &              delrho(n,2,1)*delrho(n,2,1) +
120     &              delrho(n,3,1)*delrho(n,3,1)) * 4
121            drho = dsqrt(drho2)
122            s = drho / (s_pfac*rh**F43)
123            s2 = s * s
124            s3 = s2 * s
125            s4 = s3 * s
126            s5 = s4 * s
127            s6 = s5 * s
128
129            axr43 = Ax * rh**F43
130            fs = (1 + a*s2 + b*s4 + c*s6)**m1
131            sx = axr43 * fs
132            Ex = Ex + 0.5d0 * sx * qwght(n) * fac
133            if (ldew) func(n) = func(n) + 0.5d0 * sx * fac
134
135            dfds = (2*a*s + 4*b*s3 + 6*c*s5)/(15.d0*fs**14)
136c     GC: dex/drho
137            Amat(n,1) = Amat(n,1) + F43*sx/rh * (1d0-dfds/fs*s)
138c     GC: 2*dex/d(grho2)
139            if(drho2.gt.tol_rho) then
140               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + sx/fs*dfds*s/drho2
141            endif
142
143c     xxxx missing
144#ifdef SECOND_DERIV
145            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 0
146            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + 0
147            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + 0
148#endif
149
150 150        continue
151c           beta
152            if (rho(n,3).lt.tol_rho) goto 100
153            rh = rho(n,3) * 2
154            drho2 = (delrho(n,1,2)*delrho(n,1,2) +
155     &              delrho(n,2,2)*delrho(n,2,2) +
156     &              delrho(n,3,2)*delrho(n,3,2)) * 4
157            drho = dsqrt(drho2)
158            s = drho / (s_pfac*rh**F43)
159            s2 = s * s
160            s3 = s2 * s
161            s4 = s3 * s
162            s5 = s4 * s
163            s6 = s5 * s
164
165            axr43 = Ax * rh**F43
166            fs = (1 + a*s2 + b*s4 + c*s6)**m1
167            sx = axr43 * fs
168            Ex = Ex + 0.5d0 * sx * qwght(n) * fac
169            if (ldew) func(n) = func(n) + 0.5d0 * sx * fac
170
171            dfds = (2*a*s + 4*b*s3 + 6*c*s5)/(15.d0*fs**14)
172c     GC: dex/drho
173            Amat(n,2) = Amat(n,2) + F43*sx/rh * (1d0-dfds/fs*s)
174c     GC: 2*dex/d(grho2)
175            if(drho2.gt.tol_rho) then
176               Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + sx/fs*dfds*s/drho2
177            endif
178
179c     xxxx missing
180#ifdef SECOND_DERIV
181            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 0
182            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + 0
183            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + 0
184#endif
185
186 100      continue
187      endif
188c
189      return
190      end
191#ifndef SECOND_DERIV
192#define SECOND_DERIV
193c
194c     Compile source again for the 2nd derivative case
195c
196#include "xc_xpw86.F"
197#endif
198