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     D. Mejia-Rodriguez and S.B. Trickey
30c     Phys. Rev. A 96, 052512 (2017)
31
32      Subroutine xc_cscanl(tol_rho, cfac, rho, delrho, laprho, Amat,
33     &                     Cmat, Lmat, nq, ipol, Ec, qwght, ldew, func)
34c
35      implicit none
36c
37#include "errquit.fh"
38#include "dft2drv.fh"
39c
40c     Input and other parameters
41c
42      integer ipol, nq
43      double precision dummy(1)
44      double precision cfac(*)
45      logical ldew
46      double precision func(*)
47      double precision fac
48      double precision tol_rho
49c
50c     Threshold parameters
51c
52      double precision thr1,thr2
53      parameter (thr1=0.996d0,thr2=1.004d0)
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)
66c
67c     Charge Density Laplacian
68c
69      double precision laprho(nq,ipol)
70c
71c     Quadrature Weights
72c
73      double precision qwght(nq)
74c
75c     Sampling Matrices for the XC Potential
76c
77      double precision Amat(nq,ipol), Cmat(nq,*), Lmat(nq,ipol)
78c
79c     Intermediate derivatives, results, etc.
80c
81      integer n
82      double precision ntot,n13,n83,n53,tautot
83      double precision dn2,p,rs,rs12,drsdn,dpdg,dpdn
84      double precision epsc,depscdrs,depscdzeta
85      double precision zeta,omz,opz
86      double precision phi,phi2,phi3,dphidzeta
87      double precision t2,dt2dp,dt2drs,dt2dzeta
88      double precision BETA,dBETAdrs
89      double precision A,dAdrs,dAdzeta
90      double precision At2,Gaux,Gat2
91      double precision w1fac,expw1,w1,dw1drs,dw1dzeta
92      double precision arg1,darg1dp,darg1drs,darg1dzeta
93      double precision H1,dH1dp,dH1drs,dH1dzeta
94      double precision Ec1,dEc1drs,dEc1dp,dEc1dzeta,dEc1dn,dEc1dg
95      double precision epsc0,depsc0drs,depsc0dn
96      double precision dx,ddxdzeta
97      double precision gc,dgcdzeta
98      double precision ginf,dginfdp
99      double precision w0fac,expw0,w0,dw0drs
100      double precision arg0,darg0dp,darg0drs
101      double precision H0,dH0dp,dH0drs
102      double precision Ec0,dEc0dp,dEc0drs,dEc0dzeta,dEc0dn,dEc0dg
103      double precision ds,ddsdzeta
104      double precision tueg,tvw
105      double precision alpha,dalphadzeta,dalphadn,dalphadg,dalphadt
106      double precision oma,oma2
107      double precision fca,dfcada,dfcadg,dfcadzeta,dfcadn,dfcadt
108      double precision exp5,exp6
109      double precision vcpol,vcn
110      double precision tuega, tuegb, pa, pb, qa, qb
111      double precision fsa, fsb, dfsdpa, dfsdpb, dfsdqa, dfsdqb
112      double precision dtdna, dtdnb, dtdga, dtdgb, dtdla, dtdlb, dtuegdn
113
114      double precision GAMMA,BETAzero,pi,p14a,p14b,p14f,ckf2
115      parameter (GAMMA = 0.03109069086965489503494086371273d0)
116      parameter (BETAzero = 0.06672455060314922d0)
117      parameter (p14a=0.1d0,p14b=0.1778d0)
118
119      double precision b1c,b2c,b3c,c1c,c2c,dc,dxc,xi
120      parameter (b1c=0.0285764d0,b2c=0.0889d0,b3c=0.125541d0)
121      parameter (c1c=0.64d0,c2c=1.5d0,dc=0.7d0,dxc=2.3631d0)
122      parameter (xi=0.12802585262625815d0)
123
124      double precision F4,F13,F23,F43,F53,F83,F14,F8,F5,F16,two53
125      parameter (F4=4d0,F13=1d0/3d0,F23=F13+F13,F43=F13+1d0)
126      parameter (F53=F23+1d0,F83=F53+1d0,F14=0.25d0)
127      parameter (F8=8d0,F5=5d0,F16=1d0/6d0)
128      parameter (two53=2d0**F53)
129c
130c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
131c
132      Pi = dacos(-1d0)
133      p14f = (3d0/(4d0*Pi))**F13
134      ckf2 = (3d0*Pi*Pi)**F23
135
136      do 20 n = 1, nq
137c
138         ntot=rho(n,1)
139         if (ntot.le.tol_rho) goto 20
140
141         n13=ntot**F13
142         n53=ntot**F53
143         n83=ntot**F83
144c
145         if (ipol.eq.1) then
146           zeta = 0d0
147         else
148           zeta = (rho(n,2) - rho(n,3))/ntot
149           if (zeta.lt.-1d0) zeta=-1d0
150           if (zeta.gt. 1d0) zeta= 1d0
151         endif
152c
153         if (ipol.eq.1) then
154           dn2=delrho(n,1,1)**2 + delrho(n,2,1)**2 + delrho(n,3,1)**2
155         else
156           dn2=(delrho(n,1,1)+delrho(n,1,2))**2 +
157     &         (delrho(n,2,1)+delrho(n,2,2))**2 +
158     &         (delrho(n,3,1)+delrho(n,3,2))**2
159         end if
160
161         if (ipol.eq.1) then
162           pa = dn2/(4d0*ckf2*n83)
163           qa = laprho(n,1)/(4d0*ckf2*n53)
164           tuega = 0.3d0*ckf2*n53
165           call ts_pc(tol_rho, rho(n,1), delrho(n,1:3,1), laprho(n,1),
166     &                dfsdpa, dfsdqa, fsa, 1d0)
167           tautot = tuega*fsa
168           dtdna = F53*tuega/rho(n,1)*fsa -
169     &             tuega*(F83*dfsdpa*pa/rho(n,1) +
170     &                    F53*dfsdqa*qa/rho(n,1))
171           dtdga = 2d0*tuega*dfsdpa/(4d0*ckf2*n83)
172           dtdla = tuega*dfsdqa/(4d0*ckf2*n53)
173         else
174           tuega = 0.5d0*0.3d0*two53*ckf2*rho(n,2)**F53
175           tuegb = 0.5d0*0.3d0*two53*ckf2*rho(n,3)**F53
176           call ts_pc(tol_rho, rho(n,2), delrho(n,1:3,1), laprho(n,1),
177     &                dfsdpa, dfsdqa, fsa, 2d0)
178           call ts_pc(tol_rho, rho(n,3), delrho(n,1:3,2), laprho(n,2),
179     &                dfsdpb, dfsdqb, fsb, 2d0)
180           tautot=(tuega*fsa + tuegb*fsb)
181           if (rho(n,2).gt.tol_rho) then
182             pa = (delrho(n,1,1)**2 +
183     &             delrho(n,2,1)**2 +
184     &             delrho(n,3,1)**2)/(ckf2*(2d0*rho(n,2))**F83)
185             qa = laprho(n,1)/(2d0*ckf2*(2d0*rho(n,2))**F53)
186             dtdna = F53*tuega/rho(n,2)*fsa -
187     &               tuega*(F83*dfsdpa*pa/rho(n,2) +
188     &                      F53*dfsdqa*qa/rho(n,2))
189             dtdga = tuega*dfsdpa/(ckf2*(2d0*rho(n,2))**F83)
190             dtdla = tuega*dfsdqa/(2d0*ckf2*(2d0*rho(n,2))**F53)
191           else
192             dtdna = 0d0
193             dtdga = 0d0
194             dtdla = 0d0
195           endif
196           if (rho(n,3).gt.tol_rho) then
197             pb = (delrho(n,1,2)**2 +
198     &             delrho(n,2,2)**2 +
199     &             delrho(n,3,2)**2)/(ckf2*(2d0*rho(n,3))**F83)
200             qb = laprho(n,2)/(2d0*ckf2*(2d0*rho(n,3))**F53)
201             dtdnb = F53*tuegb/rho(n,3)*fsb -
202     &               tuegb*(F83*dfsdpb*pb/rho(n,3) +
203     &                      F53*dfsdqb*qb/rho(n,3))
204             dtdgb = tuegb*dfsdpb/(ckf2*(2d0*rho(n,3))**F83)
205             dtdlb = tuegb*dfsdqb/(2d0*ckf2*(2d0*rho(n,3))**F53)
206           else
207             dtdnb = 0d0
208             dtdgb = 0d0
209             dtdlb = 0d0
210           endif
211         endif
212
213         p=dn2/(F4*ckf2*n83)
214         dpdg = 1d0/(F4*ckf2*n83)
215         dpdn = -F83*p/ntot
216
217         rs=p14f/n13
218         drsdn=-F13*rs/ntot
219         rs12=dsqrt(rs)
220c
221         call lsdac(tol_rho,rs,zeta,epsc,depscdrs,depscdzeta,dummy,
222     &              dummy,dummy)
223c
224         opz = 1d0 + zeta
225         omz = 1d0 - zeta
226         phi = 0.5d0*(opz**F23 + omz**F23)
227         phi2 = phi*phi
228         phi3 = phi2*phi
229c
230         BETA = BETAzero*(1d0 + p14a*rs)/(1d0 + p14b*rs)
231         dBETAdrs = BETAzero*(p14a - p14b)/(1d0 + p14b*rs)**2
232c
233         w1fac = epsc/(GAMMA*phi3)
234         expw1 = dexp(-w1fac)
235         w1 = expw1 - 1d0
236         dw1drs = -expw1*depscdrs/(GAMMA*phi3)
237
238         A = BETA/GAMMA/w1
239         dAdrs = dBETAdrs/GAMMA/w1 - A*dw1drs/w1
240
241         t2 = ckf2*p/(4d0*phi2*rs*2d0**F23)
242         dt2dp = ckf2/(4d0*phi2*rs*2d0**F23)
243         dt2drs = -t2/rs
244
245         At2 = A*t2
246         Gaux = 1d0 + 4d0*At2
247         GAt2 = 1d0/dsqrt(dsqrt(Gaux))
248
249         arg1 = 1d0 + w1*(1d0 - GAt2)
250         darg1drs = dw1drs*(1d0 - GAt2) +
251     &              Gat2*w1*(t2*dAdrs + A*dt2drs)/Gaux
252         darg1dp = Gat2*w1*A*dt2dp/Gaux
253
254
255         H1 = GAMMA*phi3*dlog(arg1)
256         dH1dp = GAMMA*phi3*darg1dp/arg1
257         dH1drs = GAMMA*phi3*darg1drs/arg1
258
259         Ec1 = epsc + H1
260         dEc1drs = depscdrs + dH1drs
261         dEc1dp = dH1dp
262c
263c        --------------------------------------------------------------
264c
265         epsc0 = -b1c/(1d0 + b2c*rs12 + b3c*rs)
266         depsc0drs = b1c*(b3c + 0.5d0*b2c/rs12)/
267     &               (1d0 + b2c*rs12 + b3c*rs)**2
268
269         dx = 0.5d0*(opz**F43 + omz**F43)
270         gc = (1d0 - dxc*(dx - 1d0))*(1d0 - zeta**12)
271
272         w0fac = epsc0/b1c
273         expw0 = dexp(-w0fac)
274         w0 = expw0 - 1d0
275         dw0drs = -depsc0drs*expw0/b1c
276
277         ginf = (1d0/(1d0 + 4d0*xi*p))**F14
278         dginfdp = -xi*ginf/(1d0 + 4d0*xi*p)
279
280         arg0 = 1d0 + w0*(1d0 - ginf)
281         darg0drs = dw0drs*(1d0 - ginf)
282         darg0dp = -w0*dginfdp
283
284         H0 = b1c*dlog(arg0)
285         dH0dp = b1c*darg0dp/arg0
286         dH0drs = b1c*darg0drs/arg0
287
288         Ec0 = (epsc0 + H0)*gc
289         dEc0drs = gc*(depsc0drs + dH0drs)
290         dEc0dp = gc*dH0dp
291c
292c        --------------------------------------------------------------
293c
294         ds = 0.5d0*(opz**F53 + omz**F53)
295         tueg = 0.3d0*ckf2*ds*ntot**F53
296         tvw = 0.125d0*dn2/ntot
297
298         dtuegdn = F53*tueg/ntot
299
300         alpha = (tautot - tvw)/tueg
301         if (alpha.lt.0d0) alpha=0d0
302         oma = 1d0 - alpha
303         oma2 = oma*oma
304
305         if (alpha.ge.thr1) then
306           exp5 = 0d0
307         else
308           exp5 = dexp(-c1c*alpha/oma)
309         end if
310         if (alpha.le.thr2) then
311           exp6 = 0d0
312         else
313           exp6 = dexp(c2c/oma)
314         end if
315         fca = exp5 - dc*exp6
316
317         if (alpha.ge.thr1.and.alpha.le.thr2) then
318           dfcada = 0d0
319         else
320           dfcada = -(c1c*exp5 + dc*exp6*c2c)/oma2
321         end if
322c
323c        --------------------------------------------------------------
324c
325         Ec = Ec + ntot*(Ec1 + fca*(Ec0-Ec1))*qwght(n)
326         if (ldew) func(n) = func(n) + ntot*(Ec1 + fca*(Ec0-Ec1))
327c
328         dalphadn = F13*(F8*tvw - F5*tautot)/(tueg*ntot)
329         dalphadg = -0.125d0/(tueg*ntot)
330         dalphadt = 1d0/tueg
331
332         dEc1dn = dEc1drs*drsdn + dEc1dp*dpdn
333         dEc0dn = dEc0drs*drsdn + dEc0dp*dpdn
334         dfcadn = dfcada*dalphadn
335         dfcadt = dfcada*dalphadt
336
337         vcn = Ec1 + fca*(Ec0-Ec1) +
338     &         ntot*(dEc1dn + fca*(dEc0dn-dEc1dn) + dfcadn*(Ec0-Ec1))
339
340         Amat(n,1) = Amat(n,1) + vcn + ntot*dfcadt*(Ec0-Ec1)*dtdna
341
342         dEc1dg = dEc1dp*dpdg
343         dEc0dg = dEc0dp*dpdg
344         dfcadg = dfcada*dalphadg
345
346         Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + ntot*(dEc1dg +
347     &   fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1) +
348     &   dfcadt*(Ec0-Ec1)*dtdga)
349         Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + ntot*(dEc1dg +
350     &   fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1))*2d0
351
352         Lmat(n,1) = Lmat(n,1) + ntot*dfcadt*dtdla*(Ec0-Ec1)
353
354         if (ipol.eq.2) then
355           Amat(n,2) = Amat(n,2) + vcn + ntot*dfcadt*(Ec0-Ec1)*dtdnb
356
357           Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + ntot*(dEc1dg +
358     &     fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1) +
359     &     dfcadt*(Ec0-Ec1)*dtdgb)
360
361           Lmat(n,2) = Lmat(n,2) + ntot*dfcadt*dtdlb*(Ec0-Ec1)
362
363           if (omz.lt.tol_rho) then
364             dphidzeta = 0.5d0*F23*(opz**F23/opz)
365           else if (opz.lt.tol_rho) then
366             dphidzeta = -0.5d0*F23*(omz**F23/omz)
367           else
368             dphidzeta = 0.5d0*F23*(opz**F23/opz - omz**F23/omz)
369           end if
370
371           dt2dzeta = -2d0*t2*dphidzeta/phi
372           dw1dzeta = (3d0*w1fac*dphidzeta/phi -
373     &                 depscdzeta/(GAMMA*phi3))*expw1
374           dAdzeta = -A*dw1dzeta/w1
375           darg1dzeta = dw1dzeta*(1d0 - Gat2) +
376     &                  Gat2*w1*(t2*dAdzeta + A*dt2dzeta)/Gaux
377           dH1dzeta = 3d0*H1*dphidzeta/phi + GAMMA*phi3*darg1dzeta/arg1
378           dEc1dzeta = depscdzeta + dH1dzeta
379
380           ddxdzeta = 0.5d0*F43*(opz**F13 - omz**F13)
381           dgcdzeta = -dxc*ddxdzeta*(1d0 - zeta**12) -
382     &                12d0*zeta**11*(1d0 - dxc*(dx - 1d0))
383           dEc0dzeta = dgcdzeta*(epsc0 + H0)
384
385           ddsdzeta = 0.5d0*F53*(opz**F23 - omz**F23)
386           dalphadzeta = -alpha*ddsdzeta/ds
387           dfcadzeta = dfcada*dalphadzeta
388
389           vcpol = dEc1dzeta + dfcadzeta*(Ec0-Ec1) +
390     &             fca*(dEc0dzeta-dEc1dzeta)
391
392           Amat(n,1) = Amat(n,1) + omz*vcpol
393           Amat(n,2) = Amat(n,2) - opz*vcpol
394
395         end if
396
39720    continue
398      end
399
400      Subroutine xc_cscanl_d2()
401      implicit none
402      call errquit(' xc_cscan: d2 not coded ',0,0)
403      return
404      end
405