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