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