1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Integral GKS scheme: The order of the integrals in makeCoul reflects 8!> the standard order by MOPAC 9!> \par History 10!> Teodoro Laino [tlaino] - 04.2009 : Adapted size arrays to d-orbitals and 11!> get rid of the alternative ordering. Using the 12!> CP2K one. 13!> Teodoro Laino [tlaino] - 04.2009 : Skip nullification (speed-up) 14!> Teodoro Laino [tlaino] - 04.2009 : Speed-up due to fortran arrays order 15!> optimization and collection of common pieces of 16!> code 17! ************************************************************************************************** 18MODULE semi_empirical_int_gks 19 20 USE dg_rho0_types, ONLY: dg_rho0_type 21 USE dg_types, ONLY: dg_get,& 22 dg_type 23 USE kinds, ONLY: dp 24 USE mathconstants, ONLY: fourpi,& 25 oorootpi 26 USE multipole_types, ONLY: do_multipole_none 27 USE pw_grid_types, ONLY: pw_grid_type 28 USE pw_pool_types, ONLY: pw_pool_type 29 USE semi_empirical_int_arrays, ONLY: indexb,& 30 rij_threshold 31 USE semi_empirical_mpole_types, ONLY: semi_empirical_mpole_type 32 USE semi_empirical_types, ONLY: se_int_control_type,& 33 semi_empirical_type,& 34 setup_se_int_control_type 35#include "./base/base_uses.f90" 36 37 IMPLICIT NONE 38 PRIVATE 39 40 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_gks' 41 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 42 43 PUBLIC :: corecore_gks, rotnuc_gks, drotnuc_gks, rotint_gks, drotint_gks 44 45CONTAINS 46 47! ************************************************************************************************** 48!> \brief Computes the electron-nuclei integrals 49!> \param sepi ... 50!> \param sepj ... 51!> \param rij ... 52!> \param e1b ... 53!> \param e2a ... 54!> \param se_int_control ... 55! ************************************************************************************************** 56 SUBROUTINE rotnuc_gks(sepi, sepj, rij, e1b, e2a, se_int_control) 57 TYPE(semi_empirical_type), POINTER :: sepi, sepj 58 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 59 REAL(dp), DIMENSION(45), INTENT(OUT), OPTIONAL :: e1b, e2a 60 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 61 62 CHARACTER(len=*), PARAMETER :: routineN = 'rotnuc_gks', routineP = moduleN//':'//routineN 63 64 INTEGER :: i, mu, nu 65 REAL(KIND=dp), DIMENSION(3) :: rab 66 REAL(kind=dp), DIMENSION(45, 45) :: Coul 67 68 rab = -rij 69 70 IF (se_int_control%do_ewald_gks) THEN 71 IF (DOT_PRODUCT(rij, rij) > rij_threshold) THEN 72 CALL makeCoulE(rab, sepi, sepj, Coul, se_int_control) 73 ELSE 74 CALL makeCoulE0(sepi, Coul, se_int_control) 75 END IF 76 ELSE 77 CALL makeCoul(rab, sepi, sepj, Coul, se_int_control) 78 END IF 79 80 i = 0 81 DO mu = 1, sepi%natorb 82 DO nu = 1, mu 83 i = i + 1 84 e1b(i) = -Coul(i, 1)*sepj%zeff 85 END DO 86 END DO 87 88 i = 0 89 DO mu = 1, sepj%natorb 90 DO nu = 1, mu 91 i = i + 1 92 e2a(i) = -Coul(1, i)*sepi%zeff 93 END DO 94 END DO 95 96 END SUBROUTINE rotnuc_gks 97 98! ************************************************************************************************** 99!> \brief Computes the electron-electron integrals 100!> \param sepi ... 101!> \param sepj ... 102!> \param rij ... 103!> \param w ... 104!> \param se_int_control ... 105! ************************************************************************************************** 106 SUBROUTINE rotint_gks(sepi, sepj, rij, w, se_int_control) 107 TYPE(semi_empirical_type), POINTER :: sepi, sepj 108 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 109 REAL(dp), DIMENSION(2025), INTENT(OUT), OPTIONAL :: w 110 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 111 112 CHARACTER(len=*), PARAMETER :: routineN = 'rotint_gks', routineP = moduleN//':'//routineN 113 114 INTEGER :: i, ind1, ind2, lam, mu, nu, sig 115 REAL(KIND=dp), DIMENSION(3) :: rab 116 REAL(kind=dp), DIMENSION(45, 45) :: Coul 117 118 rab = -rij 119 120 IF (se_int_control%do_ewald_gks) THEN 121 IF (DOT_PRODUCT(rij, rij) > rij_threshold) THEN 122 CALL makeCoulE(rab, sepi, sepj, Coul, se_int_control) 123 ELSE 124 CALL makeCoulE0(sepi, Coul, se_int_control) 125 END IF 126 ELSE 127 CALL makeCoul(rab, sepi, sepj, Coul, se_int_control) 128 END IF 129 130 i = 0 131 ind1 = 0 132 DO mu = 1, sepi%natorb 133 DO nu = 1, mu 134 ind1 = ind1 + 1 135 ind2 = 0 136 DO lam = 1, sepj%natorb 137 DO sig = 1, lam 138 i = i + 1 139 ind2 = ind2 + 1 140 w(i) = Coul(ind1, ind2) 141 END DO 142 END DO 143 END DO 144 END DO 145 146 END SUBROUTINE rotint_gks 147 148! ************************************************************************************************** 149!> \brief Computes the derivatives of the electron-nuclei integrals 150!> \param sepi ... 151!> \param sepj ... 152!> \param rij ... 153!> \param de1b ... 154!> \param de2a ... 155!> \param se_int_control ... 156! ************************************************************************************************** 157 SUBROUTINE drotnuc_gks(sepi, sepj, rij, de1b, de2a, se_int_control) 158 TYPE(semi_empirical_type), POINTER :: sepi, sepj 159 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 160 REAL(dp), DIMENSION(3, 45), INTENT(OUT), OPTIONAL :: de1b, de2a 161 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 162 163 CHARACTER(len=*), PARAMETER :: routineN = 'drotnuc_gks', routineP = moduleN//':'//routineN 164 165 INTEGER :: i, mu, nu 166 REAL(KIND=dp), DIMENSION(3) :: rab 167 REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul 168 169 rab = -rij 170 171 IF (se_int_control%do_ewald_gks) THEN 172 CALL makedCoulE(rab, sepi, sepj, dCoul, se_int_control) 173 ELSE 174 CALL makedCoul(rab, sepi, sepj, dCoul, se_int_control) 175 END IF 176 177 i = 0 178 DO mu = 1, sepi%natorb 179 DO nu = 1, mu 180 i = i + 1 181 de1b(1, i) = dCoul(1, i, 1)*sepj%zeff 182 de1b(2, i) = dCoul(2, i, 1)*sepj%zeff 183 de1b(3, i) = dCoul(3, i, 1)*sepj%zeff 184 END DO 185 END DO 186 187 i = 0 188 DO mu = 1, sepj%natorb 189 DO nu = 1, mu 190 i = i + 1 191 de2a(1, i) = dCoul(1, 1, i)*sepi%zeff 192 de2a(2, i) = dCoul(2, 1, i)*sepi%zeff 193 de2a(3, i) = dCoul(3, 1, i)*sepi%zeff 194 END DO 195 END DO 196 197 END SUBROUTINE drotnuc_gks 198 199! ************************************************************************************************** 200!> \brief Computes the derivatives of the electron-electron integrals 201!> \param sepi ... 202!> \param sepj ... 203!> \param rij ... 204!> \param dw ... 205!> \param se_int_control ... 206! ************************************************************************************************** 207 SUBROUTINE drotint_gks(sepi, sepj, rij, dw, se_int_control) 208 TYPE(semi_empirical_type), POINTER :: sepi, sepj 209 REAL(dp), DIMENSION(3), INTENT(IN) :: rij 210 REAL(dp), DIMENSION(3, 2025), INTENT(OUT) :: dw 211 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 212 213 CHARACTER(len=*), PARAMETER :: routineN = 'drotint_gks', routineP = moduleN//':'//routineN 214 215 INTEGER :: i, ind1, ind2, lam, mu, nu, sig 216 REAL(KIND=dp), DIMENSION(3) :: rab 217 REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul 218 219 rab = -rij 220 221 IF (se_int_control%do_ewald_gks) THEN 222 CALL makedCoulE(rab, sepi, sepj, dCoul, se_int_control) 223 ELSE 224 CALL makedCoul(rab, sepi, sepj, dCoul, se_int_control) 225 END IF 226 227 i = 0 228 ind1 = 0 229 DO mu = 1, sepi%natorb 230 DO nu = 1, mu 231 ind1 = ind1 + 1 232 ind2 = 0 233 DO lam = 1, sepj%natorb 234 DO sig = 1, lam 235 i = i + 1 236 ind2 = ind2 + 1 237 dw(1, i) = -dCoul(1, ind1, ind2) 238 dw(2, i) = -dCoul(2, ind1, ind2) 239 dw(3, i) = -dCoul(3, ind1, ind2) 240 END DO 241 END DO 242 END DO 243 END DO 244 245 END SUBROUTINE drotint_gks 246 247! ************************************************************************************************** 248!> \brief Computes the primitives of the integrals (non-periodic case) 249!> \param RAB ... 250!> \param sepi ... 251!> \param sepj ... 252!> \param Coul ... 253!> \param se_int_control ... 254! ************************************************************************************************** 255 SUBROUTINE makeCoul(RAB, sepi, sepj, Coul, se_int_control) 256 REAL(kind=dp), DIMENSION(3) :: RAB 257 TYPE(semi_empirical_type), POINTER :: sepi, sepj 258 REAL(kind=dp), DIMENSION(45, 45), INTENT(OUT) :: Coul 259 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 260 261 CHARACTER(len=*), PARAMETER :: routineN = 'makeCoul', routineP = moduleN//':'//routineN 262 263 INTEGER :: iA, iB, imA, imB, jA, jB, k1, k2, k3, k4 264 LOGICAL :: shortrange 265 REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, d1f(3), d2, & 266 d2f(3, 3), d3, d3f(3, 3, 3), d4, & 267 d4f(3, 3, 3, 3), f, rr, w, w0, w1, w2, & 268 w3, w4, w5 269 REAL(kind=dp), DIMENSION(3) :: v 270 REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B 271 REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B 272 REAL(kind=dp), DIMENSION(45) :: M0A, M0B 273 274 shortrange = se_int_control%shortrange 275 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 276 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 277 278 v(1) = RAB(1) 279 v(2) = RAB(2) 280 v(3) = RAB(3) 281 rr = SQRT(DOT_PRODUCT(v, v)) 282 283 a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB) 284 w0 = a2*rr 285 w = EXP(-w0) 286 w1 = (1.0_dp + 0.5_dp*w0) 287 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2) 288 w3 = (w2 + w0**3/6.0_dp) 289 w4 = (w3 + w0**4/30.0_dp) 290 w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp) 291 292 IF (shortrange) THEN 293 f = (-w*w1)/rr 294 d1 = -1.0_dp*(-w*w2)/rr**3 295 d2 = 3.0_dp*(-w*w3)/rr**5 296 d3 = -15.0_dp*(-w*w4)/rr**7 297 d4 = 105.0_dp*(-w*w5)/rr**9 298 ELSE 299 f = (1.0_dp - w*w1)/rr 300 d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3 301 d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5 302 d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7 303 d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9 304 ENDIF 305 306 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4) 307 308 imA = 0 309 DO iA = 1, sepi%natorb 310 DO jA = 1, iA 311 imA = imA + 1 312 313 imB = 0 314 DO iB = 1, sepj%natorb 315 DO jB = 1, iB 316 imB = imB + 1 317 318 w = M0A(imA)*M0B(imB)*f 319 DO k1 = 1, 3 320 w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1) 321 ENDDO 322 DO k2 = 1, 3 323 DO k1 = 1, 3 324 w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2) 325 ENDDO 326 ENDDO 327 DO k3 = 1, 3 328 DO k2 = 1, 3 329 DO k1 = 1, 3 330 w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3) 331 ENDDO 332 ENDDO 333 ENDDO 334 335 DO k4 = 1, 3 336 DO k3 = 1, 3 337 DO k2 = 1, 3 338 DO k1 = 1, 3 339 w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4) 340 ENDDO 341 ENDDO 342 ENDDO 343 ENDDO 344 345 Coul(imA, imB) = w 346 ENDDO 347 ENDDO 348 ENDDO 349 ENDDO 350 351 END SUBROUTINE makeCoul 352 353! ************************************************************************************************** 354!> \brief Computes the derivatives of the primitives of the integrals 355!> (non-periodic case) 356!> \param RAB ... 357!> \param sepi ... 358!> \param sepj ... 359!> \param dCoul ... 360!> \param se_int_control ... 361! ************************************************************************************************** 362 SUBROUTINE makedCoul(RAB, sepi, sepj, dCoul, se_int_control) 363 REAL(kind=dp), DIMENSION(3) :: RAB 364 TYPE(semi_empirical_type), POINTER :: sepi, sepj 365 REAL(kind=dp), DIMENSION(3, 45, 45), INTENT(OUT) :: dCoul 366 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 367 368 CHARACTER(len=*), PARAMETER :: routineN = 'makedCoul', routineP = moduleN//':'//routineN 369 370 INTEGER :: iA, iB, imA, imB, jA, jB, k1, k2, k3, k4 371 LOGICAL :: shortrange 372 REAL(kind=dp) :: a2, ACOULA, ACOULB, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), d4, & 373 d4f(3, 3, 3, 3), d5, d5f(3, 3, 3, 3, 3), f, rr, tmp, w, w0, w1, w2, w3, w4, w5, w6 374 REAL(kind=dp), DIMENSION(3) :: v, wv 375 REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B 376 REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B 377 REAL(kind=dp), DIMENSION(45) :: M0A, M0B 378 379 shortrange = se_int_control%shortrange 380 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 381 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 382 383 v(1) = RAB(1) 384 v(2) = RAB(2) 385 v(3) = RAB(3) 386 rr = SQRT(DOT_PRODUCT(v, v)) 387 388 a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB) 389 w0 = a2*rr 390 w = EXP(-w0) 391 w1 = (1.0_dp + 0.5_dp*w0) 392 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2) 393 w3 = (w2 + w0**3/6.0_dp) 394 w4 = (w3 + w0**4/30.0_dp) 395 w5 = (w3 + 4.0_dp*w0**4/105.0_dp + w0**5/210.0_dp) 396 w6 = (w3 + 15.0_dp*w0**4/378.0_dp + 2.0_dp*w0**5/315.0_dp + w0**6/1890.0_dp) 397 398 IF (shortrange) THEN 399 f = (-w*w1)/rr 400 d1 = -1.0_dp*(-w*w2)/rr**3 401 d2 = 3.0_dp*(-w*w3)/rr**5 402 d3 = -15.0_dp*(-w*w4)/rr**7 403 d4 = 105.0_dp*(-w*w5)/rr**9 404 d5 = -945.0_dp*(-w*w6)/rr**11 405 ELSE 406 f = (1.0_dp - w*w1)/rr 407 d1 = -1.0_dp*(1.0_dp - w*w2)/rr**3 408 d2 = 3.0_dp*(1.0_dp - w*w3)/rr**5 409 d3 = -15.0_dp*(1.0_dp - w*w4)/rr**7 410 d4 = 105.0_dp*(1.0_dp - w*w5)/rr**9 411 d5 = -945.0_dp*(1.0_dp - w*w6)/rr**11 412 ENDIF 413 414 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5) 415 416 imA = 0 417 DO iA = 1, sepi%natorb 418 DO jA = 1, iA 419 imA = imA + 1 420 421 imB = 0 422 DO iB = 1, sepj%natorb 423 DO jB = 1, iB 424 imB = imB + 1 425 426 tmp = M0A(imA)*M0B(imB) 427 wv(1) = tmp*d1f(1) 428 wv(2) = tmp*d1f(2) 429 wv(3) = tmp*d1f(3) 430 DO k1 = 1, 3 431 tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB) 432 wv(1) = wv(1) + tmp*d2f(1, k1) 433 wv(2) = wv(2) + tmp*d2f(2, k1) 434 wv(3) = wv(3) + tmp*d2f(3, k1) 435 ENDDO 436 DO k2 = 1, 3 437 DO k1 = 1, 3 438 tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB) 439 wv(1) = wv(1) + tmp*d3f(1, k1, k2) 440 wv(2) = wv(2) + tmp*d3f(2, k1, k2) 441 wv(3) = wv(3) + tmp*d3f(3, k1, k2) 442 ENDDO 443 ENDDO 444 DO k3 = 1, 3 445 DO k2 = 1, 3 446 DO k1 = 1, 3 447 tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB) 448 wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3) 449 wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3) 450 wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3) 451 ENDDO 452 ENDDO 453 ENDDO 454 455 DO k4 = 1, 3 456 DO k3 = 1, 3 457 DO k2 = 1, 3 458 DO k1 = 1, 3 459 tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB) 460 wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4) 461 wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4) 462 wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4) 463 ENDDO 464 ENDDO 465 ENDDO 466 ENDDO 467 468 dCoul(1, imA, imB) = wv(1) 469 dCoul(2, imA, imB) = wv(2) 470 dCoul(3, imA, imB) = wv(3) 471 ENDDO 472 ENDDO 473 ENDDO 474 ENDDO 475 476 END SUBROUTINE makedCoul 477 478! ************************************************************************************************** 479!> \brief Computes nuclei-nuclei interactions 480!> \param sepi ... 481!> \param sepj ... 482!> \param rijv ... 483!> \param enuc ... 484!> \param denuc ... 485!> \param se_int_control ... 486! ************************************************************************************************** 487 SUBROUTINE corecore_gks(sepi, sepj, rijv, enuc, denuc, se_int_control) 488 TYPE(semi_empirical_type), POINTER :: sepi, sepj 489 REAL(dp), DIMENSION(3), INTENT(IN) :: rijv 490 REAL(dp), INTENT(OUT), OPTIONAL :: enuc 491 REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: denuc 492 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 493 494 CHARACTER(len=*), PARAMETER :: routineN = 'corecore_gks', routineP = moduleN//':'//routineN 495 496 LOGICAL :: l_denuc, l_enuc 497 REAL(dp) :: alpi, alpj, dscale, rij, scale, zz 498 REAL(kind=dp), DIMENSION(3, 45, 45) :: dCoul, dCoulE 499 REAL(kind=dp), DIMENSION(45, 45) :: Coul, CoulE 500 TYPE(se_int_control_type) :: se_int_control_off 501 502 rij = DOT_PRODUCT(rijv, rijv) 503 504 l_enuc = PRESENT(enuc) 505 l_denuc = PRESENT(denuc) 506 IF ((rij > rij_threshold) .AND. (l_enuc .OR. l_denuc)) THEN 507 508 rij = SQRT(rij) 509 510 IF (se_int_control%shortrange) THEN 511 CALL setup_se_int_control_type(se_int_control_off, shortrange=.FALSE., do_ewald_r3=.FALSE., & 512 do_ewald_gks=.FALSE., integral_screening=se_int_control%integral_screening, & 513 max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.) 514 CALL makeCoul(rijv, sepi, sepj, Coul, se_int_control_off) 515 IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoul, se_int_control_off) 516 IF (se_int_control%do_ewald_gks) THEN 517 CALL makeCoulE(rijv, sepi, sepj, CoulE, se_int_control) 518 IF (l_denuc) CALL makedCoulE(rijv, sepi, sepj, dCoulE, se_int_control) 519 ELSE 520 CALL makeCoul(rijv, sepi, sepj, CoulE, se_int_control) 521 IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoulE, se_int_control) 522 END IF 523 ELSE 524 CALL makeCoul(rijv, sepi, sepj, Coul, se_int_control) 525 CoulE = Coul 526 IF (l_denuc) CALL makedCoul(rijv, sepi, sepj, dCoul, se_int_control) 527 IF (l_denuc) dCoulE = dCoul 528 END IF 529 530 scale = 0.0_dp 531 dscale = 0.0_dp 532 zz = sepi%zeff*sepj%zeff 533 alpi = sepi%alp 534 alpj = sepj%alp 535 scale = EXP(-alpi*rij) + EXP(-alpj*rij) 536 IF (l_enuc) THEN 537 enuc = zz*CoulE(1, 1) + scale*zz*Coul(1, 1) 538 END IF 539 IF (l_denuc) THEN 540 dscale = -alpi*EXP(-alpi*rij) - alpj*EXP(-alpj*rij) 541 denuc(1) = zz*dCoulE(1, 1, 1) + dscale*(rijv(1)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(1, 1, 1) 542 denuc(2) = zz*dCoulE(2, 1, 1) + dscale*(rijv(2)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(2, 1, 1) 543 denuc(3) = zz*dCoulE(3, 1, 1) + dscale*(rijv(3)/rij)*zz*Coul(1, 1) + scale*zz*dCoul(3, 1, 1) 544 END IF 545 546 ELSE 547 548 IF (se_int_control%do_ewald_gks) THEN 549 zz = sepi%zeff*sepi%zeff 550 CALL makeCoulE0(sepi, CoulE, se_int_control) 551 IF (l_enuc) THEN 552 enuc = zz*CoulE(1, 1) 553 END IF 554 IF (l_denuc) THEN 555 denuc = 0._dp 556 END IF 557 END IF 558 559 ENDIF 560 END SUBROUTINE corecore_gks 561 562! ************************************************************************************************** 563!> \brief Computes the primitives of the integrals (periodic case) 564!> \param RAB ... 565!> \param sepi ... 566!> \param sepj ... 567!> \param Coul ... 568!> \param se_int_control ... 569! ************************************************************************************************** 570 SUBROUTINE makeCoulE(RAB, sepi, sepj, Coul, se_int_control) 571 REAL(KIND=dp), DIMENSION(3) :: RAB 572 TYPE(semi_empirical_type), POINTER :: sepi, sepj 573 REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT) :: Coul 574 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 575 576 CHARACTER(len=*), PARAMETER :: routineN = 'makeCoulE', routineP = moduleN//':'//routineN 577 578 INTEGER :: gpt, imA, imB, k1, k2, k3, k4, lp, mp, np 579 INTEGER, DIMENSION(:, :), POINTER :: bds 580 REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), & 581 d4, d4f(3, 3, 3, 3), f, ff, kr, kr2, r1, r2, r3, r5, r7, r9, rr, ss, w, w0, w1, w2, w3, & 582 w4, w5 583 REAL(KIND=dp), DIMENSION(3) :: kk, v 584 REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2A, M2B 585 REAL(KIND=dp), DIMENSION(3, 45) :: M1A, M1B 586 REAL(KIND=dp), DIMENSION(45) :: M0A, M0B 587 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0 588 TYPE(dg_rho0_type), POINTER :: dg_rho0 589 TYPE(dg_type), POINTER :: dg 590 TYPE(pw_grid_type), POINTER :: pw_grid 591 TYPE(pw_pool_type), POINTER :: pw_pool 592 593 alpha = se_int_control%ewald_gks%alpha 594 pw_pool => se_int_control%ewald_gks%pw_pool 595 dg => se_int_control%ewald_gks%dg 596 CALL dg_get(dg, dg_rho0=dg_rho0) 597 rho0 => dg_rho0%density%pw%cr3d 598 pw_grid => pw_pool%pw_grid 599 bds => pw_grid%bounds 600 601 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 602 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 603 604 v(1) = RAB(1) 605 v(2) = RAB(2) 606 v(3) = RAB(3) 607 rr = SQRT(DOT_PRODUCT(v, v)) 608 609 r1 = 1.0_dp/rr 610 r2 = r1*r1 611 r3 = r2*r1 612 r5 = r3*r2 613 r7 = r5*r2 614 r9 = r7*r2 615 616 a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB) 617 618 w0 = a2*rr 619 w = EXP(-w0) 620 w1 = (1.0_dp + 0.5_dp*w0) 621 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2) 622 w3 = (w2 + w0**3/6.0_dp) 623 w4 = (w3 + w0**4/30.0_dp) 624 w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp) 625 626 f = (1.0_dp - w*w1)*r1 627 d1 = -1.0_dp*(1.0_dp - w*w2)*r3 628 d2 = 3.0_dp*(1.0_dp - w*w3)*r5 629 d3 = -15.0_dp*(1.0_dp - w*w4)*r7 630 d4 = 105.0_dp*(1.0_dp - w*w5)*r9 631 632 kr = alpha*rr 633 kr2 = kr*kr 634 w0 = 1.0_dp - erfc(kr) 635 w1 = 2.0_dp*oorootpi*EXP(-kr2) 636 w2 = w1*kr 637 638 f = f - w0*r1 639 d1 = d1 + (-w2 + w0)*r3 640 d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5 641 d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7 642 d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9 643 644 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, v=v, d1=d1, d2=d2, d3=d3, d4=d4) 645 646 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 647 DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2 648 649 w = M0A(imA)*M0B(imB)*f 650 DO k1 = 1, 3 651 w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1) 652 ENDDO 653 DO k2 = 1, 3 654 DO k1 = 1, 3 655 w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2) 656 ENDDO 657 ENDDO 658 DO k3 = 1, 3 659 DO k2 = 1, 3 660 DO k1 = 1, 3 661 w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3) 662 ENDDO 663 ENDDO 664 ENDDO 665 666 DO k4 = 1, 3 667 DO k3 = 1, 3 668 DO k2 = 1, 3 669 DO k1 = 1, 3 670 w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4) 671 ENDDO 672 ENDDO 673 ENDDO 674 ENDDO 675 676 Coul(imA, imB) = w 677 ENDDO 678 ENDDO 679 680 v(1) = RAB(1) 681 v(2) = RAB(2) 682 v(3) = RAB(3) 683 684 f = 0.0_dp 685 d1f = 0.0_dp 686 d2f = 0.0_dp 687 d3f = 0.0_dp 688 d4f = 0.0_dp 689 690 DO gpt = 1, pw_grid%ngpts_cut_local 691 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt)) 692 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt)) 693 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt)) 694 695 lp = lp + bds(1, 1) 696 mp = mp + bds(1, 2) 697 np = np + bds(1, 3) 698 699 IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE 700 kk(:) = pw_grid%g(:, gpt) 701 ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt) 702 703 kr = DOT_PRODUCT(kk, v) 704 cc = COS(kr) 705 ss = SIN(kr) 706 707 f = f + cc*ff 708 DO k1 = 1, 3 709 d1f(k1) = d1f(k1) - kk(k1)*ss*ff 710 ENDDO 711 DO k2 = 1, 3 712 DO k1 = 1, 3 713 d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff 714 ENDDO 715 ENDDO 716 DO k3 = 1, 3 717 DO k2 = 1, 3 718 DO k1 = 1, 3 719 d3f(k1, k2, k3) = d3f(k1, k2, k3) + kk(k1)*kk(k2)*kk(k3)*ss*ff 720 ENDDO 721 ENDDO 722 ENDDO 723 DO k4 = 1, 3 724 DO k3 = 1, 3 725 DO k2 = 1, 3 726 DO k1 = 1, 3 727 d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff 728 ENDDO 729 ENDDO 730 ENDDO 731 ENDDO 732 733 ENDDO 734 735 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 736 DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2 737 738 w = M0A(imA)*M0B(imB)*f 739 DO k1 = 1, 3 740 w = w + (M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB))*d1f(k1) 741 ENDDO 742 DO k2 = 1, 3 743 DO k1 = 1, 3 744 w = w + (M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB))*d2f(k1, k2) 745 ENDDO 746 ENDDO 747 DO k3 = 1, 3 748 DO k2 = 1, 3 749 DO k1 = 1, 3 750 w = w + (-M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB))*d3f(k1, k2, k3) 751 ENDDO 752 ENDDO 753 ENDDO 754 755 DO k4 = 1, 3 756 DO k3 = 1, 3 757 DO k2 = 1, 3 758 DO k1 = 1, 3 759 w = w + M2A(k1, k2, imA)*M2B(k3, k4, imB)*d4f(k1, k2, k3, k4) 760 ENDDO 761 ENDDO 762 ENDDO 763 ENDDO 764 765 Coul(imA, imB) = Coul(imA, imB) + w 766 767 ENDDO 768 ENDDO 769 770 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 771 DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2 772 w = -M0A(imA)*M0B(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2) 773 Coul(imA, imB) = Coul(imA, imB) + w 774 ENDDO 775 ENDDO 776 777 END SUBROUTINE makeCoulE 778 779! ************************************************************************************************** 780!> \brief Computes the derivatives of the primitives of the integrals 781!> (periodic case) 782!> \param RAB ... 783!> \param sepi ... 784!> \param sepj ... 785!> \param dCoul ... 786!> \param se_int_control ... 787! ************************************************************************************************** 788 SUBROUTINE makedCoulE(RAB, sepi, sepj, dCoul, se_int_control) 789 REAL(KIND=dp), DIMENSION(3) :: RAB 790 TYPE(semi_empirical_type), POINTER :: sepi, sepj 791 REAL(KIND=dp), DIMENSION(3, 45, 45), INTENT(OUT) :: dCoul 792 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 793 794 CHARACTER(len=*), PARAMETER :: routineN = 'makedCoulE', routineP = moduleN//':'//routineN 795 796 INTEGER :: gpt, imA, imB, k1, k2, k3, k4, k5, lp, & 797 mp, np 798 INTEGER, DIMENSION(:, :), POINTER :: bds 799 REAL(KIND=dp) :: a2, ACOULA, ACOULB, alpha, cc, d1, d1f(3), d2, d2f(3, 3), d3, d3f(3, 3, 3), & 800 d4, d4f(3, 3, 3, 3), d5, d5f(3, 3, 3, 3, 3), f, ff, kr, kr2, r1, r11, r2, r3, r5, r7, r9, & 801 rr, ss, tmp, w, w0, w1, w2, w3, w4, w5, w6 802 REAL(KIND=dp), DIMENSION(3) :: kk, v, wv 803 REAL(kind=dp), DIMENSION(3, 3, 45) :: M2A, M2B 804 REAL(kind=dp), DIMENSION(3, 45) :: M1A, M1B 805 REAL(kind=dp), DIMENSION(45) :: M0A, M0B 806 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0 807 TYPE(dg_rho0_type), POINTER :: dg_rho0 808 TYPE(dg_type), POINTER :: dg 809 TYPE(pw_grid_type), POINTER :: pw_grid 810 TYPE(pw_pool_type), POINTER :: pw_pool 811 812 alpha = se_int_control%ewald_gks%alpha 813 pw_pool => se_int_control%ewald_gks%pw_pool 814 dg => se_int_control%ewald_gks%dg 815 CALL dg_get(dg, dg_rho0=dg_rho0) 816 rho0 => dg_rho0%density%pw%cr3d 817 pw_grid => pw_pool%pw_grid 818 bds => pw_grid%bounds 819 820 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A, ACOULA) 821 CALL get_se_slater_multipole(sepj, M0B, M1B, M2B, ACOULB) 822 823 v(1) = RAB(1) 824 v(2) = RAB(2) 825 v(3) = RAB(3) 826 rr = SQRT(DOT_PRODUCT(v, v)) 827 828 a2 = 0.5_dp*(1.0_dp/ACOULA + 1.0_dp/ACOULB) 829 830 r1 = 1.0_dp/rr 831 r2 = r1*r1 832 r3 = r2*r1 833 r5 = r3*r2 834 r7 = r5*r2 835 r9 = r7*r2 836 r11 = r9*r2 837 838 w0 = a2*rr 839 w = EXP(-w0) 840 w1 = (1.0_dp + 0.5_dp*w0) 841 w2 = (w1 + 0.5_dp*w0 + 0.5_dp*w0**2) 842 w3 = (w2 + w0**3/6.0_dp) 843 w4 = (w3 + w0**4/30.0_dp) 844 w5 = (w3 + 8.0_dp*w0**4/210.0_dp + w0**5/210.0_dp) 845 w6 = (w3 + 5.0_dp*w0**4/126.0_dp + 2.0_dp*w0**5/315.0_dp + w0**6/1890.0_dp) 846 847 f = (1.0_dp - w*w1)*r1 848 d1 = -1.0_dp*(1.0_dp - w*w2)*r3 849 d2 = 3.0_dp*(1.0_dp - w*w3)*r5 850 d3 = -15.0_dp*(1.0_dp - w*w4)*r7 851 d4 = 105.0_dp*(1.0_dp - w*w5)*r9 852 d5 = -945.0_dp*(1.0_dp - w*w6)*r11 853 854 kr = alpha*rr 855 kr2 = kr*kr 856 w0 = 1.0_dp - erfc(kr) 857 w1 = 2.0_dp*oorootpi*EXP(-kr2) 858 w2 = w1*kr 859 860 f = f - w0*r1 861 d1 = d1 + (-w2 + w0)*r3 862 d2 = d2 + (w2*(3.0_dp + kr2*2.0_dp) - 3.0_dp*w0)*r5 863 d3 = d3 + (-w2*(15.0_dp + kr2*(10.0_dp + kr2*4.0_dp)) + 15.0_dp*w0)*r7 864 d4 = d4 + (w2*(105.0_dp + kr2*(70.0_dp + kr2*(28.0_dp + kr2*8.0_dp))) - 105.0_dp*w0)*r9 865 d5 = d5 + (-w2*(945.0_dp + kr2*(630.0_dp + kr2*(252.0_dp + kr2*(72.0_dp + kr2*16.0_dp)))) + 945.0_dp*w0)*r11 866 867 CALL build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5) 868 869 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 870 DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2 871 872 tmp = M0A(imA)*M0B(imB) 873 wv(1) = tmp*d1f(1) 874 wv(2) = tmp*d1f(2) 875 wv(3) = tmp*d1f(3) 876 877 DO k1 = 1, 3 878 tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB) 879 wv(1) = wv(1) + tmp*d2f(1, k1) 880 wv(2) = wv(2) + tmp*d2f(2, k1) 881 wv(3) = wv(3) + tmp*d2f(3, k1) 882 ENDDO 883 DO k2 = 1, 3 884 DO k1 = 1, 3 885 tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB) 886 wv(1) = wv(1) + tmp*d3f(1, k1, k2) 887 wv(2) = wv(2) + tmp*d3f(2, k1, k2) 888 wv(3) = wv(3) + tmp*d3f(3, k1, k2) 889 ENDDO 890 ENDDO 891 DO k3 = 1, 3 892 DO k2 = 1, 3 893 DO k1 = 1, 3 894 tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB) 895 wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3) 896 wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3) 897 wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3) 898 ENDDO 899 ENDDO 900 ENDDO 901 902 DO k4 = 1, 3 903 DO k3 = 1, 3 904 DO k2 = 1, 3 905 DO k1 = 1, 3 906 tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB) 907 wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4) 908 wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4) 909 wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4) 910 ENDDO 911 ENDDO 912 ENDDO 913 ENDDO 914 915 dCoul(1, imA, imB) = wv(1) 916 dCoul(2, imA, imB) = wv(2) 917 dCoul(3, imA, imB) = wv(3) 918 ENDDO 919 ENDDO 920 921 v(1) = RAB(1) 922 v(2) = RAB(2) 923 v(3) = RAB(3) 924 925 f = 0.0_dp 926 d1f = 0.0_dp 927 d2f = 0.0_dp 928 d3f = 0.0_dp 929 d4f = 0.0_dp 930 d5f = 0.0_dp 931 932 DO gpt = 1, pw_grid%ngpts_cut_local 933 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt)) 934 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt)) 935 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt)) 936 937 lp = lp + bds(1, 1) 938 mp = mp + bds(1, 2) 939 np = np + bds(1, 3) 940 941 IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE 942 kk(:) = pw_grid%g(:, gpt) 943 ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt) 944 945 kr = DOT_PRODUCT(kk, v) 946 cc = COS(kr) 947 ss = SIN(kr) 948 949 f = f + cc*ff 950 DO k1 = 1, 3 951 d1f(k1) = d1f(k1) - kk(k1)*ss*ff 952 ENDDO 953 DO k2 = 1, 3 954 DO k1 = 1, 3 955 d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*cc*ff 956 ENDDO 957 ENDDO 958 DO k3 = 1, 3 959 DO k2 = 1, 3 960 DO k1 = 1, 3 961 d3f(k1, k2, k3) = d3f(k1, k2, k3) + kk(k1)*kk(k2)*kk(k3)*ss*ff 962 ENDDO 963 ENDDO 964 ENDDO 965 DO k4 = 1, 3 966 DO k3 = 1, 3 967 DO k2 = 1, 3 968 DO k1 = 1, 3 969 d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*cc*ff 970 ENDDO 971 ENDDO 972 ENDDO 973 ENDDO 974 DO k5 = 1, 3 975 DO k4 = 1, 3 976 DO k3 = 1, 3 977 DO k2 = 1, 3 978 DO k1 = 1, 3 979 d5f(k1, k2, k3, k4, k5) = d5f(k1, k2, k3, k4, k5) - kk(k1)*kk(k2)*kk(k3)*kk(k4)*kk(k5)*ss*ff 980 ENDDO 981 ENDDO 982 ENDDO 983 ENDDO 984 ENDDO 985 ENDDO 986 987 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 988 DO imB = 1, (sepj%natorb*(sepj%natorb + 1))/2 989 tmp = M0A(imA)*M0B(imB) 990 wv(1) = tmp*d1f(1) 991 wv(2) = tmp*d1f(2) 992 wv(3) = tmp*d1f(3) 993 DO k1 = 1, 3 994 tmp = M1A(k1, imA)*M0B(imB) - M0A(imA)*M1B(k1, imB) 995 wv(1) = wv(1) + tmp*d2f(1, k1) 996 wv(2) = wv(2) + tmp*d2f(2, k1) 997 wv(3) = wv(3) + tmp*d2f(3, k1) 998 ENDDO 999 DO k2 = 1, 3 1000 DO k1 = 1, 3 1001 tmp = M2A(k1, k2, imA)*M0B(imB) - M1A(k1, imA)*M1B(k2, imB) + M0A(imA)*M2B(k1, k2, imB) 1002 wv(1) = wv(1) + tmp*d3f(1, k1, k2) 1003 wv(2) = wv(2) + tmp*d3f(2, k1, k2) 1004 wv(3) = wv(3) + tmp*d3f(3, k1, k2) 1005 ENDDO 1006 ENDDO 1007 DO k3 = 1, 3 1008 DO k2 = 1, 3 1009 DO k1 = 1, 3 1010 tmp = -M2A(k1, k2, imA)*M1B(k3, imB) + M1A(k1, imA)*M2B(k2, k3, imB) 1011 wv(1) = wv(1) + tmp*d4f(1, k1, k2, k3) 1012 wv(2) = wv(2) + tmp*d4f(2, k1, k2, k3) 1013 wv(3) = wv(3) + tmp*d4f(3, k1, k2, k3) 1014 ENDDO 1015 ENDDO 1016 ENDDO 1017 1018 DO k4 = 1, 3 1019 DO k3 = 1, 3 1020 DO k2 = 1, 3 1021 DO k1 = 1, 3 1022 tmp = M2A(k1, k2, imA)*M2B(k3, k4, imB) 1023 wv(1) = wv(1) + tmp*d5f(1, k1, k2, k3, k4) 1024 wv(2) = wv(2) + tmp*d5f(2, k1, k2, k3, k4) 1025 wv(3) = wv(3) + tmp*d5f(3, k1, k2, k3, k4) 1026 ENDDO 1027 ENDDO 1028 ENDDO 1029 ENDDO 1030 1031 dCoul(1, imA, imB) = dCoul(1, imA, imB) + wv(1) 1032 dCoul(2, imA, imB) = dCoul(2, imA, imB) + wv(2) 1033 dCoul(3, imA, imB) = dCoul(3, imA, imB) + wv(3) 1034 ENDDO 1035 ENDDO 1036 1037 END SUBROUTINE makedCoulE 1038 1039! ************************************************************************************************** 1040!> \brief Builds the tensor for the evaluation of the integrals with the 1041!> cartesian multipoles 1042!> \param d1f ... 1043!> \param d2f ... 1044!> \param d3f ... 1045!> \param d4f ... 1046!> \param d5f ... 1047!> \param v ... 1048!> \param d1 ... 1049!> \param d2 ... 1050!> \param d3 ... 1051!> \param d4 ... 1052!> \param d5 ... 1053! ************************************************************************************************** 1054 SUBROUTINE build_d_tensor_gks(d1f, d2f, d3f, d4f, d5f, v, d1, d2, d3, d4, d5) 1055 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: d1f 1056 REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: d2f 1057 REAL(KIND=dp), DIMENSION(3, 3, 3), INTENT(OUT) :: d3f 1058 REAL(KIND=dp), DIMENSION(3, 3, 3, 3), INTENT(OUT) :: d4f 1059 REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3), & 1060 INTENT(OUT), OPTIONAL :: d5f 1061 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: v 1062 REAL(KIND=dp), INTENT(IN) :: d1, d2, d3, d4 1063 REAL(KIND=dp), INTENT(IN), OPTIONAL :: d5 1064 1065 INTEGER :: k1, k2, k3, k4, k5 1066 REAL(KIND=dp) :: w 1067 1068 d1f = 0.0_dp 1069 d2f = 0.0_dp 1070 d3f = 0.0_dp 1071 d4f = 0.0_dp 1072 DO k1 = 1, 3 1073 d1f(k1) = d1f(k1) + v(k1)*d1 1074 ENDDO 1075 DO k1 = 1, 3 1076 DO k2 = 1, 3 1077 d2f(k2, k1) = d2f(k2, k1) + v(k1)*v(k2)*d2 1078 ENDDO 1079 d2f(k1, k1) = d2f(k1, k1) + d1 1080 ENDDO 1081 DO k1 = 1, 3 1082 DO k2 = 1, 3 1083 DO k3 = 1, 3 1084 d3f(k3, k2, k1) = d3f(k3, k2, k1) + v(k1)*v(k2)*v(k3)*d3 1085 ENDDO 1086 w = v(k1)*d2 1087 d3f(k1, k2, k2) = d3f(k1, k2, k2) + w 1088 d3f(k2, k1, k2) = d3f(k2, k1, k2) + w 1089 d3f(k2, k2, k1) = d3f(k2, k2, k1) + w 1090 ENDDO 1091 ENDDO 1092 DO k1 = 1, 3 1093 DO k2 = 1, 3 1094 DO k3 = 1, 3 1095 DO k4 = 1, 3 1096 d4f(k4, k3, k2, k1) = d4f(k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*d4 1097 ENDDO 1098 w = v(k1)*v(k2)*d3 1099 d4f(k1, k2, k3, k3) = d4f(k1, k2, k3, k3) + w 1100 d4f(k1, k3, k2, k3) = d4f(k1, k3, k2, k3) + w 1101 d4f(k3, k1, k2, k3) = d4f(k3, k1, k2, k3) + w 1102 d4f(k1, k3, k3, k2) = d4f(k1, k3, k3, k2) + w 1103 d4f(k3, k1, k3, k2) = d4f(k3, k1, k3, k2) + w 1104 d4f(k3, k3, k1, k2) = d4f(k3, k3, k1, k2) + w 1105 ENDDO 1106 d4f(k1, k1, k2, k2) = d4f(k1, k1, k2, k2) + d2 1107 d4f(k1, k2, k1, k2) = d4f(k1, k2, k1, k2) + d2 1108 d4f(k1, k2, k2, k1) = d4f(k1, k2, k2, k1) + d2 1109 ENDDO 1110 ENDDO 1111 IF (PRESENT(d5f) .AND. PRESENT(d5)) THEN 1112 d5f = 0.0_dp 1113 1114 DO k1 = 1, 3 1115 DO k2 = 1, 3 1116 DO k3 = 1, 3 1117 DO k4 = 1, 3 1118 DO k5 = 1, 3 1119 d5f(k5, k4, k3, k2, k1) = d5f(k5, k4, k3, k2, k1) + v(k1)*v(k2)*v(k3)*v(k4)*v(k5)*d5 1120 ENDDO 1121 w = v(k1)*v(k2)*v(k3)*d4 1122 d5f(k1, k2, k3, k4, k4) = d5f(k1, k2, k3, k4, k4) + w 1123 d5f(k1, k2, k4, k3, k4) = d5f(k1, k2, k4, k3, k4) + w 1124 d5f(k1, k4, k2, k3, k4) = d5f(k1, k4, k2, k3, k4) + w 1125 d5f(k4, k1, k2, k3, k4) = d5f(k4, k1, k2, k3, k4) + w 1126 d5f(k1, k2, k4, k4, k3) = d5f(k1, k2, k4, k4, k3) + w 1127 d5f(k1, k4, k2, k4, k3) = d5f(k1, k4, k2, k4, k3) + w 1128 d5f(k4, k1, k2, k4, k3) = d5f(k4, k1, k2, k4, k3) + w 1129 d5f(k1, k4, k4, k2, k3) = d5f(k1, k4, k4, k2, k3) + w 1130 d5f(k4, k1, k4, k2, k3) = d5f(k4, k1, k4, k2, k3) + w 1131 d5f(k4, k4, k1, k2, k3) = d5f(k4, k4, k1, k2, k3) + w 1132 ENDDO 1133 w = v(k1)*d3 1134 d5f(k1, k2, k2, k3, k3) = d5f(k1, k2, k2, k3, k3) + w 1135 d5f(k1, k2, k3, k2, k3) = d5f(k1, k2, k3, k2, k3) + w 1136 d5f(k1, k2, k3, k3, k2) = d5f(k1, k2, k3, k3, k2) + w 1137 d5f(k2, k1, k2, k3, k3) = d5f(k2, k1, k2, k3, k3) + w 1138 d5f(k2, k1, k3, k2, k3) = d5f(k2, k1, k3, k2, k3) + w 1139 d5f(k2, k1, k3, k3, k2) = d5f(k2, k1, k3, k3, k2) + w 1140 d5f(k2, k2, k1, k3, k3) = d5f(k2, k2, k1, k3, k3) + w 1141 d5f(k2, k3, k1, k2, k3) = d5f(k2, k3, k1, k2, k3) + w 1142 d5f(k2, k3, k1, k3, k2) = d5f(k2, k3, k1, k3, k2) + w 1143 d5f(k2, k2, k3, k1, k3) = d5f(k2, k2, k3, k1, k3) + w 1144 d5f(k2, k3, k2, k1, k3) = d5f(k2, k3, k2, k1, k3) + w 1145 d5f(k2, k3, k3, k1, k2) = d5f(k2, k3, k3, k1, k2) + w 1146 d5f(k2, k2, k3, k3, k1) = d5f(k2, k2, k3, k3, k1) + w 1147 d5f(k2, k3, k2, k3, k1) = d5f(k2, k3, k2, k3, k1) + w 1148 d5f(k2, k3, k3, k2, k1) = d5f(k2, k3, k3, k2, k1) + w 1149 ENDDO 1150 ENDDO 1151 ENDDO 1152 END IF 1153 END SUBROUTINE build_d_tensor_gks 1154 1155! ************************************************************************************************** 1156!> \brief ... 1157!> \param sepi ... 1158!> \param Coul ... 1159!> \param se_int_control ... 1160! ************************************************************************************************** 1161 SUBROUTINE makeCoulE0(sepi, Coul, se_int_control) 1162 TYPE(semi_empirical_type), POINTER :: sepi 1163 REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT) :: Coul 1164 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1165 1166 CHARACTER(len=*), PARAMETER :: routineN = 'makeCoulE0', routineP = moduleN//':'//routineN 1167 1168 INTEGER :: gpt, imA, imB, k1, k2, k3, k4, lp, mp, np 1169 INTEGER, DIMENSION(:, :), POINTER :: bds 1170 REAL(KIND=dp) :: alpha, d2f(3, 3), d4f(3, 3, 3, 3), f, & 1171 ff, w 1172 REAL(KIND=dp), DIMENSION(3) :: kk 1173 REAL(KIND=dp), DIMENSION(3, 3, 45) :: M2A 1174 REAL(KIND=dp), DIMENSION(3, 45) :: M1A 1175 REAL(KIND=dp), DIMENSION(45) :: M0A 1176 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0 1177 TYPE(dg_rho0_type), POINTER :: dg_rho0 1178 TYPE(dg_type), POINTER :: dg 1179 TYPE(pw_grid_type), POINTER :: pw_grid 1180 TYPE(pw_pool_type), POINTER :: pw_pool 1181 1182 alpha = se_int_control%ewald_gks%alpha 1183 pw_pool => se_int_control%ewald_gks%pw_pool 1184 dg => se_int_control%ewald_gks%dg 1185 CALL dg_get(dg, dg_rho0=dg_rho0) 1186 rho0 => dg_rho0%density%pw%cr3d 1187 pw_grid => pw_pool%pw_grid 1188 bds => pw_grid%bounds 1189 1190 CALL get_se_slater_multipole(sepi, M0A, M1A, M2A) 1191 1192 f = 0.0_dp 1193 d2f = 0.0_dp 1194 d4f = 0.0_dp 1195 1196 DO gpt = 1, pw_grid%ngpts_cut_local 1197 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt)) 1198 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt)) 1199 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt)) 1200 1201 lp = lp + bds(1, 1) 1202 mp = mp + bds(1, 2) 1203 np = np + bds(1, 3) 1204 1205 IF (pw_grid%gsq(gpt) == 0.0_dp) CYCLE 1206 kk(:) = pw_grid%g(:, gpt) 1207 ff = 2.0_dp*fourpi*rho0(lp, mp, np)**2*pw_grid%vol/pw_grid%gsq(gpt) 1208 1209 f = f + ff 1210 DO k2 = 1, 3 1211 DO k1 = 1, 3 1212 d2f(k1, k2) = d2f(k1, k2) - kk(k1)*kk(k2)*ff 1213 ENDDO 1214 ENDDO 1215 DO k4 = 1, 3 1216 DO k3 = 1, 3 1217 DO k2 = 1, 3 1218 DO k1 = 1, 3 1219 d4f(k1, k2, k3, k4) = d4f(k1, k2, k3, k4) + kk(k1)*kk(k2)*kk(k3)*kk(k4)*ff 1220 ENDDO 1221 ENDDO 1222 ENDDO 1223 ENDDO 1224 1225 ENDDO 1226 1227 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 1228 DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2 1229 1230 w = M0A(imA)*M0A(imB)*f 1231 DO k2 = 1, 3 1232 DO k1 = 1, 3 1233 w = w + (M2A(k1, k2, imA)*M0A(imB) - M1A(k1, imA)*M1A(k2, imB) + M0A(imA)*M2A(k1, k2, imB))*d2f(k1, k2) 1234 ENDDO 1235 ENDDO 1236 1237 DO k4 = 1, 3 1238 DO k3 = 1, 3 1239 DO k2 = 1, 3 1240 DO k1 = 1, 3 1241 w = w + M2A(k1, k2, imA)*M2A(k3, k4, imB)*d4f(k1, k2, k3, k4) 1242 ENDDO 1243 ENDDO 1244 ENDDO 1245 ENDDO 1246 1247 Coul(imA, imB) = w 1248 1249 ENDDO 1250 ENDDO 1251 1252 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 1253 DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2 1254 w = -M0A(imA)*M0A(imB)*0.25_dp*fourpi/(pw_grid%vol*alpha**2) 1255 Coul(imA, imB) = Coul(imA, imB) + w 1256 ENDDO 1257 ENDDO 1258 1259 DO imA = 1, (sepi%natorb*(sepi%natorb + 1))/2 1260 DO imB = 1, (sepi%natorb*(sepi%natorb + 1))/2 1261 1262 w = M0A(imA)*M0A(imB) 1263 Coul(imA, imB) = Coul(imA, imB) - 2.0_dp*alpha*oorootpi*w 1264 1265 w = 0.0_dp 1266 DO k1 = 1, 3 1267 w = w + M1A(k1, imA)*M1A(k1, imB) 1268 w = w - M0A(imA)*M2A(k1, k1, imB) 1269 w = w - M2A(k1, k1, imA)*M0A(imB) 1270 ENDDO 1271 Coul(imA, imB) = Coul(imA, imB) - 4.0_dp*alpha**3*oorootpi*w/3.0_dp 1272 1273 w = 0.0_dp 1274 DO k2 = 1, 3 1275 DO k1 = 1, 3 1276 w = w + 2.0_dp*M2A(k1, k2, imA)*M2A(k1, k2, imB) 1277 w = w + M2A(k1, k1, imA)*M2A(k2, k2, imB) 1278 ENDDO 1279 ENDDO 1280 Coul(imA, imB) = Coul(imA, imB) - 8.0_dp*alpha**5*oorootpi*w/5.0_dp 1281 ENDDO 1282 ENDDO 1283 END SUBROUTINE makeCoulE0 1284 1285! ************************************************************************************************** 1286!> \brief Retrieves the multipole for the Slater integral evaluation 1287!> \param sepi ... 1288!> \param M0 ... 1289!> \param M1 ... 1290!> \param M2 ... 1291!> \param ACOUL ... 1292! ************************************************************************************************** 1293 SUBROUTINE get_se_slater_multipole(sepi, M0, M1, M2, ACOUL) 1294 TYPE(semi_empirical_type), POINTER :: sepi 1295 REAL(kind=dp), DIMENSION(45), INTENT(OUT) :: M0 1296 REAL(kind=dp), DIMENSION(3, 45), INTENT(OUT) :: M1 1297 REAL(kind=dp), DIMENSION(3, 3, 45), INTENT(OUT) :: M2 1298 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: ACOUL 1299 1300 INTEGER :: i, j, jint, size_1c_int 1301 TYPE(semi_empirical_mpole_type), POINTER :: mpole 1302 1303 NULLIFY (mpole) 1304 size_1c_int = SIZE(sepi%w_mpole) 1305 DO jint = 1, size_1c_int 1306 mpole => sepi%w_mpole(jint)%mpole 1307 i = mpole%indi 1308 j = mpole%indj 1309 M0(indexb(i, j)) = -mpole%cs 1310 1311 M1(1, indexb(i, j)) = -mpole%ds(1) 1312 M1(2, indexb(i, j)) = -mpole%ds(2) 1313 M1(3, indexb(i, j)) = -mpole%ds(3) 1314 1315 M2(1, 1, indexb(i, j)) = -mpole%qq(1, 1)/3.0_dp 1316 M2(2, 1, indexb(i, j)) = -mpole%qq(2, 1)/3.0_dp 1317 M2(3, 1, indexb(i, j)) = -mpole%qq(3, 1)/3.0_dp 1318 1319 M2(1, 2, indexb(i, j)) = -mpole%qq(1, 2)/3.0_dp 1320 M2(2, 2, indexb(i, j)) = -mpole%qq(2, 2)/3.0_dp 1321 M2(3, 2, indexb(i, j)) = -mpole%qq(3, 2)/3.0_dp 1322 1323 M2(1, 3, indexb(i, j)) = -mpole%qq(1, 3)/3.0_dp 1324 M2(2, 3, indexb(i, j)) = -mpole%qq(2, 3)/3.0_dp 1325 M2(3, 3, indexb(i, j)) = -mpole%qq(3, 3)/3.0_dp 1326 ENDDO 1327 IF (PRESENT(ACOUL)) ACOUL = sepi%acoul 1328 END SUBROUTINE get_se_slater_multipole 1329 1330END MODULE semi_empirical_int_gks 1331