1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5!!****** cp2k/ai_overlap3 [1.0] * 6!! 7!! NAME 8!! ai_overlap3 9!! 10!! FUNCTION 11!! Calculation of three-center overlap integrals over Cartesian 12!! Gaussian-type functions. 13!! 14!! AUTHOR 15!! Matthias Krack (26.06.2001) 16!! 17!! LITERATURE 18!! S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986) 19!! 20!****************************************************************************** 21 22MODULE ai_overlap3 23 24! ************************************************************************************************** 25 26! ax,ay,az : Angular momentum index numbers of orbital a. 27! bx,by,bz : Angular momentum index numbers of orbital b. 28! coset : Cartesian orbital set pointer. 29! dab : Distance between the atomic centers a and b. 30! dac : Distance between the atomic centers a and c. 31! dbc : Distance between the atomic centers b and c. 32! l{a,b,c} : Angular momentum quantum number of shell a, b or c. 33! l{a,b}_max : Maximum angular momentum quantum number of shell a, b or c. 34! ncoset : Number of Cartesian orbitals up to l. 35! rab : Distance vector between the atomic centers a and b. 36! rac : Distance vector between the atomic centers a and c. 37! rbc : Distance vector between the atomic centers b and c. 38! rpgf{a,b,c}: Radius of the primitive Gaussian-type function a or b. 39! zet{a,b,c} : Exponents of the Gaussian-type functions a or b. 40! zetg : Reciprocal of the sum of the exponents of orbital a, b and c. 41! zetp : Reciprocal of the sum of the exponents of orbital a and b. 42 43! ************************************************************************************************** 44 45 USE kinds, ONLY: dp 46 USE mathconstants, ONLY: pi 47 USE orbital_pointers, ONLY: coset,& 48 ncoset 49#include "../base/base_uses.f90" 50 51 IMPLICIT NONE 52 53 PRIVATE 54 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_overlap3' 56 57! *** Public subroutines *** 58 59 PUBLIC :: overlap3 60 61!!*** 62! ************************************************************************************************** 63 64CONTAINS 65 66! *************************************************************************************************** 67!> \brief Calculation of three-center overlap integrals [a|b|c] over primitive 68!> Cartesian Gaussian functions 69!> \param la_max_set ... 70!> \param npgfa ... 71!> \param zeta ... 72!> \param rpgfa ... 73!> \param la_min_set ... 74!> \param lb_max_set ... 75!> \param npgfb ... 76!> \param zetb ... 77!> \param rpgfb ... 78!> \param lb_min_set ... 79!> \param lc_max_set ... 80!> \param npgfc ... 81!> \param zetc ... 82!> \param rpgfc ... 83!> \param lc_min_set ... 84!> \param rab ... 85!> \param dab ... 86!> \param rac ... 87!> \param dac ... 88!> \param rbc ... 89!> \param dbc ... 90!> \param sabc integrals [a|b|c] 91!> \param sdabc derivative [da/dAi|b|c] 92!> \param sabdc derivative [a|b|dc/dCi] 93!> \param int_abc_ext the extremal value of sabc, i.e., MAXVAL(ABS(sabc)) 94!> \par History 95!> 05.2014 created (Dorothea Golze) 96!> \author Dorothea Golze 97!> \note overlap3 essentially uses the setup of overlap3_old 98! ************************************************************************************************** 99 100 SUBROUTINE overlap3(la_max_set, npgfa, zeta, rpgfa, la_min_set, & 101 lb_max_set, npgfb, zetb, rpgfb, lb_min_set, & 102 lc_max_set, npgfc, zetc, rpgfc, lc_min_set, & 103 rab, dab, rac, dac, rbc, dbc, sabc, & 104 sdabc, sabdc, int_abc_ext) 105 106 INTEGER, INTENT(IN) :: la_max_set, npgfa 107 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta, rpgfa 108 INTEGER, INTENT(IN) :: la_min_set, lb_max_set, npgfb 109 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb, rpgfb 110 INTEGER, INTENT(IN) :: lb_min_set, lc_max_set, npgfc 111 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc, rpgfc 112 INTEGER, INTENT(IN) :: lc_min_set 113 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rab 114 REAL(KIND=dp), INTENT(IN) :: dab 115 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rac 116 REAL(KIND=dp), INTENT(IN) :: dac 117 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rbc 118 REAL(KIND=dp), INTENT(IN) :: dbc 119 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: sabc 120 REAL(KIND=dp), DIMENSION(:, :, :, :), & 121 INTENT(INOUT), OPTIONAL :: sdabc, sabdc 122 REAL(dp), INTENT(OUT), OPTIONAL :: int_abc_ext 123 124 CHARACTER(len=*), PARAMETER :: routineN = 'overlap3', routineP = moduleN//':'//routineN 125 126 INTEGER :: ax, ay, az, bx, by, bz, coa, coax, coay, coaz, coc, cocx, cocy, cocz, cx, cy, cz, & 127 handle, i, ipgf, j, jpgf, k, kpgf, l, la, la_max, la_min, la_start, lai, lb, lb_max, & 128 lb_min, lbi, lc, lc_max, lc_min, lci, na, nb, nc, nda, ndc 129 REAL(KIND=dp) :: f0, f1, f2, f3, fcx, fcy, fcz, fx, fy, & 130 fz, rcp2, zetg, zetp 131 REAL(KIND=dp), DIMENSION(3) :: rag, rbg, rcg, rcp 132 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: s 133 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc 134 135! --------------------------------------------------------------------------- 136 137 CALL timeset(routineN, handle) 138 139 NULLIFY (s, sda, sdc) 140 141 lai = 0 142 lbi = 0 143 lci = 0 144 145 IF (PRESENT(sdabc)) lai = 1 146 IF (PRESENT(sabdc)) lci = 1 147 148 la_max = la_max_set + lai 149 la_min = MAX(0, la_min_set - lai) 150 lb_max = lb_max_set 151 lb_min = lb_min_set 152 lc_max = lc_max_set + lci 153 lc_min = MAX(0, lc_min_set - lci) 154 155 ALLOCATE (s(ncoset(la_max), ncoset(lb_max), ncoset(lc_max))) 156 s = 0._dp 157 IF (PRESENT(sdabc)) THEN 158 ALLOCATE (sda(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3)) 159 sda = 0._dp 160 ENDIF 161 IF (PRESENT(sabdc)) THEN 162 ALLOCATE (sdc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max), 3)) 163 sdc = 0._dp 164 ENDIF 165 IF (PRESENT(int_abc_ext)) THEN 166 int_abc_ext = 0.0_dp 167 ENDIF 168 169! *** Loop over all pairs of primitive Gaussian-type functions *** 170 171 na = 0 172 nda = 0 173 DO ipgf = 1, npgfa 174 175 nb = 0 176 DO jpgf = 1, npgfb 177 178 ! *** Screening *** 179 IF (rpgfa(ipgf) + rpgfb(jpgf) < dab) THEN 180 nb = nb + ncoset(lb_max_set) 181 CYCLE 182 END IF 183 184 nc = 0 185 ndc = 0 186 DO kpgf = 1, npgfc 187 188 ! *** Screening *** 189 IF ((rpgfb(jpgf) + rpgfc(kpgf) < dbc) .OR. & 190 (rpgfa(ipgf) + rpgfc(kpgf) < dac)) THEN 191 nc = nc + ncoset(lc_max_set) 192 ndc = ndc + ncoset(lc_max_set) 193 CYCLE 194 END IF 195 196 ! *** Calculate some prefactors *** 197 zetg = 1.0_dp/(zeta(ipgf) + zetb(jpgf) + zetc(kpgf)) 198 zetp = 1.0_dp/(zeta(ipgf) + zetb(jpgf)) 199 f0 = (pi*zetg)**1.5_dp 200 f1 = zetb(jpgf)*zetp 201 f2 = 0.5_dp*zetg 202 rcp(:) = f1*rab(:) - rac(:) 203 rcp2 = rcp(1)*rcp(1) + rcp(2)*rcp(2) + rcp(3)*rcp(3) 204 205 ! *** Calculate the basic three-center overlap integral [s|s|s] *** 206 s(1, 1, 1) = f0*EXP(-(zeta(ipgf)*f1*dab*dab + zetc(kpgf)*zetg*rcp2/zetp)) 207 208! *** Recurrence steps: [s|s|s] -> [a|s|s] *** 209 210 IF (la_max > 0) THEN 211 212! *** Vertical recurrence steps: [s|s|s] -> [a|s|s] *** 213 214 rag(:) = zetg*(zetb(jpgf)*rab(:) + zetc(kpgf)*rac(:)) 215 216! *** [p|s|s] = (Gi - Ai)*[s|s|s] (i = x,y,z) *** 217 218 s(2, 1, 1) = rag(1)*s(1, 1, 1) 219 s(3, 1, 1) = rag(2)*s(1, 1, 1) 220 s(4, 1, 1) = rag(3)*s(1, 1, 1) 221 222! *** [a|s|s] = (Gi - Ai)*[a-1i|s|s] + f2*Ni(a-1i)*[a-2i|s|s] *** 223 224 DO la = 2, la_max 225 226! *** Increase the angular momentum component z of function a *** 227 228 s(coset(0, 0, la), 1, 1) = rag(3)*s(coset(0, 0, la - 1), 1, 1) + & 229 f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1, 1) 230 231! *** Increase the angular momentum component y of function a *** 232 233 az = la - 1 234 s(coset(0, 1, az), 1, 1) = rag(2)*s(coset(0, 0, az), 1, 1) 235 236 DO ay = 2, la 237 az = la - ay 238 s(coset(0, ay, az), 1, 1) = rag(2)*s(coset(0, ay - 1, az), 1, 1) + & 239 f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1, 1) 240 END DO 241 242! *** Increase the angular momentum component x of function a *** 243 244 DO ay = 0, la - 1 245 az = la - 1 - ay 246 s(coset(1, ay, az), 1, 1) = rag(1)*s(coset(0, ay, az), 1, 1) 247 END DO 248 249 DO ax = 2, la 250 f3 = f2*REAL(ax - 1, dp) 251 DO ay = 0, la - ax 252 az = la - ax - ay 253 s(coset(ax, ay, az), 1, 1) = rag(1)*s(coset(ax - 1, ay, az), 1, 1) + & 254 f3*s(coset(ax - 2, ay, az), 1, 1) 255 END DO 256 END DO 257 258 END DO 259 260! *** Recurrence steps: [a|s|s] -> [a|s|b] *** 261 262 IF (lb_max > 0) THEN 263 264! *** Horizontal recurrence steps *** 265 266 rbg(:) = rag(:) - rab(:) 267 268! *** [a|s|p] = [a+1i|s|s] - (Bi - Ai)*[a|s|s] *** 269 270 IF (lb_max == 1) THEN 271 la_start = la_min 272 ELSE 273 la_start = MAX(0, la_min - 1) 274 END IF 275 276 DO la = la_start, la_max - 1 277 DO ax = 0, la 278 DO ay = 0, la - ax 279 az = la - ax - ay 280 coa = coset(ax, ay, az) 281 coax = coset(ax + 1, ay, az) 282 coay = coset(ax, ay + 1, az) 283 coaz = coset(ax, ay, az + 1) 284 s(coset(ax, ay, az), 2, 1) = s(coax, 1, 1) - rab(1)*s(coa, 1, 1) 285 s(coset(ax, ay, az), 3, 1) = s(coay, 1, 1) - rab(2)*s(coa, 1, 1) 286 s(coset(ax, ay, az), 4, 1) = s(coaz, 1, 1) - rab(3)*s(coa, 1, 1) 287 END DO 288 END DO 289 END DO 290 291! *** Vertical recurrence step *** 292 293! *** [a|s|p] = (Gi - Bi)*[a|s|s] + f2*Ni(a)*[a-1i|s|s] *** 294 295 DO ax = 0, la_max 296 fx = f2*REAL(ax, dp) 297 DO ay = 0, la_max - ax 298 fy = f2*REAL(ay, dp) 299 az = la_max - ax - ay 300 fz = f2*REAL(az, dp) 301 coa = coset(ax, ay, az) 302 IF (ax == 0) THEN 303 s(coa, 2, 1) = rbg(1)*s(coa, 1, 1) 304 ELSE 305 s(coa, 2, 1) = rbg(1)*s(coa, 1, 1) + fx*s(coset(ax - 1, ay, az), 1, 1) 306 END IF 307 IF (ay == 0) THEN 308 s(coa, 3, 1) = rbg(2)*s(coa, 1, 1) 309 ELSE 310 s(coa, 3, 1) = rbg(2)*s(coa, 1, 1) + fy*s(coset(ax, ay - 1, az), 1, 1) 311 END IF 312 IF (az == 0) THEN 313 s(coa, 4, 1) = rbg(3)*s(coa, 1, 1) 314 ELSE 315 s(coa, 4, 1) = rbg(3)*s(coa, 1, 1) + fz*s(coset(ax, ay, az - 1), 1, 1) 316 END IF 317 END DO 318 END DO 319 320! *** Recurrence steps: [a|s|p] -> [a|s|b] *** 321 322 DO lb = 2, lb_max 323 324! *** Horizontal recurrence steps *** 325 326! *** [a|s|b] = [a+1i|s|b-1i] - (Bi - Ai)*[a|s|b-1i] *** 327 328 IF (lb == lb_max) THEN 329 la_start = la_min 330 ELSE 331 la_start = MAX(0, la_min - 1) 332 END IF 333 334 DO la = la_start, la_max - 1 335 DO ax = 0, la 336 DO ay = 0, la - ax 337 az = la - ax - ay 338 339 coa = coset(ax, ay, az) 340 coax = coset(ax + 1, ay, az) 341 coay = coset(ax, ay + 1, az) 342 coaz = coset(ax, ay, az + 1) 343 344! *** Shift of angular momentum component z from a to b *** 345 346 s(coa, coset(0, 0, lb), 1) = & 347 s(coaz, coset(0, 0, lb - 1), 1) - & 348 rab(3)*s(coa, coset(0, 0, lb - 1), 1) 349 350! *** Shift of angular momentum component y from a to b *** 351 352 DO by = 1, lb 353 bz = lb - by 354 s(coa, coset(0, by, bz), 1) = & 355 s(coay, coset(0, by - 1, bz), 1) - & 356 rab(2)*s(coa, coset(0, by - 1, bz), 1) 357 END DO 358 359! *** Shift of angular momentum component x from a to b *** 360 361 DO bx = 1, lb 362 DO by = 0, lb - bx 363 bz = lb - bx - by 364 s(coa, coset(bx, by, bz), 1) = & 365 s(coax, coset(bx - 1, by, bz), 1) - & 366 rab(1)*s(coa, coset(bx - 1, by, bz), 1) 367 END DO 368 END DO 369 370 END DO 371 END DO 372 END DO 373 374! *** Vertical recurrence step *** 375 376! *** [a|s|b] = (Gi - Bi)*[a|s|b-1i] + *** 377! *** f2*Ni(a)*[a-1i|s|b-1i] + *** 378! *** f2*Ni(b-1i)*[a|s|b-2i] *** 379 380 DO ax = 0, la_max 381 fx = f2*REAL(ax, dp) 382 DO ay = 0, la_max - ax 383 fy = f2*REAL(ay, dp) 384 az = la_max - ax - ay 385 fz = f2*REAL(az, dp) 386 387 coa = coset(ax, ay, az) 388 389 f3 = f2*REAL(lb - 1, dp) 390 391! *** Shift of angular momentum component z from a to b *** 392 393 IF (az == 0) THEN 394 s(coa, coset(0, 0, lb), 1) = & 395 rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + & 396 f3*s(coa, coset(0, 0, lb - 2), 1) 397 ELSE 398 coaz = coset(ax, ay, az - 1) 399 s(coa, coset(0, 0, lb), 1) = & 400 rbg(3)*s(coa, coset(0, 0, lb - 1), 1) + & 401 fz*s(coaz, coset(0, 0, lb - 1), 1) + & 402 f3*s(coa, coset(0, 0, lb - 2), 1) 403 END IF 404 405! *** Shift of angular momentum component y from a to b *** 406 407 IF (ay == 0) THEN 408 bz = lb - 1 409 s(coa, coset(0, 1, bz), 1) = & 410 rbg(2)*s(coa, coset(0, 0, bz), 1) 411 DO by = 2, lb 412 bz = lb - by 413 f3 = f2*REAL(by - 1, dp) 414 s(coa, coset(0, by, bz), 1) = & 415 rbg(2)*s(coa, coset(0, by - 1, bz), 1) + & 416 f3*s(coa, coset(0, by - 2, bz), 1) 417 END DO 418 ELSE 419 coay = coset(ax, ay - 1, az) 420 bz = lb - 1 421 s(coa, coset(0, 1, bz), 1) = & 422 rbg(2)*s(coa, coset(0, 0, bz), 1) + & 423 fy*s(coay, coset(0, 0, bz), 1) 424 DO by = 2, lb 425 bz = lb - by 426 f3 = f2*REAL(by - 1, dp) 427 s(coa, coset(0, by, bz), 1) = & 428 rbg(2)*s(coa, coset(0, by - 1, bz), 1) + & 429 fy*s(coay, coset(0, by - 1, bz), 1) + & 430 f3*s(coa, coset(0, by - 2, bz), 1) 431 END DO 432 END IF 433 434! *** Shift of angular momentum component x from a to b *** 435 436 IF (ax == 0) THEN 437 DO by = 0, lb - 1 438 bz = lb - 1 - by 439 s(coa, coset(1, by, bz), 1) = & 440 rbg(1)*s(coa, coset(0, by, bz), 1) 441 END DO 442 DO bx = 2, lb 443 f3 = f2*REAL(bx - 1, dp) 444 DO by = 0, lb - bx 445 bz = lb - bx - by 446 s(coa, coset(bx, by, bz), 1) = & 447 rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + & 448 f3*s(coa, coset(bx - 2, by, bz), 1) 449 END DO 450 END DO 451 ELSE 452 coax = coset(ax - 1, ay, az) 453 DO by = 0, lb - 1 454 bz = lb - 1 - by 455 s(coa, coset(1, by, bz), 1) = & 456 rbg(1)*s(coa, coset(0, by, bz), 1) + & 457 fx*s(coax, coset(0, by, bz), 1) 458 END DO 459 DO bx = 2, lb 460 f3 = f2*REAL(bx - 1, dp) 461 DO by = 0, lb - bx 462 bz = lb - bx - by 463 s(coa, coset(bx, by, bz), 1) = & 464 rbg(1)*s(coa, coset(bx - 1, by, bz), 1) + & 465 fx*s(coax, coset(bx - 1, by, bz), 1) + & 466 f3*s(coa, coset(bx - 2, by, bz), 1) 467 END DO 468 END DO 469 END IF 470 471 END DO 472 END DO 473 474 END DO 475 476 END IF 477 478 ELSE 479 480 IF (lb_max > 0) THEN 481 482! *** Vertical recurrence steps: [s|s|s] -> [s|s|b] *** 483 484 rbg(:) = -zetg*(zeta(ipgf)*rab(:) - zetc(kpgf)*rbc(:)) 485 486! *** [s|s|p] = (Gi - Bi)*[s|s|s] *** 487 488 s(1, 2, 1) = rbg(1)*s(1, 1, 1) 489 s(1, 3, 1) = rbg(2)*s(1, 1, 1) 490 s(1, 4, 1) = rbg(3)*s(1, 1, 1) 491 492! *** [s|s|b] = (Gi - Bi)*[s|s|b-1i] + f2*Ni(b-1i)*[s|s|b-2i] *** 493 494 DO lb = 2, lb_max 495 496! *** Increase the angular momentum component z of function b *** 497 498 s(1, coset(0, 0, lb), 1) = rbg(3)*s(1, coset(0, 0, lb - 1), 1) + & 499 f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2), 1) 500 501! *** Increase the angular momentum component y of function b *** 502 503 bz = lb - 1 504 s(1, coset(0, 1, bz), 1) = rbg(2)*s(1, coset(0, 0, bz), 1) 505 506 DO by = 2, lb 507 bz = lb - by 508 s(1, coset(0, by, bz), 1) = & 509 rbg(2)*s(1, coset(0, by - 1, bz), 1) + & 510 f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz), 1) 511 END DO 512 513! *** Increase the angular momentum component x of function b *** 514 515 DO by = 0, lb - 1 516 bz = lb - 1 - by 517 s(1, coset(1, by, bz), 1) = rbg(1)*s(1, coset(0, by, bz), 1) 518 END DO 519 520 DO bx = 2, lb 521 f3 = f2*REAL(bx - 1, dp) 522 DO by = 0, lb - bx 523 bz = lb - bx - by 524 s(1, coset(bx, by, bz), 1) = rbg(1)*s(1, coset(bx - 1, by, bz), 1) + & 525 f3*s(1, coset(bx - 2, by, bz), 1) 526 END DO 527 END DO 528 529 END DO 530 531 END IF 532 533 END IF 534 535! *** Recurrence steps: [a|s|b] -> [a|c|b] *** 536 537 IF (lc_max > 0) THEN 538 539! *** Vertical recurrence steps: [s|s|s] -> [s|c|s] *** 540 541 rcg(:) = -zetg*(zeta(ipgf)*rac(:) + zetb(jpgf)*rbc(:)) 542 543! *** [s|p|s] = (Gi - Ci)*[s|s|s] (i = x,y,z) *** 544 545 s(1, 1, 2) = rcg(1)*s(1, 1, 1) 546 s(1, 1, 3) = rcg(2)*s(1, 1, 1) 547 s(1, 1, 4) = rcg(3)*s(1, 1, 1) 548 549! *** [s|c|s] = (Gi - Ci)*[s|c-1i|s] + f2*Ni(c-1i)*[s|c-2i|s] *** 550 551 DO lc = 2, lc_max 552 553! *** Increase the angular momentum component z of function c *** 554 555 s(1, 1, coset(0, 0, lc)) = rcg(3)*s(1, 1, coset(0, 0, lc - 1)) + & 556 f2*REAL(lc - 1, dp)*s(1, 1, coset(0, 0, lc - 2)) 557 558! *** Increase the angular momentum component y of function c *** 559 560 cz = lc - 1 561 s(1, 1, coset(0, 1, cz)) = rcg(2)*s(1, 1, coset(0, 0, cz)) 562 563 DO cy = 2, lc 564 cz = lc - cy 565 s(1, 1, coset(0, cy, cz)) = rcg(2)*s(1, 1, coset(0, cy - 1, cz)) + & 566 f2*REAL(cy - 1, dp)*s(1, 1, coset(0, cy - 2, cz)) 567 END DO 568 569! *** Increase the angular momentum component x of function c *** 570 571 DO cy = 0, lc - 1 572 cz = lc - 1 - cy 573 s(1, 1, coset(1, cy, cz)) = rcg(1)*s(1, 1, coset(0, cy, cz)) 574 END DO 575 576 DO cx = 2, lc 577 f3 = f2*REAL(cx - 1, dp) 578 DO cy = 0, lc - cx 579 cz = lc - cx - cy 580 s(1, 1, coset(cx, cy, cz)) = rcg(1)*s(1, 1, coset(cx - 1, cy, cz)) + & 581 f3*s(1, 1, coset(cx - 2, cy, cz)) 582 END DO 583 END DO 584 585 END DO 586 587! *** Recurrence steps: [s|c|s] -> [a|c|b] *** 588 589 DO lc = 1, lc_max 590 591 DO cx = 0, lc 592 DO cy = 0, lc - cx 593 cz = lc - cx - cy 594 595 coc = coset(cx, cy, cz) 596 cocx = coset(MAX(0, cx - 1), cy, cz) 597 cocy = coset(cx, MAX(0, cy - 1), cz) 598 cocz = coset(cx, cy, MAX(0, cz - 1)) 599 600 fcx = f2*REAL(cx, dp) 601 fcy = f2*REAL(cy, dp) 602 fcz = f2*REAL(cz, dp) 603 604! *** Recurrence steps: [s|c|s] -> [a|c|s] *** 605 606 IF (la_max > 0) THEN 607 608! *** Vertical recurrence steps: [s|c|s] -> [a|c|s] *** 609 610 rag(:) = rcg(:) + rac(:) 611 612! *** [p|c|s] = (Gi - Ai)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] *** 613 614 s(2, 1, coc) = rag(1)*s(1, 1, coc) + fcx*s(1, 1, cocx) 615 s(3, 1, coc) = rag(2)*s(1, 1, coc) + fcy*s(1, 1, cocy) 616 s(4, 1, coc) = rag(3)*s(1, 1, coc) + fcz*s(1, 1, cocz) 617 618! *** [a|c|s] = (Gi - Ai)*[a-1i|c|s] + *** 619! *** f2*Ni(a-1i)*[a-2i|c|s] + *** 620! *** f2*Ni(c)*[a-1i|c-1i|s] *** 621 622 DO la = 2, la_max 623 624! *** Increase the angular momentum component z of a *** 625 626 s(coset(0, 0, la), 1, coc) = & 627 rag(3)*s(coset(0, 0, la - 1), 1, coc) + & 628 f2*REAL(la - 1, dp)*s(coset(0, 0, la - 2), 1, coc) + & 629 fcz*s(coset(0, 0, la - 1), 1, cocz) 630 631! *** Increase the angular momentum component y of a *** 632 633 az = la - 1 634 s(coset(0, 1, az), 1, coc) = & 635 rag(2)*s(coset(0, 0, az), 1, coc) + & 636 fcy*s(coset(0, 0, az), 1, cocy) 637 638 DO ay = 2, la 639 az = la - ay 640 s(coset(0, ay, az), 1, coc) = & 641 rag(2)*s(coset(0, ay - 1, az), 1, coc) + & 642 f2*REAL(ay - 1, dp)*s(coset(0, ay - 2, az), 1, coc) + & 643 fcy*s(coset(0, ay - 1, az), 1, cocy) 644 END DO 645 646! *** Increase the angular momentum component x of a *** 647 648 DO ay = 0, la - 1 649 az = la - 1 - ay 650 s(coset(1, ay, az), 1, coc) = & 651 rag(1)*s(coset(0, ay, az), 1, coc) + & 652 fcx*s(coset(0, ay, az), 1, cocx) 653 END DO 654 655 DO ax = 2, la 656 f3 = f2*REAL(ax - 1, dp) 657 DO ay = 0, la - ax 658 az = la - ax - ay 659 s(coset(ax, ay, az), 1, coc) = & 660 rag(1)*s(coset(ax - 1, ay, az), 1, coc) + & 661 f3*s(coset(ax - 2, ay, az), 1, coc) + & 662 fcx*s(coset(ax - 1, ay, az), 1, cocx) 663 END DO 664 END DO 665 666 END DO 667 668! *** Recurrence steps: [a|c|s] -> [a|c|b] *** 669 670 IF (lb_max > 0) THEN 671 672! *** Horizontal recurrence steps *** 673 674 rbg(:) = rag(:) - rab(:) 675 676! *** [a|c|p] = [a+1i|c|s] - (Bi - Ai)*[a|c|s] *** 677 678 IF (lb_max == 1) THEN 679 la_start = la_min 680 ELSE 681 la_start = MAX(0, la_min - 1) 682 END IF 683 684 DO la = la_start, la_max - 1 685 DO ax = 0, la 686 DO ay = 0, la - ax 687 az = la - ax - ay 688 coa = coset(ax, ay, az) 689 coax = coset(ax + 1, ay, az) 690 coay = coset(ax, ay + 1, az) 691 coaz = coset(ax, ay, az + 1) 692 s(coa, 2, coc) = s(coax, 1, coc) - rab(1)*s(coa, 1, coc) 693 s(coa, 3, coc) = s(coay, 1, coc) - rab(2)*s(coa, 1, coc) 694 s(coa, 4, coc) = s(coaz, 1, coc) - rab(3)*s(coa, 1, coc) 695 END DO 696 END DO 697 END DO 698 699! *** Vertical recurrence step *** 700 701! *** [a|c|p] = (Gi - Bi)*[a|c|s] + *** 702! f2*Ni(a)*[a-1i|c|s] + *** 703! f2*Ni(c)*[a|c-1i|s] *** 704 705 DO ax = 0, la_max 706 fx = f2*REAL(ax, dp) 707 DO ay = 0, la_max - ax 708 fy = f2*REAL(ay, dp) 709 az = la_max - ax - ay 710 fz = f2*REAL(az, dp) 711 coa = coset(ax, ay, az) 712 IF (ax == 0) THEN 713 s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + & 714 fcx*s(coa, 1, cocx) 715 ELSE 716 s(coa, 2, coc) = rbg(1)*s(coa, 1, coc) + & 717 fx*s(coset(ax - 1, ay, az), 1, coc) + & 718 fcx*s(coa, 1, cocx) 719 END IF 720 IF (ay == 0) THEN 721 s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + & 722 fcy*s(coa, 1, cocy) 723 ELSE 724 s(coa, 3, coc) = rbg(2)*s(coa, 1, coc) + & 725 fy*s(coset(ax, ay - 1, az), 1, coc) + & 726 fcy*s(coa, 1, cocy) 727 END IF 728 IF (az == 0) THEN 729 s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + & 730 fcz*s(coa, 1, cocz) 731 ELSE 732 s(coa, 4, coc) = rbg(3)*s(coa, 1, coc) + & 733 fz*s(coset(ax, ay, az - 1), 1, coc) + & 734 fcz*s(coa, 1, cocz) 735 END IF 736 END DO 737 END DO 738 739! *** Recurrence steps: [a|c|p] -> [a|c|b] *** 740 741 DO lb = 2, lb_max 742 743! *** Horizontal recurrence steps *** 744 745! *** [a|c|b] = [a+1i|c|b-1i] - (Bi - Ai)*[a|c|b-1i] *** 746 747 IF (lb == lb_max) THEN 748 la_start = la_min 749 ELSE 750 la_start = MAX(0, la_min - 1) 751 END IF 752 753 DO la = la_start, la_max - 1 754 DO ax = 0, la 755 DO ay = 0, la - ax 756 az = la - ax - ay 757 758 coa = coset(ax, ay, az) 759 coax = coset(ax + 1, ay, az) 760 coay = coset(ax, ay + 1, az) 761 coaz = coset(ax, ay, az + 1) 762 763! *** Shift of angular momentum *** 764! *** component z from a to b *** 765 766 s(coa, coset(0, 0, lb), coc) = & 767 s(coaz, coset(0, 0, lb - 1), coc) - & 768 rab(3)*s(coa, coset(0, 0, lb - 1), coc) 769 770! *** Shift of angular momentum *** 771! *** component y from a to b *** 772 773 DO by = 1, lb 774 bz = lb - by 775 s(coa, coset(0, by, bz), coc) = & 776 s(coay, coset(0, by - 1, bz), coc) - & 777 rab(2)*s(coa, coset(0, by - 1, bz), coc) 778 END DO 779 780! *** Shift of angular momentum *** 781! *** component x from a to b *** 782 783 DO bx = 1, lb 784 DO by = 0, lb - bx 785 bz = lb - bx - by 786 s(coa, coset(bx, by, bz), coc) = & 787 s(coax, coset(bx - 1, by, bz), coc) - & 788 rab(1)*s(coa, coset(bx - 1, by, bz), coc) 789 END DO 790 END DO 791 792 END DO 793 END DO 794 END DO 795 796! *** Vertical recurrence step *** 797 798! *** [a|c|b] = (Gi - Bi)*[a|c|b-1i] + *** 799! *** f2*Ni(a)*[a-1i|c|b-1i] + *** 800! *** f2*Ni(b-1i)*[a|c|b-2i] + *** 801! *** f2*Ni(c)*[a|c-1i|b-1i] *** 802 803 DO ax = 0, la_max 804 fx = f2*REAL(ax, dp) 805 DO ay = 0, la_max - ax 806 fy = f2*REAL(ay, dp) 807 az = la_max - ax - ay 808 fz = f2*REAL(az, dp) 809 810 coa = coset(ax, ay, az) 811 coax = coset(MAX(0, ax - 1), ay, az) 812 coay = coset(ax, MAX(0, ay - 1), az) 813 coaz = coset(ax, ay, MAX(0, az - 1)) 814 815 f3 = f2*REAL(lb - 1, dp) 816 817! *** Shift of angular momentum *** 818! *** component z from a to b *** 819 820 IF (az == 0) THEN 821 s(coa, coset(0, 0, lb), coc) = & 822 rbg(3)*s(coa, coset(0, 0, lb - 1), coc) + & 823 f3*s(coa, coset(0, 0, lb - 2), coc) + & 824 fcz*s(coa, coset(0, 0, lb - 1), cocz) 825 ELSE 826 s(coa, coset(0, 0, lb), coc) = & 827 rbg(3)*s(coa, coset(0, 0, lb - 1), coc) + & 828 fz*s(coaz, coset(0, 0, lb - 1), coc) + & 829 f3*s(coa, coset(0, 0, lb - 2), coc) + & 830 fcz*s(coa, coset(0, 0, lb - 1), cocz) 831 END IF 832 833! *** Shift of angular momentum *** 834! *** component y from a to b *** 835 836 IF (ay == 0) THEN 837 bz = lb - 1 838 s(coa, coset(0, 1, bz), coc) = & 839 rbg(2)*s(coa, coset(0, 0, bz), coc) + & 840 fcy*s(coa, coset(0, 0, bz), cocy) 841 DO by = 2, lb 842 bz = lb - by 843 f3 = f2*REAL(by - 1, dp) 844 s(coa, coset(0, by, bz), coc) = & 845 rbg(2)*s(coa, coset(0, by - 1, bz), coc) + & 846 f3*s(coa, coset(0, by - 2, bz), coc) + & 847 fcy*s(coa, coset(0, by - 1, bz), cocy) 848 END DO 849 ELSE 850 bz = lb - 1 851 s(coa, coset(0, 1, bz), coc) = & 852 rbg(2)*s(coa, coset(0, 0, bz), coc) + & 853 fy*s(coay, coset(0, 0, bz), coc) + & 854 fcy*s(coa, coset(0, 0, bz), cocy) 855 DO by = 2, lb 856 bz = lb - by 857 f3 = f2*REAL(by - 1, dp) 858 s(coa, coset(0, by, bz), coc) = & 859 rbg(2)*s(coa, coset(0, by - 1, bz), coc) + & 860 fy*s(coay, coset(0, by - 1, bz), coc) + & 861 f3*s(coa, coset(0, by - 2, bz), coc) + & 862 fcy*s(coa, coset(0, by - 1, bz), cocy) 863 END DO 864 END IF 865 866! *** Shift of angular momentum *** 867! *** component x from a to b *** 868 869 IF (ax == 0) THEN 870 DO by = 0, lb - 1 871 bz = lb - 1 - by 872 s(coa, coset(1, by, bz), coc) = & 873 rbg(1)*s(coa, coset(0, by, bz), coc) + & 874 fcx*s(coa, coset(0, by, bz), cocx) 875 END DO 876 DO bx = 2, lb 877 f3 = f2*REAL(bx - 1, dp) 878 DO by = 0, lb - bx 879 bz = lb - bx - by 880 s(coa, coset(bx, by, bz), coc) = & 881 rbg(1)*s(coa, coset(bx - 1, by, bz), coc) + & 882 f3*s(coa, coset(bx - 2, by, bz), coc) + & 883 fcx*s(coa, coset(bx - 1, by, bz), cocx) 884 END DO 885 END DO 886 ELSE 887 DO by = 0, lb - 1 888 bz = lb - 1 - by 889 s(coa, coset(1, by, bz), coc) = & 890 rbg(1)*s(coa, coset(0, by, bz), coc) + & 891 fx*s(coax, coset(0, by, bz), coc) + & 892 fcx*s(coa, coset(0, by, bz), cocx) 893 END DO 894 DO bx = 2, lb 895 f3 = f2*REAL(bx - 1, dp) 896 DO by = 0, lb - bx 897 bz = lb - bx - by 898 s(coa, coset(bx, by, bz), coc) = & 899 rbg(1)*s(coa, coset(bx - 1, by, bz), coc) + & 900 fx*s(coax, coset(bx - 1, by, bz), coc) + & 901 f3*s(coa, coset(bx - 2, by, bz), coc) + & 902 fcx*s(coa, coset(bx - 1, by, bz), cocx) 903 END DO 904 END DO 905 END IF 906 907 END DO 908 END DO 909 910 END DO 911 912 END IF 913 914 ELSE 915 916 IF (lb_max > 0) THEN 917 918! *** Vertical recurrence steps: [s|c|s] -> [s|c|b] *** 919 920 rbg(:) = rcg(:) + rbc(:) 921 922! *** [s|c|p] = (Gi - Bi)*[s|c|s] + f2*Ni(c)*[s|c-1i|s] *** 923 924 s(1, 2, coc) = rbg(1)*s(1, 1, coc) + fcx*s(1, 1, cocx) 925 s(1, 3, coc) = rbg(2)*s(1, 1, coc) + fcy*s(1, 1, cocy) 926 s(1, 4, coc) = rbg(3)*s(1, 1, coc) + fcz*s(1, 1, cocz) 927 928! *** [s|c|b] = (Gi - Bi)*[s|c|b-1i] + *** 929! *** f2*Ni(b-1i)*[s|c|b-2i] *** 930! *** f2*Ni(c)*[s|c-1i|b-1i] *** 931 932 DO lb = 2, lb_max 933 934! *** Increase the angular momentum component z of b *** 935 936 s(1, coset(0, 0, lb), coc) = & 937 rbg(3)*s(1, coset(0, 0, lb - 1), coc) + & 938 f2*REAL(lb - 1, dp)*s(1, coset(0, 0, lb - 2), coc) + & 939 fcz*s(1, coset(0, 0, lb - 1), cocz) 940 941! *** Increase the angular momentum component y of b *** 942 943 bz = lb - 1 944 s(1, coset(0, 1, bz), coc) = & 945 rbg(2)*s(1, coset(0, 0, bz), coc) + & 946 fcy*s(1, coset(0, 0, bz), cocy) 947 948 DO by = 2, lb 949 bz = lb - by 950 s(1, coset(0, by, bz), coc) = & 951 rbg(2)*s(1, coset(0, by - 1, bz), coc) + & 952 f2*REAL(by - 1, dp)*s(1, coset(0, by - 2, bz), coc) + & 953 fcy*s(1, coset(0, by - 1, bz), cocy) 954 END DO 955 956! *** Increase the angular momentum component x of b *** 957 958 DO by = 0, lb - 1 959 bz = lb - 1 - by 960 s(1, coset(1, by, bz), coc) = & 961 rbg(1)*s(1, coset(0, by, bz), coc) + & 962 fcx*s(1, coset(0, by, bz), cocx) 963 END DO 964 965 DO bx = 2, lb 966 f3 = f2*REAL(bx - 1, dp) 967 DO by = 0, lb - bx 968 bz = lb - bx - by 969 s(1, coset(bx, by, bz), coc) = & 970 rbg(1)*s(1, coset(bx - 1, by, bz), coc) + & 971 f3*s(1, coset(bx - 2, by, bz), coc) + & 972 fcx*s(1, coset(bx - 1, by, bz), cocx) 973 END DO 974 END DO 975 976 END DO 977 978 END IF 979 980 END IF 981 982 END DO 983 END DO 984 985 END DO 986 987 END IF 988 989! *** Store integrals 990 991 IF (PRESENT(int_abc_ext)) THEN 992 DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set) 993 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set) 994 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set) 995 sabc(na + i, nb + j, nc + k) = s(i, j, k) 996 int_abc_ext = MAX(int_abc_ext, ABS(s(i, j, k))) 997 END DO 998 END DO 999 END DO 1000 ELSE 1001 DO k = ncoset(lc_min_set - 1) + 1, ncoset(lc_max_set) 1002 DO j = ncoset(lb_min_set - 1) + 1, ncoset(lb_max_set) 1003 DO i = ncoset(la_min_set - 1) + 1, ncoset(la_max_set) 1004 sabc(na + i, nb + j, nc + k) = s(i, j, k) 1005 END DO 1006 END DO 1007 END DO 1008 ENDIF 1009 1010! *** Calculate the requested derivatives with respect to *** 1011! *** the nuclear coordinates of the atomic center a and c *** 1012 1013 IF (PRESENT(sdabc) .OR. PRESENT(sabdc)) THEN 1014 CALL derivatives_overlap3(la_max_set, la_min_set, lb_max_set, lb_min_set, & 1015 lc_max_set, lc_min_set, zeta(ipgf), zetc(kpgf), & 1016 s, sda, sdc) 1017 ENDIF 1018 1019! *** Store the first derivatives of the primitive overlap integrals *** 1020 1021 IF (PRESENT(sdabc)) THEN 1022 DO k = 1, 3 1023 DO l = 1, ncoset(lc_max_set) 1024 DO j = 1, ncoset(lb_max_set) 1025 DO i = 1, ncoset(la_max_set) 1026 sdabc(nda + i, nb + j, nc + l, k) = sda(i, j, l, k) 1027 END DO 1028 END DO 1029 ENDDO 1030 END DO 1031 END IF 1032 1033 IF (PRESENT(sabdc)) THEN 1034 DO k = 1, 3 1035 DO l = 1, ncoset(lc_max_set) 1036 DO j = 1, ncoset(lb_max_set) 1037 DO i = 1, ncoset(la_max_set) 1038 sabdc(na + i, nb + j, ndc + l, k) = sdc(i, j, l, k) 1039 END DO 1040 END DO 1041 ENDDO 1042 END DO 1043 END IF 1044 1045 nc = nc + ncoset(lc_max_set) 1046 ndc = ndc + ncoset(lc_max_set) 1047 END DO 1048 1049 nb = nb + ncoset(lb_max) 1050 END DO 1051 1052 na = na + ncoset(la_max_set) 1053 nda = nda + ncoset(la_max_set) 1054 END DO 1055 1056 DEALLOCATE (s) 1057 IF (PRESENT(sdabc)) THEN 1058 DEALLOCATE (sda) 1059 ENDIF 1060 IF (PRESENT(sabdc)) THEN 1061 DEALLOCATE (sdc) 1062 ENDIF 1063 1064 CALL timestop(handle) 1065 1066 END SUBROUTINE overlap3 1067 1068! ************************************************************************************************** 1069!> \brief Calculates the derivatives of the three-center overlap integral [a|b|c] 1070!> with respect to the nuclear coordinates of the atomic center a and c 1071!> \param la_max_set ... 1072!> \param la_min_set ... 1073!> \param lb_max_set ... 1074!> \param lb_min_set ... 1075!> \param lc_max_set ... 1076!> \param lc_min_set ... 1077!> \param zeta ... 1078!> \param zetc ... 1079!> \param s integrals [a|b|c] 1080!> \param sda derivative [da/dAi|b|c] 1081!> \param sdc derivative [a|b|dc/dCi] 1082! ************************************************************************************************** 1083 SUBROUTINE derivatives_overlap3(la_max_set, la_min_set, lb_max_set, lb_min_set, & 1084 lc_max_set, lc_min_set, zeta, zetc, s, sda, sdc) 1085 1086 INTEGER, INTENT(IN) :: la_max_set, la_min_set, lb_max_set, & 1087 lb_min_set, lc_max_set, lc_min_set 1088 REAL(KIND=dp), INTENT(IN) :: zeta, zetc 1089 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: s 1090 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sda, sdc 1091 1092 CHARACTER(len=*), PARAMETER :: routineN = 'derivatives_overlap3', & 1093 routineP = moduleN//':'//routineN 1094 1095 INTEGER :: ax, ay, az, bx, by, bz, coa, coamx, coamy, coamz, coapx, coapy, coapz, cob, coc, & 1096 cocmx, cocmy, cocmz, cocpx, cocpy, cocpz, cx, cy, cz, devx, devy, devz, handle, la, lb, lc 1097 REAL(KIND=dp) :: fax, fay, faz, fcx, fcy, fcz, fexpa, & 1098 fexpc 1099 1100 CALL timeset(routineN, handle) 1101 1102 fexpa = 2.0_dp*zeta 1103 fexpc = 2.0_dp*zetc 1104 1105! derivative with respec to x,y,z 1106 1107 devx = 1 1108 devy = 2 1109 devz = 3 1110 1111! *** [da/dAi|b|c] = 2*zeta*[a+1i|b|c] - Ni(a)[a-1i|b|c] *** 1112! *** [a|b|dc/dCi] = 2*zetc*[a|b|c+1i] - Ni(c)[a|b|c-1i] *** 1113 1114 DO la = la_min_set, la_max_set 1115 DO ax = 0, la 1116 fax = REAL(ax, dp) 1117 DO ay = 0, la - ax 1118 fay = REAL(ay, dp) 1119 az = la - ax - ay 1120 faz = REAL(az, dp) 1121 coa = coset(ax, ay, az) 1122 coamx = coset(ax - 1, ay, az) 1123 coamy = coset(ax, ay - 1, az) 1124 coamz = coset(ax, ay, az - 1) 1125 coapx = coset(ax + 1, ay, az) 1126 coapy = coset(ax, ay + 1, az) 1127 coapz = coset(ax, ay, az + 1) 1128 DO lb = lb_min_set, lb_max_set 1129 DO bx = 0, lb 1130 DO by = 0, lb - bx 1131 bz = lb - bx - by 1132 cob = coset(bx, by, bz) 1133 DO lc = lc_min_set, lc_max_set 1134 DO cx = 0, lc 1135 fcx = REAL(cx, dp) 1136 DO cy = 0, lc - cx 1137 fcy = REAL(cy, dp) 1138 cz = lc - cx - cy 1139 fcz = REAL(cz, dp) 1140 coc = coset(cx, cy, cz) 1141 cocmx = coset(cx - 1, cy, cz) 1142 cocmy = coset(cx, cy - 1, cz) 1143 cocmz = coset(cx, cy, cz - 1) 1144 cocpx = coset(cx + 1, cy, cz) 1145 cocpy = coset(cx, cy + 1, cz) 1146 cocpz = coset(cx, cy, cz + 1) 1147 IF (ASSOCIATED(sda)) THEN 1148 sda(coa, cob, coc, devx) = fexpa*s(coapx, cob, coc) - & 1149 fax*s(coamx, cob, coc) 1150 sda(coa, cob, coc, devy) = fexpa*s(coapy, cob, coc) - & 1151 fay*s(coamy, cob, coc) 1152 sda(coa, cob, coc, devz) = fexpa*s(coapz, cob, coc) - & 1153 faz*s(coamz, cob, coc) 1154 ENDIF 1155 IF (ASSOCIATED(sdc)) THEN 1156 sdc(coa, cob, coc, devx) = fexpc*s(coa, cob, cocpx) - & 1157 fcx*s(coa, cob, cocmx) 1158 sdc(coa, cob, coc, devy) = fexpc*s(coa, cob, cocpy) - & 1159 fcy*s(coa, cob, cocmy) 1160 sdc(coa, cob, coc, devz) = fexpc*s(coa, cob, cocpz) - & 1161 fcz*s(coa, cob, cocmz) 1162 ENDIF 1163 ENDDO 1164 ENDDO 1165 ENDDO 1166 END DO 1167 END DO 1168 END DO 1169 END DO 1170 END DO 1171 END DO 1172 1173 CALL timestop(handle) 1174 1175 END SUBROUTINE derivatives_overlap3 1176 1177END MODULE ai_overlap3 1178