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