1#if defined(FUJITSU_VPP) 2!ocl scalar 3#endif 4#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 5 Subroutine xc_becke88(tol_rho, fac, lfac, nlfac, rho, delrho, 6 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 7#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 8 Subroutine xc_becke88_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 9 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 10 & qwght,ldew,func) 11#else 12 Subroutine xc_becke88_d3(tol_rho, fac, lfac, nlfac, rho, delrho, 13 & Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3, 14 & nq, ipol, Ex, qwght, ldew, func) 15#endif 16c 17C$Id$ 18c 19 implicit none 20c 21#include "dft2drv.fh" 22#include "dft3drv.fh" 23c 24 double precision tol_rho, fac, Ex 25 integer nq, ipol 26 logical lfac, nlfac,ldew 27 double precision func(*) ! value of the functional [output] 28c 29c Charge Density 30c 31 double precision rho(nq,ipol*(ipol+1)/2) 32c 33c Charge Density Gradient 34c 35 double precision delrho(nq,3,ipol) 36c 37c Quadrature Weights 38c 39 double precision qwght(nq) 40c 41c Sampling Matrices for the XC Potential 42c 43 double precision Amat(nq,ipol), Cmat(nq,*) 44c 45#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 46c Second Derivatives of the Exchange Energy Functional 47 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 48#endif 49#ifdef THIRD_DERIV 50c Third Derivatives of the Exchange Energy Functional 51 double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3) 52 53#endif 54c 55 double precision BETA 56 Parameter (BETA = 0.0042D0) 57c 58c References: 59c 60c Becke, Phys. Rev. A 38, 3098 (1988) 61c Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993) 62c 63c*************************************************************************** 64c 65 integer n 66 double precision arcsinh, darcsinh, d2arcsinh 67 double precision C, rho13, rho43, gamma, x, g, gdenom, dg, 68 & dgdenom, t 69 double precision gdenom2 70#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 71 double precision rhom23, d2g, d2gdenom 72 double precision gdenom3 73#endif 74#ifdef THIRD_DERIV 75 double precision rhom53, d3g, d3gdenom 76 double precision gdenom4 77#endif 78c 79 arcsinh(x)=log(x+dsqrt(1d0+x*x)) 80 darcsinh(x)=1d0/dsqrt(1d0+x*x) 81 d2arcsinh(x) = -x/dsqrt(1d0+x*x)**3 82c 83c Uniform electron gas constant 84c 85 C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0) 86c 87 if (ipol.eq.1) then 88c 89c ======> SPIN-RESTRICTED <====== 90c 91 do 10 n = 1, nq 92 if (rho(n,1).lt.tol_rho) goto 10 93c 94c Spin alpha: 95c 96 rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0) 97 rho43 = rho13**4 98 gamma = delrho(n,1,1)*delrho(n,1,1) + 99 & delrho(n,2,1)*delrho(n,2,1) + 100 & delrho(n,3,1)*delrho(n,3,1) 101 if (dsqrt(gamma).gt.tol_rho)then 102 gamma = 0.25d0 * gamma 103 x = dsqrt(gamma) / rho43 104 else 105 x = 0d0 106 endif 107c 108 gdenom = 1d0 + 6d0*BETA*x*arcsinh(x) 109 gdenom2 = gdenom*gdenom 110 g = -BETA*x*x / gdenom 111 dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 112c dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2 113 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2 114c 115 if (lfac) then 116 Ex = Ex + 2d0*rho43*C*qwght(n)*fac 117 if(ldew)func(n) = func(n) + 2.d0*rho43*C*fac 118 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*C*fac 119 endif 120c 121 if (nlfac) then 122 Ex = Ex + 2d0*rho43*g*qwght(n)*fac 123 if(ldew)func(n) = func(n) + 2.d0*rho43*g*fac 124 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac 125 endif 126c 127 if (x.gt.tol_rho) then 128 t = 0.5d0 * dg / sqrt(gamma) * fac 129 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t 130c Cmat(n,1,1) = Cmat(n,1,1) + t*delrho(n,1,1) 131c Cmat(n,2,1) = Cmat(n,2,1) + t*delrho(n,2,1) 132c Cmat(n,3,1) = Cmat(n,3,1) + t*delrho(n,3,1) 133 endif 134c 135#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 136 if(lfac) g = g + C ! Add local contribution back to g 137 rhom23 = rho13 / (0.5d0*rho(n,1)) 138 d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0)) 139 gdenom3 = gdenom2*gdenom 140 d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2 141 & + BETA*x*x*d2gdenom/gdenom2 142 & - 2d0*BETA*x*x*(dgdenom)**2/gdenom3 143c d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2 144c & + BETA*x*x*d2gdenom/gdenom**2 145c & - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3 146c 147 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 148 & + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 149 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 150 & - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac 151 if (x.gt.tol_rho) then 152 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 153 & - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 154 endif 155#endif 156#ifdef THIRD_DERIV 157 rhom53 = rhom23 / (0.5d0*rho(n,1)) 158c 159 d3gdenom = 6.0d0*BETA* 160 1 d2arcsinh(x)*( 3.0d0 161 2 - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x)) 162c 163 gdenom4 = gdenom3*gdenom 164c 165 d3g = 6.0d0*BETA*dgdenom/gdenom2 166 1 - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3 167 2 + 6.0d0*BETA*x*d2gdenom/gdenom2 168 3 + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4 169 4 - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3 170 5 + BETA*x*x*d3gdenom/gdenom2 171c 172 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 173 1 + (8.0d0/27.0d0)*rhom53*(-g + x*dg 174 2 - 18.0d0*x*x*d2g 175 3 - 8.0d0*x*x*x*d3g)*fac 176c 177 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 178 1 + (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g 179 2 + 4.0d0*x*x*x*d3g)*fac 180c 181 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) 182 1 - (8.0d0/3.0d0)*(rhom23/rho(n,1)**3)/dsqrt(gamma) 183 2 *d3g*fac 184c 185 if (x.gt.tol_rho) then 186 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) 187 1 + (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg 188 2 - 3.0d0*x*d2g 189 3 + x*x*d3g)*fac 190 endif 191#endif 192c 193 10 continue 194c 195 else 196c 197c ======> SPIN-UNRESTRICTED <====== 198c 199 do 20 n = 1, nq 200 if (rho(n,1).lt.tol_rho) goto 20 201 if (rho(n,2).lt.tol_rho) goto 25 202c 203c Spin alpha: 204c 205 rho13 = rho(n,2)**(1.d0/3.d0) 206 rho43 = rho13*rho(n,2) 207 gamma = delrho(n,1,1)*delrho(n,1,1) + 208 & delrho(n,2,1)*delrho(n,2,1) + 209 & delrho(n,3,1)*delrho(n,3,1) 210 if (dsqrt(gamma).gt.tol_rho)then 211 x = dsqrt(gamma) / rho43 212 else 213 x = 0d0 214 endif 215c 216 gdenom = 1d0 + 6d0*BETA*x*arcsinh(x) 217 g = -BETA*x*x / gdenom 218 dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 219 gdenom2 = gdenom*gdenom 220 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2 221c 222 if (lfac) then 223 Ex = Ex + rho43*C*qwght(n)*fac 224 if (ldew)func(n) = func(n) + rho43*C*fac 225 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*C*fac 226 endif 227c 228 if (nlfac) then 229 Ex = Ex + rho43*g*qwght(n)*fac 230 if (ldew)func(n) = func(n) + rho43*g*fac 231 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac 232 endif 233c 234 if (x.gt.tol_rho) then 235 t = dg / sqrt(gamma) * fac 236 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t * 0.5d0 237c Cmat(n,1,1) = Cmat(n,1,1) + t*delrho(n,1,1) 238c Cmat(n,2,1) = Cmat(n,2,1) + t*delrho(n,2,1) 239c Cmat(n,3,1) = Cmat(n,3,1) + t*delrho(n,3,1) 240 endif 241c 242#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 243 if (lfac) g = g + C ! Add local contribution back to g 244 rhom23 = rho13 / rho(n,2) 245 d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0)) 246 gdenom3 = gdenom2*gdenom 247 d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2 248 & + BETA*x*x*d2gdenom/gdenom2 249 & - 2d0*BETA*x*x*(dgdenom)**2/gdenom3 250 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 251 & + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 252 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 253 & - (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac 254 if (x.gt.tol_rho) then 255 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 256 & - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 257 endif 258#endif 259#ifdef THIRD_DERIV 260 rhom53 = rhom23 / rho(n,2) 261 d3gdenom = 6.0d0*BETA* 262 1 d2arcsinh(x)*( 3.0d0 263 2 - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x)) 264 gdenom4 = gdenom3*gdenom 265 d3g = 6.0d0*BETA*dgdenom/gdenom2 266 1 - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3 267 2 + 6.0d0*BETA*x*d2gdenom/gdenom2 268 3 + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4 269 4 - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3 270 5 + BETA*x*x*d3gdenom/gdenom2 271c 272 Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) 273 1 + (8.0d0/27.0d0)*rhom53*(-g + x*dg 274 2 - 18.0d0*x*x*d2g 275 3 - 8.0d0*x*x*x*d3g)*fac 276c 277 Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) 278 1 + (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g 279 2 + 4.0d0*x*x*x*d3g)*fac 280c 281 Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) 282 1 - (1.0d0/3.0d0)*(rhom23/rho(n,2)**3)/dsqrt(gamma) 283 2 *d3g*fac 284c 285 if (x.gt.tol_rho) then 286 Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) 287 1 + (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg 288 2 - 3.0d0*x*d2g 289 3 + x*x*d3g)*fac 290 endif 291#endif 292c 293 25 continue 294c 295c Spin beta: 296c 297 if (rho(n,3).lt.tol_rho) goto 20 298c 299 rho13 = rho(n,3)**(1.d0/3.d0) 300 rho43 = rho13*rho(n,3) 301 gamma = delrho(n,1,2)*delrho(n,1,2) + 302 & delrho(n,2,2)*delrho(n,2,2) + 303 & delrho(n,3,2)*delrho(n,3,2) 304 if (dsqrt(gamma).gt.tol_rho)then 305 x = dsqrt(gamma) / rho43 306 else 307 x = 0d0 308 endif 309c 310 gdenom = 1d0 + 6d0*BETA*x*arcsinh(x) 311 g = -BETA*x*x / gdenom 312 dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 313 gdenom2 = gdenom*gdenom 314 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom2 315c 316 if (lfac) then 317 Ex = Ex + rho43*C*qwght(n)*fac 318 if (ldew)func(n) = func(n) + rho43*C*fac 319 Amat(n,2) = Amat(n,2) + (4d0/3d0)*rho13*C*fac 320 endif 321c 322 if (nlfac) then 323 Ex = Ex + rho43*g*qwght(n)*fac 324 if (ldew)func(n) = func(n) +rho43*g*fac 325 Amat(n,2) = Amat(n,2) + (4d0/3d0)*rho13*(g-x*dg)*fac 326 endif 327c 328 if (x.gt.tol_rho) then 329 t = dg / sqrt(gamma) * fac 330 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + t * 0.5d0 331c Cmat(n,1,2) = Cmat(n,1,2) + t*delrho(n,1,2) 332c Cmat(n,2,2) = Cmat(n,2,2) + t*delrho(n,2,2) 333c Cmat(n,3,2) = Cmat(n,3,2) + t*delrho(n,3,2) 334 endif 335c 336#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 337 if(lfac) g = g + C ! Add local contribution back to g 338 rhom23 = rho13 / rho(n,3) 339 d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0)) 340 gdenom3 = gdenom2*gdenom 341 d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom2 342 & + BETA*x*x*d2gdenom/gdenom2 343 & - 2d0*BETA*x*x*(dgdenom)**2/gdenom3 344 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 345 & + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 346 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 347 & - (2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac 348 if (x.gt.tol_rho) then 349 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 350 & - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 351 endif 352#endif 353#ifdef THIRD_DERIV 354 rhom53 = rhom23 / rho(n,3) 355c 356 d3gdenom = 6.0d0*BETA* 357 1 d2arcsinh(x)*( 3.0d0 358 2 - (2.0d0*x*x - 1.0d0)/(1.0d0 + x*x)) 359c 360 gdenom4 = gdenom3*gdenom 361c 362 d3g = 6.0d0*BETA*dgdenom/gdenom2 363 1 - 12.0d0*BETA*x*dgdenom*dgdenom/gdenom3 364 2 + 6.0d0*BETA*x*d2gdenom/gdenom2 365 3 + 6.0d0*BETA*x*x*dgdenom*dgdenom*dgdenom/gdenom4 366 4 - 6.0d0*BETA*x*x*dgdenom*d2gdenom/gdenom3 367 5 + BETA*x*x*d3gdenom/gdenom2 368c 369 Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) 370 1 + (8.0d0/27.0d0)*rhom53*(-g + x*dg 371 2 - 18.0d0*x*x*d2g 372 3 - 8.0d0*x*x*x*d3g)*fac 373c 374 Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) 375 1 + (2.0d0/9.0d0)*(rhom23/gamma)*( 7.0d0*x*x*d2g 376 2 + 4.0d0*x*x*x*d3g)*fac 377c 378 Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) 379 1 - (1.0d0/3.0d0)*(rhom23/rho(n,3)**3)/dsqrt(gamma) 380 2 *d3g*fac 381c 382 if (x.gt.tol_rho) then 383 Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) 384 1 + (1.0d0/8.0d0)*gamma**(-2.5d0)*( 3.0d0*dg 385 2 - 3.0d0*x*d2g 386 3 + x*x*d3g)*fac 387 endif 388#endif 389c 390 20 continue 391c 392 endif 393c 394 return 395 end 396#ifndef SECOND_DERIV 397#define SECOND_DERIV 398c 399c Compile source again for the 2nd derivative case 400c 401#include "xc_becke88.F" 402#endif 403#ifndef THIRD_DERIV 404#define THIRD_DERIV 405c 406c Compile source again for the 3rd derivative case 407c 408#include "xc_becke88.F" 409#endif 410