1c Copyright 2019 (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" 15#include "dft3drv.fh" 16c Nearly correct asymptotic potential (NCAP) functional 17c (Exchange part only) 18c GGA 19C utilizes ingredients: 20c rho - density 21c delrho - gradient of density 22c 23c Written by: 24c Daniel Mejia-Rodriguez 25c QTP, Department of Physics, University of Florida 26c 27c References: 28c J. Carmona-Espindola, J.L. Gazquez, A. Vela, S.B. Trickey 29c JCTC 15, 303 (2019). 30c DOI: 10.1021/acs.jctc.8b00998 31c 32c 33#if !defined SECOND_DERIV && !defined THIRD_DERIV 34 Subroutine xc_xncap(tol_rho, fac, rho, delrho, 35 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 36#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 37 Subroutine xc_xncap_d2(tol_rho, fac, rho, delrho, 38 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 39 & qwght,ldew,func) 40#else 41 Subroutine xc_xncap_d3(tol_rho, fac, rho, delrho, 42 & Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 43 & nq, ipol, Ex, qwght,ldew,func) 44#endif 45c 46 implicit none 47c 48 double precision fac, Ex 49 integer nq, ipol 50 logical ldew 51 double precision func(*) ! value of the functional [output] 52c 53c Charge Density & Its Cube Root 54c 55 double precision rho(nq,ipol*(ipol+1)/2) 56c 57c Charge Density Gradient 58c 59 double precision delrho(nq,3,ipol) 60c 61c Quadrature Weights 62c 63 double precision qwght(nq) 64c 65c Sampling Matrices for the XC Potential & Energy 66c 67 double precision amat(nq,ipol), cmat(nq,*) 68c 69#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 70 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 71#endif 72#ifdef THIRD_DERIV 73 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 74#endif 75c 76 double precision tol_rho, pi, mu, alpha, beta, zeta 77 double precision ckf, Ax 78 double precision F43, F13, F23, F49 79c 80 parameter(mu=0.2195149727645171d0, beta=0.01808569669d0) 81 parameter(zeta=0.30412141859531383d0) 82 83 parameter (F43=4.d0/3.d0, F13=1.d0/3.d0, F23=2.0d0/3.0d0) 84 parameter (F49=4d0/9d0) 85c 86 integer n 87 double precision rrho, rho43, rho13, gamma, gam12, s, s2, d1s 88 double precision arcsinh, darcsinh, d2arcsinh, d3arcsinh 89 double precision Fx, dFx, d2Fx, d3Fx 90 double precision tansin, dtansin 91 double precision factor, dfactor 92 double precision asymp, dasymp 93 double precision denom, ddenom 94 95#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 96 double precision d2denom, d2asymp, d2tansin, d2factor, d2s, rhom23 97#endif 98#if defined(THIRD_DERIV) 99 double precision d3denom, d3asymp, d3tansin, d3factor, d3s, rhom53 100#endif 101 double precision ops, omz, logops 102 103 arcsinh(s)=dlog(s+dsqrt(1d0+s*s)) 104 darcsinh(s)=1d0/dsqrt(1d0+s*s) 105 d2arcsinh(s) = -s/dsqrt(1d0+s*s)**3 106 d3arcsinh(s) = (2d0*s*s - 1d0)/dsqrt(1d0+s*s)**5 107c 108 pi = acos(-1.d0) 109 ckf = (3d0*pi*pi)**F13 110 Ax = -3d0/(4d0*pi)*ckf 111 alpha = 4d0*pi/3d0*beta/mu 112c 113 if (ipol.eq.1 )then 114c 115c ======> SPIN-RESTRICTED <====== 116c 117#ifdef IFCV81 118CDEC$ NOSWP 119#endif 120 do 10 n = 1, nq 121 if (rho(n,1).lt.tol_rho) goto 10 122 rho43 = rho(n,1)**F43 123 rrho = 1d0/rho(n,1) 124 rho13 = rho43*rrho 125 gamma = delrho(n,1,1)*delrho(n,1,1) + 126 & delrho(n,2,1)*delrho(n,2,1) + 127 & delrho(n,3,1)*delrho(n,3,1) 128 if (dsqrt(gamma).le.tol_rho**2) goto 10 129 130 s = dsqrt(gamma)/(2d0*ckf*rho43) 131 s2 = s*s 132 133 if (s.lt.1d-8) then 134 Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s + 135 & mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2 136 dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 + 137 & 4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s 138 d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s + 139 & 12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2 140 d3Fx = 6d0*alpha*zeta*mu + 141 & 24d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s 142 else 143 ops = 1d0 + s 144 omz = 1d0 - zeta 145 logops = dlog(ops) 146 147 asymp = 1d0 + alpha*(omz*s*logops + zeta*s) 148 dasymp = alpha*( (s + zeta)/ops + omz*logops ) 149 150 if (s.gt.30d0) then 151 tansin = arcsinh(s) 152 dtansin = darcsinh(s) 153 else 154 tansin = dtanh(s)*arcsinh(s) 155 dtansin = arcsinh(s)/dcosh(s)**2 + dtanh(s)*darcsinh(s) 156 endif 157 158 denom = 1d0 + beta*tansin 159 ddenom = beta*dtansin 160 161 factor = tansin/denom 162 dfactor = (dtansin - factor*ddenom)/denom 163 164 Fx = 1d0 + mu*factor*asymp 165 dFx = mu*(dfactor*asymp + factor*dasymp) 166 167#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 168 if (s.gt.30d0) then 169 d2tansin = d2arcsinh(s) 170 else 171 d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 + 172 & dtanh(s)*d2arcsinh(s) - 173 & 2d0*tansin/dcosh(s)**2 174 endif 175 d2asymp = alpha*omz*(2d0 + s)/ops**2 176 d2denom = beta*d2tansin 177 d2factor = (d2tansin - 2d0*dfactor*ddenom - 178 & factor*d2denom)/denom 179 d2Fx = mu*(d2factor*asymp + 2d0*dfactor*dasymp + 180 & factor*d2asymp) 181 182#endif 183#if defined(THIRD_DERIV) 184 if (s.gt.30d0) then 185 d3tansin = d3arcsinh(s) 186 else 187 d3tansin = 3d0*d2arcsinh(s)/dcosh(s)**2 - 188 & 4d0*darcsinh(s)*dtanh(s)/dcosh(s)**2 + 189 & dtanh(s)*d3arcsinh(s) - 190 & 2d0*dtansin/dcosh(s)**2 + 191 & 4d0*tansin*dtanh(s)/dcosh(s)**2 192 endif 193 d3asymp = -alpha*omz*(3d0 + s)/ops**3 194 d3denom = beta*d3tansin 195 d3factor = (d3tansin - 3d0*d2factor*ddenom - 196 & 3d0*dfactor*d2denom - factor*d3denom)/denom 197 d3Fx = mu*(d3factor*asymp + 3d0*d2factor*dasymp + 198 & 3d0*dfactor*d2asymp + factor*d3asymp) 199#endif 200 endif 201 202 Ex = Ex + Ax*rho43*Fx*qwght(n)*fac 203 if (ldew) func(n) = func(n) + Ax*rho43*Fx*fac 204 205 d1s = 0.5d0*s/gamma 206 Amat(n,1) = Amat(n,1) + F43*Ax*rho13*(Fx - s*dFx)*fac 207 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*Ax*rho43*dFx*d1s*fac 208 209#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 210 211 d2s = -0.5d0*d1s/gamma 212 rhom23 = rho13*rrho 213 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*F49*Ax*rhom23* 214 & (Fx - s*dFx + 4d0*s2*d2Fx)*fac 215 216 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) - 217 & (F43*Ax*rho13*d2Fx*d1s*s)*4d0*fac 218 219 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + Ax*rho43* 220 & (d2Fx*d1s**2 + dFx*d2s)*8d0*fac 221#endif 222#ifdef THIRD_DERIV 223 rhom53 = rhom23*rrho 224 d3s = -1.5d0*d2s/gamma 225 226 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) - 4d0*fac* 227 & F23*F49*Ax*rhom53*(Fx - s*dFx + 18d0*s2*d2Fx +8d0*s2*s*d3Fx) 228 229 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + 230 & F49*Ax*rhom23*d1s*s*(7d0*d2Fx + 4d0*d3Fx*s)*8d0*fac 231 232 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) - 233 & F43*Ax*rho13*(d2Fx*d1s**2 + d3Fx*d1s**2*s + 234 & d2Fx*d2s*s)*16d0*fac 235 236 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + 237 & Ax*rho43*(d3Fx*d1s**3 + 3d0*d2Fx*d1s*d2s + 238 & dFx*d3s)*32d0*fac 239 240#endif 241 10 continue 242c 243 else 244c 245c ======> SPIN-UNRESTRICTED <====== 246c 247#ifdef IFCV81 248CDEC$ NOSWP 249#endif 250 do 20 n = 1, nq 251 if (rho(n,1).lt.tol_rho) goto 20 252c 253c Alpha 254c 255 if (rho(n,2).lt.tol_rho) goto 25 256 rho43 = (2d0*rho(n,2))**F43 257 rrho = 0.5d0/rho(n,2) 258 rho13 = rho43*rrho 259c 260 gamma = delrho(n,1,1)*delrho(n,1,1) + 261 & delrho(n,2,1)*delrho(n,2,1) + 262 & delrho(n,3,1)*delrho(n,3,1) 263 gam12 = 2d0*dsqrt(gamma) 264 if (gam12.le.tol_rho**2) goto 25 265c 266 s = gam12/(2d0*ckf*rho43) 267 s2 = s*s 268 269 d1s = 0.5d0*s/gamma 270 271 if (s.lt.1d-8) then 272 Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s + 273 & mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2 274 dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 + 275 & 4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s 276 d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s + 277 & 12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2 278 d3Fx = 6d0*alpha*zeta*mu + 279 & 24d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s 280 else 281 ops = 1d0 + s 282 omz = 1d0 - zeta 283 logops = dlog(ops) 284 285 asymp = 1d0 + alpha*(omz*s*logops + zeta*s) 286 dasymp = alpha*( (s + zeta)/ops + omz*logops ) 287 288 if (s.gt.30d0) then 289 tansin = arcsinh(s) 290 dtansin = darcsinh(s) 291 else 292 tansin = dtanh(s)*arcsinh(s) 293 dtansin = arcsinh(s)/dcosh(s)**2 + dtanh(s)*darcsinh(s) 294 endif 295 296 denom = 1d0 + beta*tansin 297 ddenom = beta*dtansin 298 299 factor = tansin/denom 300 dfactor = (dtansin - factor*ddenom)/denom 301 302 Fx = 1d0 + mu*factor*asymp 303 dFx = mu*(dfactor*asymp + factor*dasymp) 304#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 305 if (s.gt.30d0) then 306 d2tansin = d2arcsinh(s) 307 else 308 d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 + 309 & dtanh(s)*d2arcsinh(s) - 310 & 2d0*tansin/dcosh(s)**2 311 endif 312 d2asymp = alpha*omz*(2d0 + s)/ops**2 313 d2denom = beta*d2tansin 314 d2factor = (d2tansin - 2d0*dfactor*ddenom - 315 & factor*d2denom)/denom 316 d2Fx = mu*(d2factor*asymp + 2d0*dfactor*dasymp + 317 & factor*d2asymp) 318#endif 319#if defined(THIRD_DERIV) 320 if (s.gt.30d0) then 321 d3tansin = d3arcsinh(s) 322 else 323 d3tansin = 3d0*d2arcsinh(s)/dcosh(s)**2 - 324 & 4d0*darcsinh(s)*dtanh(s)/dcosh(s)**2 + 325 & dtanh(s)*d3arcsinh(s) - 326 & 2d0*dtansin/dcosh(s)**2 + 327 & 4d0*tansin*dtanh(s)/dcosh(s)**2 328 endif 329 d3asymp = -alpha*omz*(3d0 + s)/ops**3 330 d3denom = beta*d3tansin 331 d3factor = (d3tansin - 3d0*d2factor*ddenom - 332 & 3d0*dfactor*d2denom - factor*d3denom)/denom 333 d3Fx = mu*(d3factor*asymp + 3d0*d2factor*dasymp + 334 & 3d0*dfactor*d2asymp + factor*d3asymp) 335#endif 336 endif 337 338 Ex = Ex + 0.5d0*Ax*rho43*Fx*qwght(n)*fac 339 if (ldew) func(n) = func(n) + 0.5d0*Ax*rho43*Fx*fac 340 341 Amat(n,1) = Amat(n,1) + F43*Ax*rho13*(Fx - s*dFx)*fac 342 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 343 & 0.5d0*Ax*rho43*dFx*d1s*fac 344 345#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 346 d2s = -0.5d0*d1s/gamma 347 rhom23 = rho13*rrho 348 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*F49*Ax*rhom23* 349 & (Fx - s*dFx + 4d0*s2*d2Fx)*fac 350 351 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) - 352 & (F43*Ax*rho13*d2Fx*d1s*s)*fac 353 354 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + fac*Ax*rho43* 355 & (d2Fx*d1s**2 + dFx*d2s)*0.5d0 356#endif 357#ifdef THIRD_DERIV 358 rhom53 = rhom23*rrho 359 d3s = -1.5d0*d2s/gamma 360 361 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) - 4d0*fac* 362 & F23*F49*Ax*rhom53*(Fx - s*dFx + 18d0*s2*d2Fx +8d0*s2*s*d3Fx) 363 364 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) + 365 & F49*Ax*rhom23*d1s*s*(7d0*d2Fx + 4d0*d3Fx*s)*2d0*fac 366 367 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) - 368 & F43*Ax*rho13*(d2Fx*d1s**2 + d3Fx*d1s**2*s + 369 & d2Fx*d2s*s)*fac 370 371 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) + 372 & Ax*rho43*(d3Fx*d1s**3 + 3d0*d2Fx*d1s*d2s + 373 & dFx*d3s)*0.5d0*fac 374 375#endif 376c 377c Beta 378c 379 25 continue 380 if (rho(n,3).lt.tol_rho) goto 20 381 rho43 = (2d0*rho(n,3))**F43 382 rrho = 0.5d0/rho(n,3) 383 rho13 = rho43*rrho 384c 385 gamma = delrho(n,1,2)*delrho(n,1,2) + 386 & delrho(n,2,2)*delrho(n,2,2) + 387 & delrho(n,3,2)*delrho(n,3,2) 388 gam12 = 2d0*dsqrt(gamma) 389 if (gam12.le.tol_rho**2) goto 20 390c 391 s = gam12/(2d0*ckf*rho43) 392 s2 = s*s 393 394 d1s = 0.5d0*s/gamma 395 396 if (s.lt.1d-8) then 397 Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s + 398 & mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2 399 dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 + 400 & 4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s 401 d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s + 402 & 12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2 403 d3Fx = 6d0*alpha*zeta*mu + 404 & 24d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s 405 else 406 ops = 1d0 + s 407 omz = 1d0 - zeta 408 logops = dlog(ops) 409 410 asymp = 1d0 + alpha*(omz*s*logops + zeta*s) 411 dasymp = alpha*( (s + zeta)/ops + omz*logops ) 412 413 if (s.gt.30d0) then 414 tansin = arcsinh(s) 415 dtansin = darcsinh(s) 416 else 417 tansin = dtanh(s)*arcsinh(s) 418 dtansin = arcsinh(s)/dcosh(s)**2 + dtanh(s)*darcsinh(s) 419 endif 420 421 denom = 1d0 + beta*tansin 422 ddenom = beta*dtansin 423 424 factor = tansin/denom 425 dfactor = (dtansin - factor*ddenom)/denom 426 427 Fx = 1d0 + mu*factor*asymp 428 dFx = mu*(dfactor*asymp + factor*dasymp) 429 430#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 431 if (s.gt.30d0) then 432 d2tansin = d2arcsinh(s) 433 else 434 d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 + 435 & dtanh(s)*d2arcsinh(s) - 436 & 2d0*tansin/dcosh(s)**2 437 endif 438 d2asymp = alpha*omz*(2d0 + s)/ops**2 439 d2denom = beta*d2tansin 440 d2factor = (d2tansin - 2d0*dfactor*ddenom - 441 & factor*d2denom)/denom 442 d2Fx = mu*(d2factor*asymp + 2d0*dfactor*dasymp + 443 & factor*d2asymp) 444#endif 445#if defined(THIRD_DERIV) 446 447 if (s.gt.30d0) then 448 d3tansin = d3arcsinh(s) 449 else 450 d3tansin = 3d0*d2arcsinh(s)/dcosh(s)**2 - 451 & 4d0*darcsinh(s)*dtanh(s)/dcosh(s)**2 + 452 & dtanh(s)*d3arcsinh(s) - 453 & 2d0*dtansin/dcosh(s)**2 + 454 & 4d0*tansin*dtanh(s)/dcosh(s)**2 455 endif 456 d3asymp = -alpha*omz*(3d0 + s)/ops**3 457 d3denom = beta*d3tansin 458 d3factor = (d3tansin - 3d0*d2factor*ddenom - 459 & 3d0*dfactor*d2denom - factor*d3denom)/denom 460 d3Fx = mu*(d3factor*asymp + 3d0*d2factor*dasymp + 461 & 3d0*dfactor*d2asymp + factor*d3asymp) 462#endif 463 endif 464 465 Ex = Ex + 0.5d0*Ax*rho43*Fx*qwght(n)*fac 466 if (ldew) func(n) = func(n) + 0.5d0*Ax*rho43*Fx*fac 467 468 Amat(n,2) = Amat(n,2) + F43*Ax*rho13*(Fx - s*dFx)*fac 469 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 0.5d0*Ax*rho43*dFx*d1s*fac 470 471#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 472 d2s = -0.5d0*d1s/gamma 473 rhom23 = rho13*rrho 474 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 2d0*F49*Ax*rhom23* 475 & (Fx - s*dFx + 4d0*s2*d2Fx)*fac 476 477 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) - 478 & (F43*Ax*rho13*d2Fx*d1s*s)*fac 479 480 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + fac*Ax*rho43* 481 & (d2Fx*d1s**2 + dFx*d2s)*0.5d0 482#endif 483#ifdef THIRD_DERIV 484 rhom53 = rhom23*rrho 485 d3s = -1.5d0*d2s/gamma 486 487 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) - 4d0*fac* 488 & F23*F49*Ax*rhom53*(Fx - s*dFx + 18d0*s2*d2Fx +8d0*s2*s*d3Fx) 489 490 Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) + 491 & F49*Ax*rhom23*d1s*s*(7d0*d2Fx + 4d0*d3Fx*s)*2d0*fac 492 493 Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) - 494 & F43*Ax*rho13*(d2Fx*d1s**2 + d3Fx*d1s**2*s + 495 & d2Fx*d2s*s)*fac 496 497 Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) + 498 & Ax*rho43*(d3Fx*d1s**3 + 3d0*d2Fx*d1s*d2s + 499 & dFx*d3s)*0.5d0*fac 500#endif 501c 502 20 continue 503 endif 504c 505 return 506 end 507#ifndef SECOND_DERIV 508#define SECOND_DERIV 509c 510c Compile source again for the 2nd derivative case 511c 512#include "xc_xncap.F" 513#endif 514#ifndef THIRD_DERIV 515#define THIRD_DERIV 516c 517c Compile source again for the 3rd derivative case 518c 519#include "xc_xncap.F" 520#endif 521c $Id$ 522