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
14#include "dft2drv.fh"
15c     Strongly constrained and appropriately normed (SCAN)
16c     functional (Exchange part only)
17c           META GGA
18C         utilizes ingredients:
19c                              rho   -  density
20c                              delrho - gradient of density
21c                              tau - K.S kinetic energy density
22c
23c     Written by:
24c     Daniel Mejia-Rodriguez
25c     QTP, Department of Physics, University of Florida
26c
27c     References:
28c     D. Mejia-Rodriguez, S.B. Trickey
29c     Phys. Rev. A 96, 052512 (2017)
30
31      Subroutine xc_xscanl(tol_rho, fac,  rho, delrho, laprho,
32     &                    Amat, Cmat, Lmat, nq, ipol, Ex,
33     &                    qwght, ldew, func)
34      implicit none
35c
36      double precision fac, Ex
37      integer nq, ipol
38      logical ldew
39      double precision func(*)  ! value of the functional [output]
40c
41c     Charge Density & Its Cube Root
42c
43      double precision rho(nq,ipol*(ipol+1)/2)
44c
45c     Charge Density Gradient
46c
47      double precision delrho(nq,3,ipol)
48c
49c     Charge Density Laplacian
50c
51      double precision laprho(nq,ipol)
52c
53c     Quadrature Weights
54c
55      double precision qwght(nq)
56c
57c     Sampling Matrices for the XC Potential & Energy
58c
59      double precision Amat(nq,ipol), Cmat(nq,*), Lmat(nq,ipol)
60c
61c     Laplacian of the density
62c
63      double precision tol_rho
64c
65      integer ispin,cmatpos
66c
67      if (ipol.eq.1 )then
68c
69c     SPIN-RESTRICTED
70c     Ex = Ex[n]
71c
72         call xc_xscanl_cs(tol_rho, fac,  rho, delrho, laprho,
73     &                    Amat, Cmat, Lmat, nq, Ex, 1d0,
74     &                    qwght, ldew, func)
75      else
76c
77c     SPIN-UNRESTRICTED
78c     Ex = Ex[2n_up]/2 + Ex[2n_down]/2
79
80         do ispin=1,2
81            if (ispin.eq.1) cmatpos=D1_GAA
82            if (ispin.eq.2) cmatpos=D1_GBB
83            call xc_xscanl_cs(tol_rho, fac,
84     R           rho(1,ispin+1), delrho(1,1,ispin), laprho(1,ispin),
85     &           Amat(1,ispin), Cmat(1,cmatpos), Lmat(1,ispin),
86     &           nq, Ex, 2d0,
87     &           qwght, ldew, func)
88         enddo
89
90      endif
91      return
92      end
93      Subroutine xc_xscanl_cs(tol_rho, fac,  rho, delrho, laprho,
94     &                     Amat, Cmat, Lmat, nq, Ex, facttwo,
95     &                     qwght, ldew, func)
96      implicit none
97c
98      double precision fac, Ex
99      integer nq
100      logical ldew
101      double precision func(*)  ! value of the functional [output]
102c
103c     Charge Density & Its Cube Root
104c
105      double precision rho(*)
106c
107c     Charge Density Gradient
108c
109      double precision delrho(nq,3)
110c
111c     Charge Density Laplacian
112c
113      double precision laprho(nq)
114c
115c     Quadrature Weights
116c
117      double precision qwght(nq)
118c
119c     Sampling Matrices for the XC Potential & Energy
120c
121      double precision Amat(nq), Cmat(nq), Lmat(nq)
122c
123      double precision facttwo, afact2 ! 2 for o.s. 1 for c.s.
124c
125      integer n
126      double precision tol_rho, pi
127      double precision rhoval, rrho, rho13, rho43, rho53, rho83
128      double precision g2, lapval
129      double precision p, p14, a, q
130c
131      double precision thr1, thr2
132      double precision F13, F23, F43, F53, F83
133      double precision Ax, Pconst
134      double precision rH0, rK1, rA1, rC1, rC2, rD, rMu
135      double precision rB1, rB2, rB3, rB4
136
137      double precision oma, oma2
138      double precision exp1, exp2, exp3, exp4, exp5
139      double precision x1, x2, x
140      double precision H, Hden, Hnum
141      double precision G
142      double precision Fa, Fs, Fx
143
144c     functional derivatives below FFFFFFFFFFFF
145
146      double precision derivr, derivg, derivl
147      double precision dFada
148      double precision dGdp, dHdp, dHdx, dHda
149      double precision dxdp, dx1dp, dx2dp, dxda
150      double precision dadp, dadq, dpdg, dpdr, dqdr, dqdl
151      double precision dfsdp, dfsdq
152      double precision dFxda, dFxdp, dFxdr, dFxdg, dFxdl
153
154c     functional derivatives above FFFFFFFFFFFF
155
156      parameter (F13=1.0d0/3.0d0)
157      parameter (F23=2.0d0*F13)
158      parameter (F43=1.0d0+F13)
159      parameter (F53=1.0d0+F23)
160      parameter (F83=1.0d0+F53)
161
162      parameter (thr1=0.996d0, thr2=1.004d0)
163
164      parameter (rH0=1.174d0)
165      parameter (rK1=0.065d0)
166      parameter (rA1=4.9479d0)
167      parameter (rC1=0.667d0)
168      parameter (rC2=0.8d0)
169      parameter (rD=1.24d0)
170      parameter (rMu=10.0d0/81.0d0)
171      parameter (rB3=0.5d0)
172
173      rB2=dsqrt(5913.d0/405000.d0)
174      rB1=(511.d0/13500.d0)/(2.d0*rB2)
175      rB4=rMu*rMu/rK1-1606.d0/18225.d0-rB1*rB1
176c
177      pi=acos(-1d0)
178      Ax = (-0.75d0)*(3d0/pi)**F13
179      Pconst = (3.d0*pi**2)**F23
180      afact2=1d0/facttwo
181c
182      do n = 1, nq
183         if (rho(n).ge.tol_rho) then
184c
185            call ts_pc(tol_rho, rho(n), delrho(n,1:3), laprho(n),
186     &                 dfsdp, dfsdq, fs, facttwo)
187
188            rhoval=rho(n)*facttwo
189            rho43 = rhoval**F43
190            rrho  = 1d0/rhoval
191            rho13 = rho43*rrho
192            rho83 = rho43*rho43
193            rho53 = rho43*rho13
194
195            g2 = delrho(n,1)*delrho(n,1) +
196     &           delrho(n,2)*delrho(n,2) +
197     &           delrho(n,3)*delrho(n,3)
198
199            g2 = g2 *facttwo*facttwo
200
201            lapval = laprho(n)*facttwo
202c
203c     Evaluate the Fx
204c
205            p   = g2/(4d0*Pconst*rho83)
206            p14 = dsqrt(dsqrt(p))
207
208            q = lapval/(4d0*Pconst*rho53)
209c
210            a=fs-F53*p
211            if(a.lt.0d0)  a=0d0
212            oma = 1d0 - a
213            oma2 = oma*oma
214
215            exp1 = dexp(-rB4/rMu*p)
216            exp2 = dexp(-rB3*oma2)
217            x1 = rMu*p*(1d0 + rB4/rMu*p*exp1)
218            x2 = rB1*p + rB2*oma*exp2
219
220            x = x1 + x2*x2
221
222            Hden = rK1 + x
223            Hnum = hden + rK1*x
224            H = Hnum/Hden
225
226            if (p14.lt.0.002d0) then
227              exp3 = 0d0
228            else
229              exp3 = dexp(-rA1/p14)
230            endif
231            G = 1d0 - exp3
232
233            if (a.ge.thr1) then
234              exp4 = 0d0
235            else
236              exp4 = dexp(-rC1*a/oma)
237            end if
238
239            if (a.le.thr2) then
240              exp5 = 0d0
241            else
242              exp5 = dexp(rC2/oma)
243            end if
244
245            Fa = exp4 - rD*exp5
246
247            Fx = G*(H + Fa*(rH0 - H))
248
249            Ex = Ex + Fx*Ax*rho43*qwght(n)*fac*afact2
250            if (ldew)  func(n)= func(n) + Fx*Ax*rho43*fac*afact2
251
252c     functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF
253
254            dpdr = -F83*p*rrho
255            dpdg = 1d0/(4d0*Pconst*rho83)
256
257            dqdr = -F53*q*rrho
258            dqdl = 1.d0/(4d0*Pconst*rho53)
259
260            dadp = dfsdp - F53
261            dadq = dfsdq
262
263            if (p14.lt.0.001d0) then
264              dGdp = 0d0
265            else
266              dGdp = -0.25d0*rA1*exp3/(p*p14)
267            end if
268
269            dx1dp = rMu + rB4*p*exp1*(2d0 - p*rB4/rMu)
270            dx2dp = rB1
271            dxdp = dx1dp + 2d0*x2*dx2dp
272            dxda = 2d0*rB2*exp2*x2*(2d0*rB3*oma2 - 1d0)
273
274            dHdx = (rK1/Hden)**2
275            dHdp = dHdx*dxdp
276            dHda = dHdx*dxda
277
278            if ((a.ge.thr1).and.(a.le.thr2)) then
279              dFada = 0d0
280            else
281              dFada = -(rC1*exp4 + rD*exp5*rC2)/oma2
282            end if
283
284            dFxdp = dGdp*(H + Fa*(rH0 - H)) + G*dHdp*(1d0 - Fa)
285            dFxda = G*(dHda + dFada*(rH0 - H) - Fa*dHda)
286
287            dFxdr = dFxda*(dadp*dpdr + dadq*dqdr) + dFxdp*dpdr
288            dFxdg = dFxda*dadp*dpdg + dFxdp*dpdg
289            dFxdl = dFxda*dadq*dqdl
290
291            derivr = F43*Ax*rho13*Fx + Ax*rho43*dFxdr
292            derivg = Ax*rho43*dFxdg
293            derivl = Ax*rho43*dFxdl
294
295            Amat(n) = Amat(n) + derivr*fac
296c
297c     4x factor comes from gamma_aa = gamma_total/4
298c
299            Cmat(n)=  Cmat(n) + 2.0d0*derivg*fac
300            Lmat(n)=  Lmat(n) + derivl*fac
301         endif
302      enddo
303      return
304      end
305
306      Subroutine xc_xscanl_d2()
307      call errquit(' xscanl: d2 not coded ',0,0)
308      return
309      end
310