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     J. Sun, A. Ruzsinszky, J.P. Perdew
29c     PRL 115, 036402 (2015)
30c     DOI: 10.1103/PhysRevLett.115036402
31
32      Subroutine xc_xscan(whichfx, tol_rho, fac,  rho, delrho,
33     &                    Amat, Cmat, nq, ipol, Ex,
34     &                    qwght, ldew, func, tau,Mmat)
35      implicit none
36c
37      character*(*) whichfx
38c
39      double precision fac, Ex
40      integer nq, ipol
41      logical ldew
42      double precision func(*)  ! value of the functional [output]
43c
44c     Charge Density & Its Cube Root
45c
46      double precision rho(nq,ipol*(ipol+1)/2)
47c
48c     Charge Density Gradient
49c
50      double precision delrho(nq,3,ipol)
51c
52c     Quadrature Weights
53c
54      double precision qwght(nq)
55c
56c     Sampling Matrices for the XC Potential & Energy
57c
58      double precision Amat(nq,ipol), Cmat(nq,*)
59c
60c     kinetic energy density   or  tau
61c
62      double precision tau(nq,ipol), Mmat(nq,*)
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_xscan_cs(whichfx, tol_rho, fac,  rho, delrho,
73     &                    Amat, Cmat, nq, Ex, 1d0,
74     &                    qwght, ldew, func, tau,Mmat)
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_xscan_cs(whichfx, tol_rho, fac,
84     R           rho(1,ispin+1), delrho(1,1,ispin),
85     &           Amat(1,ispin), Cmat(1,cmatpos),
86     &           nq, Ex, 2d0,
87     &           qwght, ldew, func,
88     T           tau(1,ispin),Mmat(1,ispin))
89         enddo
90
91      endif
92      return
93      end
94      Subroutine xc_xscan_cs(whichfx, tol_rho, fac,  rho, delrho,
95     &                     Amat, Cmat, nq, Ex, facttwo,
96     &                     qwght, ldew, func, tau,Mmat)
97      implicit none
98c
99      character*(*) whichfx
100      double precision fac, Ex
101      integer nq
102      logical ldew
103      double precision func(*)  ! value of the functional [output]
104c
105c     Charge Density & Its Cube Root
106c
107      double precision rho(*)
108c
109c     Charge Density Gradient
110c
111      double precision delrho(nq,3)
112c
113c     Quadrature Weights
114c
115      double precision qwght(nq)
116c
117c     Sampling Matrices for the XC Potential & Energy
118c
119      double precision Amat(nq), Cmat(nq)
120c
121c     kinetic energy density   or  tau
122c
123      double precision tau(nq,*), Mmat(nq)
124c
125      double precision facttwo, afact2 ! 2 for o.s. 1 for c.s.
126c
127      integer n, ifx
128      double precision tol_rho, pi
129      double precision rhoval, rrho, rho13, rho43, rho83
130      double precision tauN, tauW, tauU
131      double precision p, p14, a, g2
132
133      double precision F13, F23, F43, F53, F83, F18
134      double precision Ax, Pconst, thr1, thr2
135      double precision rH0, rK1, rA1, rC1, rC2, rD, rMu
136      double precision rB1, rB2, rB3, rB4
137
138      double precision oma, oma2
139      double precision exp1, exp2, exp3, exp4, exp5
140      double precision rtemp(0:7), rRegu(0:7), rRegu1(7)
141      double precision regalpha, dregalpha
142      double precision x1, x2, x
143      double precision H, Hden, Hnum
144      double precision G
145      double precision Fx, Fa
146
147c     functional derivatives below FFFFFFFFFFFF
148
149      double precision derivr, derivg, derivt
150      double precision dFada
151      double precision dGdp, dHdp, dHdx, dHda
152      double precision dxdp, dx1dp, dx2dp, dxda
153      double precision dpdg, dpdr
154      double precision dFxda, dFxdp, dFxdr, dFxdg, dFxdt
155      double precision dadp,dadg,dadt,dadr
156
157c     functional derivatives above FFFFFFFFFFFF
158
159      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
160      parameter (F83=8.d0/3.d0, F23=2.d0/3.d0)
161      parameter (F18=1.d0/8.d0, F53=5.d0/3.d0)
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      parameter (rRegu = (/ 1d0, -0.667d0, -0.4445555d0,
174     &                    -6.63086601049291d-1, 1.45129704448975d0,
175     &                    -8.87998041596655d-1, 2.34528941478571d-1,
176     &                    -2.31858433223407d-2/) )
177      parameter (rRegu1 = (/ -0.667d0, -2d0*0.4445555d0,
178     &             -3d0*6.63086601049291d-1, 4d0*1.45129704448975d0,
179     &             -5d0*8.87998041596655d-1, 6d0*2.34528941478571d-1,
180     &             -7d0*2.31858433223407d-2/) )
181
182      rB2=dsqrt(5913.d0/405000.d0)
183      rB1=(511.d0/13500.d0)/(2.d0*rB2)
184      rB4=rMu*rMu/rK1-1606.d0/18225.d0-rB1*rB1
185c
186      pi=acos(-1d0)
187      Ax = (-0.75d0)*(3d0/pi)**F13
188      Pconst = (3.d0*pi**2)**F23
189      afact2=1d0/facttwo
190c
191      do n = 1, nq
192         if (rho(n).ge.tol_rho) then
193
194            rhoval=rho(n)*facttwo
195            rho43 = rhoval**F43  ! rho^4/3
196            rrho  = 1d0/rhoval   ! reciprocal of rho
197            rho13 = rho43*rrho
198            rho83 = rho43*rho43
199
200
201            g2 = delrho(n,1)*delrho(n,1) +
202     &           delrho(n,2)*delrho(n,2) +
203     &           delrho(n,3)*delrho(n,3)
204
205            g2 = g2 *facttwo*facttwo
206
207            tauN = tau(n,1)*facttwo
208            tauW = F18*g2*rrho
209            tauU = 0.3d0*Pconst*rhoval**F53
210c
211c     Evaluate the Fx
212c
213            p   = g2/(4d0*Pconst*rho83)
214            p14 = dsqrt(dsqrt(p))
215c
216            if (whichfx.eq.'orig') then
217              a=(tauN-tauW)/tauU
218              if(a.lt.0d0)  a=0d0
219              regalpha = a
220              dregalpha = 1d0
221            else if (whichfx.eq.'regu') then
222              a=(tauN-tauW)/(tauU + 1d-4*facttwo**f53)
223              if(a.lt.0d0)  a=0d0
224              regalpha = a**3/(a**2 + 1d-3)
225              dregalpha = a/(a**2+1d-3)*(3d0*a - 2d0*regalpha)
226            endif
227
228            oma = 1d0 - regalpha
229            oma2 = oma*oma
230
231            exp1 = dexp(-rB4/rMu*p)
232            exp2 = dexp(-rB3*oma2)
233            x1 = rMu*p*(1d0 + rB4/rMu*p*exp1)
234            x2 = rB1*p + rB2*oma*exp2
235
236            x = x1 + x2*x2
237
238            Hden = rK1 + x
239            Hnum = hden + rK1*x
240            H = Hnum/Hden
241
242            if (p14.lt.0.002d0) then
243              exp3 = 0d0
244            else
245              exp3 = dexp(-rA1/p14)
246            endif
247            G = 1d0 - exp3
248
249c
250c Switching function
251c
252            if (whichfx.eq.'orig') then
253              if (regalpha.ge.thr1) then
254                exp4 = 0d0
255              else
256                exp4 = dexp(-rC1*regalpha/oma)
257              end if
258              if (regalpha.le.thr2) then
259                exp5 = 0d0
260              else
261                exp5 = dexp(rC2/oma)
262              end if
263              Fa = exp4 - rD*exp5
264            else if (whichfx.eq.'regu') then
265              if (regalpha.lt.2.5d0) then
266                rtemp(0) = 1d0
267                do ifx=1,7
268                  rtemp(ifx) = rtemp(ifx-1)*regalpha
269                enddo
270                Fa = dot_product(rRegu,rtemp)
271              else
272                exp5 = dexp(rC2/oma)
273                Fa = -rD*exp5
274              end if
275            end if
276
277            Fx = G*(H + Fa*(rH0 - H))
278
279            Ex = Ex + Fx*Ax*rho43*qwght(n)*fac*afact2
280            if (ldew)  func(n)= func(n) + Fx*Ax*rho43*fac*afact2
281
282c     functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF
283
284            dpdr = -F83*p*rrho
285            dpdg = 1d0/(4d0*Pconst*rho83)
286
287            if (whichfx.eq.'orig') then
288              dadg = -F18*rrho/tauU
289              dadt = 1d0/tauU
290              dadr = F13*(8d0*tauW - 5.0*tauN)*rrho/tauU
291            else if (whichfx.eq.'regu') then
292              dadg = -F18*rrho/(tauU + 1d-4*facttwo**f53)
293              dadt = 1d0/(tauU + 1d-4*facttwo**f53)
294              dadr = 5d0/3d0*(p - a)*tauU*rrho/(tauU +1d-4*facttwo**f53)
295            endif
296
297            if (p14.lt.0.001d0) then
298              dGdp = 0d0
299            else
300              dGdp = -0.25d0*rA1*exp3/(p*p14)
301            end if
302
303            dx1dp = rMu + rB4*p*exp1*(2d0 - p*rB4/rMu)
304            dx2dp = rB1
305            dxdp = dx1dp + 2d0*x2*dx2dp
306            dxda = 2d0*rB2*exp2*x2*(2d0*rB3*oma2 - 1d0)*dregalpha
307
308            dHdx = (rK1/Hden)**2
309            dHdp = dHdx*dxdp
310            dHda = dHdx*dxda
311
312c
313c Switching function
314c
315            if (whichfx.eq.'orig') then
316              if ((regalpha.ge.thr1).and.(regalpha.le.thr2)) then
317                dFada = 0d0
318              else
319                dFada = -(rC1*exp4 + rD*exp5*rC2)/oma2
320              end if
321            else if (whichfx.eq.'regu') then
322              if (regalpha.lt.2.5d0) then
323                dFada = dot_product(rRegu1,rtemp(0:6))*dregalpha
324              else
325                dFada = -rD*exp5*rC2/oma2*dregalpha
326              endif
327            endif
328
329            dFxdp = dGdp*(H + Fa*(rH0 - H)) + G*dHdp*(1d0 - Fa)
330            dFxda = G*(dHda + dFada*(rH0 - H) - Fa*dHda)
331
332            dFxdr = dFxda*dadr + dFxdp*dpdr
333            dFxdg = dFxda*dadg + dFxdp*dpdg
334            dFxdt = dFxda*dadt
335
336            derivr = F43*Ax*rho13*Fx + Ax*rho43*dFxdr
337            derivg = Ax*rho43*dFxdg
338            derivt = Ax*rho43*dFxdt
339
340            Amat(n) = Amat(n) + derivr*fac
341c
342c     4x factor comes from gamma_aa = gamma_total/4
343c
344            Cmat(n)=  Cmat(n) + 2.0d0*derivg*fac
345            Mmat(n)=  Mmat(n) + 0.5d0*derivt*fac
346         endif
347      enddo
348      return
349      end
350
351      Subroutine xc_xscan_d2()
352      call errquit(' xscan: d2 not coded ',0,0)
353      return
354      end
355
356
357