1c    M06 suite correlation functional
2c           META GGA
3C         utilizes ingredients:
4c                              rho   -  density
5c                              delrho - gradient of density
6c                              tau (tauN)- K.S kinetic energy density
7c                              ijzy - 1  M06-L
8c                              ijzy - 2  M06-HF
9c                              ijzy - 3  M06
10c                              ijzy - 4  M06-2X
11c                              ijzy - 5 revM06-L
12c                              ijzy - 6 revM06
13
14      Subroutine xc_cm06(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
15     &                     nq, ipol, Ec, qwght, ldew, func,
16     &                     tau, Amat, Cmat, Mmat,ijzy)
17
18
19c
20c$Id$
21c
22c
23c     [a]   Zhao, Y. and  Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101;
24c     [b]   Zhao, Y. and  Truhlar, D. G. J. Phys. Chem. A (2006),110(49),13126-13130.
25c     [c]   Grafenstein, J., Izotov, D. and Cremer, D. J. Chem. Phys. 2007, 127, 214103.
26
27
28      implicit none
29c
30#include "errquit.fh"
31c
32c
33c     Input and other parameters
34c
35      integer ipol, nq
36
37      double precision cfac
38      logical lcfac, nlcfac
39
40      logical lfac, nlfac
41      double precision fac
42      double precision tol_rho
43
44c
45c     Threshold parameters
46c
47      double precision F1, F2, F3, F4,COpp
48      Data COpp/0.0031d0/,F1/1.0d0/,F2/2.0d0/,
49     & F3/3.0d0/,F4/4.0d0/
50c
51c     Correlation energy
52c
53      double precision Ec
54c
55c     Charge Density
56c
57      double precision rho(nq,ipol*(ipol+1)/2)
58c
59c     Charge Density Gradient
60c
61      double precision delrho(nq,3,ipol), gammaval, gam12
62
63c
64c     Kinetic Energy Density
65c
66      double precision tau(nq,ipol), tauN
67
68c
69c     Quadrature Weights
70c
71      double precision qwght(nq)
72c
73      logical ldew
74      double precision func(*)
75c
76c     Sampling Matrices for the XC Potential
77c
78      double precision Amat(nq,ipol), Cmat(nq,*)
79      double precision Mmat(nq,*)
80
81      integer n, ijzy
82
83c    call to the m05css subroutine
84      double precision PA,GAA,TA,FA,FPA,FGA,FTA,EUA,EUEGA,ChiA,EUPA
85     &,ChiAP,ChiAG
86      double precision PB,GBB,TB,FB,FPB,FGB,FTB,EUB,EUEGB,ChiB,EUPB
87     &,ChiBP,ChiBG
88c
89      double precision  sop, sopp0, sopp1,sopp2, sopp3, sopp4
90      double precision Pi, F6, F43, Pi34, F13,
91     &RS,RSP,Zeta,dZdA,dZdB,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ
92      double precision P, EUEG, U, W
93      double precision dUdChiA,dUdChiB,dUdPA,dUdPB,dUdGA,dUdGB,
94     &dWdU,dWdPA,dWdPB, dWdGA,dWdGB,EUEGPA,EUEGPB
95
96c
97c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
98c
99      sop=1.0d0
100      sopp0= 0d0
101      sopp1= 0d0
102      sopp2= 0d0
103      sopp3= 0d0
104      sopp4= 0d0
105      if (ijzy.eq.1) then
106C     Parameters for M06-L Correlation
107         sopp0= 6.042374D-01
108         sopp1= 1.776783D+02
109         sopp2= -2.513252D+02
110         sopp3= 7.635173D+01
111         sopp4= -1.255699D+01
112      elseif (ijzy.eq.2) then
113C     Parameters for M06-HF Correlation
114         sopp0= 1.674634D+00
115         sopp1= 5.732017D+01
116         sopp2= 5.955416D+01
117         sopp3= -2.311007D+02
118         sopp4= 1.255199D+02
119      elseif (ijzy.eq.3) then
120C     Parameters for M06 Correlation
121         sopp0= 3.741539D+00
122         sopp1= 2.187098D+02
123         sopp2= -4.531252D+02
124         sopp3= 2.936479D+02
125         sopp4= -6.287470D+01
126      elseif (ijzy.eq.4) then
127C     Parameters for M06-2X Correlation
128         sopp0= 8.833596D-01
129         sopp1= 3.357972D+01
130         sopp2= -7.043548D+01
131         sopp3= 4.978271D+01
132         sopp4= -1.852891D+01
133      elseif (ijzy.eq.5) then
134C     Parameters for revM06-L Correlation
135         sopp0= 3.44360696D-01
136         sopp1= -5.57080242D-01
137         sopp2= -2.009821162D+00
138         sopp3= -1.857641887D+00
139         sopp4= -1.076639864D+00
140      elseif (ijzy.eq.6) then
141C     Parameters for revM06 Correlation
142         sopp0= 1.222401598D+00
143         sopp1= 6.613907336D-01
144         sopp2= -1.884581043D+00
145         sopp3= -2.780360568D+00
146         sopp4= -3.068579344D+00
147      else
148         call errquit("xc_cm06: illegal value of ijzy",ijzy,UERR)
149      endif
150
151      call xc_cvs98(tol_rho, 1.0d0, lfac, nlcfac,
152     &           rho, delrho,  nq, ipol,
153     &           Ec, qwght, ldew,func,tau,Amat,Cmat,Mmat,ijzy+1)
154
155
156
157      Pi = F4*ATan(F1)
158      F6=6.0d0
159      F43 = F4 / F3
160      Pi34 = F3 / (F4*Pi)
161      F13 = F1 / F3
162
163      do 20 n = 1, nq
164       if (rho(n,1).lt.Tol_Rho) goto 20
165       if (ipol.eq.1) then
166c
167c    get the density, gradient, and tau for the alpha spin from the total
168c
169         PA = rho(n,1)/F2
170         GAA = (    delrho(n,1,1)*delrho(n,1,1) +
171     &                 delrho(n,2,1)*delrho(n,2,1) +
172     &                 delrho(n,3,1)*delrho(n,3,1))/F4
173         if(sqrt(gaa).lt.tol_rho) goto 20
174c  In the m05css subroutine, we use 2*TA as the tau, so we do not divide
175c  the tau by 2 here
176
177         TA = tau(n,1)
178         if(ta.lt.tol_rho) goto 20
179
180         Call m06css(Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
181     &                ChiA,EUPA,ChiAP,ChiAG,ijzy)
182         PB = PA
183         GBB = GAA
184         TB = TA
185         FB = FA
186         FPB = FPA
187         FGB = FGA
188         FTB = FTA
189         EUB = EUA
190         ChiB = ChiA
191         EUPB = EUPA
192         ChiBP = ChiAP
193         ChiBG = ChiAG
194
195         Ec = Ec + 2.d0*FA*qwght(n)            !factor of 2 account for both spin
196         if(ldew) func(n)=func(n)+ FA*2d0
197         Amat(n,1)=Amat(n,1)+ FPA
198         Cmat(n,1)=  Cmat(n,1) + FGA
199         Mmat(n,1)=  Mmat(n,1) + FTA
200c         write (*,*) "PA,GAA,TA",PA,GAA,TA
201c         write (*,*) "FPA,FGA,FTA",FPA,FGA,FTA
202c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
203      else  ! ipol=2
204c
205c        ======> SPIN-UNRESTRICTED <======
206c
207         ChiA  = 0.0d0
208         ChiAP = 0.0d0
209         ChiAG = 0.0d0
210         ChiB  = 0.0d0
211         ChiBP = 0.0d0
212         ChiBG = 0.0d0
213         EUA   = 0.0d0
214         EUB   = 0.0d0
215         EUPA  = 0.0d0
216         EUPB  = 0.0d0
217c
218c       alpha
219c
220
221         PA = rho(n,2)
222         if (PA.le.Tol_Rho) go to 25
223         GAA =   delrho(n,1,1)*delrho(n,1,1) +
224     &           delrho(n,2,1)*delrho(n,2,1) +
225     &          delrho(n,3,1)*delrho(n,3,1)
226c
227c  In the m05css subroutine, we use 2*TA as the tau
228c
229         TA = 2*tau(n,1)
230
231         if(sqrt(gaa).lt.tol_rho) goto 25
232         if(ta.lt.tol_rho) goto 25
233
234         Call m06css(Tol_Rho,PA,GAA,TA,FA,FPA,FGA,FTA,EUA,
235     &                ChiA,EUPA,ChiAP,ChiAG,ijzy)
236         Ec = Ec + FA*qwght(n)
237         if(ldew) func(n)=func(n)+ FA
238         Amat(n,1)=Amat(n,1)+ FPA
239         Cmat(n,1)=  Cmat(n,1) + FGA
240         Mmat(n,1)=  Mmat(n,1) + FTA
241c
242c  In the m05css subroutine, we use 2*TB as the tau,
243c
244c
245c       Beta
246c
247 25      continue
248         PB = rho(n,3)
249         if (PB.le.Tol_Rho) go to 30
250         GBB =   delrho(n,1,2)*delrho(n,1,2) +
251     &           delrho(n,2,2)*delrho(n,2,2) +
252     &          delrho(n,3,2)*delrho(n,3,2)
253
254         TB = 2*tau(n,2)
255
256         if(sqrt(gbb).lt.tol_rho) goto 30
257         if(tb.lt.tol_rho) goto 30
258
259         Call m06css(Tol_Rho,PB,GBB,TB,FB,FPB,FGB,FTB,EUB,
260     &                ChiB,EUPB,ChiBP,ChiBG,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,3)=  Cmat(n,3) + FGB
265         Mmat(n,2)=  Mmat(n,2) + FTB
266      endif
267
268 30   continue
269      P = PA + PB
270
271      If((PA.gt.Tol_Rho).and.(PB.gt.Tol_Rho)) then
272          RS = (Pi34/P) ** F13
273          RSP = -RS/(F3*P)
274          Zeta = (PA-PB)/P
275          dZdA = (F1-Zeta)/P
276          dZdB = (-F1-Zeta)/P
277          Call lsdac(tol_rho,
278     R         RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,
279     $      d2LdZZ)
280          EUEG = P*PotLC - EUA - EUB
281          U = COpp*(ChiA+ChiB)/(F1 + COpp*(ChiA+ChiB))
282          W = sopp0+U*(sopp1+U*(sopp2+U*(sopp3+U*sopp4)))
283          Ec = Ec + sop*EUEG*W*qwght(n)
284          if(ldew) func(n)=func(n)+ sop*EUEG*W
285          dUdChiA =COpp/(F1 + COpp*(ChiA+ChiB))**2
286          dUdChiB =COpp/(F1 + COpp*(ChiA+ChiB))**2
287          dUdPA= dUdChiA*ChiAP
288          dUdPB= dUdChiB*ChiBP
289          dUdGA= dUdChiA*ChiAG
290          dUdGB= dUdChiB*ChiBG
291          dWdU =sopp1+U*(F2*sopp2+U*(F3*sopp3+U*F4*sopp4))
292          dWdPA= dWdU*dUdPA
293          dWdPB= dWdU*dUdPB
294          dWdGA= dWdU*dUdGA
295          dWdGB= dWdU*dUdGB
296          EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA
297          EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB
298          if (ipol.eq.1) then
299            Amat(n,1) = Amat(n,1) + sop*(EUEGPA*W + EUEG*dWdPA)
300            Cmat(n,1)=  Cmat(n,1) + sop*(EUEG*dWdGA)
301          else
302            Amat(n,1) = Amat(n,1) + sop*(EUEGPA*W + EUEG*dWdPA)
303            Amat(n,2) = Amat(n,2) + sop*(EUEGPB*W + EUEG*dWdPB)
304            Cmat(n,1) = Cmat(n,1) + sop*EUEG*dWdGA
305            Cmat(n,3) = Cmat(n,3) + sop*(EUEG*dWdGB)
306          endif
307      endIf
308c      write (*,*) "PA, PB, GAA, GBB,ipol",PA, PB, GAA, GBB,ipol
309c      write (*,*) "FA, FB,FGA, FGB",FA, FB,FGA, FGB
310c      Stop
31120    continue
312      end
313
314      Subroutine xc_cm06_d2()
315      implicit none
316      call errquit(' xc06: d2 not coded ',0,0)
317      return
318      end
319
320
321
322
323      Subroutine m06css(Tol_Rho,PX,GX,TX,F,FP,FG,FT,EUEG,Chi,EUEGP,
324     &                   ChiP,ChiG,ijzy)
325      Implicit none
326c
327#include "errquit.fh"
328C
329C     Compute the same-spin part of the m05 correlation functional for one grid
330C     point and one spin-case.
331C
332C
333      integer ijzy
334      double precision PX, GX, TX, F, FP, FG, FT, Tol_Rho
335      double precision EUEG, Chi, EUEGP, ChiP, ChiG
336      double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11
337      double precision ss, sss0,sss1, sss2, sss3, sss4, Css
338      double precision Pi, Pi34, F13, F23, F43, F53, F83, F113
339      double precision RS, FDUEG, D, Fscc, RSP, dFsccP, dFsccG
340      double precision E, W, U, dFsccT, dUdChi, dWdU, dWdP, dWdG
341      double precision d2LdSS,d2LdSZ,d2LdZZ,PotLC,dLdS,dLdZ
342
343
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     $  Css/0.06d0/
348C
349c      Tol_Rho=1.0D-7
350c      write (*,*) Tol_Rho
351      sss0=  0d0
352      sss1=  0d0
353      sss2=  0d0
354      sss3=  0d0
355      sss4=  0d0
356      ss=1.0
357      if (ijzy.eq.1) then
358C     Parameters for M06-L Correlation
359         sss0=  5.349466D-01
360         sss1=  5.396620D-01
361         sss2=  -3.161217D+01
362         sss3=  5.149592D+01
363         sss4=  -2.919613D+01
364      elseif (ijzy.eq.2) then
365C     Parameters for M06-HF Correlation
366         sss0=  1.023254D-01
367         sss1=  -2.453783D+00
368         sss2=  2.913180D+01
369         sss3=  -3.494358D+01
370         sss4=  2.315955D+01
371      elseif (ijzy.eq.3) then
372C     Parameters for M06 Correlation
373         sss0=  5.094055D-01
374         sss1=  -1.491085D+00
375         sss2=  1.723922D+01
376         sss3=  -3.859018D+01
377         sss4=  2.845044D+01
378      elseif (ijzy.eq.4) then
379C     Parameters for M06-2X Correlation
380         sss0=  3.097855D-01
381         sss1=  -5.528642D+00
382         sss2=  1.347420D+01
383         sss3=  -3.213623D+01
384         sss4=  2.846742D+01
385      elseif (ijzy.eq.5) then
386C     Parameters for revM06-L Correlation
387         sss0=  1.227659748D+00
388         sss1=  8.55201283D-01
389         sss2=  -3.113346677D+00
390         sss3=  -2.239678026D+00
391         sss4=  3.54638962D-01
392      elseif (ijzy.eq.6) then
393C     Parameters for revM06 Correlation
394         sss0=  9.017224575D-01
395         sss1=  2.079991827D-01
396         sss2=  -1.823747562D+00
397         sss3=  -1.384430429D+00
398         sss4=  -4.423253381D-01
399      else
400         call errquit("m06css: illegal value of ijzy",ijzy,UERR)
401      endif
402
403      If ((PX.le.Tol_Rho))  then
404        EUEG = Zero
405        Chi = Zero
406        EUEGP = Zero
407        ChiP = Zero
408        ChiG = Zero
409        PX = Zero
410        GX = Zero
411        TX = Zero
412        F  = Zero
413        FP = Zero
414        FG = Zero
415        FT = Zero
416      else
417        Pi = F4*ATan(F1)
418        Pi34 = F3 / (F4*Pi)
419        F13 = F1 / F3
420        F23 = F2 / F3
421        F43 = F2 * F23
422        F53 = F5 / F3
423        F83 = F8 / F3
424        F113 = F11 / F3
425        FDUEG = (F3/F5)*(F6*Pi*Pi)**F23
426        RS = (Pi34/PX) ** F13
427        Call lsdac(tol_rho,
428     R       RS,F1,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ)
429        EUEG = PX*PotLC
430        D = TX - Pt25*GX/PX
431C        DUEG = FDUEG*PX**F53
432        Chi = GX/PX**F83
433        U = Css*Chi/(F1 + Css*Chi)
434        W = sss0+U*(sss1+U*(sss2+U*(sss3+U*sss4)))
435        Fscc=D/TX
436        E = Fscc*W*EUEG
437        F = E*ss
438        RSP = -RS/(F3*Px)
439        ChiG = F1/PX**F83
440        ChiP = -F83*Chi/PX
441        dFsccP=Pt25*GX/(TX*PX**2)
442        dFsccG=-Pt25/(TX*PX)
443        dFsccT=Pt25*GX/(PX*TX**2)
444        dUdChi=Css/((F1+Css*Chi)**2)
445        dWdU=sss1+U*(F2*sss2+U*(F3*sss3+U*F4*sss4))
446        dWdP=dWdU*dUdChi*ChiP
447        dWdG=dWdU*dUdChi*ChiG
448        EUEGP = PotLC + PX*dLdS*RSP
449        FP = ss*(dFsccP*W*EUEG
450     $                 + Fscc*dWdP*EUEG
451     $                 + Fscc*W*EUEGP)
452        FG = ss*(dFsccG*W*EUEG
453     $                 + Fscc*dWdG*EUEG)
454
455        FT = ss*(dFsccT*W*EUEG)
456       Endif
457
458       Return
459       End
460
461
462