LCOV - code coverage report
Current view: top level - src - debug_os_integrals.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 117 270 43.3 %
Date: 2024-03-28 07:31:50 Functions: 3 6 50.0 %

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

Generated by: LCOV version 1.15