1c Copyright 2018 (C) Orbital-free DFT group at University of Florida
2c Licensed under the Educational Community License, Version 2.0
3c (the "License"); you may not use this file except in compliance with
4c the License. You may obtain a copy of the License at
5c
6c    http://www.osedu.org/licenses/ECL-2.0
7c
8c Unless required by applicable law or agreed to in writing,
9c software distributed under the License is distributed on an "AS IS"
10c BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
11c or implied. See the License for the specific language governing
12c permissions and limitations under the License.
13c
14c ---------------------------------------------------------------------
15c
16c     Strongly constrained and appropriately normed (SCAN)
17c     functional (Correlation part only)
18c           META GGA
19C         utilizes ingredients:
20c                              rho   -  density
21c                              delrho - gradient of density
22c                              tau - K.S kinetic energy density
23c
24c     Written by:
25c     Daniel Mejia-Rodriguez
26c     QTP, Department of Physics, University of Florida
27c
28c     References:
29c     J. Sun, A. Ruzsinszky, J.P. Perdew
30c     PRL 115, 036402 (2015)
31c     DOI: 10.1103/PhysRevLett.115036402
32c
33c     A.P. Bartok and J.R. Yates
34c     JCP 150, 161101 (2019)
35c     DOI: 10.1063/1.5094646
36
37      Subroutine xc_cscan(whichfc, tol_rho, cfac, rho, delrho, Amat,
38     &                    Cmat, nq, ipol, Ec, qwght, ldew, func,
39     &                    tau, Mmat)
40c
41      implicit none
42c
43#include "errquit.fh"
44#include "dft2drv.fh"
45c
46c     Input and other parameters
47c
48      character(*) whichfc
49      integer ipol, nq
50      double precision dummy(1)
51      double precision cfac
52      logical ldew
53      double precision func(*)
54      double precision fac
55      double precision tol_rho
56c
57c     Threshold parameters
58c
59      double precision thr1,thr2
60      parameter (thr1=0.996d0,thr2=1.004d0)
61c
62c     Correlation energy
63c
64      double precision Ec
65c
66c     Charge Density
67c
68      double precision rho(nq,ipol*(ipol+1)/2)
69c
70c     Charge Density Gradient
71c
72      double precision delrho(nq,3,ipol), gammaval, gam12
73c
74c     Kinetic Energy Density
75c
76      double precision tau(nq,ipol), tauN
77c
78c     Quadrature Weights
79c
80      double precision qwght(nq)
81c
82c     Sampling Matrices for the XC Potential
83c
84      double precision Amat(nq,ipol), Cmat(nq,*)
85      double precision Mmat(nq,*)
86c
87c     Intermediate derivatives, results, etc.
88c
89      integer n, ifc
90      double precision ntot,n13,n83,tautot
91      double precision dn2,p,rs,rs12,drsdn,dpdg,dpdn
92      double precision epsc,depscdrs,depscdzeta
93      double precision zeta,omz,opz
94      double precision phi,phi2,phi3,dphidzeta
95      double precision t2,dt2dp,dt2drs,dt2dzeta
96      double precision BETA,dBETAdrs
97      double precision A,dAdrs,dAdzeta
98      double precision At2,Gaux,Gat2
99      double precision w1fac,expw1,w1,dw1drs,dw1dzeta
100      double precision arg1,darg1dp,darg1drs,darg1dzeta
101      double precision H1,dH1dp,dH1drs,dH1dzeta
102      double precision Ec1,dEc1drs,dEc1dp,dEc1dzeta,dEc1dn,dEc1dg
103      double precision epsc0,depsc0drs,depsc0dn
104      double precision dx,ddxdzeta
105      double precision gc,dgcdzeta
106      double precision ginf,dginfdp
107      double precision w0fac,expw0,w0,dw0drs
108      double precision arg0,darg0dp,darg0drs
109      double precision H0,dH0dp,dH0drs
110      double precision Ec0,dEc0dp,dEc0drs,dEc0dzeta,dEc0dn,dEc0dg
111      double precision ds,ddsdzeta
112      double precision tueg,tvw
113      double precision alpha,dalphadzeta,dalphadn,dalphadg,dalphadt
114      double precision oma,oma2
115      double precision fca,dfcada,dfcadg,dfcadzeta,dfcadn
116      double precision exp5,exp6
117      double precision vcpol,vcn
118
119      double precision GAMMA,BETAzero,pi,p14a,p14b,p14f,ckf,ckf2
120      parameter (GAMMA = 0.03109069086965489503494086371273d0)
121      parameter (BETAzero = 0.06672455060314922d0)
122      parameter (p14a=0.1d0,p14b=0.1778d0)
123
124      double precision b1c,b2c,b3c,c1c,c2c,dc,dxc,xi
125      parameter (b1c=0.0285764d0,b2c=0.0889d0,b3c=0.125541d0)
126      parameter (c1c=0.64d0,c2c=1.5d0,dc=0.7d0,dxc=2.3631d0)
127      parameter (xi=0.12802585262625815d0)
128
129      double precision F4,F13,F23,F43,F53,F83,F14,F8,F5,F16
130      parameter (F4=4d0,F13=1d0/3d0,F23=F13+F13,F43=F13+1d0)
131      parameter (F53=F23+1d0,F83=F53+1d0,F14=0.25d0)
132      parameter (F8=8d0,F5=5d0,F16=1d0/6d0)
133
134      double precision rRegu(0:7), rRegu1(7), rtemp(0:7)
135      double precision regalpha, dregalpha
136      parameter ( rRegu = (/ 1d0, -0.64d0, -0.4352d0,
137     &            -1.535685604549d0, 3.061560252175d0,
138     &            -1.915710236206d0, 5.16884468372d-1,
139     &            -5.1848879792d-2 /) )
140      parameter ( rRegu1 = (/ -0.64d0, -2d0*0.4352d0,
141     &            -3d0*1.535685604549d0, 4d0*3.061560252175d0,
142     &            -5d0*1.915710236206d0, 6d0*5.16884468372d-1,
143     &            -7d0*5.1848879792d-2 /) )
144c
145c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
146c
147      Pi = dacos(-1d0)
148      p14f = (3d0/(4d0*Pi))**F13
149      ckf = (3d0*Pi*Pi)**F13
150      ckf2 = ckf*ckf
151
152      do 20 n = 1, nq
153c
154         ntot=rho(n,1)
155         if (ntot.le.tol_rho) goto 20
156
157         n13=ntot**F13
158         n83=ntot**F83
159c
160         if (ipol.eq.1) then
161           zeta = 0d0
162         else
163           zeta = (rho(n,2) - rho(n,3))/ntot
164           if (zeta.lt.-1d0) zeta=-1d0
165           if (zeta.gt. 1d0) zeta= 1d0
166         endif
167c
168         if (ipol.eq.1) then
169           dn2=delrho(n,1,1)**2 + delrho(n,2,1)**2 + delrho(n,3,1)**2
170         else
171           dn2=(delrho(n,1,1)+delrho(n,1,2))**2 +
172     &         (delrho(n,2,1)+delrho(n,2,2))**2 +
173     &         (delrho(n,3,1)+delrho(n,3,2))**2
174         end if
175
176         if (ipol.eq.1) then
177           tautot=tau(n,1)
178         else
179           tautot=tau(n,1)+tau(n,2)
180         endif
181
182         p=dn2/(F4*ckf2*n83)
183         dpdg = 1d0/(F4*ckf2*n83)
184         dpdn = -F83*p/ntot
185
186         rs=p14f/n13
187         drsdn=-F13*rs/ntot
188         rs12=dsqrt(rs)
189c
190         call lsdac(tol_rho,rs,zeta,epsc,depscdrs,depscdzeta,dummy,
191     &              dummy,dummy)
192c
193         opz = 1d0 + zeta
194         omz = 1d0 - zeta
195         phi = 0.5d0*(opz**F23 + omz**F23)
196         phi2 = phi*phi
197         phi3 = phi2*phi
198c
199         BETA = BETAzero*(1d0 + p14a*rs)/(1d0 + p14b*rs)
200         dBETAdrs = BETAzero*(p14a - p14b)/(1d0 + p14b*rs)**2
201c
202         w1fac = epsc/(GAMMA*phi3)
203         expw1 = dexp(-w1fac)
204         w1 = expw1 - 1d0
205         dw1drs = -expw1*depscdrs/(GAMMA*phi3)
206
207         A = BETA/GAMMA/w1
208         dAdrs = dBETAdrs/GAMMA/w1 - A*dw1drs/w1
209
210         t2 = ckf2*p/(4d0*phi2*rs*2d0**F23)
211         dt2dp = ckf2/(4d0*phi2*rs*2d0**F23)
212         dt2drs = -t2/rs
213
214         At2 = A*t2
215         Gaux = 1d0 + 4d0*At2
216         GAt2 = 1d0/dsqrt(dsqrt(Gaux))
217
218         arg1 = 1d0 + w1*(1d0 - GAt2)
219         darg1drs = dw1drs*(1d0 - GAt2) +
220     &              Gat2*w1*(t2*dAdrs + A*dt2drs)/Gaux
221         darg1dp = Gat2*w1*A*dt2dp/Gaux
222
223
224         H1 = GAMMA*phi3*dlog(arg1)
225         dH1dp = GAMMA*phi3*darg1dp/arg1
226         dH1drs = GAMMA*phi3*darg1drs/arg1
227
228         Ec1 = epsc + H1
229         dEc1drs = depscdrs + dH1drs
230         dEc1dp = dH1dp
231c
232c        --------------------------------------------------------------
233c
234         epsc0 = -b1c/(1d0 + b2c*rs12 + b3c*rs)
235         depsc0drs = b1c*(b3c + 0.5d0*b2c/rs12)/
236     &               (1d0 + b2c*rs12 + b3c*rs)**2
237
238         dx = 0.5d0*(opz**F43 + omz**F43)
239         gc = (1d0 - dxc*(dx - 1d0))*(1d0 - zeta**12)
240
241         w0fac = epsc0/b1c
242         expw0 = dexp(-w0fac)
243         w0 = expw0 - 1d0
244         dw0drs = -depsc0drs*expw0/b1c
245
246         ginf = (1d0/(1d0 + 4d0*xi*p))**F14
247         dginfdp = -xi*ginf/(1d0 + 4d0*xi*p)
248
249         arg0 = 1d0 + w0*(1d0 - ginf)
250         darg0drs = dw0drs*(1d0 - ginf)
251         darg0dp = -w0*dginfdp
252
253         H0 = b1c*dlog(arg0)
254         dH0dp = b1c*darg0dp/arg0
255         dH0drs = b1c*darg0drs/arg0
256
257         Ec0 = (epsc0 + H0)*gc
258         dEc0drs = gc*(depsc0drs + dH0drs)
259         dEc0dp = gc*dH0dp
260c
261c        --------------------------------------------------------------
262c
263         ds = 0.5d0*(opz**F53 + omz**F53)
264         tueg = 0.3d0*ckf*ckf*ds*ntot**F53
265         tvw = 0.125d0*dn2/ntot
266         if (whichfc.eq.'orig') then
267           alpha = (tautot - tvw)/tueg
268           if (alpha.lt.0d0) alpha=0d0
269           regalpha = alpha
270           dregalpha = 1d0
271         else if (whichfc.eq.'regu') then
272           alpha = (tautot - tvw)/(tueg + 1d-4*ds)
273           if (alpha.lt.0d0) alpha=0d0
274           regalpha = alpha**3/(alpha**2 + 1d-3)
275           dregalpha = alpha/(alpha**2 + 1d-3) *
276     &                 (3d0*alpha - 2d0*regalpha)
277         endif
278         oma = 1d0 - regalpha
279         oma2 = oma*oma
280
281         if (whichfc.eq.'orig') then
282           if (regalpha.ge.thr1) then
283             exp5 = 0d0
284           else
285             exp5 = dexp(-c1c*regalpha/oma)
286           end if
287           if (regalpha.le.thr2) then
288             exp6 = 0d0
289           else
290             exp6 = dexp(c2c/oma)
291           end if
292           fca = exp5 - dc*exp6
293           if (regalpha.ge.thr1.and.regalpha.le.thr2) then
294             dfcada = 0d0
295           else
296             dfcada = -(c1c*exp5 + dc*exp6*c2c)/oma2
297           end if
298         else if(whichfc.eq.'regu') then
299           if (regalpha.lt.2.5d0) then
300             rtemp(0) = 1d0
301             do ifc=1,7
302               rtemp(ifc) = rtemp(ifc-1)*regalpha
303             enddo
304             fca = dot_product(rRegu,rtemp)
305             dfcada = dot_product(rRegu1, rtemp(0:6))*dregalpha
306           else
307             exp6 = dexp(c2c/oma)
308             fca = -dc*exp6
309             dfcada = -dc*exp6*c2c/oma2*dregalpha
310           endif
311         endif
312c
313c        --------------------------------------------------------------
314c
315         Ec = Ec + cfac*ntot*(Ec1 + fca*(Ec0-Ec1))*qwght(n)
316         if (ldew) func(n) = func(n) + cfac*ntot*(Ec1 + fca*(Ec0-Ec1))
317c
318         if (whichfc.eq.'orig') then
319           dalphadn = f53*(p/ds-alpha)/ntot
320           dalphadg = -0.125d0/(tueg*ntot)
321           dalphadt = 1d0/tueg
322         else if (whichfc.eq.'regu') then
323           dalphadn = 5d0/3d0*(p/ds-alpha)*tueg/((tueg+1d-4*ds)*ntot)
324           dalphadg = -0.125d0/((tueg+1d-4*ds)*ntot)
325           dalphadt = 1d0/(tueg+1d-4*ds)
326         endif
327
328         dEc1dn = dEc1drs*drsdn + dEc1dp*dpdn
329         dEc0dn = dEc0drs*drsdn + dEc0dp*dpdn
330         dfcadn = dfcada*dalphadn
331         vcn = Ec1 + fca*(Ec0-Ec1) +
332     &         ntot*(dEc1dn + fca*(dEc0dn-dEc1dn) + dfcadn*(Ec0-Ec1))
333
334         Amat(n,1) = Amat(n,1) + cfac*vcn
335
336         dEc1dg = dEc1dp*dpdg
337         dEc0dg = dEc0dp*dpdg
338         dfcadg = dfcada*dalphadg
339
340         Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + cfac*ntot*(dEc1dg +
341     &   fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1))
342         Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + cfac*ntot*(dEc1dg +
343     &   fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1))*2d0
344
345         Mmat(n,1) = Mmat(n,1) + 0.5d0*cfac*
346     &                           ntot*dfcada*dalphadt*(Ec0-Ec1)
347
348         if (ipol.eq.2) then
349           Amat(n,2) = Amat(n,2) + cfac*vcn
350
351           Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + cfac*ntot*(dEc1dg +
352     &     fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1))
353
354           Mmat(n,2) = Mmat(n,2) + 0.5d0*cfac*
355     &                             ntot*dfcada*dalphadt*(Ec0-Ec1)
356
357           if (omz.lt.tol_rho) then
358             dphidzeta = 0.5d0*F23*(opz**F23/opz)
359           else if (opz.lt.tol_rho) then
360             dphidzeta = -0.5d0*F23*(omz**F23/omz)
361           else
362             dphidzeta = 0.5d0*F23*(opz**F23/opz - omz**F23/omz)
363           end if
364
365           dt2dzeta = -2d0*t2*dphidzeta/phi
366           dw1dzeta = (3d0*w1fac*dphidzeta/phi -
367     &                 depscdzeta/(GAMMA*phi3))*expw1
368           dAdzeta = -A*dw1dzeta/w1
369           darg1dzeta = dw1dzeta*(1d0 - Gat2) +
370     &                  Gat2*w1*(t2*dAdzeta + A*dt2dzeta)/Gaux
371           dH1dzeta = 3d0*H1*dphidzeta/phi + GAMMA*phi3*darg1dzeta/arg1
372           dEc1dzeta = depscdzeta + dH1dzeta
373
374           ddxdzeta = 0.5d0*F43*(opz**F13 - omz**F13)
375           dgcdzeta = -dxc*ddxdzeta*(1d0 - zeta**12) -
376     &                12d0*zeta**11*(1d0 - dxc*(dx - 1d0))
377           dEc0dzeta = dgcdzeta*(epsc0 + H0)
378
379           ddsdzeta = 0.5d0*F53*(opz**F23 - omz**F23)
380
381           dalphadzeta = -alpha*ddsdzeta/ds
382
383           dfcadzeta = dfcada*dalphadzeta
384
385           vcpol = dEc1dzeta + dfcadzeta*(Ec0-Ec1) +
386     &             fca*(dEc0dzeta-dEc1dzeta)
387
388           Amat(n,1) = Amat(n,1) + cfac*omz*vcpol
389           Amat(n,2) = Amat(n,2) - cfac*opz*vcpol
390
391         end if
392
39320    continue
394      end
395
396      Subroutine xc_cscan_d2()
397      implicit none
398      call errquit(' xc_cscan: d2 not coded ',0,0)
399      return
400      end
401