1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Debugs Obara-Saika integral matrices 8!> \par History 9!> created [07.2014] 10!> \authors Dorothea Golze 11! ************************************************************************************************** 12MODULE debug_os_integrals 13 14 USE ai_overlap, ONLY: overlap 15 USE ai_overlap3, ONLY: overlap3 16 USE ai_overlap3_debug, ONLY: init_os_overlap3,& 17 os_overlap3 18 USE ai_overlap_aabb, ONLY: overlap_aabb 19 USE ai_overlap_debug, ONLY: init_os_overlap2,& 20 os_overlap2 21 USE kinds, ONLY: dp 22 USE orbital_pointers, ONLY: coset,& 23 indco,& 24 ncoset 25#include "./base/base_uses.f90" 26 27 IMPLICIT NONE 28 29 PRIVATE 30 31! ************************************************************************************************** 32 33 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'debug_os_integrals' 34 35 PUBLIC :: overlap_ab_test, overlap_abc_test, overlap_aabb_test 36 37! ************************************************************************************************** 38 39CONTAINS 40 41! *************************************************************************************************** 42!> \brief recursive test routines for integral (a,b) for only two exponents 43! ************************************************************************************************** 44 SUBROUTINE overlap_ab_test_simple() 45 46 CHARACTER(LEN=*), PARAMETER :: routineN = 'overlap_ab_test_simple', & 47 routineP = moduleN//':'//routineN 48 49 INTEGER :: ia1, iax, iay, iaz, ib1, ibx, iby, ibz, & 50 la_max, la_min, lb_max, lb_min, lds, & 51 ma, maxl, mb 52 INTEGER, DIMENSION(3) :: na, nb 53 REAL(KIND=dp) :: dab, dmax, res1, xa, xb 54 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab 55 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: swork 56 REAL(KIND=dp), DIMENSION(1) :: rpgfa, rpgfb, xa_work, xb_work 57 REAL(KIND=dp), DIMENSION(3) :: A, B, rab 58 59 xa = 0.783300000000_dp ! exponents 60 xb = 1.239648746700_dp 61 62 A = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions 63 B = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr 64 65 la_min = 0 66 lb_min = 0 67 68 la_max = 3 69 lb_max = 4 70 71 !--------------------------------------- 72 ALLOCATE (sab(ncoset(la_max), ncoset(lb_max))) 73 74 maxl = MAX(la_max, lb_max) 75 lds = ncoset(maxl) 76 ALLOCATE (swork(lds, lds, 1)) 77 sab = 0._dp 78 rab(:) = B(:) - A(:) 79 dab = SQRT(DOT_PRODUCT(rab, rab)) 80 xa_work(1) = xa 81 xb_work(1) = xb 82 rpgfa = 20._dp 83 rpgfb = 20._dp 84 CALL overlap(la_max_set=la_max, la_min_set=la_min, npgfa=1, rpgfa=rpgfa, zeta=xa_work, & 85 lb_max_set=lb_max, lb_min_set=lb_min, npgfb=1, rpgfb=rpgfb, zetb=xb_work, & 86 rab=rab, dab=dab, sab=sab, da_max_set=0, return_derivatives=.FALSE., s=swork, lds=lds) 87 !--------------------------------------- 88 89 CALL init_os_overlap2(xa, xb, A, B) 90 91 dmax = 0._dp 92 DO ma = la_min, la_max 93 DO mb = lb_min, lb_max 94 DO iax = 0, ma 95 DO iay = 0, ma - iax 96 iaz = ma - iax - iay 97 na(1) = iax; na(2) = iay; na(3) = iaz 98 ia1 = coset(iax, iay, iaz) 99 DO ibx = 0, mb 100 DO iby = 0, mb - ibx 101 ibz = mb - ibx - iby 102 nb(1) = ibx; nb(2) = iby; nb(3) = ibz 103 ib1 = coset(ibx, iby, ibz) 104 res1 = os_overlap2(na, nb) 105 ! write(*,*) "la, lb,na, nb, res1", ma, mb, na, nb, res1 106 ! write(*,*) "sab ia1, ib1", ia1, ib1, sab(ia1,ib1) 107 dmax = MAX(dmax, ABS(res1 - sab(ia1, ib1))) 108 END DO 109 END DO 110 END DO 111 END DO 112 END DO 113 END DO 114 115 DEALLOCATE (sab, swork) 116 117 END SUBROUTINE overlap_ab_test_simple 118 119! *************************************************************************************************** 120!> \brief recursive test routines for integral (a,b) 121!> \param la_max ... 122!> \param la_min ... 123!> \param npgfa ... 124!> \param zeta ... 125!> \param lb_max ... 126!> \param lb_min ... 127!> \param npgfb ... 128!> \param zetb ... 129!> \param ra ... 130!> \param rb ... 131!> \param sab ... 132!> \param dmax ... 133! ************************************************************************************************** 134 SUBROUTINE overlap_ab_test(la_max, la_min, npgfa, zeta, lb_max, lb_min, npgfb, zetb, & 135 ra, rb, sab, dmax) 136 137 INTEGER, INTENT(IN) :: la_max, la_min, npgfa 138 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 139 INTEGER, INTENT(IN) :: lb_max, lb_min, npgfb 140 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 141 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb 142 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sab 143 REAL(KIND=dp), INTENT(INOUT) :: dmax 144 145 CHARACTER(LEN=*), PARAMETER :: routineN = 'overlap_ab_test', & 146 routineP = moduleN//':'//routineN 147 148 INTEGER :: coa, cob, ia1, iax, iay, iaz, ib1, ibx, & 149 iby, ibz, ipgf, jpgf, ma, mb 150 INTEGER, DIMENSION(3) :: na, nb 151 REAL(KIND=dp) :: res1, res2, xa, xb 152 REAL(KIND=dp), DIMENSION(3) :: A, B 153 154 coa = 0 155 DO ipgf = 1, npgfa 156 cob = 0 157 DO jpgf = 1, npgfb 158 xa = zeta(ipgf) !exponents 159 xb = zetb(jpgf) 160 A = ra !positions 161 B = rb 162 CALL init_os_overlap2(xa, xb, A, B) 163 DO ma = la_min, la_max 164 DO mb = lb_min, lb_max 165 DO iax = 0, ma 166 DO iay = 0, ma - iax 167 iaz = ma - iax - iay 168 na(1) = iax; na(2) = iay; na(3) = iaz 169 ia1 = coset(iax, iay, iaz) 170 DO ibx = 0, mb 171 DO iby = 0, mb - ibx 172 ibz = mb - ibx - iby 173 nb(1) = ibx; nb(2) = iby; nb(3) = ibz 174 ib1 = coset(ibx, iby, ibz) 175 res1 = os_overlap2(na, nb) 176 res2 = sab(coa + ia1, cob + ib1) 177 dmax = MAX(dmax, ABS(res1 - res2)) 178 END DO 179 END DO 180 END DO 181 END DO 182 END DO 183 END DO 184 cob = cob + ncoset(lb_max) 185 END DO 186 coa = coa + ncoset(la_max) 187 END DO 188 !WRITE(*,*) "dmax overlap_ab_test", dmax 189 190 END SUBROUTINE overlap_ab_test 191 192! *************************************************************************************************** 193!> \brief recursive test routines for integral (a,b,c) for only three exponents 194! ************************************************************************************************** 195 SUBROUTINE overlap_abc_test_simple() 196 197 CHARACTER(LEN=*), PARAMETER :: routineN = 'overlap_abc_test_simple', & 198 routineP = moduleN//':'//routineN 199 200 INTEGER :: ia1, iax, iay, iaz, ib1, ibx, iby, ibz, & 201 ic1, icx, icy, icz, la_max, la_min, & 202 lb_max, lb_min, lc_max, lc_min, ma, & 203 mb, mc 204 INTEGER, DIMENSION(3) :: na, nb, nc 205 REAL(KIND=dp) :: dab, dac, dbc, dmax, res1, xa, xb, xc 206 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sabc 207 REAL(KIND=dp), DIMENSION(1) :: rpgfa, rpgfb, rpgfc, xa_work, xb_work, & 208 xc_work 209 REAL(KIND=dp), DIMENSION(3) :: A, B, C, rab, rac, rbc 210 211 xa = 0.783300000000_dp ! exponents 212 xb = 1.239648746700_dp 213 xc = 0.548370000000_dp 214 215 A = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions 216 B = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr 217 C = (/0.032380000000_dp, 1.23470000000_dp, 0.11137400000_dp/) !* bohr 218 219 la_min = 0 220 lb_min = 0 221 lc_min = 0 222 223 la_max = 0 224 lb_max = 0 225 lc_max = 1 226 227 !--------------------------------------- 228 rab(:) = B(:) - A(:) 229 dab = SQRT(DOT_PRODUCT(rab, rab)) 230 rac(:) = C(:) - A(:) 231 dac = SQRT(DOT_PRODUCT(rac, rac)) 232 rbc(:) = C(:) - B(:) 233 dbc = SQRT(DOT_PRODUCT(rbc, rbc)) 234 ALLOCATE (sabc(ncoset(la_max), ncoset(lb_max), ncoset(lc_max))) 235 xa_work(1) = xa 236 xb_work(1) = xb 237 xc_work(1) = xc 238 rpgfa = 20._dp 239 rpgfb = 20._dp 240 rpgfc = 20._dp 241 sabc = 0._dp 242 243 CALL overlap3(la_max_set=la_max, npgfa=1, zeta=xa_work, rpgfa=rpgfa, la_min_set=la_min, & 244 lb_max_set=lb_max, npgfb=1, zetb=xb_work, rpgfb=rpgfb, lb_min_set=lb_min, & 245 lc_max_set=lc_max, npgfc=1, zetc=xc_work, rpgfc=rpgfc, lc_min_set=lc_min, & 246 rab=rab, dab=dab, rac=rac, dac=dac, rbc=rbc, dbc=dbc, sabc=sabc) 247 248 !--------------------------------------- 249 250 CALL init_os_overlap3(xa, xb, xc, A, B, C) 251 252 dmax = 0._dp 253 DO ma = la_min, la_max 254 DO mc = lc_min, lc_max 255 DO mb = lb_min, lb_max 256 DO iax = 0, ma 257 DO iay = 0, ma - iax 258 iaz = ma - iax - iay 259 na(1) = iax; na(2) = iay; na(3) = iaz 260 ia1 = coset(iax, iay, iaz) 261 DO icx = 0, mc 262 DO icy = 0, mc - icx 263 icz = mc - icx - icy 264 nc(1) = icx; nc(2) = icy; nc(3) = icz 265 ic1 = coset(icx, icy, icz) 266 DO ibx = 0, mb 267 DO iby = 0, mb - ibx 268 ibz = mb - ibx - iby 269 nb(1) = ibx; nb(2) = iby; nb(3) = ibz 270 ib1 = coset(ibx, iby, ibz) 271 res1 = os_overlap3(na, nc, nb) 272 !write(*,*) "la, lc, lb,na,nc, nb, res1", ma, mc, mb, na, nc, nb, res1 273 !write(*,*) "sabc ia1, ib1, ic1", ia1, ib1, ic1, sabc(ia1,ib1,ic1) 274 dmax = MAX(dmax, ABS(res1 - sabc(ia1, ib1, ic1))) 275 END DO 276 END DO 277 END DO 278 END DO 279 END DO 280 END DO 281 END DO 282 END DO 283 END DO 284 285 DEALLOCATE (sabc) 286 287 END SUBROUTINE overlap_abc_test_simple 288 289! *************************************************************************************************** 290!> \brief recursive test routines for integral (a,b,c) 291!> \param la_max ... 292!> \param npgfa ... 293!> \param zeta ... 294!> \param la_min ... 295!> \param lb_max ... 296!> \param npgfb ... 297!> \param zetb ... 298!> \param lb_min ... 299!> \param lc_max ... 300!> \param npgfc ... 301!> \param zetc ... 302!> \param lc_min ... 303!> \param ra ... 304!> \param rb ... 305!> \param rc ... 306!> \param sabc ... 307!> \param dmax ... 308! ************************************************************************************************** 309 SUBROUTINE overlap_abc_test(la_max, npgfa, zeta, la_min, & 310 lb_max, npgfb, zetb, lb_min, & 311 lc_max, npgfc, zetc, lc_min, & 312 ra, rb, rc, sabc, dmax) 313 314 INTEGER, INTENT(IN) :: la_max, npgfa 315 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta 316 INTEGER, INTENT(IN) :: la_min, lb_max, npgfb 317 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb 318 INTEGER, INTENT(IN) :: lb_min, lc_max, npgfc 319 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetc 320 INTEGER, INTENT(IN) :: lc_min 321 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb, rc 322 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sabc 323 REAL(KIND=dp), INTENT(INOUT) :: dmax 324 325 CHARACTER(LEN=*), PARAMETER :: routineN = 'overlap_abc_test', & 326 routineP = moduleN//':'//routineN 327 328 INTEGER :: coa, cob, coc, ia1, iax, iay, iaz, ib1, & 329 ibx, iby, ibz, ic1, icx, icy, icz, & 330 ipgf, jpgf, kpgf, ma, mb, mc 331 INTEGER, DIMENSION(3) :: na, nb, nc 332 REAL(KIND=dp) :: res1, res2, xa, xb, xc 333 REAL(KIND=dp), DIMENSION(3) :: A, B, C 334 335 coa = 0 336 DO ipgf = 1, npgfa 337 cob = 0 338 DO jpgf = 1, npgfb 339 coc = 0 340 DO kpgf = 1, npgfc 341 342 xa = zeta(ipgf) ! exponents 343 xb = zetb(jpgf) 344 xc = zetc(kpgf) 345 346 A = Ra !positions 347 B = Rb 348 C = Rc 349 350 CALL init_os_overlap3(xa, xb, xc, A, B, C) 351 352 DO ma = la_min, la_max 353 DO mc = lc_min, lc_max 354 DO mb = lb_min, lb_max 355 DO iax = 0, ma 356 DO iay = 0, ma - iax 357 iaz = ma - iax - iay 358 na(1) = iax; na(2) = iay; na(3) = iaz 359 ia1 = coset(iax, iay, iaz) 360 DO icx = 0, mc 361 DO icy = 0, mc - icx 362 icz = mc - icx - icy 363 nc(1) = icx; nc(2) = icy; nc(3) = icz 364 ic1 = coset(icx, icy, icz) 365 DO ibx = 0, mb 366 DO iby = 0, mb - ibx 367 ibz = mb - ibx - iby 368 nb(1) = ibx; nb(2) = iby; nb(3) = ibz 369 ib1 = coset(ibx, iby, ibz) 370 res1 = os_overlap3(na, nc, nb) 371 res2 = sabc(coa + ia1, cob + ib1, coc + ic1) 372 dmax = MAX(dmax, ABS(res1 - res2)) 373 !IF(dmax > 1.E-10) WRITE(*,*) "dmax in loop", dmax 374 END DO 375 END DO 376 END DO 377 END DO 378 END DO 379 END DO 380 END DO 381 END DO 382 END DO 383 coc = coc + ncoset(lc_max) 384 END DO 385 cob = cob + ncoset(lb_max) 386 END DO 387 coa = coa + ncoset(la_max) 388 END DO 389 !WRITE(*,*) "dmax abc", dmax 390 391 END SUBROUTINE overlap_abc_test 392 393! *************************************************************************************************** 394!> \brief recursive test routines for integral (aa,bb) for only four exponents 395! ************************************************************************************************** 396 SUBROUTINE overlap_aabb_test_simple() 397 398 CHARACTER(LEN=*), PARAMETER :: routineN = 'overlap_aabb_test_simple', & 399 routineP = moduleN//':'//routineN 400 401 INTEGER :: i, iax, iay, iaz, ibx, iby, ibz, j, k, l, la_max, la_max1, la_max2, la_min, & 402 la_min1, la_min2, lb_max, lb_max1, lb_max2, lb_min, lb_min1, lb_min2, lds, ma, maxl, mb 403 INTEGER, DIMENSION(3) :: na, naa, nb, nbb 404 REAL(KIND=dp) :: dab, dmax, res1, xa, xa1, xa2, xb, xb1, & 405 xb2 406 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: swork 407 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: saabb 408 REAL(KIND=dp), DIMENSION(1) :: rpgfa1, rpgfa2, rpgfb1, rpgfb2, & 409 xa_work1, xa_work2, xb_work1, xb_work2 410 REAL(KIND=dp), DIMENSION(3) :: A, B, rab 411 412 xa1 = 0.783300000000_dp ! exponents 413 xb1 = 1.239648746700_dp 414 xa2 = 0.3400000000_dp ! exponents 415 xb2 = 2.76_dp 416 417 A = (/0.329309000000_dp, 0.28408240000_dp, 0.28408240000_dp/) !* bohr !positions 418 B = (/0.983983000000_dp, 0.00453720000_dp, 0.00432740000_dp/) !* bohr 419 420 la_min1 = 0 421 la_min2 = 0 422 lb_min1 = 3 423 lb_min2 = 1 424 425 la_max1 = 1 426 la_max2 = 2 427 lb_max1 = 3 428 lb_max2 = 4 429 430 !--------------------------------------- 431 ALLOCATE (saabb(ncoset(la_max1), ncoset(la_max2), ncoset(lb_max1), ncoset(lb_max2))) 432 433 maxl = MAX(la_max1 + la_max2, lb_max1 + lb_max2) 434 lds = ncoset(maxl) 435 ALLOCATE (swork(lds, lds)) 436 saabb = 0._dp 437 rab(:) = B(:) - A(:) 438 dab = SQRT(DOT_PRODUCT(rab, rab)) 439 xa_work1(1) = xa1 440 xa_work2(1) = xa2 441 xb_work1(1) = xb1 442 xb_work2(1) = xb2 443 rpgfa1 = 20._dp 444 rpgfa2 = 20._dp 445 rpgfb1 = 20._dp 446 rpgfb2 = 20._dp 447 CALL overlap_aabb(la_max_set1=la_max1, la_min_set1=la_min1, npgfa1=1, rpgfa1=rpgfa1, zeta1=xa_work1, & 448 la_max_set2=la_max2, la_min_set2=la_min2, npgfa2=1, rpgfa2=rpgfa2, zeta2=xa_work2, & 449 lb_max_set1=lb_max1, lb_min_set1=lb_min1, npgfb1=1, rpgfb1=rpgfb1, zetb1=xb_work1, & 450 lb_max_set2=lb_max2, lb_min_set2=lb_min2, npgfb2=1, rpgfb2=rpgfb2, zetb2=xb_work2, & 451 asets_equal=.FALSE., bsets_equal=.FALSE., rab=rab, dab=dab, saabb=saabb, s=swork, lds=lds) 452 !--------------------------------------- 453 454 xa = xa1 + xa2 455 xb = xb1 + xb2 456 la_min = la_min1 + la_min2 457 la_max = la_max1 + la_max2 458 lb_min = lb_min1 + lb_min2 459 lb_max = lb_max1 + lb_max2 460 461 CALL init_os_overlap2(xa, xb, A, B) 462 463 dmax = 0._dp 464 DO ma = la_min, la_max 465 DO mb = lb_min, lb_max 466 DO iax = 0, ma 467 DO iay = 0, ma - iax 468 iaz = ma - iax - iay 469 na(1) = iax; na(2) = iay; na(3) = iaz 470 DO ibx = 0, mb 471 DO iby = 0, mb - ibx 472 ibz = mb - ibx - iby 473 nb(1) = ibx; nb(2) = iby; nb(3) = ibz 474 res1 = os_overlap2(na, nb) 475 DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1) 476 DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2) 477 naa = indco(1:3, i) + indco(1:3, j) 478 DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1) 479 DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2) 480 nbb = indco(1:3, k) + indco(1:3, l) 481 IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN 482 dmax = MAX(dmax, ABS(res1 - saabb(i, j, k, l))) 483 END IF 484 END DO 485 END DO 486 END DO 487 END DO 488 END DO 489 END DO 490 END DO 491 END DO 492 END DO 493 END DO 494 495 DEALLOCATE (saabb, swork) 496 497 END SUBROUTINE overlap_aabb_test_simple 498 499! *************************************************************************************************** 500!> \brief recursive test routines for integral (aa,bb) 501!> \param la_max1 ... 502!> \param la_min1 ... 503!> \param npgfa1 ... 504!> \param zeta1 ... 505!> \param la_max2 ... 506!> \param la_min2 ... 507!> \param npgfa2 ... 508!> \param zeta2 ... 509!> \param lb_max1 ... 510!> \param lb_min1 ... 511!> \param npgfb1 ... 512!> \param zetb1 ... 513!> \param lb_max2 ... 514!> \param lb_min2 ... 515!> \param npgfb2 ... 516!> \param zetb2 ... 517!> \param ra ... 518!> \param rb ... 519!> \param saabb ... 520!> \param dmax ... 521! ************************************************************************************************** 522 SUBROUTINE overlap_aabb_test(la_max1, la_min1, npgfa1, zeta1, & 523 la_max2, la_min2, npgfa2, zeta2, & 524 lb_max1, lb_min1, npgfb1, zetb1, & 525 lb_max2, lb_min2, npgfb2, zetb2, & 526 ra, rb, saabb, dmax) 527 528 INTEGER, INTENT(IN) :: la_max1, la_min1, npgfa1 529 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta1 530 INTEGER, INTENT(IN) :: la_max2, la_min2, npgfa2 531 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zeta2 532 INTEGER, INTENT(IN) :: lb_max1, lb_min1, npgfb1 533 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb1 534 INTEGER, INTENT(IN) :: lb_max2, lb_min2, npgfb2 535 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: zetb2 536 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb 537 REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN) :: saabb 538 REAL(KIND=dp), INTENT(INOUT) :: dmax 539 540 CHARACTER(LEN=*), PARAMETER :: routineN = 'overlap_aabb_test', & 541 routineP = moduleN//':'//routineN 542 543 INTEGER :: coa1, coa2, cob1, cob2, i, iax, iay, & 544 iaz, ibx, iby, ibz, ipgf, j, jpgf, k, & 545 kpgf, l, la_max, la_min, lb_max, & 546 lb_min, lpgf, ma, mb 547 INTEGER, DIMENSION(3) :: na, naa, nb, nbb 548 REAL(KIND=dp) :: res1, xa, xb 549 REAL(KIND=dp), DIMENSION(3) :: A, B 550 551 coa1 = 0 552 DO ipgf = 1, npgfa1 553 coa2 = 0 554 DO jpgf = 1, npgfa2 555 cob1 = 0 556 DO kpgf = 1, npgfb1 557 cob2 = 0 558 DO lpgf = 1, npgfb2 559 560 xa = zeta1(ipgf) + zeta2(jpgf) ! exponents 561 xb = zetb1(kpgf) + zetb2(lpgf) ! exponents 562 la_max = la_max1 + la_max2 563 lb_max = lb_max1 + lb_max2 564 la_min = la_min1 + la_min2 565 lb_min = lb_min1 + lb_min2 566 567 A = ra !positions 568 B = rb 569 570 CALL init_os_overlap2(xa, xb, A, B) 571 572 DO ma = la_min, la_max 573 DO mb = lb_min, lb_max 574 DO iax = 0, ma 575 DO iay = 0, ma - iax 576 iaz = ma - iax - iay 577 na(1) = iax; na(2) = iay; na(3) = iaz 578 DO ibx = 0, mb 579 DO iby = 0, mb - ibx 580 ibz = mb - ibx - iby 581 nb(1) = ibx; nb(2) = iby; nb(3) = ibz 582 res1 = os_overlap2(na, nb) 583 DO i = ncoset(la_min1 - 1) + 1, ncoset(la_max1) 584 DO j = ncoset(la_min2 - 1) + 1, ncoset(la_max2) 585 naa = indco(1:3, i) + indco(1:3, j) 586 DO k = ncoset(lb_min1 - 1) + 1, ncoset(lb_max1) 587 DO l = ncoset(lb_min2 - 1) + 1, ncoset(lb_max2) 588 nbb = indco(1:3, k) + indco(1:3, l) 589 IF (ALL(na == naa) .AND. ALL(nb == nbb)) THEN 590 dmax = MAX(dmax, ABS(res1 - saabb(coa1 + i, coa2 + j, cob1 + k, cob2 + l))) 591 ENDIF 592 END DO 593 END DO 594 END DO 595 END DO 596 END DO 597 END DO 598 END DO 599 END DO 600 END DO 601 END DO 602 cob2 = cob2 + ncoset(lb_max2) 603 END DO 604 cob1 = cob1 + ncoset(lb_max1) 605 END DO 606 coa2 = coa2 + ncoset(la_max2) 607 END DO 608 coa1 = coa1 + ncoset(la_max1) 609 END DO 610 611 !WRITE(*,*) "dmax aabb", dmax 612 613 END SUBROUTINE overlap_aabb_test 614 615END MODULE debug_os_integrals 616