1#include "dft2drv.fh"
2c    VSXC correlation functional
3c           META GGA
4C         utilizes ingredients:
5c                              rho   -  density
6c                              delrho - gradient of density
7c                              tau (tauN)- K.S kinetic energy density
8c                              ijzy - 1  VS98
9c                              ijzy - 2  M06-L
10c                              ijzy - 3  M06-HF
11c                              ijzy - 4  M06
12c                              ijzy - 5  M06-2X
13c                              ijzy - 6  revM06
14c                              ijzy - 7  revM06-L
15c
16      Subroutine xc_cvs98(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
17     &                     nq, ipol, Ec, qwght, ldew, func,
18     &                     tau, Amat, Cmat, Mmat, ijzy)
19
20
21c
22c$Id$
23c
24c  Reference
25c   [a] T. V. Voorhis and G. E. Scuseria, J. Chem. Phys. 109, 400 (1998).
26c   [b] Y. Zhao and D. G. Truhlar, J. Chem. Phys. 125, 194101 (2006).
27
28c
29      implicit none
30c
31#include "errquit.fh"
32c
33c
34c     Input and other parameters
35c
36      integer ipol, nq
37
38      double precision cfac
39      logical lcfac, nlcfac
40
41      logical lfac, nlfac
42      double precision fac
43      double precision tol_rho
44c
45      logical ldew
46      double precision func(*)
47c
48c     Threshold parameters
49c
50      double precision DTol,F1, F2, F3, F4, gab, cf
51      Data F1/1.0d0/,F2/2.0d0/,
52     & F3/3.0d0/,F4/4.0d0/,gab/0.00304966d0/,
53     & cf/9.115599720d0/
54c
55c     Correlation energy
56c
57      double precision Ec
58c
59c     Charge Density
60c
61      double precision rho(nq,ipol*(ipol+1)/2)
62c
63c     Charge Density Gradient
64c
65      double precision delrho(nq,3,ipol)
66
67c
68c     Kinetic Energy Density
69c
70      double precision tau(nq,ipol), tauN
71
72c
73c     Quadrature Weights
74c
75      double precision qwght(nq)
76c
77c     Sampling Matrices for the XC Potential
78c
79      double precision Amat(nq,ipol), Cmat(nq,*)
80      double precision Mmat(nq,*)
81
82      integer n, ijzy
83
84c    call to the vs98css subroutine
85      double precision PA,GAA,TA,FA,FPA,FGA,FTA,EUA,EUEGA,ChiA,EUPA
86     &,ChiAP,ChiAG,ZA,ZAP,ZAT
87      double precision PB,GBB,TB,FB,FPB,FGB,FTB,EUB,EUEGB,ChiB,EUPB
88     &,ChiBP,ChiBG,ZB,ZBP,ZBT
89c
90      double precision Pi, F43, F13, Pi34, F6, PotLC,
91     &RS,RSP,Zeta,dZdA,dZdB,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ
92      double precision P, EUEG, ZAB, XAB, kab, xk, zk
93      double precision dgdx,dgdz,dgdPA,dgdGA,dgdTA,dgdPB,dgdGB,dgdTB
94      double precision EUEGPA,EUEGPB,gcab
95      double precision r7, r8, r9, r10, r11, r12
96
97
98c
99c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
100c
101c      DTol=1.0d-7
102      dtol=tol_rho
103C     Parameters for VS98
104      if (ijzy.eq.1) then
105              r7=   7.035010d-01
106              r8=   7.694574d-03
107              r9=   5.152765d-02
108              r10=   3.394308d-05
109              r11=  -1.269420d-03
110              r12=   1.296118d-03
111C     Parameters for M06-L
112      elseif (ijzy.eq.2) then
113              r7=      3.957626D-01
114              r8=      -5.614546D-01
115              r9=      1.403963D-02
116              r10=     9.831442D-04
117              r11=     -3.577176D-03
118              r12=     0.000000D+00
119C     Parameters for M06-HF
120      elseif (ijzy.eq.3) then
121              r7=    -6.746338D-01
122              r8=    -1.534002D-01
123              r9=    -9.021521D-02
124              r10=   -1.292037D-03
125              r11=   -2.352983D-04
126              r12=   0.000000D+00
127
128C     Parameters for M06
129      elseif (ijzy.eq.4) then
130               r7= -2.741539D+00
131               r8= -6.720113D-01
132               r9= -7.932688D-02
133               r10=1.918681D-03
134               r11=-2.032902D-03
135               r12=0.000000D+00
136
137C     Parameters for M06-2X
138      elseif (ijzy.eq.5) then
139              r7=  1.166404D-01
140              r8=  -9.120847D-02
141              r9=  -6.726189D-02
142              r10= 6.720580D-05
143              r11= 8.448011D-04
144              r12= 0.000000D+00
145C     Parameters for revM06-L
146      elseif (ijzy.eq.6) then
147              r7=  4.007146D-01
148              r8=  1.5796569D-02
149              r9=  -3.2680984D-02
150              r10= 0.000000D+00
151              r11= 0.000000D+00
152              r12= 1.260132D-03
153C     Parameters for revM06
154      elseif (ijzy.eq.7) then
155              r7=  -3.390666720D-01
156              r8=  3.790156384D-03
157              r9=  -2.762485975D-02
158              r10= 0.000000D+00
159              r11= 0.000000D+00
160              r12= 4.076285162D-04
161      else
162         call errquit("xc_cvs98: illegal value of ijzy",ijzy,UERR)
163      endif
164      Pi = F4*ATan(F1)
165      F6=6.0d0
166      F43 = F4 / F3
167      Pi34 = F3 / (F4*Pi)
168      F13 = F1 / F3
169
170
171      do 20 n = 1, nq
172       if (rho(n,1).lt.DTol) goto 20
173       if (ipol.eq.1) then
174c
175c    get the density, gradient, and tau for the alpha spin from the total
176c
177         PA = rho(n,1)/F2
178         GAA = (    delrho(n,1,1)*delrho(n,1,1) +
179     &                 delrho(n,2,1)*delrho(n,2,1) +
180     &                 delrho(n,3,1)*delrho(n,3,1))/4.0d0
181c  In the bc95css subroutine, we use 2*TA as the tau, so we do not divide
182c  the tau by 2 here
183
184         TA = tau(n,1)
185
186         Call vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
187     &                ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,ijzy)
188         PB = PA
189         GBB = GAA
190         TB = TA
191         FB = FA
192         FPB = FPA
193         FGB = FGA
194         FTB = FTA
195         EUB = EUA
196         ZB = ZA
197         ChiB = ChiA
198         EUPB = EUPA
199         ChiBP = ChiAP
200         ChiBG = ChiAG
201         ZBP = ZAP
202         ZBT = ZAT
203
204         Ec = Ec + 2.0d0*FA*qwght(n)            !factor of 2 account for both spin
205         if(ldew) func(n)=func(n)+ 2.0d0*FA
206         Amat(n,1)=Amat(n,1)+ FPA
207         Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + FGA
208         Mmat(n,1)=  Mmat(n,1) + FTA
209c if 0
210c       write (0,'(A,3F20.6)') " PA,EUA",PA,EUA
211c endif
212
213
214c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
215      else  ! ipol=2
216c
217c        ======> SPIN-UNRESTRICTED <======
218c
219c
220c       alpha
221c
222
223         PA = rho(n,2)
224         if (PA.le.DTol) go to 25
225         GAA =   delrho(n,1,1)*delrho(n,1,1) +
226     &           delrho(n,2,1)*delrho(n,2,1) +
227     &          delrho(n,3,1)*delrho(n,3,1)
228c
229c  In the bc95css subroutine, we use 2*TA as the tau
230c
231         TA = tau(n,1)*2.0d0
232
233         Call vs98ss(tol_rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,ZA,
234     &                ChiA,EUPA,ChiAP,ChiAG,ZAP,ZAT,ijzy)
235         Ec = Ec + FA*qwght(n)
236         if(ldew) func(n)=func(n)+ FA
237         Amat(n,1)=Amat(n,1)+ FPA
238         Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + FGA
239c      2*0.5=1.0 for Mmat
240         Mmat(n,1)=  Mmat(n,1) + FTA
241#if 0
242      write (0,'(A,3F20.6)') "AAmat Cmat Mmat",FPA,FGA,FTA
243#endif
244
245c
246c  In the vs98ss subroutine, we use 2*TA as the tau,
247c
248c
249c       Beta
250c
251 25       continue
252         PB = rho(n,3)
253         GBB =   delrho(n,1,2)*delrho(n,1,2) +
254     &           delrho(n,2,2)*delrho(n,2,2) +
255     &          delrho(n,3,2)*delrho(n,3,2)
256
257         TB = tau(n,2)*2.0d0
258
259         Call vs98ss(tol_rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,ZB,
260     &                ChiB,EUPB,ChiBP,ChiBG,ZBP,ZBT,ijzy)
261         Ec = Ec + FB*qwght(n)
262         if(ldew) func(n)=func(n)+ FB
263         Amat(n,2)= Amat(n,2)+ FPB
264         Cmat(n,D1_GBB)=  Cmat(n,D1_GBB) + FGB
265         Mmat(n,2)=  Mmat(n,2) + FTB
266
267#if 0
268      write (0,'(A,3F20.6)') "BAmat Cmat Mmat",FPB,FGB,FTB
269#endif
270      endif
271 30   continue
272      P = rho(n,1)
273      If(PA.gt.DTol.and.PB.gt.DTol) then
274          RS = (Pi34/P) ** F13
275          RSP = -RS/(F3*P)
276          Zeta = (PA-PB)/P
277          dZdA = (F1-Zeta)/P
278          dZdB = (-F1-Zeta)/P
279          Call lsdac(tol_rho,
280     A         RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,
281     $      d2LdZZ)
282          EUEG = P*PotLC - EUA - EUB
283          ZAB = ZA + ZB
284          XAB = ChiA+ChiB
285          kab = F1 + gab*(XAB+ZAB)
286          xk = XAB/kab
287          zk = ZAB/kab
288       call gvt4(gcab,dgdx,dgdz,xk,zk,kab,gab,gab,r7,r8,r9,r10,r11,r12)
289          Ec = Ec + gcab*EUEG*qwght(n)
290          if(ldew) func(n)=func(n)+ gcab*EUEG
291          dgdPA = dgdx*ChiAP + dgdz*ZAP
292          dgdGA = dgdx*ChiAG
293          dgdTA = dgdz*ZAT
294          dgdPB = dgdx*ChiBP + dgdz*ZBP
295          dgdGB = dgdx*ChiBG
296          dgdTB = dgdz*ZBT
297          EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA
298          EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB
299          if (ipol.eq.1) then
300           Amat(n,1) = Amat(n,1) + (EUEGPA*gcab + EUEG*dgdPA)
301           Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + EUEG*dgdGA
302           Mmat(n,1)=  Mmat(n,1) + EUEG*dgdTA
303          else
304            Amat(n,1) = Amat(n,1) + (EUEGPA*gcab + EUEG*dgdPA)
305            Amat(n,2) = Amat(n,2) + (EUEGPB*gcab + EUEG*dgdPB)
306            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + EUEG*dgdGA
307            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + EUEG*dgdGB
308            Mmat(n,1)=  Mmat(n,1) + EUEG*dgdTA
309            Mmat(n,2)=  Mmat(n,2) + EUEG*dgdTB
310          endif
311      endIf
312c      write (*,*) "Amat(n,1),Cmat(n,1),Mmat(n,1)",Amat(n,1),Cmat(n,1)
313c     & ,Mmat(n,1),ipol
314c      stop
31520    continue
316      end
317
318      Subroutine xc_cvs98_d2()
319      call errquit(' cvs98: d2 not coded ',0,0)
320      return
321      end
322
323      Subroutine vs98ss(tol_rho,PX,GX,TX,F,FP,FG,FT,EUEG,Z,Chi,EUEGP,
324     &                   ChiP,ChiG,ZP,ZT,ijzy)
325      Implicit none
326c
327#include "errquit.fh"
328C
329C     Compute the same-spin part of the vs98 correlation functional for one grid
330C     point and one spin-case.
331C
332
333      integer ijzy
334      double precision tol_rho
335      double precision r13, r14, r15, r16, r17, r18
336      double precision PX, GX, TX, F, FP, FG, FT, DTol, Z, ZP, ZT
337      double precision EUEG, Chi, EUEGP, ChiP, ChiG, cf, gcc
338      double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11
339      double precision Pi, Pi34, F13, F23, F43, F53, F83, F113
340      double precision RS, D, RSP, PotLC, DX, DZ, dgdP, dgdG, dgdT
341      double precision E,DP, DG, DT, rhoo, rho43, rho53, rho83
342      double precision rrho, F4o3, rho13, kc, xk, zk, gc, dgdx, dgdz
343      double precision d2LdSS, d2LdSZ, d2LdZZ, dLdS, dLdZ
344
345      Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/,
346     $  F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/,
347     $  gcc/0.00515088d0/,cf/9.115599720d0/
348
349
350      F4o3 = 4.0d0/3.0d0
351C     Parameters for VS98
352      if (ijzy.eq.1) then
353              r13=   3.270912d-01
354              r14=  -3.228915d-02
355              r15=  -2.942406d-02
356              r16=   2.134222d-03
357              r17=  -5.451559d-03
358              r18=   1.577575d-02
359C     Parameters for M06-L
360      elseif (ijzy.eq.2) then
361              r13=   4.650534D-01
362              r14=   1.617589D-01
363              r15=   1.833657D-01
364              r16=   4.692100D-04
365              r17=  -4.990573D-03
366              r18=   0.000000D+00
367C     Parameters for M06-HF
368      elseif (ijzy.eq.3) then
369              r13=   8.976746D-01
370              r14=  -2.345830D-01
371              r15=   2.368173D-01
372              r16=  -9.913890D-04
373              r17=  -1.146165D-02
374              r18=   0.000000D+00
375C     Parameters for M06
376      elseif (ijzy.eq.4) then
377               r13=  4.905945D-01
378               r14= -1.437348D-01
379               r15=  2.357824D-01
380               r16=  1.871015D-03
381               r17= -3.788963D-03
382               r18=  0.000000D+00
383C     Parameters for M06-2X
384      elseif (ijzy.eq.5) then
385              r13=  6.902145D-01
386              r14=  9.847204D-02
387              r15=  2.214797D-01
388              r16= -1.968264D-03
389              r17= -6.775479D-03
390              r18=  0.000000D+00
391C     Parameters for revM06-L
392      elseif (ijzy.eq.6) then
393              r13=  -5.38821292D-01
394              r14=  -2.829603D-02
395              r15=  2.3889696D-02
396              r16=  0.000000D+00
397              r17=  0.000000D+00
398              r18=  -2.437902D-03
399C     Parameters for revM06
400      elseif (ijzy.eq.7) then
401              r13=  -1.467095900D-01
402              r14=  -1.832187007D-04
403              r15=  8.484372430D-02
404              r16=  0.000000D+00
405              r17=  0.000000D+00
406              r18=  2.280677172D-04
407      else
408        call errquit("vs98ss: illegal value of ijzy",ijzy,UERR)
409      endif
410
411      dtol=tol_rho
412      If(PX.le.DTol.or.gx.le.dtol.or.tx.le.dtol) then
413        EUEG = Zero
414        Chi = Zero
415        EUEGP = Zero
416        ChiP = Zero
417        ChiG = Zero
418        PX = Zero
419        GX = Zero
420        TX = Zero
421        F  = Zero
422        FP = Zero
423        FG = Zero
424        FT = Zero
425        Z  = Zero
426        ZP = Zero
427        ZT = Zero
428      else
429        Pi = F4*ATan(F1)
430        Pi34 = F3 / (F4*Pi)
431        F13 = F1 / F3
432        F23 = F2 / F3
433        F43 = F2 * F23
434        F53 = F5 / F3
435        F83 = F8 / F3
436        F113 = F11 / F3
437        rhoo = PX
438        rrho = 1.0d0/rhoo
439        rho43 = rhoo**F4o3
440        rho13 = rho43*rrho
441        rho53 = rhoo**F53
442        rho83 = rho53*rhoo
443
444        RS = (Pi34/PX) ** F13
445        Call lsdac(tol_rho,
446     A       RS,F1,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
447        EUEG = PX*PotLC
448        Chi = GX/rho83
449        Z = (TX/rho53) - cf
450        kc = F1 + gcc*(Chi + Z)
451        xk = Chi/kc
452        zk = Z/kc
453        D = F1 - Chi/(F4*(Z + cf))
454        call gvt4(gc,dgdx,dgdz,xk,zk,kc,gcc,gcc,r13,r14,r15,r16,r17,r18)
455        E = D*EUEG*gc
456c         write (*,*) "Chi, Z, gc", CHi, Z, gc
457        F = E
458c
459        RSP = -RS/(F3*Px)
460        ChiG = F1/PX**F83
461        ChiP = -F83*Chi/PX
462        ZP = -F53 * TX/rho83
463        ZT =  F1/rho53
464        DZ = Chi/(F4*(Z + cf)*(Z + cf))
465        DX = -F1/(F4*(Z + cf))
466        DP = DZ*ZP + DX*ChiP
467        DG = DX*ChiG
468        DT = DZ*ZT
469        dgdP = dgdx*ChiP + dgdz*ZP
470        dgdG = dgdx*ChiG
471        dgdT = dgdz*ZT
472        EUEGP = PotLC + PX*dLdS*RSP
473        FP = DP*EUEG*gc + D*EUEGP*gc + D*EUEG*dgdP
474        FG = DG*EUEG*gc + D*EUEG*dgdG
475        FT = DT*EUEG*gc + D*EUEG*dgdT
476       Endif
477       Return
478       End
479
480
481