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