1c   VS98  exchange functional
2c           META GGA
3C         utilizes ingredients:
4c                              rho   -  density
5c                              delrho - gradient of density
6c                              tau - K.S kinetic energy density
7c                              tauU - uniform-gas KE density
8c                              ijzy - 1  VS98
9c                              ijzy - 2  M06-L
10c                              ijzy - 3  M06-HF
11c                              ijzy - 4  M06
12c                              ijzy - 6  revM06-L
13c                              ijzy - 7  revM06
14c
15c     References:
16c
17c     [a] T. V. Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998).
18c     [b] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006).
19
20
21
22      Subroutine xc_xvs98(tol_rho, fac,lfac,nlfac, rho, delrho,
23     &                     Amat, Cmat, nq, ipol, Ex,
24     &                     qwght, ldew, func, tau, Mmat, ijzy)
25
26
27c
28c$Id$
29c
30      implicit none
31c
32#include "errquit.fh"
33c
34      double precision fac, Ex
35      integer nq, ipol
36      logical lfac, nlfac,ldew
37      double precision func(*)  ! value of the functional [output]
38c
39c     Charge Density & Its Cube Root
40c
41      double precision rho(nq,ipol*(ipol+1)/2)
42c
43c     Charge Density Gradient
44c
45      double precision delrho(nq,3,ipol)
46c
47c     Quadrature Weights
48c
49      double precision qwght(nq)
50c
51c     Sampling Matrices for the XC Potential & Energy
52c
53      double precision Amat(nq,ipol), Cmat(nq,*)
54
55
56      double precision tol_rho, pi
57
58c
59      integer n, ijzy
60      double precision rrho, rho43, rho13, rhoo, rho53, rho83
61      double precision  Gamma
62c
63c     kinetic energy density or tau
64c
65      double precision tau(nq,ipol), Mmat(nq,*)
66
67      double precision tauN,tauu,DTol
68      double precision Tiny, f13, f43, f53, f83, f113
69      double precision gx, gg, x, z, kx,xk,zk
70      double precision One, Two, Three, Four, Five, Six, Seven, Eight
71      double precision Nine, F10, F11
72      double precision cf, Axlsda, r1, r2, r3, r4, r5, r6
73
74c      functional derivatives below FFFFFFFFFFFF
75
76       double precision dxdr, dxdg, dzdr, dzdt, dgdx, dgdz
77
78c      functional derivatives above FFFFFFFFFFFF
79
80
81cedo       parameter( pi = 3.1415926535897932384626433832795d0 )
82
83       parameter (cf = 9.115599720d0, Axlsda = -0.9305257363491d0 )
84       parameter (gg  = 0.00186726d0)
85       parameter (f13=1.d0/3.d0,f43=4.0d0/3.0d0,f53=5.0d0/3.0d0)
86       parameter (f83=8.d0/3.0d0, F113=11.0d0/3.d0)
87       parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0,
88     &             Five=5.0d0,Six=6.0d0, Seven=7.0d0,
89     &             Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0)
90      pi=acos(-1d0)
91
92
93      if (ijzy.eq.1) then
94c
95c     Parameters for VS98
96c
97        r1=  -9.800683d-01
98        r2=  -3.556788d-03
99        r3=   6.250326d-03
100        r4=  -2.354518d-05
101        r5=  -1.282732d-04
102        r6=   3.574822d-04
103      elseif (ijzy.eq.2) then
104c
105c     Parameters for M06-L
106c
107        r1 =   6.012244D-01*Axlsda
108        r2 =   4.748822D-03*Axlsda
109        r3 =  -8.635108D-03*Axlsda
110        r4 =  -9.308062D-06*Axlsda
111        r5 =   4.482811D-05*Axlsda
112        r6 =   0.000000D+00
113      elseif (ijzy.eq.3) then
114c
115c     Parameters for M06-HF
116c
117        r1 =   -1.179732D-01*Axlsda
118        r2 =   -2.500000D-03*Axlsda
119        r3 =   -1.180065D-02*Axlsda
120        r4 =   0.000000D+00
121        r5 =   0.000000D+00
122        r6 =   0.000000D+00
123      elseif (ijzy.eq.4) then
124c
125c     Parameters for M06
126c
127        r1 =   1.422057D-01*Axlsda
128        r2 =   7.370319D-04*Axlsda
129        r3 =   -1.601373D-02*Axlsda
130        r4 =   0.000000D+00
131        r5 =   0.000000D+00
132        r6 =   0.000000D+00
133      elseif (ijzy.eq.6) then
134c
135c     Parameters for revM06-L
136c
137        r1 =   -4.23227252D-01*Axlsda
138        r2 =   0.000000D+00*Axlsda
139        r3 =   3.724234D-03*Axlsda
140        r4 =   0.000000D+00*Axlsda
141        r5 =   0.000000D+00*Axlsda
142        r6 =   0.000000D+00
143      elseif (ijzy.eq.7) then
144c
145c     Parameters for revM06
146c
147        r1 =   -5.523940140D-02*Axlsda
148        r2 =   0.000000D+00*Axlsda
149        r3 =  -3.782631233D-03*Axlsda
150        r4 =   0.000000D+00*Axlsda
151        r5 =   0.000000D+00*Axlsda
152        r6 =   0.000000D+00
153      else
154        call errquit("xc_xvs98: illegal value of ijzy",ijzy,UERR)
155      endif
156
157      DTol = tol_rho
158
159
160c
161
162
163c
164      if (ipol.eq.1 )then
165c
166c        ======> SPIN-RESTRICTED <======
167c                     or
168c                SPIN-UNPOLARIZED
169c
170c
171         do 10 n = 1, nq
172            rhoo = rho(n,1)/Two   ! rho_sigma
173            if (rhoo.lt.DTol) goto 10
174            rho43 = rhoo**f43
175            rrho = 1d0/rhoo       ! reciprocal of rho
176            rho13 = rho43*rrho
177            rho53 = rhoo**f53
178            rho83 = rho53*rhoo
179c
180
181            tauN = tau(n,1)
182            if(taun.lt.dtol) goto 10
183            tauu=tauN
184            gamma = delrho(n,1,1)*delrho(n,1,1) +
185     &           delrho(n,2,1)*delrho(n,2,1) +
186     &           delrho(n,3,1)*delrho(n,3,1)
187         if(sqrt(gamma).lt.dtol) goto 10
188            gamma = gamma/Four
189            x = gamma/rho83
190            dxdr = -f83*x*rrho
191            dxdg = One/rho83
192            z = tauu/rho53 - cf
193            dzdr = -f53 * tauu/rho83
194            dzdt = One/rho53
195            kx = One + gg*x + gg*z
196            xk = x/kx
197            zk = z/kx
198            call gvt4(gx,dgdx,dgdz,xk,zk,kx,gg,gg,r1,r2,r3,r4,r5,r6)
199
200            Ex = Ex + Two*rho43*gx*qwght(n)
201            if(ldew) func(n)=func(n)+ Two*rho43*gx
202c
203c     functional derivatives
204c
205            Amat(n,1) = Amat(n,1) + f43*rho13*gx +
206     &                  rho43*(dgdx*dxdr + dgdz*dzdr)
207            Cmat(n,1)=  Cmat(n,1) + rho43*(dgdx*dxdg)
208            Mmat(n,1)=  Mmat(n,1) + rho43*(dgdz*dzdt)
209
21010      continue
211
212
213c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
214      else  ! ipol=2
215c
216c        ======> SPIN-UNRESTRICTED <======
217
218
219         do 20 n = 1, nq
220           if (rho(n,1).lt.DTol) goto 20
221c
222c     Alpha            ALPHA               ALPHA
223c
224            if (rho(n,2).lt.DTol) goto 25
225             rhoo = rho(n,2)
226             rho43 = rhoo**f43
227             rrho = 1.0d0/rhoo       ! reciprocal of rho
228             rho13 = rho43*rrho
229             rho53 = rhoo**f53
230             rho83 = rho53*rhoo
231
232c
233
234             tauN = tau(n,1)
235             tauu = tauN*Two
236
237            Gamma = delrho(n,1,1)*delrho(n,1,1) +
238     &              delrho(n,2,1)*delrho(n,2,1) +
239     &              delrho(n,3,1)*delrho(n,3,1)
240
241            x = gamma/rho83
242            dxdr = -f83*x*rrho
243            dxdg = One/rho83
244            z = tauu/rho53 - cf
245            dzdr = -f53 * tauu/rho83
246            dzdt = One/rho53
247            kx = One + gg*x + gg*z
248            xk = x/kx
249            zk = z/kx
250            call gvt4(gx,dgdx,dgdz,xk,zk,kx,gg,gg,r1,r2,r3,r4,r5,r6)
251            Ex = Ex + rho43*gx*qwght(n)
252            if(ldew) func(n)=func(n)+ rho43*gx
253c
254c     functional derivatives
255c
256            Amat(n,1) = Amat(n,1) + f43*rho13*gx +
257     &                  rho43*(dgdx*dxdr + dgdz*dzdr)
258            Cmat(n,1)=  Cmat(n,1) + rho43*(dgdx*dxdg)
259            Mmat(n,1)=  Mmat(n,1) + rho43*(dgdz*dzdt)
260
261
262c
263c     Beta               BETA           BETA
264c
265
26625         continue
267
268c
269c     Beta
270c
271            if (rho(n,3).lt.DTol) goto 20
272             rhoo = rho(n,3)
273             rho43 = rhoo**f43
274             rrho = 1.0d0/rhoo       ! reciprocal of rho
275             rho13 = rho43*rrho
276             rho53 = rhoo**f53
277             rho83 = rho53*rhoo
278
279c
280
281             tauN = tau(n,2)
282             tauu = Two*tauN
283
284
285            Gamma = delrho(n,1,2)*delrho(n,1,2) +
286     &              delrho(n,2,2)*delrho(n,2,2) +
287     &              delrho(n,3,2)*delrho(n,3,2)
288
289            x = gamma/rho83
290            dxdr = -f83*x*rrho
291            dxdg = One/rho83
292            z = tauu/rho53 - cf
293            dzdr = -f53 * tauu/rho83
294            dzdt = One/rho53
295            kx = One + gg*x + gg*z
296            xk = x/kx
297            zk = z/kx
298            call gvt4(gx,dgdx,dgdz,xk,zk,kx,gg,gg,r1,r2,r3,r4,r5,r6)
299
300            Ex = Ex + rho43*gx*qwght(n)
301            if(ldew) func(n)=func(n)+ rho43*gx
302c
303c     functional derivatives
304c
305            Amat(n,2) = Amat(n,2) + f43*rho13*gx +
306     &                  rho43*(dgdx*dxdr + dgdz*dzdr)
307            Cmat(n,3)=  Cmat(n,3) + rho43*(dgdx*dxdg)
308            Mmat(n,2)=  Mmat(n,2) + rho43*(dgdz*dzdt)
309
310c
31120      continue
312      endif
313c
314      return
315      end
316
317
318      Subroutine xc_xvs98_d2()
319      implicit none
320      call errquit(' xvs98: d2 not coded ',0,0)
321      return
322      end
323
324      Subroutine gvt4(g,dgdx,dgdz,xk,zk,k,c,ct,r1,r2,r3,r4,r5,r6)
325      Implicit none
326c
327c     Evaluate the GVT4 form in the VS98 functional
328c
329c
330c    some working variables
331      double precision g,dgdx,dgdz,xk,zk,k,c,ct,r1,r2,r3,r4,r5,r6
332      double precision sxk,szk,sk,sc,sct,sr1,sr2,sr3,sr4,sr5,sr6,sk2
333      double precision One, Two, Three, Four, Six
334      Data One/1.0d0/, Two/2.0d0/, Three/3.0d0/, Four/4.0d0/, Six/6.0d0/
335C
336      sxk = xk
337      szk = zk
338      sk = k
339      sc = c
340      sct = ct
341      sr1 = r1
342      sr2 = r2
343      sr3 = r3
344      sr4 = r4
345      sr5 = r5
346      sr6 = r6
347      sk2 = sk*sk
348      g =  (sr1 + sr2*sxk + sr3*szk
349     $  +sr4*sxk*sxk + sr5*szk*sxk + sr6*szk*szk)/sk
350      dgdx =   (-sr1*sc
351     $  +sr2*(One-Two*sc*sxk)
352     $  -Two*sr3*szk*sc
353     $  +sr4*(Two*sxk-Three*sxk*sxk*sc)
354     $  +sr5*(szk -Three*szk*sxk*sc)
355     $  -Three*sr6*szk*szk*sc )/sk2
356      dgdz =   (-sr1*sct
357     $  -Two*sr2*sxk*sct
358     $  +sr3*(One-Two*szk*sct)
359     $  -Three*sr4*sxk*sxk*sct
360     $  +sr5*(sxk-Three*sxk*szk*sct)
361     $  +sr6*(Two*szk-Three*szk*szk*sct))/sk2
362
363      return
364      end
365
366
367