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 Contraction of integrals over primitive Cartesian Gaussians based on the contraction
      10             : !>        matrix sphi which is part of the gto_basis_set_type
      11             : !> \par History
      12             : !>      -added libxsmm_abc_contract routine, A. Bussy (04.2020)
      13             : !> \author Dorothea Golze (05.2016)
      14             : ! **************************************************************************************************
      15             : MODULE ai_contraction_sphi
      16             : 
      17             :    USE kinds, ONLY: dp
      18             : #if(__LIBXSMM)
      19             :    USE libxsmm, ONLY: LIBXSMM_PREFETCH_NONE, &
      20             :                       libxsmm_blasint_kind, &
      21             :                       libxsmm_dgemm, &
      22             :                       libxsmm_dispatch, &
      23             :                       libxsmm_available, &
      24             :                       libxsmm_dmmcall, &
      25             :                       libxsmm_dmmfunction
      26             : #endif
      27             : 
      28             : #include "../base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             : 
      32             :    PRIVATE
      33             : 
      34             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_contraction_sphi'
      35             : 
      36             :    PUBLIC :: ab_contract, abc_contract, abcd_contract, libxsmm_abc_contract
      37             : 
      38             : CONTAINS
      39             : 
      40             : ! **************************************************************************************************
      41             : !> \brief contract overlap integrals (a,b) and transfer to spherical Gaussians
      42             : !> \param abint contracted, normalized integrals of spherical Gaussians
      43             : !> \param sab uncontracted, unnormalized integrals of primitive Cartesian Gaussians
      44             : !> \param sphi_a contraction matrix for center a
      45             : !> \param sphi_b contraction matrix for center b
      46             : !> \param ncoa number of cartesian orbitals on a
      47             : !> \param ncob number of cartesian orbitals on b
      48             : !> \param nsgfa number of spherical Gaussian functions on a
      49             : !> \param nsgfb number of spherical Gaussian functions on b
      50             : ! **************************************************************************************************
      51     1108678 :    SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
      52             : 
      53             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT)      :: abint
      54             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: sab, sphi_a, sphi_b
      55             :       INTEGER, INTENT(IN)                                :: ncoa, ncob, nsgfa, nsgfb
      56             : 
      57             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'ab_contract'
      58             : 
      59             :       INTEGER                                            :: handle, m1, m2, msphia, msphib, nn
      60     1108678 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: cpp
      61             : 
      62     1108678 :       CALL timeset(routineN, handle)
      63             : 
      64     1108678 :       msphia = SIZE(sphi_a, 1)
      65     1108678 :       msphib = SIZE(sphi_b, 1)
      66             : 
      67     1108678 :       m1 = SIZE(sab, 1)
      68     1108678 :       m2 = SIZE(sab, 2)
      69             : 
      70     1108678 :       nn = SIZE(abint, 1)
      71             : 
      72     4434712 :       ALLOCATE (cpp(nsgfa, m2))
      73             : 
      74    29894247 :       CALL dgemm("T", "N", nsgfa, m2, ncoa, 1._dp, sphi_a, msphia, sab, m1, 0.0_dp, cpp, nsgfa)
      75             :       CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp, nsgfa, sphi_b, msphib, 0.0_dp, &
      76    66497896 :                  abint, nn)
      77             : 
      78     1108678 :       DEALLOCATE (cpp)
      79             : 
      80     1108678 :       CALL timestop(handle)
      81             : 
      82     1108678 :    END SUBROUTINE ab_contract
      83             : 
      84             : ! **************************************************************************************************
      85             : !> \brief contract three-center overlap integrals (a,b,c) and transfer
      86             : !>        to spherical Gaussians
      87             : !> \param abcint contracted, normalized integrals of spherical Gaussians
      88             : !> \param sabc uncontracted, unnormalized integrals of primitive Cartesian Gaussians
      89             : !> \param sphi_a contraction matrix for center a
      90             : !> \param sphi_b contraction matrix for center b
      91             : !> \param sphi_c contraction matrix for center c
      92             : !> \param ncoa number of cartesian orbitals on a
      93             : !> \param ncob number of cartesian orbitals on b
      94             : !> \param ncoc number of cartesian orbitals on c
      95             : !> \param nsgfa number of spherical Gaussian functions on a
      96             : !> \param nsgfb number of spherical Gaussian functions on b
      97             : !> \param nsgfc number of spherical Gaussian functions on c
      98             : ! **************************************************************************************************
      99       96316 :    SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
     100             :                            nsgfa, nsgfb, nsgfc)
     101             : 
     102             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: abcint, sabc
     103             :       REAL(KIND=dp), DIMENSION(:, :)                     :: sphi_a, sphi_b, sphi_c
     104             :       INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
     105             : 
     106             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'abc_contract'
     107             : 
     108             :       INTEGER                                            :: handle, i, m1, m2, m3, msphia, msphib, &
     109             :                                                             msphic
     110       96316 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: cpc, cpp
     111             : 
     112       96316 :       CALL timeset(routineN, handle)
     113             : 
     114       96316 :       CPASSERT(SIZE(abcint, 1) == nsgfa)
     115       96316 :       CPASSERT(SIZE(abcint, 2) == nsgfb)
     116             : 
     117       96316 :       msphia = SIZE(sphi_a, 1)
     118       96316 :       msphib = SIZE(sphi_b, 1)
     119       96316 :       msphic = SIZE(sphi_c, 1)
     120             : 
     121       96316 :       m1 = SIZE(sabc, 1)
     122       96316 :       m2 = SIZE(sabc, 2)
     123       96316 :       m3 = SIZE(sabc, 3)
     124             : 
     125      866844 :       ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
     126             : 
     127       96316 :       CALL dgemm("T", "N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, cpp, nsgfa)
     128     1076476 :       DO i = 1, m3
     129             :          CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, cpp(:, :, i), nsgfa, sphi_b, msphib, &
     130     1076476 :                     0.0_dp, cpc(:, :, i), nsgfa)
     131             :       END DO
     132             :       CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, 1._dp, cpc, nsgfa*nsgfb, sphi_c, msphic, 0.0_dp, &
     133     1125848 :                  abcint, nsgfa*nsgfb)
     134             : 
     135       96316 :       DEALLOCATE (cpp, cpc)
     136             : 
     137       96316 :       CALL timestop(handle)
     138             : 
     139       96316 :    END SUBROUTINE abc_contract
     140             : 
     141             : ! **************************************************************************************************
     142             : !> \brief contract four-center overlap integrals (a,b,c,d) and transfer
     143             : !>        to spherical Gaussians
     144             : !> \param abcdint contracted, normalized integrals of spherical Gaussians
     145             : !> \param sabcd uncontracted, unnormalized integrals of primitive Cartesian Gaussians
     146             : !> \param sphi_a contraction matrix for center a
     147             : !> \param sphi_b contraction matrix for center b
     148             : !> \param sphi_c contraction matrix for center c
     149             : !> \param sphi_d contraction matrix for center d
     150             : !> \param ncoa number of cartesian orbitals on a
     151             : !> \param ncob number of cartesian orbitals on b
     152             : !> \param ncoc number of cartesian orbitals on c
     153             : !> \param ncod number of cartesian orbitals on d
     154             : !> \param nsgfa number of spherical Gaussian functions on a
     155             : !> \param nsgfb number of spherical Gaussian functions on b
     156             : !> \param nsgfc number of spherical Gaussian functions on c
     157             : !> \param nsgfd number of spherical Gaussian functions on d
     158             : ! **************************************************************************************************
     159           9 :    SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, &
     160             :                             ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
     161             : 
     162             :       REAL(KIND=dp), DIMENSION(:, :, :, :), &
     163             :          INTENT(INOUT)                                   :: abcdint
     164             :       REAL(KIND=dp), DIMENSION(:, :, :, :), INTENT(IN)   :: sabcd
     165             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: sphi_a, sphi_b, sphi_c, sphi_d
     166             :       INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, &
     167             :                                                             nsgfc, nsgfd
     168             : 
     169             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'abcd_contract'
     170             : 
     171             :       INTEGER                                            :: handle, isgfc, isgfd, m1, m2, m3, m4, &
     172             :                                                             msphia, msphib, msphic, msphid
     173           9 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: temp_cccc, work_cpcc
     174           9 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: temp_cpcc, work_cppc
     175           9 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: cpcc, cppc, cppp
     176             : 
     177           9 :       CALL timeset(routineN, handle)
     178             : 
     179           9 :       msphia = SIZE(sphi_a, 1)
     180           9 :       msphib = SIZE(sphi_b, 1)
     181           9 :       msphic = SIZE(sphi_c, 1)
     182           9 :       msphid = SIZE(sphi_d, 1)
     183             : 
     184           9 :       m1 = SIZE(sabcd, 1)
     185           9 :       m2 = SIZE(sabcd, 2)
     186           9 :       m3 = SIZE(sabcd, 3)
     187           9 :       m4 = SIZE(sabcd, 4)
     188             : 
     189             :       ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
     190         144 :                 cpcc(nsgfa, m2, nsgfc, nsgfd))
     191             : 
     192          81 :       ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
     193       35541 :       work_cppc = 0._dp
     194        5085 :       temp_cpcc = 0._dp
     195             : 
     196          63 :       ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
     197        1269 :       work_cpcc = 0._dp
     198         189 :       temp_cccc = 0._dp
     199             : 
     200             :       CALL dgemm("T", "N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
     201           9 :                  0.0_dp, cppp, nsgfa)
     202             :       CALL dgemm("N", "N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
     203           9 :                  sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
     204             : 
     205          45 :       DO isgfd = 1, nsgfd
     206      142164 :          work_cppc(:, :, :) = cppc(:, :, :, isgfd)
     207             :          CALL dgemm("N", "N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
     208          36 :                     sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
     209       20340 :          cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
     210         189 :          DO isgfc = 1, nsgfc
     211       20304 :             work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
     212             :             CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
     213         144 :                        msphib, 0.0_dp, temp_cccc, nsgfa)
     214        3060 :             abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
     215             :          END DO
     216             :       END DO
     217             : 
     218           9 :       DEALLOCATE (cpcc, cppc, cppp)
     219           9 :       DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
     220             : 
     221           9 :       CALL timestop(handle)
     222             : 
     223           9 :    END SUBROUTINE abcd_contract
     224             : 
     225             : ! **************************************************************************************************
     226             : !> \brief 3-center contraction routine from primitive cartesain Gaussians to spherical Gaussian
     227             : !>        functions. Exploits LIBXSMM for performance, falls back to BLAS if LIBXSMM not available.
     228             : !>        Requires pre-allocation of work buffers and pre-transposition of the sphi_a array. Requires
     229             : !>        the LIBXSMM library to be initialized somewhere before this routine is called.
     230             : !> \param abcint contracted integrals
     231             : !> \param sabc uncontracted integrals
     232             : !> \param tsphi_a assumed to have dimensions nsgfa x ncoa
     233             : !> \param sphi_b assumed to have dimensions ncob x nsgfb
     234             : !> \param sphi_c assumed to have dimensions ncoc x nsgfc
     235             : !> \param ncoa ...
     236             : !> \param ncob ...
     237             : !> \param ncoc ...
     238             : !> \param nsgfa ...
     239             : !> \param nsgfb ...
     240             : !> \param nsgfc ...
     241             : !> \param cpp_buffer ...
     242             : !> \param ccp_buffer ...
     243             : !> \note tested from version 1.9.0 of libxsmm
     244             : ! **************************************************************************************************
     245     2480509 :    SUBROUTINE libxsmm_abc_contract(abcint, sabc, tsphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
     246     2480509 :                                    nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer)
     247             : 
     248             :       REAL(dp), DIMENSION(*)                             :: abcint, sabc
     249             :       REAL(dp), DIMENSION(:, :)                          :: tsphi_a, sphi_b, sphi_c
     250             :       INTEGER, INTENT(IN)                                :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
     251             :       REAL(dp), DIMENSION(nsgfa*ncob)                    :: cpp_buffer
     252             :       REAL(dp), DIMENSION(nsgfa*nsgfb*ncoc)              :: ccp_buffer
     253             : 
     254             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'libxsmm_abc_contract'
     255             : 
     256             :       INTEGER                                            :: handle, i
     257             :       LOGICAL                                            :: libxsmm_kernels_available
     258             : #if(__LIBXSMM)
     259             :       INTEGER(libxsmm_blasint_kind)                      :: m, n, k
     260             :       TYPE(libxsmm_dmmfunction)                          :: xmm1, xmm2
     261             : #endif
     262             : 
     263     2480509 :       CALL timeset(routineN, handle)
     264     2480509 :       libxsmm_kernels_available = .FALSE.
     265             : 
     266             : #if(__LIBXSMM)
     267             :       !We make use of libxsmm code dispatching feature and call the same kernel multiple times
     268             :       !We loop over the last index of the matrix and call libxsmm each time
     269             : 
     270             :       !pre-fetch the kernels
     271     2480509 :       m = nsgfa; n = ncob; k = ncoa
     272     2480509 :       CALL libxsmm_dispatch(xmm1, m, n, k, beta=0.0_dp, prefetch=LIBXSMM_PREFETCH_NONE)
     273     2480509 :       m = nsgfa; n = nsgfb; k = ncob
     274     2480509 :       CALL libxsmm_dispatch(xmm2, m, n, k, beta=0.0_dp, prefetch=LIBXSMM_PREFETCH_NONE)
     275             : 
     276     2480509 :       libxsmm_kernels_available = libxsmm_available(xmm1) .AND. libxsmm_available(xmm2)
     277             : 
     278             :       IF (libxsmm_kernels_available) THEN
     279             :          ! contractions over a and b
     280    36337251 :          DO i = 1, ncoc
     281    33856742 :             CALL libxsmm_dmmcall(xmm1, tsphi_a, sabc((i - 1)*ncoa*ncob + 1), cpp_buffer)
     282    36337251 :             CALL libxsmm_dmmcall(xmm2, cpp_buffer, sphi_b, ccp_buffer((i - 1)*nsgfa*nsgfb + 1))
     283             :          END DO
     284             :       ELSE
     285           0 :          CPWARN("libxsmm available, but kernels are not, fallback to dgemm")
     286             :       END IF
     287             : #endif
     288             : 
     289     2480509 :       IF (.NOT. libxsmm_kernels_available) THEN
     290             :          ! we follow the same flow as for libxsmm above, but use BLAS dgemm
     291             :          ! contractions over a and b
     292           0 :          DO i = 1, ncoc
     293             :             CALL dgemm("N", "N", nsgfa, ncob, ncoa, 1.0_dp, tsphi_a, nsgfa, sabc((i - 1)*ncoa*ncob + 1), &
     294           0 :                        ncoa, 0.0_dp, cpp_buffer, nsgfa)
     295             :             CALL dgemm("N", "N", nsgfa, nsgfb, ncob, 1.0_dp, cpp_buffer, nsgfa, sphi_b, ncob, 0.0_dp, &
     296           0 :                        ccp_buffer((i - 1)*nsgfa*nsgfb + 1), nsgfa)
     297             :          END DO
     298             :       END IF
     299             : 
     300             :       ! last contraction, over c, as a larger MM
     301             :       CALL dgemm("N", "N", nsgfa*nsgfb, nsgfc, ncoc, 1.0_dp, ccp_buffer, nsgfa*nsgfb, sphi_c, ncoc, &
     302     2480509 :                  0.0_dp, abcint, nsgfa*nsgfb)
     303             : 
     304     2480509 :       CALL timestop(handle)
     305             : 
     306     2480509 :    END SUBROUTINE libxsmm_abc_contract
     307             : 
     308             : END MODULE ai_contraction_sphi