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 14c --------------------------------------------------------------------- 15c 16c Strongly constrained and appropriately normed (SCAN) 17c functional (Correlation part only) 18c META GGA 19C utilizes ingredients: 20c rho - density 21c delrho - gradient of density 22c tau - K.S kinetic energy density 23c 24c Written by: 25c Daniel Mejia-Rodriguez 26c QTP, Department of Physics, University of Florida 27c 28c References: 29c J. Sun, A. Ruzsinszky, J.P. Perdew 30c PRL 115, 036402 (2015) 31c DOI: 10.1103/PhysRevLett.115036402 32c 33c A.P. Bartok and J.R. Yates 34c JCP 150, 161101 (2019) 35c DOI: 10.1063/1.5094646 36 37 Subroutine xc_cscan(whichfc, tol_rho, cfac, rho, delrho, Amat, 38 & Cmat, nq, ipol, Ec, qwght, ldew, func, 39 & tau, Mmat) 40c 41 implicit none 42c 43#include "errquit.fh" 44#include "dft2drv.fh" 45c 46c Input and other parameters 47c 48 character(*) whichfc 49 integer ipol, nq 50 double precision dummy(1) 51 double precision cfac 52 logical ldew 53 double precision func(*) 54 double precision fac 55 double precision tol_rho 56c 57c Threshold parameters 58c 59 double precision thr1,thr2 60 parameter (thr1=0.996d0,thr2=1.004d0) 61c 62c Correlation energy 63c 64 double precision Ec 65c 66c Charge Density 67c 68 double precision rho(nq,ipol*(ipol+1)/2) 69c 70c Charge Density Gradient 71c 72 double precision delrho(nq,3,ipol), gammaval, gam12 73c 74c Kinetic Energy Density 75c 76 double precision tau(nq,ipol), tauN 77c 78c Quadrature Weights 79c 80 double precision qwght(nq) 81c 82c Sampling Matrices for the XC Potential 83c 84 double precision Amat(nq,ipol), Cmat(nq,*) 85 double precision Mmat(nq,*) 86c 87c Intermediate derivatives, results, etc. 88c 89 integer n, ifc 90 double precision ntot,n13,n83,tautot 91 double precision dn2,p,rs,rs12,drsdn,dpdg,dpdn 92 double precision epsc,depscdrs,depscdzeta 93 double precision zeta,omz,opz 94 double precision phi,phi2,phi3,dphidzeta 95 double precision t2,dt2dp,dt2drs,dt2dzeta 96 double precision BETA,dBETAdrs 97 double precision A,dAdrs,dAdzeta 98 double precision At2,Gaux,Gat2 99 double precision w1fac,expw1,w1,dw1drs,dw1dzeta 100 double precision arg1,darg1dp,darg1drs,darg1dzeta 101 double precision H1,dH1dp,dH1drs,dH1dzeta 102 double precision Ec1,dEc1drs,dEc1dp,dEc1dzeta,dEc1dn,dEc1dg 103 double precision epsc0,depsc0drs,depsc0dn 104 double precision dx,ddxdzeta 105 double precision gc,dgcdzeta 106 double precision ginf,dginfdp 107 double precision w0fac,expw0,w0,dw0drs 108 double precision arg0,darg0dp,darg0drs 109 double precision H0,dH0dp,dH0drs 110 double precision Ec0,dEc0dp,dEc0drs,dEc0dzeta,dEc0dn,dEc0dg 111 double precision ds,ddsdzeta 112 double precision tueg,tvw 113 double precision alpha,dalphadzeta,dalphadn,dalphadg,dalphadt 114 double precision oma,oma2 115 double precision fca,dfcada,dfcadg,dfcadzeta,dfcadn 116 double precision exp5,exp6 117 double precision vcpol,vcn 118 119 double precision GAMMA,BETAzero,pi,p14a,p14b,p14f,ckf,ckf2 120 parameter (GAMMA = 0.03109069086965489503494086371273d0) 121 parameter (BETAzero = 0.06672455060314922d0) 122 parameter (p14a=0.1d0,p14b=0.1778d0) 123 124 double precision b1c,b2c,b3c,c1c,c2c,dc,dxc,xi 125 parameter (b1c=0.0285764d0,b2c=0.0889d0,b3c=0.125541d0) 126 parameter (c1c=0.64d0,c2c=1.5d0,dc=0.7d0,dxc=2.3631d0) 127 parameter (xi=0.12802585262625815d0) 128 129 double precision F4,F13,F23,F43,F53,F83,F14,F8,F5,F16 130 parameter (F4=4d0,F13=1d0/3d0,F23=F13+F13,F43=F13+1d0) 131 parameter (F53=F23+1d0,F83=F53+1d0,F14=0.25d0) 132 parameter (F8=8d0,F5=5d0,F16=1d0/6d0) 133 134 double precision rRegu(0:7), rRegu1(7), rtemp(0:7) 135 double precision regalpha, dregalpha 136 parameter ( rRegu = (/ 1d0, -0.64d0, -0.4352d0, 137 & -1.535685604549d0, 3.061560252175d0, 138 & -1.915710236206d0, 5.16884468372d-1, 139 & -5.1848879792d-2 /) ) 140 parameter ( rRegu1 = (/ -0.64d0, -2d0*0.4352d0, 141 & -3d0*1.535685604549d0, 4d0*3.061560252175d0, 142 & -5d0*1.915710236206d0, 6d0*5.16884468372d-1, 143 & -7d0*5.1848879792d-2 /) ) 144c 145c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 146c 147 Pi = dacos(-1d0) 148 p14f = (3d0/(4d0*Pi))**F13 149 ckf = (3d0*Pi*Pi)**F13 150 ckf2 = ckf*ckf 151 152 do 20 n = 1, nq 153c 154 ntot=rho(n,1) 155 if (ntot.le.tol_rho) goto 20 156 157 n13=ntot**F13 158 n83=ntot**F83 159c 160 if (ipol.eq.1) then 161 zeta = 0d0 162 else 163 zeta = (rho(n,2) - rho(n,3))/ntot 164 if (zeta.lt.-1d0) zeta=-1d0 165 if (zeta.gt. 1d0) zeta= 1d0 166 endif 167c 168 if (ipol.eq.1) then 169 dn2=delrho(n,1,1)**2 + delrho(n,2,1)**2 + delrho(n,3,1)**2 170 else 171 dn2=(delrho(n,1,1)+delrho(n,1,2))**2 + 172 & (delrho(n,2,1)+delrho(n,2,2))**2 + 173 & (delrho(n,3,1)+delrho(n,3,2))**2 174 end if 175 176 if (ipol.eq.1) then 177 tautot=tau(n,1) 178 else 179 tautot=tau(n,1)+tau(n,2) 180 endif 181 182 p=dn2/(F4*ckf2*n83) 183 dpdg = 1d0/(F4*ckf2*n83) 184 dpdn = -F83*p/ntot 185 186 rs=p14f/n13 187 drsdn=-F13*rs/ntot 188 rs12=dsqrt(rs) 189c 190 call lsdac(tol_rho,rs,zeta,epsc,depscdrs,depscdzeta,dummy, 191 & dummy,dummy) 192c 193 opz = 1d0 + zeta 194 omz = 1d0 - zeta 195 phi = 0.5d0*(opz**F23 + omz**F23) 196 phi2 = phi*phi 197 phi3 = phi2*phi 198c 199 BETA = BETAzero*(1d0 + p14a*rs)/(1d0 + p14b*rs) 200 dBETAdrs = BETAzero*(p14a - p14b)/(1d0 + p14b*rs)**2 201c 202 w1fac = epsc/(GAMMA*phi3) 203 expw1 = dexp(-w1fac) 204 w1 = expw1 - 1d0 205 dw1drs = -expw1*depscdrs/(GAMMA*phi3) 206 207 A = BETA/GAMMA/w1 208 dAdrs = dBETAdrs/GAMMA/w1 - A*dw1drs/w1 209 210 t2 = ckf2*p/(4d0*phi2*rs*2d0**F23) 211 dt2dp = ckf2/(4d0*phi2*rs*2d0**F23) 212 dt2drs = -t2/rs 213 214 At2 = A*t2 215 Gaux = 1d0 + 4d0*At2 216 GAt2 = 1d0/dsqrt(dsqrt(Gaux)) 217 218 arg1 = 1d0 + w1*(1d0 - GAt2) 219 darg1drs = dw1drs*(1d0 - GAt2) + 220 & Gat2*w1*(t2*dAdrs + A*dt2drs)/Gaux 221 darg1dp = Gat2*w1*A*dt2dp/Gaux 222 223 224 H1 = GAMMA*phi3*dlog(arg1) 225 dH1dp = GAMMA*phi3*darg1dp/arg1 226 dH1drs = GAMMA*phi3*darg1drs/arg1 227 228 Ec1 = epsc + H1 229 dEc1drs = depscdrs + dH1drs 230 dEc1dp = dH1dp 231c 232c -------------------------------------------------------------- 233c 234 epsc0 = -b1c/(1d0 + b2c*rs12 + b3c*rs) 235 depsc0drs = b1c*(b3c + 0.5d0*b2c/rs12)/ 236 & (1d0 + b2c*rs12 + b3c*rs)**2 237 238 dx = 0.5d0*(opz**F43 + omz**F43) 239 gc = (1d0 - dxc*(dx - 1d0))*(1d0 - zeta**12) 240 241 w0fac = epsc0/b1c 242 expw0 = dexp(-w0fac) 243 w0 = expw0 - 1d0 244 dw0drs = -depsc0drs*expw0/b1c 245 246 ginf = (1d0/(1d0 + 4d0*xi*p))**F14 247 dginfdp = -xi*ginf/(1d0 + 4d0*xi*p) 248 249 arg0 = 1d0 + w0*(1d0 - ginf) 250 darg0drs = dw0drs*(1d0 - ginf) 251 darg0dp = -w0*dginfdp 252 253 H0 = b1c*dlog(arg0) 254 dH0dp = b1c*darg0dp/arg0 255 dH0drs = b1c*darg0drs/arg0 256 257 Ec0 = (epsc0 + H0)*gc 258 dEc0drs = gc*(depsc0drs + dH0drs) 259 dEc0dp = gc*dH0dp 260c 261c -------------------------------------------------------------- 262c 263 ds = 0.5d0*(opz**F53 + omz**F53) 264 tueg = 0.3d0*ckf*ckf*ds*ntot**F53 265 tvw = 0.125d0*dn2/ntot 266 if (whichfc.eq.'orig') then 267 alpha = (tautot - tvw)/tueg 268 if (alpha.lt.0d0) alpha=0d0 269 regalpha = alpha 270 dregalpha = 1d0 271 else if (whichfc.eq.'regu') then 272 alpha = (tautot - tvw)/(tueg + 1d-4*ds) 273 if (alpha.lt.0d0) alpha=0d0 274 regalpha = alpha**3/(alpha**2 + 1d-3) 275 dregalpha = alpha/(alpha**2 + 1d-3) * 276 & (3d0*alpha - 2d0*regalpha) 277 endif 278 oma = 1d0 - regalpha 279 oma2 = oma*oma 280 281 if (whichfc.eq.'orig') then 282 if (regalpha.ge.thr1) then 283 exp5 = 0d0 284 else 285 exp5 = dexp(-c1c*regalpha/oma) 286 end if 287 if (regalpha.le.thr2) then 288 exp6 = 0d0 289 else 290 exp6 = dexp(c2c/oma) 291 end if 292 fca = exp5 - dc*exp6 293 if (regalpha.ge.thr1.and.regalpha.le.thr2) then 294 dfcada = 0d0 295 else 296 dfcada = -(c1c*exp5 + dc*exp6*c2c)/oma2 297 end if 298 else if(whichfc.eq.'regu') then 299 if (regalpha.lt.2.5d0) then 300 rtemp(0) = 1d0 301 do ifc=1,7 302 rtemp(ifc) = rtemp(ifc-1)*regalpha 303 enddo 304 fca = dot_product(rRegu,rtemp) 305 dfcada = dot_product(rRegu1, rtemp(0:6))*dregalpha 306 else 307 exp6 = dexp(c2c/oma) 308 fca = -dc*exp6 309 dfcada = -dc*exp6*c2c/oma2*dregalpha 310 endif 311 endif 312c 313c -------------------------------------------------------------- 314c 315 Ec = Ec + cfac*ntot*(Ec1 + fca*(Ec0-Ec1))*qwght(n) 316 if (ldew) func(n) = func(n) + cfac*ntot*(Ec1 + fca*(Ec0-Ec1)) 317c 318 if (whichfc.eq.'orig') then 319 dalphadn = f53*(p/ds-alpha)/ntot 320 dalphadg = -0.125d0/(tueg*ntot) 321 dalphadt = 1d0/tueg 322 else if (whichfc.eq.'regu') then 323 dalphadn = 5d0/3d0*(p/ds-alpha)*tueg/((tueg+1d-4*ds)*ntot) 324 dalphadg = -0.125d0/((tueg+1d-4*ds)*ntot) 325 dalphadt = 1d0/(tueg+1d-4*ds) 326 endif 327 328 dEc1dn = dEc1drs*drsdn + dEc1dp*dpdn 329 dEc0dn = dEc0drs*drsdn + dEc0dp*dpdn 330 dfcadn = dfcada*dalphadn 331 vcn = Ec1 + fca*(Ec0-Ec1) + 332 & ntot*(dEc1dn + fca*(dEc0dn-dEc1dn) + dfcadn*(Ec0-Ec1)) 333 334 Amat(n,1) = Amat(n,1) + cfac*vcn 335 336 dEc1dg = dEc1dp*dpdg 337 dEc0dg = dEc0dp*dpdg 338 dfcadg = dfcada*dalphadg 339 340 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + cfac*ntot*(dEc1dg + 341 & fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1)) 342 Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + cfac*ntot*(dEc1dg + 343 & fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1))*2d0 344 345 Mmat(n,1) = Mmat(n,1) + 0.5d0*cfac* 346 & ntot*dfcada*dalphadt*(Ec0-Ec1) 347 348 if (ipol.eq.2) then 349 Amat(n,2) = Amat(n,2) + cfac*vcn 350 351 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + cfac*ntot*(dEc1dg + 352 & fca*(dEc0dg-dEc1dg) + dfcadg*(Ec0-Ec1)) 353 354 Mmat(n,2) = Mmat(n,2) + 0.5d0*cfac* 355 & ntot*dfcada*dalphadt*(Ec0-Ec1) 356 357 if (omz.lt.tol_rho) then 358 dphidzeta = 0.5d0*F23*(opz**F23/opz) 359 else if (opz.lt.tol_rho) then 360 dphidzeta = -0.5d0*F23*(omz**F23/omz) 361 else 362 dphidzeta = 0.5d0*F23*(opz**F23/opz - omz**F23/omz) 363 end if 364 365 dt2dzeta = -2d0*t2*dphidzeta/phi 366 dw1dzeta = (3d0*w1fac*dphidzeta/phi - 367 & depscdzeta/(GAMMA*phi3))*expw1 368 dAdzeta = -A*dw1dzeta/w1 369 darg1dzeta = dw1dzeta*(1d0 - Gat2) + 370 & Gat2*w1*(t2*dAdzeta + A*dt2dzeta)/Gaux 371 dH1dzeta = 3d0*H1*dphidzeta/phi + GAMMA*phi3*darg1dzeta/arg1 372 dEc1dzeta = depscdzeta + dH1dzeta 373 374 ddxdzeta = 0.5d0*F43*(opz**F13 - omz**F13) 375 dgcdzeta = -dxc*ddxdzeta*(1d0 - zeta**12) - 376 & 12d0*zeta**11*(1d0 - dxc*(dx - 1d0)) 377 dEc0dzeta = dgcdzeta*(epsc0 + H0) 378 379 ddsdzeta = 0.5d0*F53*(opz**F23 - omz**F23) 380 381 dalphadzeta = -alpha*ddsdzeta/ds 382 383 dfcadzeta = dfcada*dalphadzeta 384 385 vcpol = dEc1dzeta + dfcadzeta*(Ec0-Ec1) + 386 & fca*(dEc0dzeta-dEc1dzeta) 387 388 Amat(n,1) = Amat(n,1) + cfac*omz*vcpol 389 Amat(n,2) = Amat(n,2) - cfac*opz*vcpol 390 391 end if 392 39320 continue 394 end 395 396 Subroutine xc_cscan_d2() 397 implicit none 398 call errquit(' xc_cscan: d2 not coded ',0,0) 399 return 400 end 401