1c
2c     LDA screened exchange functional
3c
4c     BGJ - 2/01
5c
6c     Reference:
7c
8c     D. M. Bylander and L. Kleinman, Phys. Rev. B 41, 7868 (1990).
9c
10#ifndef SECOND_DERIV
11      Subroutine xc_dirac_screened(tol_rho, fac, lfac, nlfac, rho,
12     &     Amat, nq, ipol, Ex, Ks, qwght, ldew, func)
13#else
14#include "dft2drv.fh"
15      Subroutine xc_dirac_screened_d2(tol_rho, fac, lfac, nlfac, rho,
16     &     Amat, Amat2, nq, ipol, Ex, Ks, qwght, ldew, func)
17#endif
18c
19C$Id$
20c
21      implicit none
22#include "errquit.fh"
23c
24c
25      double precision tol_rho, fac, Ex
26      integer nq, ipol
27      logical lfac, nlfac,ldew
28      double precision Ks       ! Thomas-Fermi screening constant [input]
29      double precision func(*)  ! value of the functional [output]
30c
31c     Charge Density
32c
33      double precision rho(nq,ipol*(ipol+1)/2)
34c
35c     Quadrature Weights
36c
37      double precision qwght(nq)
38c
39c     Sampling Matrices for the XC Potential
40c
41      double precision Amat(nq,ipol)
42#ifdef SECOND_DERIV
43      double precision Amat2(nq,NCOL_AMAT2)
44#endif
45c
46      double precision ONE3, FOUR3
47      Parameter (ONE3 = 1d0/3d0, FOUR3 = 4d0/3d0)
48      integer n
49      double precision pi, C, TPP, rhoval, rho13, rho43, z, d1z,
50     &     f, d1f
51      double precision kf, d1kf, f1, f2, d1f1, d1f2
52      double precision rz, Tworz, atan2rz, datan2rz, rdatan, lnrdatan,
53     &     z243, Eightrz2, Eightrz3, dzdp
54#ifdef SECOND_DERIV
55      double precision d2z, d2f
56c
57      call errquit('Must implement 2nd derivs in xc_dirac_screened',0,
58     &       CAPMIS_ERR)
59#endif
60c
61      pi = acos(-1d0)
62      C = -1.5d0*(0.75d0/pi)**ONE3
63      TPP = (3.d0*pi*pi)**ONE3
64c
65      if (ipol.eq.1) then
66c
67c        ======> SPIN-RESTRICTED <======
68c
69         do 10 n = 1, nq
70            if (rho(n,1).lt.tol_rho) goto 10
71            rhoval = 0.5d0*rho(n,1)
72c
73c           Spin alpha:
74c
75            rho13 = rhoval**ONE3
76            rho43 = rho13*rhoval
77c
78            kf = TPP*rho13
79            z = Ks/kf
80            rz = 1d0/z
81            Tworz = 2d0*rz
82            Eightrz2 = 4d0*Tworz*rz
83            Eightrz3 = Eightrz2*rz
84            rdatan = 1d0 + Tworz*Tworz
85            datan2rz = 1d0/rdatan
86            lnrdatan = log(rdatan)
87            atan2rz = atan(Tworz)
88            z243 = z*z/4d0 + 3d0
89            f1 = -z*z/6d0
90            f2 = 1d0 - z243*lnrdatan + 8d0*rz*atan2rz
91            d1f1 = -z/3d0
92            d1f2 = -z*0.5d0*lnrdatan - Eightrz2*atan2rz
93     &           - Eightrz3*datan2rz*(2d0 - z243)
94            f = 1d0 + f1*f2
95            dzdp = -ONE3*z/rhoval
96            d1f = (d1f1*f2 + f1*d1f2)*dzdp
97c
98            if (lfac) then
99               Ex = Ex + 2d0*rho43*C*f*qwght(n)*fac
100               if(ldew)func(n) = func(n) + 2.d0*rho43*C*f*fac
101               Amat(n,1) = Amat(n,1) + (FOUR3*rho13*f+rho43*d1f)*C*fac
102#ifdef SECOND_DERIV
103               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
104     &              + ONE3*(FOUR3*rho13/rhoval)*C*fac
105#endif
106            endif
107c
108 10      continue
109c
110      else
111c
112c        ======> SPIN-UNRESTRICTED <======
113c
114         do 20 n = 1, nq
115            if (rho(n,1).lt.tol_rho) goto 20
116            if (rho(n,2).lt.tol_rho) goto 25
117c
118c           Spin alpha:
119c
120            rhoval = rho(n,2)
121            rho13 = rhoval**ONE3
122            rho43 = rho13*rhoval
123c
124            kf = TPP*rho13
125            z = Ks/kf
126            rz = 1d0/z
127            Tworz = 2d0*rz
128            Eightrz2 = 4d0*Tworz*rz
129            Eightrz3 = Eightrz2*rz
130            rdatan = 1d0 + Tworz*Tworz
131            datan2rz = 1d0/rdatan
132            lnrdatan = log(rdatan)
133            atan2rz = atan(Tworz)
134            z243 = z*z/4d0 + 3d0
135            f1 = -z*z/6d0
136            f2 = 1d0 - z243*lnrdatan + 8d0*rz*atan2rz
137            d1f1 = -z/3d0
138            d1f2 = -z*0.5d0*lnrdatan - Eightrz2*atan2rz
139     &           - Eightrz3*datan2rz*(2d0 - z243)
140            f = 1d0 + f1*f2
141            dzdp = -ONE3*z/rhoval
142            d1f = (d1f1*f2 + f1*d1f2)*dzdp
143c
144            if (lfac) then
145               Ex = Ex + rho43*C*f*qwght(n)*fac
146               if(ldew)func(n) = func(n) + rho43*C*f*fac
147               Amat(n,1) = Amat(n,1) + (FOUR3*rho13*f+rho43*d1f)*C*fac
148#ifdef SECOND_DERIV
149               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
150     &              + ONE3*(FOUR3*rho13/rhoval)*C*fac
151#endif
152            endif
153c
154 25         continue
155c
156c           Spin beta:
157c
158            if (rho(n,3).lt.tol_rho) goto 20
159c
160            rhoval = rho(n,3)
161            rho13 = rhoval**ONE3
162            rho43 = rho13*rhoval
163c
164            kf = TPP*rho13
165            z = Ks/kf
166            rz = 1d0/z
167            Tworz = 2d0*rz
168            Eightrz2 = 4d0*Tworz*rz
169            Eightrz3 = Eightrz2*rz
170            rdatan = 1d0 + Tworz*Tworz
171            datan2rz = 1d0/rdatan
172            lnrdatan = log(rdatan)
173            atan2rz = atan(Tworz)
174            z243 = z*z/4d0 + 3d0
175            f1 = -z*z/6d0
176            f2 = 1d0 - z243*lnrdatan + 8d0*rz*atan2rz
177c            f2 = 1d0 - z243*log(1d0+rz) + 8d0*rz*atan2rz
178            d1f1 = -z/3d0
179            d1f2 = -z*0.5d0*lnrdatan - Eightrz2*atan2rz
180     &           - Eightrz3*datan2rz*(2d0 - z243)
181            f = 1d0 + f1*f2
182            dzdp = -ONE3*z/rhoval
183            d1f = (d1f1*f2 + f1*d1f2)*dzdp
184c
185            if (lfac) then
186               Ex = Ex + rho43*C*f*qwght(n)*fac
187               if(ldew)func(n) = func(n) + rho43*C*f*fac
188               Amat(n,2) = Amat(n,2) + (FOUR3*rho13*f+rho43*d1f)*C*fac
189#ifdef SECOND_DERIV
190               Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
191     &              + ONE3*(FOUR3*rho13/rhoval)*C*fac
192#endif
193            endif
194c
195 20      continue
196c
197      endif
198c
199      return
200      end
201#ifndef SECOND_DERIV
202#define SECOND_DERIV
203c
204c     Compile source again for the 2nd derivative case
205c
206#include "xc_dirac_screened.F"
207#endif
208