LCOV - code coverage report
Current view: top level - src - generic_os_integrals.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 442 446 99.1 %
Date: 2024-03-28 07:31:50 Functions: 7 7 100.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 Calculation of contracted, spherical Gaussian integrals using the (OS) integral
      10             : !>        scheme. Routines for the following two-center integrals:
      11             : !>        i)  (a|O(r12)|b) where O(r12) is the overlap, coulomb operator etc.
      12             : !>        ii) (aba) and (abb) s-overlaps
      13             : !> \par History
      14             : !>      created [06.2015]
      15             : !>      05.2019: Added truncated coulomb operator (A. Bussy)
      16             : !> \author Dorothea Golze
      17             : ! **************************************************************************************************
      18             : MODULE generic_os_integrals
      19             :    USE ai_contraction_sphi,             ONLY: ab_contract,&
      20             :                                               abc_contract,&
      21             :                                               abcd_contract
      22             :    USE ai_derivatives,                  ONLY: dabdr_noscreen
      23             :    USE ai_operator_ra2m,                ONLY: operator_ra2m
      24             :    USE ai_operators_r12,                ONLY: ab_sint_os,&
      25             :                                               cps_coulomb2,&
      26             :                                               cps_gauss2,&
      27             :                                               cps_truncated2,&
      28             :                                               cps_verf2,&
      29             :                                               cps_verfc2,&
      30             :                                               cps_vgauss2,&
      31             :                                               operator2
      32             :    USE ai_overlap,                      ONLY: overlap_aab,&
      33             :                                               overlap_ab,&
      34             :                                               overlap_abb
      35             :    USE ai_overlap_aabb,                 ONLY: overlap_aabb
      36             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      37             :                                               gto_basis_set_type
      38             :    USE constants_operator,              ONLY: operator_coulomb,&
      39             :                                               operator_gauss,&
      40             :                                               operator_truncated,&
      41             :                                               operator_verf,&
      42             :                                               operator_verfc,&
      43             :                                               operator_vgauss
      44             :    USE debug_os_integrals,              ONLY: overlap_aabb_test,&
      45             :                                               overlap_ab_test,&
      46             :                                               overlap_abc_test
      47             :    USE kinds,                           ONLY: dp
      48             :    USE orbital_pointers,                ONLY: ncoset
      49             : #include "./base/base_uses.f90"
      50             : 
      51             :    IMPLICIT NONE
      52             : 
      53             :    PRIVATE
      54             : 
      55             : ! **************************************************************************************************
      56             : 
      57             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_os_integrals'
      58             : 
      59             :    PUBLIC :: int_operators_r12_ab_os, int_overlap_ab_os, int_ra2m_ab_os, int_overlap_aba_os, &
      60             :              int_overlap_abb_os, int_overlap_aabb_os
      61             : 
      62             : CONTAINS
      63             : 
      64             : ! **************************************************************************************************
      65             : !> \brief Calcululates the two-center integrals of the type (a|O(r12)|b) using the OS scheme
      66             : !> \param r12_operator the integral operator, which depends on r12=|r1-r2|
      67             : !> \param vab integral matrix of spherical contracted Gaussian functions
      68             : !> \param dvab derivative of the integrals
      69             : !> \param rab distance vector between center A and B
      70             : !> \param fba basis at center A
      71             : !> \param fbb basis at center B
      72             : !> \param omega parameter in the operator
      73             : !> \param r_cutoff the cutoff in case of truncated coulomb operator
      74             : !> \param calculate_forces ...
      75             : ! **************************************************************************************************
      76         196 :    SUBROUTINE int_operators_r12_ab_os(r12_operator, vab, dvab, rab, fba, fbb, omega, r_cutoff, &
      77             :                                       calculate_forces)
      78             : 
      79             :       INTEGER, INTENT(IN)                                :: r12_operator
      80             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
      81             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
      82             :          OPTIONAL                                        :: dvab
      83             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
      84             :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
      85             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: omega, r_cutoff
      86             :       LOGICAL, INTENT(IN)                                :: calculate_forces
      87             : 
      88             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operators_r12_ab_os'
      89             : 
      90             :       INTEGER                                            :: handle
      91             :       REAL(KIND=dp)                                      :: my_omega, my_r_cutoff
      92             : 
      93             :       PROCEDURE(ab_sint_os), POINTER                     :: cps_operator2
      94             : 
      95         196 :       NULLIFY (cps_operator2)
      96         196 :       CALL timeset(routineN, handle)
      97             : 
      98         196 :       my_omega = 1.0_dp
      99         196 :       my_r_cutoff = 1.0_dp
     100             : 
     101         264 :       SELECT CASE (r12_operator)
     102             :       CASE (operator_coulomb)
     103          68 :          cps_operator2 => cps_coulomb2
     104             :       CASE (operator_verf)
     105          32 :          cps_operator2 => cps_verf2
     106          32 :          IF (PRESENT(omega)) my_omega = omega
     107             :       CASE (operator_verfc)
     108          32 :          cps_operator2 => cps_verfc2
     109          32 :          IF (PRESENT(omega)) my_omega = omega
     110             :       CASE (operator_vgauss)
     111          32 :          cps_operator2 => cps_vgauss2
     112          32 :          IF (PRESENT(omega)) my_omega = omega
     113             :       CASE (operator_gauss)
     114          32 :          cps_operator2 => cps_gauss2
     115          32 :          IF (PRESENT(omega)) my_omega = omega
     116             :       CASE (operator_truncated)
     117           0 :          cps_operator2 => cps_truncated2
     118           0 :          IF (PRESENT(r_cutoff)) my_r_cutoff = r_cutoff
     119             :       CASE DEFAULT
     120         196 :          CPABORT("Operator not available")
     121             :       END SELECT
     122             : 
     123             :       CALL int_operator_ab_os_low(cps_operator2, vab, dvab, rab, fba, fbb, my_omega, my_r_cutoff, &
     124         232 :                                   calculate_forces)
     125             : 
     126         196 :       CALL timestop(handle)
     127             : 
     128         196 :    END SUBROUTINE int_operators_r12_ab_os
     129             : 
     130             : ! **************************************************************************************************
     131             : !> \brief calculate integrals (a|O(r12)|b)
     132             : !> \param cps_operator2 procedure pointer for the respective operator.
     133             : !> \param vab integral matrix of spherical contracted Gaussian functions
     134             : !> \param dvab derivative of the integrals
     135             : !> \param rab distance vector between center A and B
     136             : !> \param fba basis at center A
     137             : !> \param fbb basis at center B
     138             : !> \param omega parameter in the operator
     139             : !> \param r_cutoff ...
     140             : !> \param calculate_forces ...
     141             : ! **************************************************************************************************
     142         196 :    SUBROUTINE int_operator_ab_os_low(cps_operator2, vab, dvab, rab, fba, fbb, omega, r_cutoff, &
     143             :                                      calculate_forces)
     144             : 
     145             :       PROCEDURE(ab_sint_os), POINTER                     :: cps_operator2
     146             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: vab
     147             :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     148             :          INTENT(INOUT)                                   :: dvab
     149             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     150             :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     151             :       REAL(KIND=dp), INTENT(IN)                          :: omega, r_cutoff
     152             :       LOGICAL, INTENT(IN)                                :: calculate_forces
     153             : 
     154             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'int_operator_ab_os_low'
     155             : 
     156             :       INTEGER                                            :: handle, i, iset, jset, lds, m1, m2, &
     157             :                                                             maxco, maxcoa, maxcob, maxl, maxla, &
     158             :                                                             maxlb, ncoa, ncoap, ncob, ncobp, &
     159             :                                                             nseta, nsetb, sgfa, sgfb
     160         196 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     161         196 :                                                             npgfb, nsgfa, nsgfb
     162         196 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     163             :       REAL(KIND=dp)                                      :: dab, rab2
     164         196 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: vac, vac_plus
     165         196 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: devab, vwork
     166         196 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: sphi_a, sphi_b, zeta, zetb
     167             : 
     168         196 :       CALL timeset(routineN, handle)
     169         196 :       NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
     170         196 :                first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb)
     171             : 
     172             :       ! basis ikind
     173         196 :       first_sgfa => fba%first_sgf
     174         196 :       la_max => fba%lmax
     175         196 :       la_min => fba%lmin
     176         196 :       npgfa => fba%npgf
     177         196 :       nseta = fba%nset
     178         196 :       nsgfa => fba%nsgf_set
     179         196 :       sphi_a => fba%sphi
     180         196 :       zeta => fba%zet
     181             :       ! basis jkind
     182         196 :       first_sgfb => fbb%first_sgf
     183         196 :       lb_max => fbb%lmax
     184         196 :       lb_min => fbb%lmin
     185         196 :       npgfb => fbb%npgf
     186         196 :       nsetb = fbb%nset
     187         196 :       nsgfb => fbb%nsgf_set
     188         196 :       sphi_b => fbb%sphi
     189         196 :       zetb => fbb%zet
     190             : 
     191         196 :       CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
     192         196 :       CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
     193         196 :       maxco = MAX(maxcoa, maxcob)
     194         196 :       IF (calculate_forces) THEN
     195         784 :          maxl = MAX(maxla + 1, maxlb)
     196             :       ELSE
     197         784 :          maxl = MAX(maxla, maxlb)
     198             :       END IF
     199         784 :       lds = ncoset(maxl)
     200             : 
     201         784 :       rab2 = SUM(rab*rab)
     202             :       dab = SQRT(rab2)
     203             : 
     204     5241692 :       vab = 0.0_dp
     205    14458436 :       IF (calculate_forces) dvab = 0.0_dp
     206             : 
     207        1906 :       DO iset = 1, nseta
     208             : 
     209        1710 :          ncoa = npgfa(iset)*ncoset(la_max(iset))
     210        1710 :          ncoap = npgfa(iset)*ncoset(la_max(iset) + 1)
     211        1710 :          sgfa = first_sgfa(1, iset)
     212             : 
     213       26924 :          DO jset = 1, nsetb
     214             : 
     215       25018 :             ncob = npgfb(jset)*ncoset(lb_max(jset))
     216       25018 :             ncobp = npgfb(jset)*ncoset(lb_max(jset) + 1)
     217       25018 :             sgfb = first_sgfb(1, jset)
     218       25018 :             m1 = sgfa + nsgfa(iset) - 1
     219       25018 :             m2 = sgfb + nsgfb(jset) - 1
     220             : 
     221             :             ! calculate integrals
     222       25018 :             IF (calculate_forces) THEN
     223           0 :                ALLOCATE (vwork(ncoap, ncobp, la_max(iset) + lb_max(jset) + 3), &
     224      360000 :                          vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3))
     225    25392000 :                devab = 0._dp
     226   232342880 :                vwork = 0.0_dp
     227     8456000 :                vac = 0.0_dp
     228             :                CALL operator2(cps_operator2, la_max(iset) + 1, npgfa(iset), zeta(:, iset), la_min(iset), &
     229             :                               lb_max(jset) + 1, npgfb(jset), zetb(:, jset), lb_min(jset), &
     230       24000 :                               omega, r_cutoff, rab, rab2, vac, vwork, maxder=1, vac_plus=vac_plus)
     231             :                CALL dabdr_noscreen(la_max(iset), npgfa(iset), zeta(:, iset), lb_max(jset), npgfb(jset), &
     232       24000 :                                    vac_plus, devab(:, :, 1), devab(:, :, 2), devab(:, :, 3))
     233       96000 :                DO i = 1, 3
     234             :                   CALL ab_contract(dvab(sgfa:m1, sgfb:m2, i), devab(:, :, i), sphi_a(:, sgfa:), &
     235       96000 :                                    sphi_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
     236             :                END DO
     237             : 
     238             :             ELSE
     239           0 :                ALLOCATE (vwork(ncoa, ncob, la_max(iset) + lb_max(jset) + 1), &
     240       15270 :                          vac(ncoa, ncob), vac_plus(ncoap, ncobp), devab(ncoa, ncob, 3))
     241     9661740 :                vwork = 0.0_dp
     242     1489646 :                vac = 0.0_dp
     243             :                CALL operator2(cps_operator2, la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
     244             :                               lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
     245        1018 :                               omega, r_cutoff, rab, rab2, vac, vwork)
     246             :             END IF
     247             : 
     248             :             CALL ab_contract(vab(sgfa:m1, sgfb:m2), vac(1:ncoa, 1:ncob), sphi_a(:, sgfa:), sphi_b(:, sgfb:), &
     249       25018 :                              ncoa, ncob, nsgfa(iset), nsgfb(jset))
     250       26728 :             DEALLOCATE (vwork, vac, vac_plus, devab)
     251             :          END DO
     252             :       END DO
     253             : 
     254         196 :       CALL timestop(handle)
     255             : 
     256         392 :    END SUBROUTINE int_operator_ab_os_low
     257             : 
     258             : ! **************************************************************************************************
     259             : !> \brief calculate overlap integrals (a,b)
     260             : !> \param sab integral (a,b)
     261             : !> \param dsab derivative of sab with respect to A
     262             : !> \param rab distance vector between center A and B
     263             : !> \param fba basis at center A
     264             : !> \param fbb basis at center B
     265             : !> \param calculate_forces ...
     266             : !> \param debug integrals are debugged by recursive routines if requested
     267             : !> \param dmax maximal deviation between integrals when debugging
     268             : ! **************************************************************************************************
     269        2296 :    SUBROUTINE int_overlap_ab_os(sab, dsab, rab, fba, fbb, calculate_forces, debug, dmax)
     270             : 
     271             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: sab
     272             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
     273             :          OPTIONAL                                        :: dsab
     274             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     275             :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     276             :       LOGICAL, INTENT(IN)                                :: calculate_forces, debug
     277             :       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
     278             : 
     279             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'int_overlap_ab_os'
     280             : 
     281             :       INTEGER                                            :: handle, i, iset, jset, m1, m2, maxco, &
     282             :                                                             maxcoa, maxcob, maxl, maxla, maxlb, &
     283             :                                                             ncoa, ncob, nseta, nsetb, sgfa, sgfb
     284        2296 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     285        2296 :                                                             npgfb, nsgfa, nsgfb
     286        2296 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     287             :       LOGICAL                                            :: failure
     288             :       REAL(KIND=dp)                                      :: dab, ra(3), rb(3)
     289             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: sint
     290             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dsint
     291        2296 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     292        2296 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, scon_a, scon_b, zeta, zetb
     293             : 
     294        2296 :       failure = .FALSE.
     295        2296 :       CALL timeset(routineN, handle)
     296        2296 :       NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
     297        2296 :                first_sgfa, first_sgfb, scon_a, scon_b, zeta, zetb)
     298             : 
     299             :       ! basis ikind
     300        2296 :       first_sgfa => fba%first_sgf
     301        2296 :       la_max => fba%lmax
     302        2296 :       la_min => fba%lmin
     303        2296 :       npgfa => fba%npgf
     304        2296 :       nseta = fba%nset
     305        2296 :       nsgfa => fba%nsgf_set
     306        2296 :       rpgfa => fba%pgf_radius
     307        2296 :       set_radius_a => fba%set_radius
     308        2296 :       scon_a => fba%scon
     309        2296 :       zeta => fba%zet
     310             :       ! basis jkind
     311        2296 :       first_sgfb => fbb%first_sgf
     312        2296 :       lb_max => fbb%lmax
     313        2296 :       lb_min => fbb%lmin
     314        2296 :       npgfb => fbb%npgf
     315        2296 :       nsetb = fbb%nset
     316        2296 :       nsgfb => fbb%nsgf_set
     317        2296 :       rpgfb => fbb%pgf_radius
     318        2296 :       set_radius_b => fbb%set_radius
     319        2296 :       scon_b => fbb%scon
     320        2296 :       zetb => fbb%zet
     321             : 
     322        2296 :       CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
     323        2296 :       CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
     324        2296 :       maxco = MAX(maxcoa, maxcob)
     325        2296 :       maxl = MAX(maxla, maxlb)
     326        9184 :       ALLOCATE (sint(maxco, maxco))
     327       11480 :       ALLOCATE (dsint(maxco, maxco, 3))
     328             : 
     329        9184 :       dab = SQRT(SUM(rab**2))
     330             : 
     331     4523933 :       sab = 0.0_dp
     332        2296 :       IF (calculate_forces) THEN
     333     6857394 :          IF (PRESENT(dsab)) dsab = 0.0_dp
     334             :       END IF
     335             : 
     336       11826 :       DO iset = 1, nseta
     337             : 
     338        9530 :          ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     339        9530 :          sgfa = first_sgfa(1, iset)
     340             : 
     341       75365 :          DO jset = 1, nsetb
     342             : 
     343       63539 :             IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     344             : 
     345       22981 :             ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     346       22981 :             sgfb = first_sgfb(1, jset)
     347       22981 :             m1 = sgfa + nsgfa(iset) - 1
     348       22981 :             m2 = sgfb + nsgfb(jset) - 1
     349             : 
     350       22981 :             IF (calculate_forces) THEN
     351             :                CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     352             :                                lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     353       11517 :                                rab, sint, dsint)
     354             :             ELSE
     355             :                CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     356             :                                lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     357       11464 :                                rab, sint)
     358             : 
     359             :             END IF
     360             : 
     361             :             ! debug if requested
     362       22981 :             IF (debug) THEN
     363          47 :                ra = 0.0_dp
     364          47 :                rb = rab
     365             :                CALL overlap_ab_test(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
     366             :                                     lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), &
     367          47 :                                     ra, rb, sint, dmax)
     368             :             END IF
     369             : 
     370             :             CALL ab_contract(sab(sgfa:m1, sgfb:m2), sint(1:ncoa, 1:ncob), scon_a(:, sgfa:), scon_b(:, sgfb:), &
     371       22981 :                              ncoa, ncob, nsgfa(iset), nsgfb(jset))
     372       32511 :             IF (calculate_forces) THEN
     373       46068 :                DO i = 1, 3
     374             :                   CALL ab_contract(dsab(sgfa:m1, sgfb:m2, i), dsint(1:ncoa, 1:ncob, i), scon_a(:, sgfa:), &
     375       46068 :                                    scon_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
     376             :                END DO
     377             :             END IF
     378             :          END DO
     379             :       END DO
     380             : 
     381        2296 :       DEALLOCATE (sint, dsint)
     382             : 
     383        2296 :       CALL timestop(handle)
     384             : 
     385        2296 :    END SUBROUTINE int_overlap_ab_os
     386             : 
     387             : ! **************************************************************************************************
     388             : !> \brief calculate integrals (a|(r-Ra)^(2m)|b)
     389             : !> \param sab integral (a|(r-Ra)^(2m)|b)
     390             : !> \param dsab derivative of sab with respect to A
     391             : !> \param rab distance vector between center A and B
     392             : !> \param fba fba basis at center A
     393             : !> \param fbb fbb basis at center B
     394             : !> \param m exponent m in operator (r-Ra)^(2m)
     395             : !> \param calculate_forces ...
     396             : ! **************************************************************************************************
     397          32 :    SUBROUTINE int_ra2m_ab_os(sab, dsab, rab, fba, fbb, m, calculate_forces)
     398             : 
     399             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: sab
     400             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT), &
     401             :          OPTIONAL                                        :: dsab
     402             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     403             :       TYPE(gto_basis_set_type), POINTER                  :: fba, fbb
     404             :       INTEGER, INTENT(IN)                                :: m
     405             :       LOGICAL, INTENT(IN)                                :: calculate_forces
     406             : 
     407             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'int_ra2m_ab_os'
     408             : 
     409             :       INTEGER                                            :: handle, i, iset, jset, m1, m2, maxco, &
     410             :                                                             maxcoa, maxcob, maxl, maxla, maxlb, &
     411             :                                                             ncoa, ncob, nseta, nsetb, sgfa, sgfb
     412          32 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     413          32 :                                                             npgfb, nsgfa, nsgfb
     414          32 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     415             :       LOGICAL                                            :: failure
     416             :       REAL(KIND=dp)                                      :: dab
     417             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: sint
     418             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dsint
     419          32 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: scon_a, scon_b, zeta, zetb
     420             : 
     421          32 :       failure = .FALSE.
     422          32 :       CALL timeset(routineN, handle)
     423          32 :       NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
     424          32 :                first_sgfa, first_sgfb, scon_a, scon_b, zeta, zetb)
     425             : 
     426             :       ! basis ikind
     427          32 :       first_sgfa => fba%first_sgf
     428          32 :       la_max => fba%lmax
     429          32 :       la_min => fba%lmin
     430          32 :       npgfa => fba%npgf
     431          32 :       nseta = fba%nset
     432          32 :       nsgfa => fba%nsgf_set
     433          32 :       scon_a => fba%scon
     434          32 :       zeta => fba%zet
     435             :       ! basis jkind
     436          32 :       first_sgfb => fbb%first_sgf
     437          32 :       lb_max => fbb%lmax
     438          32 :       lb_min => fbb%lmin
     439          32 :       npgfb => fbb%npgf
     440          32 :       nsetb = fbb%nset
     441          32 :       nsgfb => fbb%nsgf_set
     442          32 :       scon_b => fbb%scon
     443          32 :       zetb => fbb%zet
     444             : 
     445          32 :       CALL get_gto_basis_set(fba, maxco=maxcoa, maxl=maxla)
     446          32 :       CALL get_gto_basis_set(fbb, maxco=maxcob, maxl=maxlb)
     447          32 :       maxco = MAX(maxcoa, maxcob)
     448          32 :       maxl = MAX(maxla, maxlb)
     449         128 :       ALLOCATE (sint(maxco, maxco))
     450         160 :       ALLOCATE (dsint(maxco, maxco, 3))
     451             : 
     452         128 :       dab = SQRT(SUM(rab**2))
     453             : 
     454      963872 :       sab = 0.0_dp
     455          32 :       IF (calculate_forces) THEN
     456     2891680 :          IF (PRESENT(dsab)) dsab = 0.0_dp
     457             :       END IF
     458             : 
     459         352 :       DO iset = 1, nseta
     460             : 
     461         320 :          ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     462         320 :          sgfa = first_sgfa(1, iset)
     463             : 
     464        5152 :          DO jset = 1, nsetb
     465             : 
     466        4800 :             ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     467        4800 :             sgfb = first_sgfb(1, jset)
     468        4800 :             m1 = sgfa + nsgfa(iset) - 1
     469        4800 :             m2 = sgfb + nsgfb(jset) - 1
     470             : 
     471             :             CALL operator_ra2m(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
     472             :                                lb_max(jset), lb_min(jset), npgfb(jset), zetb(:, jset), &
     473        4800 :                                m, rab, sint, dsint, calculate_forces)
     474             : 
     475             :             CALL ab_contract(sab(sgfa:m1, sgfb:m2), sint, scon_a(:, sgfa:), scon_b(:, sgfb:), &
     476        4800 :                              ncoa, ncob, nsgfa(iset), nsgfb(jset))
     477        5120 :             IF (calculate_forces) THEN
     478       19200 :                DO i = 1, 3
     479             :                   CALL ab_contract(dsab(sgfa:m1, sgfb:m2, i), dsint(:, :, i), scon_a(:, sgfa:), &
     480       19200 :                                    scon_b(:, sgfb:), ncoa, ncob, nsgfa(iset), nsgfb(jset))
     481             :                END DO
     482             :             END IF
     483             :          END DO
     484             :       END DO
     485             : 
     486          32 :       DEALLOCATE (sint, dsint)
     487             : 
     488          32 :       CALL timestop(handle)
     489             : 
     490          32 :    END SUBROUTINE int_ra2m_ab_os
     491             : 
     492             : ! **************************************************************************************************
     493             : !> \brief calculate integrals (a,b,fa)
     494             : !> \param abaint integral (a,b,fa)
     495             : !> \param dabdaint derivative of abaint with respect to A
     496             : !> \param rab distance vector between center A and B
     497             : !> \param oba orbital basis at center A
     498             : !> \param obb orbital basis at center B
     499             : !> \param fba auxiliary basis set at center A
     500             : !> \param calculate_forces ...
     501             : !> \param debug integrals are debugged by recursive routines if requested
     502             : !> \param dmax maximal deviation between integrals when debugging
     503             : ! **************************************************************************************************
     504        1148 :    SUBROUTINE int_overlap_aba_os(abaint, dabdaint, rab, oba, obb, fba, &
     505             :                                  calculate_forces, debug, dmax)
     506             : 
     507             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: abaint
     508             :       REAL(KIND=dp), ALLOCATABLE, &
     509             :          DIMENSION(:, :, :, :), OPTIONAL                 :: dabdaint
     510             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     511             :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fba
     512             :       LOGICAL, INTENT(IN)                                :: calculate_forces, debug
     513             :       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
     514             : 
     515             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aba_os'
     516             : 
     517             :       INTEGER                                            :: handle, i, iset, jset, kaset, m1, m2, &
     518             :                                                             m3, ncoa, ncob, ncoc, nseta, nsetb, &
     519             :                                                             nsetca, sgfa, sgfb, sgfc
     520        1148 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, lca_max, &
     521        1148 :                                                             lca_min, npgfa, npgfb, npgfca, nsgfa, &
     522        1148 :                                                             nsgfb, nsgfca
     523        1148 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfca
     524             :       REAL(KIND=dp)                                      :: dab, dac, dbc
     525        1148 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: saba
     526        1148 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: sdaba
     527             :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rac, rb, rbc
     528        1148 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_ca
     529        1148 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, rpgfca, scon_a, scon_b, &
     530        1148 :                                                             scon_ca, zeta, zetb, zetca
     531             : 
     532        1148 :       CALL timeset(routineN, handle)
     533        1148 :       NULLIFY (la_max, la_min, lb_max, lb_min, lca_max, lca_min, npgfa, npgfb, &
     534        1148 :                npgfca, nsgfa, nsgfb, nsgfca)
     535        1148 :       NULLIFY (first_sgfa, first_sgfb, first_sgfca, set_radius_a, set_radius_b, &
     536        1148 :                set_radius_ca, rpgfa, rpgfb, rpgfca, scon_a, scon_b, scon_ca, &
     537        1148 :                zeta, zetb, zetca)
     538             : 
     539             :       ! basis ikind
     540        1148 :       first_sgfa => oba%first_sgf
     541        1148 :       la_max => oba%lmax
     542        1148 :       la_min => oba%lmin
     543        1148 :       npgfa => oba%npgf
     544        1148 :       nseta = oba%nset
     545        1148 :       nsgfa => oba%nsgf_set
     546        1148 :       rpgfa => oba%pgf_radius
     547        1148 :       set_radius_a => oba%set_radius
     548        1148 :       scon_a => oba%scon
     549        1148 :       zeta => oba%zet
     550             :       ! basis jkind
     551        1148 :       first_sgfb => obb%first_sgf
     552        1148 :       lb_max => obb%lmax
     553        1148 :       lb_min => obb%lmin
     554        1148 :       npgfb => obb%npgf
     555        1148 :       nsetb = obb%nset
     556        1148 :       nsgfb => obb%nsgf_set
     557        1148 :       rpgfb => obb%pgf_radius
     558        1148 :       set_radius_b => obb%set_radius
     559        1148 :       scon_b => obb%scon
     560        1148 :       zetb => obb%zet
     561             : 
     562             :       ! basis RI A
     563        1148 :       first_sgfca => fba%first_sgf
     564        1148 :       lca_max => fba%lmax
     565        1148 :       lca_min => fba%lmin
     566        1148 :       npgfca => fba%npgf
     567        1148 :       nsetca = fba%nset
     568        1148 :       nsgfca => fba%nsgf_set
     569        1148 :       rpgfca => fba%pgf_radius
     570        1148 :       set_radius_ca => fba%set_radius
     571        1148 :       scon_ca => fba%scon
     572        1148 :       zetca => fba%zet
     573             : 
     574        4592 :       dab = SQRT(SUM(rab**2))
     575             : 
     576     1044602 :       abaint = 0.0_dp
     577        1148 :       IF (calculate_forces) THEN
     578     1704876 :          IF (PRESENT(dabdaint)) dabdaint = 0.0_dp
     579             :       END IF
     580             : 
     581        2352 :       DO iset = 1, nseta
     582             : 
     583        1204 :          ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     584        1204 :          sgfa = first_sgfa(1, iset)
     585             : 
     586        3836 :          DO jset = 1, nsetb
     587             : 
     588        1484 :             IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     589             : 
     590        1484 :             ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     591        1484 :             sgfb = first_sgfb(1, jset)
     592        1484 :             m1 = sgfa + nsgfa(iset) - 1
     593        1484 :             m2 = sgfb + nsgfb(jset) - 1
     594             : 
     595             :             ! calculate integrals abaint and derivative [d(a,b,a)/dA] dabdaint if requested
     596             :             rac = 0._dp
     597        1484 :             dac = 0._dp
     598        5936 :             rbc = -rab
     599             :             dbc = dab
     600             : 
     601       14742 :             DO kaset = 1, nsetca
     602             : 
     603       12054 :                IF (set_radius_b(jset) + set_radius_ca(kaset) < dab) CYCLE
     604             : 
     605        7997 :                ncoc = npgfca(kaset)*(ncoset(lca_max(kaset)) - ncoset(lca_min(kaset) - 1))
     606        7997 :                sgfc = first_sgfca(1, kaset)
     607        7997 :                m3 = sgfc + nsgfca(kaset) - 1
     608        9481 :                IF (ncoa*ncob*ncoc > 0) THEN
     609       39985 :                   ALLOCATE (saba(ncoa, ncob, ncoc))
     610             :                   ! integrals
     611        7997 :                   IF (calculate_forces) THEN
     612       20382 :                      ALLOCATE (sdaba(ncoa, ncob, ncoc, 3))
     613             :                      CALL overlap_aab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     614             :                                       lca_max(kaset), lca_min(kaset), npgfca(kaset), rpgfca(:, kaset), zetca(:, kaset), &
     615             :                                       lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     616        3397 :                                       rab, saba=saba, daba=sdaba)
     617             : 
     618       13588 :                      DO i = 1, 3
     619             :                         CALL abc_contract(dabdaint(sgfa:m1, sgfb:m2, sgfc:m3, i), sdaba(1:ncoa, 1:ncob, 1:ncoc, i), &
     620             :                                           scon_a(:, sgfa:), scon_b(:, sgfb:), scon_ca(:, sgfc:), &
     621       13588 :                                           ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfca(kaset))
     622             :                      END DO
     623             : 
     624        3397 :                      DEALLOCATE (sdaba)
     625             :                   ELSE
     626             :                      CALL overlap_aab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     627             :                                       lca_max(kaset), lca_min(kaset), npgfca(kaset), rpgfca(:, kaset), zetca(:, kaset), &
     628             :                                       lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     629        4600 :                                       rab, saba=saba)
     630             : 
     631             :                   END IF
     632             :                   ! debug if requested
     633        7997 :                   IF (debug) THEN
     634           7 :                      ra = 0.0_dp
     635           7 :                      rb = rab
     636             :                      CALL overlap_abc_test(la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
     637             :                                            lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
     638             :                                            lca_max(kaset), npgfca(kaset), zetca(:, kaset), lca_min(kaset), &
     639           7 :                                            ra, rb, ra, saba, dmax)
     640             :                   END IF
     641             :                   CALL abc_contract(abaint(sgfa:m1, sgfb:m2, sgfc:m3), saba(1:ncoa, 1:ncob, 1:ncoc), &
     642             :                                     scon_a(:, sgfa:), scon_b(:, sgfb:), scon_ca(:, sgfc:), &
     643        7997 :                                     ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfca(kaset))
     644        7997 :                   DEALLOCATE (saba)
     645             :                END IF
     646             :             END DO
     647             :          END DO
     648             :       END DO
     649             : 
     650        1148 :       CALL timestop(handle)
     651             : 
     652        2296 :    END SUBROUTINE int_overlap_aba_os
     653             : 
     654             : ! **************************************************************************************************
     655             : !> \brief calculate integrals (a,b,fb)
     656             : !> \param abbint integral (a,b,fb)
     657             : !> \param dabbint derivative of abbint with respect to A
     658             : !> \param rab distance vector between center A and B
     659             : !> \param oba orbital basis at center A
     660             : !> \param obb orbital basis at center B
     661             : !> \param fbb auxiliary basis set at center B
     662             : !> \param calculate_forces ...
     663             : !> \param debug integrals are debugged by recursive routines if requested
     664             : !> \param dmax maximal deviation between integrals when debugging
     665             : ! **************************************************************************************************
     666        1148 :    SUBROUTINE int_overlap_abb_os(abbint, dabbint, rab, oba, obb, fbb, calculate_forces, debug, dmax)
     667             : 
     668             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: abbint
     669             :       REAL(KIND=dp), ALLOCATABLE, &
     670             :          DIMENSION(:, :, :, :), OPTIONAL                 :: dabbint
     671             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     672             :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb, fbb
     673             :       LOGICAL, INTENT(IN)                                :: calculate_forces, debug
     674             :       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
     675             : 
     676             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_abb_os'
     677             : 
     678             :       INTEGER                                            :: handle, i, iset, jset, kbset, m1, m2, &
     679             :                                                             m3, ncoa, ncob, ncoc, nseta, nsetb, &
     680             :                                                             nsetcb, sgfa, sgfb, sgfc
     681        1148 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, lcb_max, &
     682        1148 :                                                             lcb_min, npgfa, npgfb, npgfcb, nsgfa, &
     683        1148 :                                                             nsgfb, nsgfcb
     684        1148 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb, first_sgfcb
     685             :       REAL(KIND=dp)                                      :: dab, dac, dbc
     686        1148 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: sabb
     687        1148 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: sdabb
     688             :       REAL(KIND=dp), DIMENSION(3)                        :: ra, rac, rb, rbc
     689        1148 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b, set_radius_cb
     690        1148 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, rpgfcb, scon_a, scon_b, &
     691        1148 :                                                             scon_cb, zeta, zetb, zetcb
     692             : 
     693        1148 :       CALL timeset(routineN, handle)
     694        1148 :       NULLIFY (la_max, la_min, lb_max, lb_min, lcb_max, lcb_min, npgfa, npgfb, &
     695        1148 :                npgfcb, nsgfa, nsgfb, nsgfcb)
     696        1148 :       NULLIFY (first_sgfa, first_sgfb, first_sgfcb, set_radius_a, set_radius_b, &
     697        1148 :                set_radius_cb, rpgfa, rpgfb, rpgfcb, scon_a, scon_b, scon_cb, &
     698        1148 :                zeta, zetb, zetcb)
     699             : 
     700             :       ! basis ikind
     701        1148 :       first_sgfa => oba%first_sgf
     702        1148 :       la_max => oba%lmax
     703        1148 :       la_min => oba%lmin
     704        1148 :       npgfa => oba%npgf
     705        1148 :       nseta = oba%nset
     706        1148 :       nsgfa => oba%nsgf_set
     707        1148 :       rpgfa => oba%pgf_radius
     708        1148 :       set_radius_a => oba%set_radius
     709        1148 :       scon_a => oba%scon
     710        1148 :       zeta => oba%zet
     711             :       ! basis jkind
     712        1148 :       first_sgfb => obb%first_sgf
     713        1148 :       lb_max => obb%lmax
     714        1148 :       lb_min => obb%lmin
     715        1148 :       npgfb => obb%npgf
     716        1148 :       nsetb = obb%nset
     717        1148 :       nsgfb => obb%nsgf_set
     718        1148 :       rpgfb => obb%pgf_radius
     719        1148 :       set_radius_b => obb%set_radius
     720        1148 :       scon_b => obb%scon
     721        1148 :       zetb => obb%zet
     722             : 
     723             :       ! basis RI B
     724        1148 :       first_sgfcb => fbb%first_sgf
     725        1148 :       lcb_max => fbb%lmax
     726        1148 :       lcb_min => fbb%lmin
     727        1148 :       npgfcb => fbb%npgf
     728        1148 :       nsetcb = fbb%nset
     729        1148 :       nsgfcb => fbb%nsgf_set
     730        1148 :       rpgfcb => fbb%pgf_radius
     731        1148 :       set_radius_cb => fbb%set_radius
     732        1148 :       scon_cb => fbb%scon
     733        1148 :       zetcb => fbb%zet
     734             : 
     735        4592 :       dab = SQRT(SUM(rab**2))
     736             : 
     737     1048274 :       abbint = 0.0_dp
     738        1148 :       IF (calculate_forces) THEN
     739     1704876 :          IF (PRESENT(dabbint)) dabbint = 0.0_dp
     740             :       END IF
     741             : 
     742        2352 :       DO iset = 1, nseta
     743             : 
     744        1204 :          ncoa = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     745        1204 :          sgfa = first_sgfa(1, iset)
     746             : 
     747        3836 :          DO jset = 1, nsetb
     748             : 
     749        1484 :             IF (set_radius_a(iset) + set_radius_b(jset) < dab) CYCLE
     750             : 
     751        1484 :             ncob = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     752        1484 :             sgfb = first_sgfb(1, jset)
     753        1484 :             m1 = sgfa + nsgfa(iset) - 1
     754        1484 :             m2 = sgfb + nsgfb(jset) - 1
     755             : 
     756             :             ! calculate integrals abbint and derivative [d(a,b,b)/dA] dabbint if requested
     757        1484 :             rac = rab
     758        1484 :             dac = dab
     759             :             rbc = 0._dp
     760        1484 :             dbc = 0._dp
     761             : 
     762       14742 :             DO kbset = 1, nsetcb
     763             : 
     764       12054 :                IF (set_radius_a(iset) + set_radius_cb(kbset) < dab) CYCLE
     765             : 
     766        7997 :                ncoc = npgfcb(kbset)*(ncoset(lcb_max(kbset)) - ncoset(lcb_min(kbset) - 1))
     767        7997 :                sgfc = first_sgfcb(1, kbset)
     768        7997 :                m3 = sgfc + nsgfcb(kbset) - 1
     769        9481 :                IF (ncoa*ncob*ncoc > 0) THEN
     770       39985 :                   ALLOCATE (sabb(ncoa, ncob, ncoc))
     771        7997 :                   IF (calculate_forces) THEN
     772       20382 :                      ALLOCATE (sdabb(ncoa, ncob, ncoc, 3))
     773             :                      CALL overlap_abb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     774             :                                       lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     775             :                                       lcb_max(kbset), lcb_min(kbset), npgfcb(kbset), rpgfcb(:, kbset), zetcb(:, kbset), &
     776        3397 :                                       rab, sabb, sdabb)
     777             : 
     778       13588 :                      DO i = 1, 3
     779             :                         CALL abc_contract(dabbint(sgfa:m1, sgfb:m2, sgfc:m3, i), sdabb(1:ncoa, 1:ncob, 1:ncoc, i), &
     780             :                                           scon_a(:, sgfa:), scon_b(:, sgfb:), scon_cb(:, sgfc:), &
     781       13588 :                                           ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfcb(kbset))
     782             :                      END DO
     783        3397 :                      DEALLOCATE (sdabb)
     784             :                   ELSE
     785             :                      CALL overlap_abb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     786             :                                       lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     787             :                                       lcb_max(kbset), lcb_min(kbset), npgfcb(kbset), rpgfcb(:, kbset), zetcb(:, kbset), &
     788        4600 :                                       rab, sabb)
     789             :                   END IF
     790             :                   ! debug if requested
     791        7997 :                   IF (debug) THEN
     792           7 :                      ra = 0.0_dp
     793           7 :                      rb = rab
     794             :                      CALL overlap_abc_test(la_max(iset), npgfa(iset), zeta(:, iset), la_min(iset), &
     795             :                                            lb_max(jset), npgfb(jset), zetb(:, jset), lb_min(jset), &
     796             :                                            lcb_max(kbset), npgfcb(kbset), zetcb(:, kbset), lcb_min(kbset), &
     797           7 :                                            ra, rb, rb, sabb, dmax)
     798             :                   END IF
     799             : 
     800             :                   CALL abc_contract(abbint(sgfa:m1, sgfb:m2, sgfc:m3), sabb(1:ncoa, 1:ncob, 1:ncoc), &
     801             :                                     scon_a(:, sgfa:), scon_b(:, sgfb:), scon_cb(:, sgfc:), &
     802        7997 :                                     ncoa, ncob, ncoc, nsgfa(iset), nsgfb(jset), nsgfcb(kbset))
     803        7997 :                   DEALLOCATE (sabb)
     804             :                END IF
     805             :             END DO
     806             : 
     807             :          END DO
     808             :       END DO
     809             : 
     810        1148 :       CALL timestop(handle)
     811             : 
     812        2296 :    END SUBROUTINE int_overlap_abb_os
     813             : 
     814             : ! **************************************************************************************************
     815             : !> \brief calculate overlap integrals (aa,bb)
     816             : !> \param saabb integral (aa,bb)
     817             : !> \param oba orbital basis at center A
     818             : !> \param obb orbital basis at center B
     819             : !> \param rab ...
     820             : !> \param debug integrals are debugged by recursive routines if requested
     821             : !> \param dmax maximal deviation between integrals when debugging
     822             : ! **************************************************************************************************
     823           9 :    SUBROUTINE int_overlap_aabb_os(saabb, oba, obb, rab, debug, dmax)
     824             : 
     825             :       REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER      :: saabb
     826             :       TYPE(gto_basis_set_type), POINTER                  :: oba, obb
     827             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: rab
     828             :       LOGICAL, INTENT(IN)                                :: debug
     829             :       REAL(KIND=dp), INTENT(INOUT)                       :: dmax
     830             : 
     831             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'int_overlap_aabb_os'
     832             : 
     833             :       INTEGER :: handle, iset, isgfa1, jset, jsgfa2, kset, ksgfb1, lds, lset, lsgfb2, m1, m2, m3, &
     834             :          m4, maxco, maxcoa, maxcob, maxl, maxla, maxlb, ncoa1, ncoa2, ncob1, ncob2, nseta, nsetb, &
     835             :          sgfa1, sgfa2, sgfb1, sgfb2
     836           9 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     837           9 :                                                             npgfb, nsgfa, nsgfb
     838           9 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     839             :       LOGICAL                                            :: asets_equal, bsets_equal
     840             :       REAL(KIND=dp)                                      :: dab, ra(3), rb(3)
     841           9 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: swork
     842             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: sint
     843           9 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     844           9 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
     845             : 
     846           9 :       CALL timeset(routineN, handle)
     847           9 :       NULLIFY (la_max, la_min, lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, &
     848           9 :                first_sgfa, first_sgfb, set_radius_a, set_radius_b, rpgfa, rpgfb, &
     849           9 :                sphi_a, sphi_b, zeta, zetb)
     850             : 
     851             :       ! basis ikind
     852           9 :       first_sgfa => oba%first_sgf
     853           9 :       la_max => oba%lmax
     854           9 :       la_min => oba%lmin
     855           9 :       npgfa => oba%npgf
     856           9 :       nseta = oba%nset
     857           9 :       nsgfa => oba%nsgf_set
     858           9 :       rpgfa => oba%pgf_radius
     859           9 :       set_radius_a => oba%set_radius
     860           9 :       sphi_a => oba%sphi
     861           9 :       zeta => oba%zet
     862             :       ! basis jkind
     863           9 :       first_sgfb => obb%first_sgf
     864           9 :       lb_max => obb%lmax
     865           9 :       lb_min => obb%lmin
     866           9 :       npgfb => obb%npgf
     867           9 :       nsetb = obb%nset
     868           9 :       nsgfb => obb%nsgf_set
     869           9 :       rpgfb => obb%pgf_radius
     870           9 :       set_radius_b => obb%set_radius
     871           9 :       sphi_b => obb%sphi
     872           9 :       zetb => obb%zet
     873             : 
     874           9 :       CALL get_gto_basis_set(oba, maxco=maxcoa, maxl=maxla)
     875           9 :       CALL get_gto_basis_set(obb, maxco=maxcob, maxl=maxlb)
     876           9 :       maxco = MAX(maxcoa, maxcob)
     877           9 :       maxla = 2*maxla
     878           9 :       maxlb = 2*maxlb
     879           9 :       maxl = MAX(maxla, maxlb)
     880           9 :       lds = ncoset(maxl)
     881          54 :       ALLOCATE (sint(maxco, maxco, maxco, maxco))
     882          36 :       ALLOCATE (swork(lds, lds))
     883     5736789 :       sint = 0._dp
     884         999 :       swork = 0._dp
     885             : 
     886          36 :       dab = SQRT(SUM(rab**2))
     887             : 
     888          18 :       DO iset = 1, nseta
     889             : 
     890           9 :          ncoa1 = npgfa(iset)*ncoset(la_max(iset))
     891           9 :          sgfa1 = first_sgfa(1, iset)
     892           9 :          m1 = sgfa1 + nsgfa(iset) - 1
     893             : 
     894          27 :          DO jset = iset, nseta
     895             : 
     896           9 :             ncoa2 = npgfa(jset)*ncoset(la_max(jset))
     897           9 :             sgfa2 = first_sgfa(1, jset)
     898           9 :             m2 = sgfa2 + nsgfa(jset) - 1
     899             : 
     900          27 :             DO kset = 1, nsetb
     901             : 
     902           9 :                ncob1 = npgfb(kset)*ncoset(lb_max(kset))
     903           9 :                sgfb1 = first_sgfb(1, kset)
     904           9 :                m3 = sgfb1 + nsgfb(kset) - 1
     905             : 
     906          27 :                DO lset = kset, nsetb
     907             : 
     908           9 :                   ncob2 = npgfb(lset)*ncoset(lb_max(lset))
     909           9 :                   sgfb2 = first_sgfb(1, lset)
     910           9 :                   m4 = sgfb2 + nsgfb(lset) - 1
     911             : 
     912             :                   ! check if sets are identical to spare some integral evaluation
     913           9 :                   asets_equal = .FALSE.
     914           9 :                   IF (iset == jset) asets_equal = .TRUE.
     915           9 :                   bsets_equal = .FALSE.
     916           9 :                   IF (kset == lset) bsets_equal = .TRUE.
     917             :                   ! calculate integrals
     918             :                   CALL overlap_aabb(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     919             :                                     la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), &
     920             :                                     lb_max(kset), lb_min(kset), npgfb(kset), rpgfb(:, kset), zetb(:, kset), &
     921             :                                     lb_max(lset), lb_min(lset), npgfb(lset), rpgfb(:, lset), zetb(:, lset), &
     922           9 :                                     asets_equal, bsets_equal, rab, dab, sint, swork, lds)
     923             :                   ! debug if requested
     924           9 :                   IF (debug) THEN
     925           3 :                      ra = 0.0_dp
     926           3 :                      rb = rab
     927             :                      CALL overlap_aabb_test(la_max(iset), la_min(iset), npgfa(iset), zeta(:, iset), &
     928             :                                             la_max(jset), la_min(jset), npgfa(jset), zeta(:, jset), &
     929             :                                             lb_max(kset), lb_min(kset), npgfb(kset), zetb(:, kset), &
     930             :                                             lb_max(lset), lb_min(lset), npgfb(lset), zetb(:, lset), &
     931           3 :                                             ra, rb, sint, dmax)
     932             :                   END IF
     933             : 
     934             :                   CALL abcd_contract(saabb(sgfa1:m1, sgfa2:m2, sgfb1:m3, sgfb2:m4), sint, sphi_a(:, sgfa1:), &
     935             :                                      sphi_a(:, sgfa2:), sphi_b(:, sgfb1:), sphi_b(:, sgfb2:), ncoa1, ncoa2, &
     936           9 :                                      ncob1, ncob2, nsgfa(iset), nsgfa(jset), nsgfb(kset), nsgfb(lset))
     937             : 
     938             :                   ! account for the fact that some integrals are alike
     939          54 :                   DO isgfa1 = sgfa1, m1
     940         189 :                      DO jsgfa2 = sgfa2, m2
     941         756 :                         DO ksgfb1 = sgfb1, m3
     942        3024 :                            DO lsgfb2 = sgfb2, m4
     943        2304 :                               saabb(jsgfa2, isgfa1, ksgfb1, lsgfb2) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2)
     944        2304 :                               saabb(isgfa1, jsgfa2, lsgfb2, ksgfb1) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2)
     945        2880 :                               saabb(jsgfa2, isgfa1, lsgfb2, ksgfb1) = saabb(isgfa1, jsgfa2, ksgfb1, lsgfb2)
     946             :                            END DO
     947             :                         END DO
     948             :                      END DO
     949             :                   END DO
     950             : 
     951             :                END DO
     952             :             END DO
     953             :          END DO
     954             :       END DO
     955             : 
     956           9 :       DEALLOCATE (sint, swork)
     957             : 
     958           9 :       CALL timestop(handle)
     959             : 
     960           9 :    END SUBROUTINE int_overlap_aabb_os
     961             : 
     962             : END MODULE generic_os_integrals

Generated by: LCOV version 1.15