1c
2c$Id$
3c
4#include "dft2drv.fh"
5c    Tao,Perdew, Staroverov, Scuseria exchange functional
6c           META GGA
7C         utilizes ingredients:
8c                              rho   -  density
9c                              delrho - gradient of density
10c                              tau - K.S kinetic energy density
11c                              tauW - von Weiszacker kinetic energy density
12c                              tauU - uniform-gas KE density
13c     References:
14c     [a] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria
15c         PRL 91, 146401 (2003).
16c     [b] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria
17c         JCP 120, 6898 (2004).
18      Subroutine xc_xtpss03(tol_rho, fac,  rho, delrho,
19     &                     Amat, Cmat, nq, ipol, Ex,
20     &                     qwght, ldew, func, tau,Mmat)
21      implicit none
22c
23      double precision fac, Ex
24      integer nq, ipol
25      logical ldew
26      double precision func(*)  ! value of the functional [output]
27c
28c     Charge Density & Its Cube Root
29c
30      double precision rho(nq,ipol*(ipol+1)/2)
31c
32c     Charge Density Gradient
33c
34      double precision delrho(nq,3,ipol)
35c
36c     Quadrature Weights
37c
38      double precision qwght(nq)
39c
40c     Sampling Matrices for the XC Potential & Energy
41c
42      double precision Amat(nq,ipol), Cmat(nq,*)
43c
44c     kinetic energy density   or  tau
45c
46      double precision tau(nq,ipol), Mmat(nq,*)
47      double precision tol_rho
48c
49      integer ispin,cmatpos
50c
51      if (ipol.eq.1 )then
52c
53c     SPIN-RESTRICTED
54c     Ex = Ex[n]
55c
56         call xc_xtpss03_cs(tol_rho, fac,  rho, delrho,
57     &                     Amat, Cmat, nq, Ex, 1d0,
58     &                     qwght, ldew, func, tau,Mmat)
59      else
60c
61c     SPIN-UNRESTRICTED
62c     Ex = Ex[2n_up]/2 + Ex[2n_down]/2
63
64         do ispin=1,2
65            if (ispin.eq.1) cmatpos=D1_GAA
66            if (ispin.eq.2) cmatpos=D1_GBB
67            call xc_xtpss03_cs(tol_rho, fac,
68     R           rho(1,ispin+1), delrho(1,1,ispin),
69     &           Amat(1,ispin), Cmat(1,cmatpos),
70     &           nq, Ex, 2d0,
71     &           qwght, ldew, func,
72     T           tau(1,ispin),Mmat(1,ispin))
73         enddo
74
75      endif
76      return
77      end
78      Subroutine xc_xtpss03_cs(tol_rho, fac,  rho, delrho,
79     &                     Amat, Cmat, nq, Ex, facttwo,
80     &                     qwght, ldew, func, tau,Mmat)
81      implicit none
82c
83      double precision fac, Ex
84      integer nq
85      logical ldew
86      double precision func(*)  ! value of the functional [output]
87c
88c     Charge Density & Its Cube Root
89c
90      double precision rho(*)
91c
92c     Charge Density Gradient
93c
94      double precision delrho(nq,3)
95c
96c     Quadrature Weights
97c
98      double precision qwght(nq)
99c
100c     Sampling Matrices for the XC Potential & Energy
101c
102      double precision Amat(nq), Cmat(nq)
103c
104c     kinetic energy density   or  tau
105c
106      double precision tau(nq,*), Mmat(nq)
107c
108      double precision facttwo ! 2 for o.s. 1 for c.s.
109c
110      double precision tol_rho, pi
111      integer n
112      double precision rrho, rho43, rho13, gamma
113      double precision tauN, tauW, tauU
114
115      double precision  p, qtil, x,  al, mt, z
116      double precision   F83, F23, F53, F43, F13
117      double precision   G920
118      double precision  b,c,e,es
119      double precision    C1, C2, C3
120      double precision    kap, mu
121      double precision xb,xc,xd
122      double precision x1,x2,x3,x4,x5,x6,x7
123      double precision   P32, Ax
124c     functional derivatives below FFFFFFFFFFFF
125      double precision dzdn, dpdn, daldn, dqtdn
126      double precision til1, til2
127      double precision dtil2dn, dtil1dn
128      double precision ax1, bx1, dx1dn
129      double precision dx2dn
130      double precision dxddn, dxcdn, dx3dn
131      double precision dx4dn, dx5dn, dx6dn, dx7dn
132      double precision  dxdn
133      double precision xmany, dxmanydn
134      double precision dmtdn, derivn
135
136      double precision dzdg, dpdg, daldg, dqtdg
137      double precision dtil2dg
138      double precision dx1dg, dx2dg
139      double precision dxcdg, dxddg,dx3dg
140      double precision dx4dg, dx5dg, dx6dg, dx7dg
141      double precision dxmanydg, dxdg, dmtdg, derivg
142
143      double precision dzdt, daldt, dqtdt
144      double precision dx1dt, dx2dt, dx3dt
145      double precision dx5dt
146      double precision dxmanydt, dxdt, dmtdt, derivt
147      double precision afact2
148      double precision rhoval
149
150c     functional derivatives above FFFFFFFFFFFF
151
152      parameter(kap=0.8040d0, mu=0.21951d0)
153      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
154      parameter (F83=8.d0/3.d0, F23=2.d0/3.d0, F53=5.d0/3.d0)
155      parameter (G920  =9.d0/20.d0 )
156
157      parameter(b=0.40d0, c=1.59096d0, e=1.537d0)
158      parameter (C1  =  10.d0/81.d0,
159     &     C2  = 146.d0/2025.d0,
160     &     C3  = -73.d0/405.d0 )
161c
162      pi=acos(-1d0)
163      Ax = (-0.75d0)*(3d0/pi)**F13
164      P32 = (3.d0*pi**2)**F23
165      es=dsqrt(e)
166      afact2=1d0/facttwo
167c
168      do n = 1, nq
169         if (rho(n).ge.tol_rho) then
170
171c     rho43= n*e_x^unif=exchange energy per electron for uniform electron gas
172c     = n* Ax*n^(1/3)   or n*C*n^(1/3)
173
174            rhoval=rho(n)*facttwo
175            rho43 = Ax*rhoval**F43 ! Ax*n^4/3
176            rrho = 1d0/rhoval   ! reciprocal of rho
177            rho13 = rho43*rrho
178
179C     Below we just sum up the LDA contribution to the functional
180            Ex = Ex + rho43*qwght(n)*fac*afact2
181            if (ldew)  func(n)= func(n) + rho43*fac*afact2
182            Amat(n) = Amat(n) + F43*rho13*fac
183
184c
185            gamma = delrho(n,1)*delrho(n,1) +
186     &           delrho(n,2)*delrho(n,2) +
187     &           delrho(n,3)*delrho(n,3)
188            gamma=gamma*facttwo*facttwo
189            tauN = tau(n,1)*facttwo
190            tauW=0.125d0*gamma*rrho
191            tauU=0.3d0*P32*rhoval**F53
192c
193c     Evaluate the Fx, i.e. mt(x) = Fx - 1 (LDA bit already there)
194c
195            p=gamma/(rhoval**F83*P32*4.d0)
196            if(abs(p).gt.tol_rho) then
197            z=tauW/tauN
198            al=(tauN-tauW)/tauU
199c     al=dabs(al)
200            if(al.lt.0d0)  al=0d0
201
202            qtil=(G920*(al-1.d0)/((1.d0+b*al*(al-1.d0))**.5d0)) +
203     +           F23*p
204
205            xb=(c*z**2)/( (1+z**2)**2 )
206            x1=(C1 + xb)*p
207            x2=C2*qtil*qtil
208            xc=C3*qtil
209            xd=(0.5d0*(.6d0*z)**2  + .5d0*p*p)**.5d0
210            x3=xc*xd
211            x4=C1*C1*p*p/kap
212            x5=2.d0*es*C1*(.6d0*z)**2
213            x6= e*mu*p*p*p
214            x7 = (1.d0+es*p)**(-2.d0)
215
216            x=(x1+x2+x3 +x4 +x5+x6)*x7
217
218            if (abs(x).lt.tol_rho) write(0,*) ' x for fx ',x
219
220c     functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF
221
222C     Derivatives wrt n, density   below
223            dzdn=-z*rrho
224            dpdn = -p*rrho*F83
225            daldn=F53*( -p*dzdn/z**2 +dpdn*(-1.d0+1.d0/z) )
226
227            til1=al-1.d0
228            til2=(1.d0+b*al*(al-1.d0))**(-0.5d0)
229            dtil1dn=daldn
230            dtil2dn=b*daldn*(2.d0*al-1d0)*
231     &           (-.5d0)*(til2**3)
232            dqtdn = G920*(til2*dtil1dn+til1*dtil2dn)+F23*dpdn
233
234            ax1=c*p*z*z
235            bx1=(1+z*z)**(-2.d0)
236            dx1dn=(x1/p)*dpdn + 2d0*c*p*z*dzdn/((1d0+z*z)**3)*(1d0-z*z)
237            dx2dn=2.d0*C2*qtil*dqtdn
238
239            dxddn=.5d0/xd*( (3d0/5d0)**2*z*dzdn +
240     +           p*dpdn)
241            dxcdn=C3*dqtdn
242            dx3dn=xc*dxddn+xd*dxcdn
243
244            dx4dn=(2.d0*x4/p)*dpdn
245            dx5dn=(2.d0*x5/z)*dzdn
246            dx6dn=(3.d0*x6/p)*dpdn
247            dx7dn=-2.d0*es*dpdn/(1.d0+es*p)**3
248
249            xmany=x1+x2+x3 +x4 +x5+x6
250            dxmanydn= dx1dn+dx2dn+dx3dn+dx4dn+dx5dn+dx6dn
251            dxdn=x7*dxmanydn+xmany*dx7dn
252C     Derivatives wrt n, density   above
253
254C     Derivatives wrt gamma,    below
255
256            dpdg=p/gamma
257            dzdg=z/gamma
258            daldg=(al/p)*dpdg-F53*(p/(z*z))*dzdg
259
260            dtil2dg=-0.5d0*daldg*b*(2.d0*al-1d0)*til2**3.d0
261            dqtdg=G920*(til1*dtil2dg+til2*daldg)+F23*dpdg
262            dx1dg=(x1/p)*dpdg + 2d0*c*p*z*dzdg/((1d0+z*z)**3)*(1d0-z*z)
263
264            dx2dg=C2*2.d0*qtil*dqtdg
265
266            dxcdg=C3*dqtdg
267            dxddg=.5d0/xd*( (3d0/5d0)**2*z*dzdg +
268     +           p*dpdg)
269            dx3dg=xc*dxddg+xd*dxcdg
270
271            dx4dg=(2.d0*x4/p)*dpdg
272            dx5dg=(2.d0*x5/z)*dzdg
273            dx6dg=(3.d0*x6/p)*dpdg
274
275            dx7dg=-2.d0*es*dpdg*(1.d0+p*es)**(-3.d0)
276
277            dxmanydg= dx1dg+dx2dg+dx3dg+dx4dg+dx5dg+dx6dg
278            dxdg=x7*dxmanydg+xmany*dx7dg
279
280C     Derivatives wrt tau,    below
281c     ttttttttttttttttttttttttttttttttttttttttttttttttt
282            dzdt= -z/tauN
283            daldt=1.d0/tauU
284
285            dqtdt=g920*daldt*til2*(1d0-
286     -           0.5d0*b*til1*til2*til2*(2d0*al-1d0))
287
288            dx1dt=c*p*dzdt*2d0*z*(1d0-z*z)/((1.d0+z*z)**3)
289            dx2dt=2*c2*qtil*dqtdt
290            dx3dt=x3*(dqtdt/qtil +
291     &           0.5d0*(3d0/5d0)**2*z*dzdt/(xd*xd))
292            dx5dt=2d0*(x5/z)*dzdt
293
294            dxmanydt= dx1dt+dx2dt+dx3dt+dx5dt
295            dxdt=x7*dxmanydt
296c     ttttttttttttttttttttttttttttttttttttttttttttttttttt
297
298            mt = kap - kap/(1.d0 + x/kap)
299
300            Ex = Ex + mt*rho43*qwght(n)*fac*afact2
301            if (ldew)  func(n)= func(n) + mt*rho43*fac*afact2
302
303            dmtdn=dxdn/(1.d0+x/kap)**2
304            derivn=mt*F43*rho13+rho43*dmtdn
305
306            dmtdg=dxdg/(1.d0+x/kap)**2
307            derivg = rho43*dmtdg
308
309            dmtdt=dxdt/(1.d0+x/kap)**2
310            derivt = rho43*dmtdt
311            Amat(n) = Amat(n) + derivn*fac
312c
313c     4x factor comes from gamma_aa = gamma_total/4
314c
315            Cmat(n)=  Cmat(n) + 2d0*derivg*fac
316            Mmat(n)=  Mmat(n) +0.5d0*derivt*fac
317         endif
318         endif
319      enddo
320      return
321      end
322
323      Subroutine xc_xtpss03_d2()
324      call errquit(' xtpss03: d2 not coded ',0,0)
325      return
326      end
327
328
329