LCOV - code coverage report
Current view: top level - src - qs_tensors.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 94.9 % 1150 1091
Test Date: 2025-07-25 12:55:17 Functions: 90.5 % 21 19

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Utility methods to build 3-center integral tensors of various types.
      10              : ! **************************************************************************************************
      11              : MODULE qs_tensors
      12              :    USE OMP_LIB,                         ONLY: omp_get_num_threads,&
      13              :                                               omp_get_thread_num
      14              :    USE ai_contraction,                  ONLY: block_add
      15              :    USE ai_contraction_sphi,             ONLY: ab_contract,&
      16              :                                               abc_contract_xsmm
      17              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      18              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      19              :                                               gto_basis_set_p_type,&
      20              :                                               gto_basis_set_type
      21              :    USE block_p_types,                   ONLY: block_p_type
      22              :    USE cell_types,                      ONLY: cell_type,&
      23              :                                               real_to_scaled
      24              :    USE cp_array_utils,                  ONLY: cp_2d_r_p_type
      25              :    USE cp_control_types,                ONLY: dft_control_type
      26              :    USE cp_dbcsr_api,                    ONLY: dbcsr_filter,&
      27              :                                               dbcsr_finalize,&
      28              :                                               dbcsr_get_block_p,&
      29              :                                               dbcsr_get_matrix_type,&
      30              :                                               dbcsr_has_symmetry,&
      31              :                                               dbcsr_type,&
      32              :                                               dbcsr_type_antisymmetric,&
      33              :                                               dbcsr_type_no_symmetry
      34              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      35              :    USE cp_files,                        ONLY: close_file,&
      36              :                                               open_file
      37              :    USE dbt_api,                         ONLY: &
      38              :         dbt_blk_sizes, dbt_clear, dbt_copy, dbt_create, dbt_destroy, dbt_filter, dbt_get_block, &
      39              :         dbt_get_info, dbt_get_num_blocks, dbt_get_nze_total, dbt_iterator_next_block, &
      40              :         dbt_iterator_num_blocks, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
      41              :         dbt_ndims, dbt_put_block, dbt_reserve_blocks, dbt_type
      42              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      43              :    USE distribution_2d_types,           ONLY: distribution_2d_type
      44              :    USE gamma,                           ONLY: init_md_ftable
      45              :    USE hfx_compression_methods,         ONLY: hfx_add_mult_cache_elements,&
      46              :                                               hfx_add_single_cache_element,&
      47              :                                               hfx_decompress_first_cache,&
      48              :                                               hfx_flush_last_cache,&
      49              :                                               hfx_get_mult_cache_elements,&
      50              :                                               hfx_get_single_cache_element,&
      51              :                                               hfx_reset_cache_and_container
      52              :    USE hfx_types,                       ONLY: alloc_containers,&
      53              :                                               dealloc_containers,&
      54              :                                               hfx_cache_type,&
      55              :                                               hfx_compression_type,&
      56              :                                               hfx_container_type,&
      57              :                                               hfx_init_container
      58              :    USE input_constants,                 ONLY: do_potential_coulomb,&
      59              :                                               do_potential_id,&
      60              :                                               do_potential_mix_cl_trunc,&
      61              :                                               do_potential_short,&
      62              :                                               do_potential_truncated
      63              :    USE input_section_types,             ONLY: section_vals_val_get
      64              :    USE kinds,                           ONLY: dp,&
      65              :                                               int_8
      66              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      67              :                                               kpoint_type
      68              :    USE libint_2c_3c,                    ONLY: cutoff_screen_factor,&
      69              :                                               eri_2center,&
      70              :                                               eri_2center_derivs,&
      71              :                                               eri_3center,&
      72              :                                               eri_3center_derivs,&
      73              :                                               libint_potential_type
      74              :    USE libint_wrapper,                  ONLY: &
      75              :         cp_libint_cleanup_2eri, cp_libint_cleanup_2eri1, cp_libint_cleanup_3eri, &
      76              :         cp_libint_cleanup_3eri1, cp_libint_init_2eri, cp_libint_init_2eri1, cp_libint_init_3eri, &
      77              :         cp_libint_init_3eri1, cp_libint_set_contrdepth, cp_libint_t
      78              :    USE message_passing,                 ONLY: mp_para_env_type
      79              :    USE molecule_types,                  ONLY: molecule_type
      80              :    USE orbital_pointers,                ONLY: ncoset
      81              :    USE particle_types,                  ONLY: particle_type
      82              :    USE qs_environment_types,            ONLY: get_qs_env,&
      83              :                                               qs_environment_type
      84              :    USE qs_kind_types,                   ONLY: qs_kind_type
      85              :    USE qs_neighbor_list_types,          ONLY: &
      86              :         get_iterator_info, get_neighbor_list_set_p, neighbor_list_iterate, &
      87              :         neighbor_list_iterator_create, neighbor_list_iterator_p_type, &
      88              :         neighbor_list_iterator_release, neighbor_list_set_p_type, nl_sub_iterate, &
      89              :         release_neighbor_list_sets
      90              :    USE qs_neighbor_lists,               ONLY: atom2d_build,&
      91              :                                               atom2d_cleanup,&
      92              :                                               build_neighbor_lists,&
      93              :                                               local_atoms_type,&
      94              :                                               pair_radius_setup
      95              :    USE qs_tensors_types,                ONLY: &
      96              :         distribution_3d_destroy, distribution_3d_type, neighbor_list_3c_iterator_type, &
      97              :         neighbor_list_3c_type, symmetric_ij, symmetric_ijk, symmetric_jk, symmetric_none, &
      98              :         symmetrik_ik
      99              :    USE t_c_g0,                          ONLY: get_lmax_init,&
     100              :                                               init
     101              :    USE util,                            ONLY: get_limit
     102              : 
     103              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
     104              : #include "./base/base_uses.f90"
     105              : 
     106              :    IMPLICIT NONE
     107              : 
     108              :    PRIVATE
     109              : 
     110              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tensors'
     111              : 
     112              :    PUBLIC :: build_3c_neighbor_lists, &
     113              :              neighbor_list_3c_destroy, neighbor_list_3c_iterate, neighbor_list_3c_iterator_create, &
     114              :              neighbor_list_3c_iterator_destroy, get_3c_iterator_info, build_3c_integrals, &
     115              :              build_2c_neighbor_lists, build_2c_integrals, cutoff_screen_factor, &
     116              :              get_tensor_occupancy, compress_tensor, decompress_tensor, &
     117              :              build_3c_derivatives, build_2c_derivatives, calc_2c_virial, calc_3c_virial
     118              : 
     119              :    TYPE one_dim_int_array
     120              :       INTEGER, DIMENSION(:), ALLOCATABLE    :: array
     121              :    END TYPE
     122              : 
     123              :    ! cache size for integral compression
     124              :    INTEGER, PARAMETER, PRIVATE :: cache_size = 1024
     125              : 
     126              : CONTAINS
     127              : 
     128              : ! **************************************************************************************************
     129              : !> \brief Build 2-center neighborlists adapted to different operators
     130              : !>        This mainly wraps build_neighbor_lists for consistency with build_3c_neighbor_lists
     131              : !> \param ij_list 2c neighbor list for atom pairs i, j
     132              : !> \param basis_i basis object for atoms i
     133              : !> \param basis_j basis object for atoms j
     134              : !> \param potential_parameter ...
     135              : !> \param name name of 2c neighbor list
     136              : !> \param qs_env ...
     137              : !> \param sym_ij Symmetry in i, j (default .TRUE.)
     138              : !> \param molecular ...
     139              : !> \param dist_2d optionally a custom 2d distribution
     140              : !> \param pot_to_rad which radius (1 for i, 2 for j) should be adapted to incorporate potential
     141              : ! **************************************************************************************************
     142        11908 :    SUBROUTINE build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, &
     143              :                                       sym_ij, molecular, dist_2d, pot_to_rad)
     144              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     145              :          POINTER                                         :: ij_list
     146              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j
     147              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
     148              :       CHARACTER(LEN=*), INTENT(IN)                       :: name
     149              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     150              :       LOGICAL, INTENT(IN), OPTIONAL                      :: sym_ij, molecular
     151              :       TYPE(distribution_2d_type), OPTIONAL, POINTER      :: dist_2d
     152              :       INTEGER, INTENT(IN), OPTIONAL                      :: pot_to_rad
     153              : 
     154              :       INTEGER                                            :: ikind, nkind, pot_to_rad_prv
     155              :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: i_present, j_present
     156        11908 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: pair_radius
     157              :       REAL(kind=dp)                                      :: subcells
     158              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: i_radius, j_radius
     159        11908 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     160              :       TYPE(cell_type), POINTER                           :: cell
     161              :       TYPE(distribution_1d_type), POINTER                :: local_particles
     162              :       TYPE(distribution_2d_type), POINTER                :: dist_2d_prv
     163        11908 :       TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:)  :: atom2d
     164        11908 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     165        11908 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     166              : 
     167        11908 :       NULLIFY (atomic_kind_set, cell, local_particles, molecule_set, &
     168        11908 :                particle_set, dist_2d_prv)
     169              : 
     170        11908 :       IF (PRESENT(pot_to_rad)) THEN
     171         1392 :          pot_to_rad_prv = pot_to_rad
     172              :       ELSE
     173              :          pot_to_rad_prv = 1
     174              :       END IF
     175              : 
     176              :       CALL get_qs_env(qs_env, &
     177              :                       nkind=nkind, &
     178              :                       cell=cell, &
     179              :                       particle_set=particle_set, &
     180              :                       atomic_kind_set=atomic_kind_set, &
     181              :                       local_particles=local_particles, &
     182              :                       distribution_2d=dist_2d_prv, &
     183        11908 :                       molecule_set=molecule_set)
     184              : 
     185        11908 :       CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
     186              : 
     187        47632 :       ALLOCATE (i_present(nkind), j_present(nkind))
     188        47632 :       ALLOCATE (i_radius(nkind), j_radius(nkind))
     189              : 
     190        29761 :       i_present = .FALSE.
     191        29761 :       j_present = .FALSE.
     192        29761 :       i_radius = 0.0_dp
     193        29761 :       j_radius = 0.0_dp
     194              : 
     195        11908 :       IF (PRESENT(dist_2d)) dist_2d_prv => dist_2d
     196              : 
     197              :       !  Set up the radii, depending on the operator type
     198        11908 :       IF (potential_parameter%potential_type == do_potential_id) THEN
     199              : 
     200              :          !overlap => use the kind radius for both i and j
     201        25295 :          DO ikind = 1, nkind
     202        15133 :             IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
     203        15133 :                i_present(ikind) = .TRUE.
     204        15133 :                CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
     205              :             END IF
     206        25295 :             IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
     207        15133 :                j_present(ikind) = .TRUE.
     208        15133 :                CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
     209              :             END IF
     210              :          END DO
     211              : 
     212              :       ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
     213              : 
     214              :          !Coulomb operator, virtually infinite range => set j_radius to arbitrarily large number
     215         1110 :          DO ikind = 1, nkind
     216          736 :             IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
     217          736 :                i_present(ikind) = .TRUE.
     218          736 :                IF (pot_to_rad_prv == 1) i_radius(ikind) = 1000000.0_dp
     219              :             END IF
     220         1110 :             IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
     221          736 :                j_present(ikind) = .TRUE.
     222          736 :                IF (pot_to_rad_prv == 2) j_radius(ikind) = 1000000.0_dp
     223              :             END IF
     224              :          END DO !ikind
     225              : 
     226              :       ELSE IF (potential_parameter%potential_type == do_potential_truncated .OR. &
     227              :                potential_parameter%potential_type == do_potential_short .OR. &
     228              :                potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
     229              : 
     230              :          !Truncated coulomb/short range: set j_radius to r_cutoff + the kind_radii
     231         3356 :          DO ikind = 1, nkind
     232         1984 :             IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
     233         1984 :                i_present(ikind) = .TRUE.
     234         1984 :                CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
     235         1984 :                IF (pot_to_rad_prv == 1) i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
     236              :             END IF
     237         3356 :             IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
     238         1984 :                j_present(ikind) = .TRUE.
     239         1984 :                CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
     240         1984 :                IF (pot_to_rad_prv == 2) j_radius(ikind) = j_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
     241              :             END IF
     242              :          END DO
     243              : 
     244              :       ELSE
     245            0 :          CPABORT("Operator not implemented.")
     246              :       END IF
     247              : 
     248        47632 :       ALLOCATE (pair_radius(nkind, nkind))
     249        59540 :       pair_radius = 0.0_dp
     250        11908 :       CALL pair_radius_setup(i_present, j_present, i_radius, j_radius, pair_radius)
     251              : 
     252        53577 :       ALLOCATE (atom2d(nkind))
     253              : 
     254              :       CALL atom2d_build(atom2d, local_particles, dist_2d_prv, atomic_kind_set, &
     255        11908 :                         molecule_set, molecule_only=.FALSE., particle_set=particle_set)
     256              :       CALL build_neighbor_lists(ij_list, particle_set, atom2d, cell, pair_radius, subcells, &
     257        11908 :                                 symmetric=sym_ij, molecular=molecular, nlname=TRIM(name))
     258              : 
     259        11908 :       CALL atom2d_cleanup(atom2d)
     260              : 
     261        35724 :    END SUBROUTINE
     262              : 
     263              : ! **************************************************************************************************
     264              : !> \brief Build a 3-center neighbor list
     265              : !> \param ijk_list 3c neighbor list for atom triples i, j, k
     266              : !> \param basis_i basis object for atoms i
     267              : !> \param basis_j basis object for atoms j
     268              : !> \param basis_k basis object for atoms k
     269              : !> \param dist_3d 3d distribution object
     270              : !> \param potential_parameter ...
     271              : !> \param name name of 3c neighbor list
     272              : !> \param qs_env ...
     273              : !> \param sym_ij Symmetry in i, j (default .FALSE.)
     274              : !> \param sym_jk Symmetry in j, k (default .FALSE.)
     275              : !> \param sym_ik Symmetry in i, k (default .FALSE.)
     276              : !> \param molecular ??? not tested
     277              : !> \param op_pos ...
     278              : !> \param own_dist ...
     279              : ! **************************************************************************************************
     280          696 :    SUBROUTINE build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, &
     281              :                                       dist_3d, potential_parameter, name, qs_env, &
     282              :                                       sym_ij, sym_jk, sym_ik, molecular, op_pos, &
     283              :                                       own_dist)
     284              :       TYPE(neighbor_list_3c_type), INTENT(OUT)           :: ijk_list
     285              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j, basis_k
     286              :       TYPE(distribution_3d_type), INTENT(IN)             :: dist_3d
     287              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
     288              :       CHARACTER(LEN=*), INTENT(IN)                       :: name
     289              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     290              :       LOGICAL, INTENT(IN), OPTIONAL                      :: sym_ij, sym_jk, sym_ik, molecular
     291              :       INTEGER, INTENT(IN), OPTIONAL                      :: op_pos
     292              :       LOGICAL, INTENT(IN), OPTIONAL                      :: own_dist
     293              : 
     294              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_3c_neighbor_lists'
     295              : 
     296              :       INTEGER                                            :: handle, op_pos_prv, sym_level
     297              :       TYPE(libint_potential_type)                        :: pot_par_1, pot_par_2
     298              : 
     299          696 :       CALL timeset(routineN, handle)
     300              : 
     301          696 :       IF (PRESENT(op_pos)) THEN
     302          512 :          op_pos_prv = op_pos
     303              :       ELSE
     304              :          op_pos_prv = 1
     305              :       END IF
     306              : 
     307          490 :       SELECT CASE (op_pos_prv)
     308              :       CASE (1)
     309          490 :          pot_par_1 = potential_parameter
     310          490 :          pot_par_2%potential_type = do_potential_id
     311              :       CASE (2)
     312          206 :          pot_par_2 = potential_parameter
     313          512 :          pot_par_1%potential_type = do_potential_id
     314              :       END SELECT
     315              : 
     316              :       CALL build_2c_neighbor_lists(ijk_list%ij_list, basis_i, basis_j, pot_par_1, TRIM(name)//"_sub_1", &
     317              :                                    qs_env, sym_ij=.FALSE., molecular=molecular, &
     318          696 :                                    dist_2d=dist_3d%dist_2d_1, pot_to_rad=1)
     319              : 
     320              :       CALL build_2c_neighbor_lists(ijk_list%jk_list, basis_j, basis_k, pot_par_2, TRIM(name)//"_sub_2", &
     321              :                                    qs_env, sym_ij=.FALSE., molecular=molecular, &
     322          696 :                                    dist_2d=dist_3d%dist_2d_2, pot_to_rad=2)
     323              : 
     324          696 :       ijk_list%sym = symmetric_none
     325              : 
     326          696 :       sym_level = 0
     327          696 :       IF (PRESENT(sym_ij)) THEN
     328          160 :          IF (sym_ij) THEN
     329            0 :             ijk_list%sym = symmetric_ij
     330            0 :             sym_level = sym_level + 1
     331              :          END IF
     332              :       END IF
     333              : 
     334          696 :       IF (PRESENT(sym_jk)) THEN
     335          486 :          IF (sym_jk) THEN
     336          428 :             ijk_list%sym = symmetric_jk
     337          428 :             sym_level = sym_level + 1
     338              :          END IF
     339              :       END IF
     340              : 
     341          696 :       IF (PRESENT(sym_ik)) THEN
     342            0 :          IF (sym_ik) THEN
     343            0 :             ijk_list%sym = symmetrik_ik
     344            0 :             sym_level = sym_level + 1
     345              :          END IF
     346              :       END IF
     347              : 
     348          696 :       IF (sym_level >= 2) THEN
     349            0 :          ijk_list%sym = symmetric_ijk
     350              :       END IF
     351              : 
     352          696 :       ijk_list%dist_3d = dist_3d
     353          696 :       IF (PRESENT(own_dist)) THEN
     354          696 :          ijk_list%owns_dist = own_dist
     355              :       ELSE
     356            0 :          ijk_list%owns_dist = .FALSE.
     357              :       END IF
     358              : 
     359          696 :       CALL timestop(handle)
     360          696 :    END SUBROUTINE
     361              : 
     362              : ! **************************************************************************************************
     363              : !> \brief Symmetry criterion
     364              : !> \param a ...
     365              : !> \param b ...
     366              : !> \return ...
     367              : ! **************************************************************************************************
     368      4207758 :    PURE FUNCTION include_symmetric(a, b)
     369              :       INTEGER, INTENT(IN)                                :: a, b
     370              :       LOGICAL                                            :: include_symmetric
     371              : 
     372      4207758 :       IF (a > b) THEN
     373      1691358 :          include_symmetric = (MODULO(a + b, 2) /= 0)
     374              :       ELSE
     375      2516400 :          include_symmetric = (MODULO(a + b, 2) == 0)
     376              :       END IF
     377              : 
     378      4207758 :    END FUNCTION
     379              : 
     380              : ! **************************************************************************************************
     381              : !> \brief Destroy 3c neighborlist
     382              : !> \param ijk_list ...
     383              : ! **************************************************************************************************
     384          696 :    SUBROUTINE neighbor_list_3c_destroy(ijk_list)
     385              :       TYPE(neighbor_list_3c_type), INTENT(INOUT)         :: ijk_list
     386              : 
     387          696 :       CALL release_neighbor_list_sets(ijk_list%ij_list)
     388          696 :       CALL release_neighbor_list_sets(ijk_list%jk_list)
     389              : 
     390          696 :       IF (ijk_list%owns_dist) THEN
     391          696 :          CALL distribution_3d_destroy(ijk_list%dist_3d)
     392              :       END IF
     393              : 
     394          696 :    END SUBROUTINE
     395              : 
     396              : ! **************************************************************************************************
     397              : !> \brief Create a 3-center neighborlist iterator
     398              : !> \param iterator ...
     399              : !> \param ijk_nl ...
     400              : ! **************************************************************************************************
     401        69768 :    SUBROUTINE neighbor_list_3c_iterator_create(iterator, ijk_nl)
     402              :       TYPE(neighbor_list_3c_iterator_type), INTENT(OUT)  :: iterator
     403              :       TYPE(neighbor_list_3c_type), INTENT(IN)            :: ijk_nl
     404              : 
     405              :       CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterator_create'
     406              : 
     407              :       INTEGER                                            :: handle
     408              : 
     409         7752 :       CALL timeset(routineN, handle)
     410              : 
     411         7752 :       CALL neighbor_list_iterator_create(iterator%iter_ij, ijk_nl%ij_list)
     412         7752 :       CALL neighbor_list_iterator_create(iterator%iter_jk, ijk_nl%jk_list, search=.TRUE.)
     413              : 
     414         7752 :       iterator%iter_level = 0
     415         7752 :       iterator%ijk_nl = ijk_nl
     416              : 
     417        23256 :       iterator%bounds_i = 0
     418        23256 :       iterator%bounds_j = 0
     419        23256 :       iterator%bounds_k = 0
     420              : 
     421         7752 :       CALL timestop(handle)
     422         7752 :    END SUBROUTINE
     423              : 
     424              : ! **************************************************************************************************
     425              : !> \brief impose atomic upper and lower bounds
     426              : !> \param iterator ...
     427              : !> \param bounds_i ...
     428              : !> \param bounds_j ...
     429              : !> \param bounds_k ...
     430              : ! **************************************************************************************************
     431         7672 :    SUBROUTINE nl_3c_iter_set_bounds(iterator, bounds_i, bounds_j, bounds_k)
     432              :       TYPE(neighbor_list_3c_iterator_type), &
     433              :          INTENT(INOUT)                                   :: iterator
     434              :       INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL        :: bounds_i, bounds_j, bounds_k
     435              : 
     436        24190 :       IF (PRESENT(bounds_i)) iterator%bounds_i = bounds_i
     437        14260 :       IF (PRESENT(bounds_j)) iterator%bounds_j = bounds_j
     438        18292 :       IF (PRESENT(bounds_k)) iterator%bounds_k = bounds_k
     439              : 
     440         7672 :    END SUBROUTINE
     441              : 
     442              : ! **************************************************************************************************
     443              : !> \brief Destroy 3c-nl iterator
     444              : !> \param iterator ...
     445              : ! **************************************************************************************************
     446         7752 :    SUBROUTINE neighbor_list_3c_iterator_destroy(iterator)
     447              :       TYPE(neighbor_list_3c_iterator_type), &
     448              :          INTENT(INOUT)                                   :: iterator
     449              : 
     450              :       CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterator_destroy'
     451              : 
     452              :       INTEGER                                            :: handle
     453              : 
     454         7752 :       CALL timeset(routineN, handle)
     455         7752 :       CALL neighbor_list_iterator_release(iterator%iter_ij)
     456         7752 :       CALL neighbor_list_iterator_release(iterator%iter_jk)
     457         7752 :       NULLIFY (iterator%iter_ij)
     458         7752 :       NULLIFY (iterator%iter_jk)
     459              : 
     460         7752 :       CALL timestop(handle)
     461         7752 :    END SUBROUTINE
     462              : 
     463              : ! **************************************************************************************************
     464              : !> \brief Iterate 3c-nl iterator
     465              : !> \param iterator ...
     466              : !> \return 0 if successful; 1 if end was reached
     467              : ! **************************************************************************************************
     468      5706732 :    RECURSIVE FUNCTION neighbor_list_3c_iterate(iterator) RESULT(iter_stat)
     469              :       TYPE(neighbor_list_3c_iterator_type), &
     470              :          INTENT(INOUT)                                   :: iterator
     471              :       INTEGER                                            :: iter_stat
     472              : 
     473              :       INTEGER                                            :: iatom, iter_level, jatom, jatom_1, &
     474              :                                                             jatom_2, katom
     475              :       LOGICAL                                            :: skip_this
     476              : 
     477      5706732 :       iter_level = iterator%iter_level
     478              : 
     479      5706732 :       IF (iter_level == 0) THEN
     480       195772 :          iter_stat = neighbor_list_iterate(iterator%iter_ij)
     481              : 
     482       195772 :          IF (iter_stat /= 0) THEN
     483              :             RETURN
     484              :          END IF
     485              : 
     486       188020 :          CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom)
     487       188020 :          skip_this = .FALSE.
     488              :          IF ((iterator%bounds_i(1) > 0 .AND. iatom < iterator%bounds_i(1)) &
     489       188020 :              .OR. (iterator%bounds_i(2) > 0 .AND. iatom > iterator%bounds_i(2))) skip_this = .TRUE.
     490              :          IF ((iterator%bounds_j(1) > 0 .AND. jatom < iterator%bounds_j(1)) &
     491       188020 :              .OR. (iterator%bounds_j(2) > 0 .AND. jatom > iterator%bounds_j(2))) skip_this = .TRUE.
     492              : 
     493       188020 :          IF (skip_this) THEN
     494        85323 :             iter_stat = neighbor_list_3c_iterate(iterator)
     495        85323 :             RETURN
     496              :          END IF
     497              : 
     498              :       END IF
     499      5613657 :       iter_stat = nl_sub_iterate(iterator%iter_jk, iterator%iter_ij)
     500      5613657 :       IF (iter_stat /= 0) THEN
     501       102697 :          iterator%iter_level = 0
     502       102697 :          iter_stat = neighbor_list_3c_iterate(iterator)
     503       102697 :          RETURN
     504              :       ELSE
     505      5510960 :          iterator%iter_level = 1
     506              :       END IF
     507              : 
     508              :       CPASSERT(iter_stat == 0)
     509              :       CPASSERT(iterator%iter_level == 1)
     510      5510960 :       CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom_1)
     511      5510960 :       CALL get_iterator_info(iterator%iter_jk, iatom=jatom_2, jatom=katom)
     512              : 
     513      5510960 :       CPASSERT(jatom_1 == jatom_2)
     514      5510960 :       jatom = jatom_1
     515              : 
     516      5510960 :       skip_this = .FALSE.
     517              :       IF ((iterator%bounds_k(1) > 0 .AND. katom < iterator%bounds_k(1)) &
     518      5510960 :           .OR. (iterator%bounds_k(2) > 0 .AND. katom > iterator%bounds_k(2))) skip_this = .TRUE.
     519              : 
     520              :       IF (skip_this) THEN
     521       312756 :          iter_stat = neighbor_list_3c_iterate(iterator)
     522       312756 :          RETURN
     523              :       END IF
     524              : 
     525      5198204 :       SELECT CASE (iterator%ijk_nl%sym)
     526              :       CASE (symmetric_none)
     527            0 :          skip_this = .FALSE.
     528              :       CASE (symmetric_ij)
     529            0 :          skip_this = .NOT. include_symmetric(iatom, jatom)
     530              :       CASE (symmetric_jk)
     531      4207758 :          skip_this = .NOT. include_symmetric(jatom, katom)
     532              :       CASE (symmetrik_ik)
     533            0 :          skip_this = .NOT. include_symmetric(iatom, katom)
     534              :       CASE (symmetric_ijk)
     535            0 :          skip_this = .NOT. include_symmetric(iatom, jatom) .OR. .NOT. include_symmetric(jatom, katom)
     536              :       CASE DEFAULT
     537      5198204 :          CPABORT("should not happen")
     538              :       END SELECT
     539              : 
     540      4207758 :       IF (skip_this) THEN
     541      1692660 :          iter_stat = neighbor_list_3c_iterate(iterator)
     542      1692660 :          RETURN
     543              :       END IF
     544              : 
     545              :    END FUNCTION
     546              : 
     547              : ! **************************************************************************************************
     548              : !> \brief Get info of current iteration
     549              : !> \param iterator ...
     550              : !> \param ikind ...
     551              : !> \param jkind ...
     552              : !> \param kkind ...
     553              : !> \param nkind ...
     554              : !> \param iatom ...
     555              : !> \param jatom ...
     556              : !> \param katom ...
     557              : !> \param rij ...
     558              : !> \param rjk ...
     559              : !> \param rik ...
     560              : !> \param cell_j ...
     561              : !> \param cell_k ...
     562              : !> \return ...
     563              : ! **************************************************************************************************
     564      3505544 :    SUBROUTINE get_3c_iterator_info(iterator, ikind, jkind, kkind, nkind, iatom, jatom, katom, &
     565              :                                    rij, rjk, rik, cell_j, cell_k)
     566              :       TYPE(neighbor_list_3c_iterator_type), &
     567              :          INTENT(INOUT)                                   :: iterator
     568              :       INTEGER, INTENT(OUT), OPTIONAL                     :: ikind, jkind, kkind, nkind, iatom, &
     569              :                                                             jatom, katom
     570              :       REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: rij, rjk, rik
     571              :       INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL       :: cell_j, cell_k
     572              : 
     573              :       INTEGER, DIMENSION(2)                              :: atoms_1, atoms_2, kinds_1, kinds_2
     574              :       INTEGER, DIMENSION(3)                              :: cell_1, cell_2
     575              :       REAL(KIND=dp), DIMENSION(3)                        :: r_1, r_2
     576              : 
     577      3505544 :       CPASSERT(iterator%iter_level == 1)
     578              : 
     579              :       CALL get_iterator_info(iterator%iter_ij, &
     580              :                              ikind=kinds_1(1), jkind=kinds_1(2), nkind=nkind, &
     581              :                              iatom=atoms_1(1), jatom=atoms_1(2), r=r_1, &
     582      3505544 :                              cell=cell_1)
     583              : 
     584              :       CALL get_iterator_info(iterator%iter_jk, &
     585              :                              ikind=kinds_2(1), jkind=kinds_2(2), &
     586              :                              iatom=atoms_2(1), jatom=atoms_2(2), r=r_2, &
     587      3505544 :                              cell=cell_2)
     588              : 
     589      3505544 :       IF (PRESENT(ikind)) ikind = kinds_1(1)
     590      3505544 :       IF (PRESENT(jkind)) jkind = kinds_1(2)
     591      3505544 :       IF (PRESENT(kkind)) kkind = kinds_2(2)
     592      3505544 :       IF (PRESENT(iatom)) iatom = atoms_1(1)
     593      3505544 :       IF (PRESENT(jatom)) jatom = atoms_1(2)
     594      3505544 :       IF (PRESENT(katom)) katom = atoms_2(2)
     595              : 
     596      3505544 :       IF (PRESENT(rij)) rij = r_1
     597      3505544 :       IF (PRESENT(rjk)) rjk = r_2
     598     17527720 :       IF (PRESENT(rik)) rik = r_1 + r_2
     599              : 
     600      3505544 :       IF (PRESENT(cell_j)) cell_j = cell_1
     601     16848520 :       IF (PRESENT(cell_k)) cell_k = cell_1 + cell_2
     602              : 
     603      3505544 :    END SUBROUTINE
     604              : 
     605              : ! **************************************************************************************************
     606              : !> \brief Allocate blocks of a 3-center tensor based on neighborlist
     607              : !> \param t3c empty DBCSR tensor
     608              : !>            Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages)
     609              : !>            if k-points are used
     610              : !> \param nl_3c 3-center neighborlist
     611              : !> \param basis_i ...
     612              : !> \param basis_j ...
     613              : !> \param basis_k ...
     614              : !> \param qs_env ...
     615              : !> \param potential_parameter ...
     616              : !> \param op_pos ...
     617              : !> \param do_kpoints ...
     618              : !> \param do_hfx_kpoints ...
     619              : !> \param bounds_i ...
     620              : !> \param bounds_j ...
     621              : !> \param bounds_k ...
     622              : !> \param RI_range ...
     623              : !> \param img_to_RI_cell ...
     624              : !> \param cell_to_index ...
     625              : !> \param cell_sym ...
     626              : ! **************************************************************************************************
     627         2552 :    SUBROUTINE alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, op_pos, &
     628              :                              do_kpoints, do_hfx_kpoints, bounds_i, bounds_j, bounds_k, RI_range, &
     629         2552 :                              img_to_RI_cell, cell_to_index, cell_sym)
     630              :       TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT)     :: t3c
     631              :       TYPE(neighbor_list_3c_type), INTENT(INOUT)         :: nl_3c
     632              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j, basis_k
     633              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     634              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
     635              :       INTEGER, INTENT(IN), OPTIONAL                      :: op_pos
     636              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_kpoints, do_hfx_kpoints
     637              :       INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL        :: bounds_i, bounds_j, bounds_k
     638              :       REAL(dp), INTENT(IN), OPTIONAL                     :: RI_range
     639              :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: img_to_RI_cell
     640              :       INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER     :: cell_to_index
     641              :       LOGICAL, INTENT(IN), OPTIONAL                      :: cell_sym
     642              : 
     643              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'alloc_block_3c'
     644              : 
     645              :       INTEGER                                            :: handle, iatom, ikind, j_img, jatom, &
     646              :                                                             jcell, jkind, k_img, katom, kcell, &
     647              :                                                             kkind, natom, ncell_RI, nimg, op_ij, &
     648              :                                                             op_jk, op_pos_prv
     649              :       INTEGER(int_8)                                     :: a, b, nblk_per_thread
     650         2552 :       INTEGER(int_8), ALLOCATABLE, DIMENSION(:, :)       :: nblk
     651         2552 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: img_to_RI_cell_prv
     652              :       INTEGER, DIMENSION(3)                              :: blk_idx, cell_j, cell_k, &
     653              :                                                             kp_index_lbounds, kp_index_ubounds
     654              :       LOGICAL                                            :: cell_sym_prv, do_hfx_kpoints_prv, &
     655              :                                                             do_kpoints_prv
     656              :       REAL(KIND=dp)                                      :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
     657              :                                                             kind_radius_i, kind_radius_j, &
     658              :                                                             kind_radius_k
     659              :       REAL(KIND=dp), DIMENSION(3)                        :: rij, rik, rjk
     660         2552 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     661              :       TYPE(cell_type), POINTER                           :: cell
     662              :       TYPE(dft_control_type), POINTER                    :: dft_control
     663              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     664              :       TYPE(neighbor_list_3c_iterator_type)               :: nl_3c_iter
     665              :       TYPE(one_dim_int_array), ALLOCATABLE, &
     666         2552 :          DIMENSION(:, :)                                 :: alloc_i, alloc_j, alloc_k
     667         2552 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     668              : 
     669         2552 :       CALL timeset(routineN, handle)
     670         2552 :       NULLIFY (qs_kind_set, atomic_kind_set, cell)
     671              : 
     672         2552 :       IF (PRESENT(do_kpoints)) THEN
     673          464 :          do_kpoints_prv = do_kpoints
     674              :       ELSE
     675              :          do_kpoints_prv = .FALSE.
     676              :       END IF
     677         2552 :       IF (PRESENT(do_hfx_kpoints)) THEN
     678          222 :          do_hfx_kpoints_prv = do_hfx_kpoints
     679              :       ELSE
     680              :          do_hfx_kpoints_prv = .FALSE.
     681              :       END IF
     682          222 :       IF (do_hfx_kpoints_prv) THEN
     683          222 :          CPASSERT(do_kpoints_prv)
     684          222 :          CPASSERT(PRESENT(RI_range))
     685          222 :          CPASSERT(PRESENT(img_to_RI_cell))
     686              :       END IF
     687              : 
     688         2330 :       IF (PRESENT(img_to_RI_cell)) THEN
     689          666 :          ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
     690        10432 :          img_to_RI_cell_prv(:) = img_to_RI_cell
     691              :       END IF
     692              : 
     693         2552 :       IF (PRESENT(cell_sym)) THEN
     694         1978 :          cell_sym_prv = cell_sym
     695              :       ELSE
     696              :          cell_sym_prv = .FALSE.
     697              :       END IF
     698              : 
     699         2552 :       dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
     700              : 
     701         2552 :       op_ij = do_potential_id; op_jk = do_potential_id
     702              : 
     703         2552 :       IF (PRESENT(op_pos)) THEN
     704         2552 :          op_pos_prv = op_pos
     705              :       ELSE
     706              :          op_pos_prv = 1
     707              :       END IF
     708              : 
     709         2330 :       SELECT CASE (op_pos_prv)
     710              :       CASE (1)
     711         2330 :          op_ij = potential_parameter%potential_type
     712              :       CASE (2)
     713         2552 :          op_jk = potential_parameter%potential_type
     714              :       END SELECT
     715              : 
     716              :       IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
     717         2552 :           .OR. op_ij == do_potential_mix_cl_trunc) THEN
     718         1488 :          dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
     719         1488 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
     720         1064 :       ELSEIF (op_ij == do_potential_coulomb) THEN
     721          436 :          dr_ij = 1000000.0_dp
     722          436 :          dr_ik = 1000000.0_dp
     723              :       END IF
     724              : 
     725              :       IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
     726         2552 :           .OR. op_jk == do_potential_mix_cl_trunc) THEN
     727           20 :          dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
     728           20 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
     729         2532 :       ELSEIF (op_jk == do_potential_coulomb) THEN
     730            0 :          dr_jk = 1000000.0_dp
     731            0 :          dr_ik = 1000000.0_dp
     732              :       END IF
     733              : 
     734              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, natom=natom, &
     735         2552 :                       dft_control=dft_control, para_env=para_env, cell=cell)
     736              : 
     737         2552 :       IF (do_kpoints_prv) THEN
     738          234 :          CPASSERT(PRESENT(cell_to_index))
     739          234 :          CPASSERT(ASSOCIATED(cell_to_index))
     740              : !         nimg = dft_control%nimages
     741        33696 :          nimg = MAXVAL(cell_to_index)
     742          234 :          ncell_RI = nimg
     743          234 :          IF (do_hfx_kpoints_prv) THEN
     744          222 :             nimg = SIZE(t3c, 1)
     745          222 :             ncell_RI = SIZE(t3c, 2)
     746              :          END IF
     747              :       ELSE
     748              :          nimg = 1
     749              :          ncell_RI = 1
     750              :       END IF
     751              : 
     752         2552 :       IF (do_kpoints_prv) THEN
     753          936 :          kp_index_lbounds = LBOUND(cell_to_index)
     754          936 :          kp_index_ubounds = UBOUND(cell_to_index)
     755              :       END IF
     756              : 
     757              :       !Do a first loop over the nl and count the blocks present
     758        10208 :       ALLOCATE (nblk(nimg, ncell_RI))
     759        14156 :       nblk(:, :) = 0
     760              : 
     761         2552 :       CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
     762              : 
     763         2552 :       CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
     764              : 
     765      1089046 :       DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
     766              :          CALL get_3c_iterator_info(nl_3c_iter, iatom=iatom, ikind=ikind, jkind=jkind, kkind=kkind, &
     767      1086494 :                                    rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
     768              : 
     769      4345976 :          djk = NORM2(rjk)
     770      4345976 :          dij = NORM2(rij)
     771      4345976 :          dik = NORM2(rik)
     772              : 
     773      1086494 :          IF (do_kpoints_prv) THEN
     774              : 
     775       801066 :             IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
     776              :                 ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
     777              : 
     778       114438 :             jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
     779       114438 :             IF (jcell > nimg .OR. jcell < 1) CYCLE
     780              : 
     781       702612 :             IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
     782              :                 ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
     783              : 
     784        96315 :             kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
     785        96315 :             IF (kcell > nimg .OR. kcell < 1) CYCLE
     786        66939 :             IF (do_hfx_kpoints_prv) THEN
     787        63309 :                IF (dik > RI_range) CYCLE
     788              :                kcell = 1
     789              :             END IF
     790              :          ELSE
     791              :             jcell = 1; kcell = 1
     792              :          END IF
     793              : 
     794       985236 :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
     795       985236 :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
     796       985236 :          CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
     797              : 
     798       985236 :          IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
     799       985236 :          IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
     800       985236 :          IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
     801              : 
     802      1080522 :          nblk(jcell, kcell) = nblk(jcell, kcell) + 1
     803              :       END DO
     804         2552 :       CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
     805              : 
     806              :       !Do a second loop over the nl to give block indices
     807        21812 :       ALLOCATE (alloc_i(nimg, ncell_RI))
     808        19260 :       ALLOCATE (alloc_j(nimg, ncell_RI))
     809        19260 :       ALLOCATE (alloc_k(nimg, ncell_RI))
     810         5180 :       DO k_img = 1, ncell_RI
     811        14156 :          DO j_img = 1, nimg
     812        23711 :             ALLOCATE (alloc_i(j_img, k_img)%array(nblk(j_img, k_img)))
     813        23711 :             ALLOCATE (alloc_j(j_img, k_img)%array(nblk(j_img, k_img)))
     814        26339 :             ALLOCATE (alloc_k(j_img, k_img)%array(nblk(j_img, k_img)))
     815              :          END DO
     816              :       END DO
     817        14156 :       nblk(:, :) = 0
     818              : 
     819         2552 :       CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
     820              : 
     821         2552 :       CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
     822              : 
     823      1089046 :       DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
     824              :          CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
     825              :                                    iatom=iatom, jatom=jatom, katom=katom, &
     826      1086494 :                                    rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
     827              : 
     828      4345976 :          djk = NORM2(rjk)
     829      4345976 :          dij = NORM2(rij)
     830      4345976 :          dik = NORM2(rik)
     831              : 
     832      1086494 :          IF (do_kpoints_prv) THEN
     833              : 
     834       801066 :             IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
     835              :                 ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
     836              : 
     837       114438 :             jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
     838       114438 :             IF (jcell > nimg .OR. jcell < 1) CYCLE
     839              : 
     840       702612 :             IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
     841              :                 ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
     842              : 
     843        96315 :             kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
     844        96315 :             IF (kcell > nimg .OR. kcell < 1) CYCLE
     845        66939 :             IF (do_hfx_kpoints_prv) THEN
     846        63309 :                IF (dik > RI_range) CYCLE
     847         9550 :                kcell = img_to_RI_cell_prv(kcell)
     848              :             END IF
     849              :          ELSE
     850              :             jcell = 1; kcell = 1
     851              :          END IF
     852              : 
     853      3940944 :          blk_idx = [iatom, jatom, katom]
     854       985236 :          IF (do_hfx_kpoints_prv) THEN
     855         9550 :             blk_idx(3) = (kcell - 1)*natom + katom
     856         9550 :             kcell = 1
     857              :          END IF
     858              : 
     859       985236 :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
     860       985236 :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
     861       985236 :          CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
     862              : 
     863       985236 :          IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
     864       985236 :          IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
     865       985236 :          IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
     866              : 
     867       497535 :          nblk(jcell, kcell) = nblk(jcell, kcell) + 1
     868              : 
     869              :          !Note: there may be repeated indices due to periodic images => dbt_reserve_blocks takes care of it
     870       497535 :          alloc_i(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(1)
     871       497535 :          alloc_j(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(2)
     872      1080522 :          alloc_k(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(3)
     873              : 
     874              :       END DO
     875         2552 :       CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
     876              : 
     877              : !TODO: Parallelize creation of block list.
     878              : !$OMP PARALLEL DEFAULT(NONE) SHARED(t3c,nimg,nblk,alloc_i,alloc_j,alloc_k,ncell_RI,cell_sym_prv) &
     879         2552 : !$OMP PRIVATE(k_img,j_img,nblk_per_thread,A,b)
     880              :       DO k_img = 1, ncell_RI
     881              :          DO j_img = 1, nimg
     882              :             IF (cell_sym_prv .AND. j_img < k_img) CYCLE
     883              :             IF (ALLOCATED(alloc_i(j_img, k_img)%array)) THEN
     884              :                nblk_per_thread = nblk(j_img, k_img)/omp_get_num_threads() + 1
     885              :                a = omp_get_thread_num()*nblk_per_thread + 1
     886              :                b = MIN(a + nblk_per_thread, nblk(j_img, k_img))
     887              :                CALL dbt_reserve_blocks(t3c(j_img, k_img), &
     888              :                                        alloc_i(j_img, k_img)%array(a:b), &
     889              :                                        alloc_j(j_img, k_img)%array(a:b), &
     890              :                                        alloc_k(j_img, k_img)%array(a:b))
     891              :             END IF
     892              :          END DO
     893              :       END DO
     894              : !$OMP END PARALLEL
     895              : 
     896         2552 :       CALL timestop(handle)
     897              : 
     898        52448 :    END SUBROUTINE
     899              : 
     900              : ! **************************************************************************************************
     901              : !> \brief Build 3-center derivative tensors
     902              : !> \param t3c_der_i empty DBCSR tensor which will contain the 1st center derivatives
     903              : !> \param t3c_der_k empty DBCSR tensor which will contain the 3rd center derivatives
     904              : !> \param filter_eps Filter threshold for tensor blocks
     905              : !> \param qs_env ...
     906              : !> \param nl_3c 3-center neighborlist
     907              : !> \param basis_i ...
     908              : !> \param basis_j ...
     909              : !> \param basis_k ...
     910              : !> \param potential_parameter ...
     911              : !> \param der_eps neglect integrals smaller than der_eps
     912              : !> \param op_pos operator position.
     913              : !>        1: calculate (i|jk) integrals,
     914              : !>        2: calculate (ij|k) integrals
     915              : !> \param do_kpoints ...
     916              : !> this routine requires that libint has been static initialised somewhere else
     917              : !> \param do_hfx_kpoints ...
     918              : !> \param bounds_i ...
     919              : !> \param bounds_j ...
     920              : !> \param bounds_k ...
     921              : !> \param RI_range ...
     922              : !> \param img_to_RI_cell ...
     923              : ! **************************************************************************************************
     924          574 :    SUBROUTINE build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, &
     925          574 :                                    nl_3c, basis_i, basis_j, basis_k, &
     926              :                                    potential_parameter, der_eps, &
     927              :                                    op_pos, do_kpoints, do_hfx_kpoints, &
     928              :                                    bounds_i, bounds_j, bounds_k, &
     929          574 :                                    RI_range, img_to_RI_cell)
     930              : 
     931              :       TYPE(dbt_type), DIMENSION(:, :, :), INTENT(INOUT)  :: t3c_der_i, t3c_der_k
     932              :       REAL(KIND=dp), INTENT(IN)                          :: filter_eps
     933              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     934              :       TYPE(neighbor_list_3c_type), INTENT(INOUT)         :: nl_3c
     935              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j, basis_k
     936              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
     937              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: der_eps
     938              :       INTEGER, INTENT(IN), OPTIONAL                      :: op_pos
     939              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_kpoints, do_hfx_kpoints
     940              :       INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL        :: bounds_i, bounds_j, bounds_k
     941              :       REAL(dp), INTENT(IN), OPTIONAL                     :: RI_range
     942              :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: img_to_RI_cell
     943              : 
     944              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_derivatives'
     945              : 
     946              :       INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
     947              :          block_start_k, egfi, handle, handle2, i, i_img, i_xyz, iatom, ibasis, ikind, ilist, imax, &
     948              :          iset, j_img, jatom, jcell, jkind, jset, katom, kcell, kkind, kset, m_max, max_ncoj, &
     949              :          max_nset, max_nsgfi, maxli, maxlj, maxlk, mepos, natom, nbasis, ncell_RI, ncoi, ncoj, &
     950              :          ncok, nimg, nseti, nsetj, nsetk, nthread, op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, &
     951              :          unit_id
     952          574 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: img_to_RI_cell_prv
     953              :       INTEGER, DIMENSION(2)                              :: bo
     954              :       INTEGER, DIMENSION(3)                              :: blk_idx, blk_size, cell_j, cell_k, &
     955              :                                                             kp_index_lbounds, kp_index_ubounds, sp
     956          574 :       INTEGER, DIMENSION(:), POINTER                     :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
     957         1148 :                                                             lmin_k, npgfi, npgfj, npgfk, nsgfi, &
     958          574 :                                                             nsgfj, nsgfk
     959          574 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf_i, first_sgf_j, first_sgf_k
     960          574 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     961              :       LOGICAL                                            :: do_hfx_kpoints_prv, do_kpoints_prv, &
     962              :                                                             found, skip
     963              :       LOGICAL, DIMENSION(3)                              :: block_j_not_zero, block_k_not_zero, &
     964              :                                                             der_j_zero, der_k_zero
     965              :       REAL(dp), DIMENSION(3)                             :: der_ext_i, der_ext_j, der_ext_k
     966              :       REAL(KIND=dp)                                      :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
     967              :                                                             kind_radius_i, kind_radius_j, &
     968              :                                                             kind_radius_k, prefac
     969          574 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ccp_buffer, cpp_buffer, &
     970          574 :                                                             max_contraction_i, max_contraction_j, &
     971          574 :                                                             max_contraction_k
     972          574 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dijk_contr, dummy_block_t
     973          574 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: block_t_i, block_t_j, block_t_k, dijk_j, &
     974         1148 :                                                             dijk_k, tmp_ijk_i, tmp_ijk_j
     975              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rij, rik, rj, rjk, rk
     976          574 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_i, set_radius_j, set_radius_k
     977          574 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
     978         1148 :                                                             sphi_k, zeti, zetj, zetk
     979          574 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     980         1148 :       TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER     :: spi, spk, tspj
     981              :       TYPE(cp_libint_t)                                  :: lib
     982         4018 :       TYPE(dbt_type)                                     :: t3c_tmp
     983          574 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t3c_template
     984          574 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :, :)    :: t3c_der_j
     985              :       TYPE(dft_control_type), POINTER                    :: dft_control
     986              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     987              :       TYPE(kpoint_type), POINTER                         :: kpoints
     988              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     989              :       TYPE(neighbor_list_3c_iterator_type)               :: nl_3c_iter
     990          574 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     991              : 
     992          574 :       CALL timeset(routineN, handle)
     993              : 
     994          574 :       IF (PRESENT(do_kpoints)) THEN
     995           82 :          do_kpoints_prv = do_kpoints
     996              :       ELSE
     997              :          do_kpoints_prv = .FALSE.
     998              :       END IF
     999              : 
    1000          574 :       IF (PRESENT(do_hfx_kpoints)) THEN
    1001           82 :          do_hfx_kpoints_prv = do_hfx_kpoints
    1002              :       ELSE
    1003              :          do_hfx_kpoints_prv = .FALSE.
    1004              :       END IF
    1005           82 :       IF (do_hfx_kpoints_prv) THEN
    1006           82 :          CPASSERT(do_kpoints_prv)
    1007           82 :          CPASSERT(PRESENT(RI_range))
    1008           82 :          CPASSERT(PRESENT(img_to_RI_cell))
    1009              :       END IF
    1010              : 
    1011          492 :       IF (PRESENT(img_to_RI_cell)) THEN
    1012          246 :          ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
    1013         3834 :          img_to_RI_cell_prv(:) = img_to_RI_cell
    1014              :       END IF
    1015              : 
    1016          574 :       op_ij = do_potential_id; op_jk = do_potential_id
    1017              : 
    1018          574 :       IF (PRESENT(op_pos)) THEN
    1019          574 :          op_pos_prv = op_pos
    1020              :       ELSE
    1021            0 :          op_pos_prv = 1
    1022              :       END IF
    1023              : 
    1024          492 :       SELECT CASE (op_pos_prv)
    1025              :       CASE (1)
    1026          492 :          op_ij = potential_parameter%potential_type
    1027              :       CASE (2)
    1028          574 :          op_jk = potential_parameter%potential_type
    1029              :       END SELECT
    1030              : 
    1031          574 :       dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
    1032              : 
    1033              :       IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
    1034          574 :           .OR. op_ij == do_potential_mix_cl_trunc) THEN
    1035           50 :          dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
    1036           50 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
    1037          524 :       ELSEIF (op_ij == do_potential_coulomb) THEN
    1038          244 :          dr_ij = 1000000.0_dp
    1039          244 :          dr_ik = 1000000.0_dp
    1040              :       END IF
    1041              : 
    1042              :       IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
    1043          574 :           .OR. op_jk == do_potential_mix_cl_trunc) THEN
    1044            8 :          dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
    1045            8 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
    1046          566 :       ELSEIF (op_jk == do_potential_coulomb) THEN
    1047            0 :          dr_jk = 1000000.0_dp
    1048            0 :          dr_ik = 1000000.0_dp
    1049              :       END IF
    1050              : 
    1051          574 :       NULLIFY (qs_kind_set, atomic_kind_set)
    1052              : 
    1053              :       ! get stuff
    1054              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
    1055          574 :                       natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
    1056              : 
    1057          574 :       IF (do_kpoints_prv) THEN
    1058           82 :          nimg = dft_control%nimages
    1059           82 :          ncell_RI = nimg
    1060           82 :          IF (do_hfx_kpoints_prv) THEN
    1061           82 :             nimg = SIZE(t3c_der_k, 1)
    1062           82 :             ncell_RI = SIZE(t3c_der_k, 2)
    1063              :          END IF
    1064           82 :          CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
    1065              :       ELSE
    1066              :          nimg = 1
    1067              :          ncell_RI = 1
    1068              :       END IF
    1069              : 
    1070          574 :       IF (do_hfx_kpoints_prv) THEN
    1071           82 :          CPASSERT(op_pos_prv == 2)
    1072              :       ELSE
    1073         1968 :          CPASSERT(ALL(SHAPE(t3c_der_k) == [nimg, ncell_RI, 3]))
    1074              :       END IF
    1075              : 
    1076         8910 :       ALLOCATE (t3c_template(nimg, ncell_RI))
    1077         1148 :       DO j_img = 1, ncell_RI
    1078         3744 :          DO i_img = 1, nimg
    1079         3170 :             CALL dbt_create(t3c_der_i(i_img, j_img, 1), t3c_template(i_img, j_img))
    1080              :          END DO
    1081              :       END DO
    1082              : 
    1083              :       CALL alloc_block_3c(t3c_template, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, &
    1084              :                           op_pos=op_pos_prv, do_kpoints=do_kpoints, do_hfx_kpoints=do_hfx_kpoints, &
    1085              :                           bounds_i=bounds_i, bounds_j=bounds_j, bounds_k=bounds_k, &
    1086         1066 :                           RI_range=RI_range, img_to_RI_cell=img_to_RI_cell, cell_to_index=cell_to_index)
    1087         2296 :       DO i_xyz = 1, 3
    1088         4018 :          DO j_img = 1, ncell_RI
    1089        11232 :             DO i_img = 1, nimg
    1090         7788 :                CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_i(i_img, j_img, i_xyz))
    1091         9510 :                CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_k(i_img, j_img, i_xyz))
    1092              :             END DO
    1093              :          END DO
    1094              :       END DO
    1095              : 
    1096         1148 :       DO j_img = 1, ncell_RI
    1097         3744 :          DO i_img = 1, nimg
    1098         3170 :             CALL dbt_destroy(t3c_template(i_img, j_img))
    1099              :          END DO
    1100              :       END DO
    1101         3170 :       DEALLOCATE (t3c_template)
    1102              : 
    1103          574 :       IF (nl_3c%sym == symmetric_jk) THEN
    1104         9840 :          ALLOCATE (t3c_der_j(nimg, ncell_RI, 3))
    1105         1968 :          DO i_xyz = 1, 3
    1106         3444 :             DO j_img = 1, ncell_RI
    1107         4428 :                DO i_img = 1, nimg
    1108         1476 :                   CALL dbt_create(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
    1109         2952 :                   CALL dbt_copy(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
    1110              :                END DO
    1111              :             END DO
    1112              :          END DO
    1113              :       END IF
    1114              : 
    1115              :       !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
    1116          574 :       nbasis = SIZE(basis_i)
    1117          574 :       max_nsgfi = 0
    1118          574 :       max_nset = 0
    1119          574 :       maxli = 0
    1120         1638 :       DO ibasis = 1, nbasis
    1121              :          CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
    1122         1064 :                                 nset=iset, nsgf_set=nsgfi, npgf=npgfi)
    1123         1064 :          maxli = MAX(maxli, imax)
    1124         1064 :          max_nset = MAX(max_nset, iset)
    1125         5792 :          max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
    1126              :       END DO
    1127              :       max_ncoj = 0
    1128              :       maxlj = 0
    1129         1638 :       DO ibasis = 1, nbasis
    1130              :          CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
    1131         1064 :                                 nset=jset, nsgf_set=nsgfj, npgf=npgfj)
    1132         1064 :          maxlj = MAX(maxlj, imax)
    1133         1064 :          max_nset = MAX(max_nset, jset)
    1134         4484 :          max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
    1135              :       END DO
    1136              :       maxlk = 0
    1137         1638 :       DO ibasis = 1, nbasis
    1138              :          CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
    1139         1064 :                                 nset=kset, nsgf_set=nsgfk, npgf=npgfk)
    1140         1064 :          maxlk = MAX(maxlk, imax)
    1141         1638 :          max_nset = MAX(max_nset, kset)
    1142              :       END DO
    1143          574 :       m_max = maxli + maxlj + maxlk + 1
    1144              : 
    1145              :       !To minimize expensive memory ops and generally optimize contraction, pre-allocate
    1146              :       !contiguous sphi arrays (and transposed in the case of sphi_i)
    1147              : 
    1148          574 :       NULLIFY (tspj, spi, spk)
    1149        26756 :       ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
    1150              : 
    1151         1638 :       DO ibasis = 1, nbasis
    1152         7962 :          DO iset = 1, max_nset
    1153         6324 :             NULLIFY (spi(iset, ibasis)%array)
    1154         6324 :             NULLIFY (tspj(iset, ibasis)%array)
    1155         7388 :             NULLIFY (spk(iset, ibasis)%array)
    1156              :          END DO
    1157              :       END DO
    1158              : 
    1159         2296 :       DO ilist = 1, 3
    1160         5488 :          DO ibasis = 1, nbasis
    1161         3192 :             IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
    1162         3192 :             IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
    1163         3192 :             IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
    1164              : 
    1165        14608 :             DO iset = 1, basis_set%nset
    1166              : 
    1167         9694 :                ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
    1168         9694 :                sgfi = basis_set%first_sgf(1, iset)
    1169         9694 :                egfi = sgfi + basis_set%nsgf_set(iset) - 1
    1170              : 
    1171        12886 :                IF (ilist == 1) THEN
    1172        16616 :                   ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
    1173      3075498 :                   spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
    1174              : 
    1175         5540 :                ELSE IF (ilist == 2) THEN
    1176        11384 :                   ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
    1177       177774 :                   tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
    1178              : 
    1179              :                ELSE
    1180        10776 :                   ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
    1181       838246 :                   spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
    1182              :                END IF
    1183              : 
    1184              :             END DO !iset
    1185              :          END DO !ibasis
    1186              :       END DO !ilist
    1187              : 
    1188              :       !Init the truncated Coulomb operator
    1189              :       IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated .OR. &
    1190          574 :           op_ij == do_potential_mix_cl_trunc .OR. op_jk == do_potential_mix_cl_trunc) THEN
    1191              : 
    1192           58 :          IF (m_max > get_lmax_init()) THEN
    1193            8 :             IF (para_env%mepos == 0) THEN
    1194            4 :                CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
    1195              :             END IF
    1196            8 :             CALL init(m_max, unit_id, para_env%mepos, para_env)
    1197            8 :             IF (para_env%mepos == 0) THEN
    1198            4 :                CALL close_file(unit_id)
    1199              :             END IF
    1200              :          END IF
    1201              :       END IF
    1202              : 
    1203          574 :       CALL init_md_ftable(nmax=m_max)
    1204              : 
    1205          574 :       IF (do_kpoints_prv) THEN
    1206          328 :          kp_index_lbounds = LBOUND(cell_to_index)
    1207          328 :          kp_index_ubounds = UBOUND(cell_to_index)
    1208              :       END IF
    1209              : 
    1210          574 :       nthread = 1
    1211          574 : !$    nthread = omp_get_max_threads()
    1212              : 
    1213              : !$OMP PARALLEL DEFAULT(NONE) &
    1214              : !$OMP SHARED (nthread,do_kpoints_prv,kp_index_lbounds,kp_index_ubounds,maxli,maxlk,maxlj,bounds_i,&
    1215              : !$OMP         bounds_j,bounds_k,nimg,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,op_pos_prv,&
    1216              : !$OMP         potential_parameter,der_eps,tspj,spi,spk,cell_to_index,RI_range,natom,nl_3c,&
    1217              : !$OMP         t3c_der_i,t3c_der_k,t3c_der_j,ncell_RI,img_to_RI_cell_prv,do_hfx_kpoints_prv) &
    1218              : !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,cell_j,cell_k,&
    1219              : !$OMP          prefac,jcell,kcell,first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
    1220              : !$OMP          sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
    1221              : !$OMP          set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
    1222              : !$OMP          rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
    1223              : !$OMP          sgfk,dijk_j,dijk_k,ri,rj,rk,max_contraction_i,max_contraction_j,blk_idx,&
    1224              : !$OMP          max_contraction_k,iset,jset,kset,block_t_i,blk_size,dijk_contr,cpp_buffer,ccp_buffer,&
    1225              : !$OMP          block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
    1226              : !$OMP          dummy_block_t,sp,handle2,mepos,bo,block_t_k,der_ext_i,der_ext_j,der_ext_k,tmp_ijk_i,&
    1227          574 : !$OMP          block_k_not_zero,der_k_zero,skip,der_j_zero,block_t_j,block_j_not_zero,tmp_ijk_j,i)
    1228              : 
    1229              :       mepos = 0
    1230              : !$    mepos = omp_get_thread_num()
    1231              : 
    1232              :       CALL cp_libint_init_3eri1(lib, MAX(maxli, maxlj, maxlk))
    1233              :       CALL cp_libint_set_contrdepth(lib, 1)
    1234              :       CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
    1235              : 
    1236              :       !We split the provided bounds among the threads such that each threads works on a different set of atoms
    1237              :       IF (PRESENT(bounds_i)) THEN
    1238              :          bo = get_limit(bounds_i(2) - bounds_i(1) + 1, nthread, mepos)
    1239              :          bo(:) = bo(:) + bounds_i(1) - 1
    1240              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
    1241              :       ELSE IF (PRESENT(bounds_j)) THEN
    1242              :          bo = get_limit(bounds_j(2) - bounds_j(1) + 1, nthread, mepos)
    1243              :          bo(:) = bo(:) + bounds_j(1) - 1
    1244              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bo, bounds_k)
    1245              :       ELSE IF (PRESENT(bounds_k)) THEN
    1246              :          bo = get_limit(bounds_k(2) - bounds_k(1) + 1, nthread, mepos)
    1247              :          bo(:) = bo(:) + bounds_k(1) - 1
    1248              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bo)
    1249              :       ELSE
    1250              :          bo = get_limit(natom, nthread, mepos)
    1251              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
    1252              :       END IF
    1253              : 
    1254              :       skip = .FALSE.
    1255              :       IF (bo(1) > bo(2)) skip = .TRUE.
    1256              : 
    1257              :       DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
    1258              :          CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
    1259              :                                    iatom=iatom, jatom=jatom, katom=katom, &
    1260              :                                    rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
    1261              :          IF (skip) EXIT
    1262              : 
    1263              :          djk = NORM2(rjk)
    1264              :          dij = NORM2(rij)
    1265              :          dik = NORM2(rik)
    1266              : 
    1267              :          IF (do_kpoints_prv) THEN
    1268              :             prefac = 0.5_dp
    1269              :          ELSEIF (nl_3c%sym == symmetric_jk) THEN
    1270              :             IF (jatom == katom) THEN
    1271              :                prefac = 0.5_dp
    1272              :             ELSE
    1273              :                prefac = 1.0_dp
    1274              :             END IF
    1275              :          ELSE
    1276              :             prefac = 1.0_dp
    1277              :          END IF
    1278              :          IF (do_hfx_kpoints_prv) prefac = 1.0_dp
    1279              : 
    1280              :          IF (do_kpoints_prv) THEN
    1281              : 
    1282              :             IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
    1283              :                 ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
    1284              : 
    1285              :             jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
    1286              :             IF (jcell > nimg .OR. jcell < 1) CYCLE
    1287              : 
    1288              :             IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
    1289              :                 ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
    1290              : 
    1291              :             kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
    1292              :             IF (kcell > nimg .OR. kcell < 1) CYCLE
    1293              :             !In case of HFX k-points, we only consider P^k if d_ik <= RI_range
    1294              :             IF (do_hfx_kpoints_prv) THEN
    1295              :                IF (dik > RI_range) CYCLE
    1296              :                kcell = img_to_RI_cell_prv(kcell)
    1297              :             END IF
    1298              :          ELSE
    1299              :             jcell = 1; kcell = 1
    1300              :          END IF
    1301              : 
    1302              :          blk_idx = [iatom, jatom, katom]
    1303              :          IF (do_hfx_kpoints_prv) THEN
    1304              :             blk_idx(3) = (kcell - 1)*natom + katom
    1305              :             kcell = 1
    1306              :          END IF
    1307              : 
    1308              :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
    1309              :                                 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
    1310              :                                 sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
    1311              : 
    1312              :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
    1313              :                                 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
    1314              :                                 sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
    1315              : 
    1316              :          CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
    1317              :                                 npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
    1318              :                                 sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
    1319              : 
    1320              :          IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
    1321              :          IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
    1322              :          IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
    1323              : 
    1324              :          ALLOCATE (max_contraction_i(nseti))
    1325              :          max_contraction_i = 0.0_dp
    1326              :          DO iset = 1, nseti
    1327              :             sgfi = first_sgf_i(1, iset)
    1328              :             max_contraction_i(iset) = MAXVAL((/(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/))
    1329              :          END DO
    1330              : 
    1331              :          ALLOCATE (max_contraction_j(nsetj))
    1332              :          max_contraction_j = 0.0_dp
    1333              :          DO jset = 1, nsetj
    1334              :             sgfj = first_sgf_j(1, jset)
    1335              :             max_contraction_j(jset) = MAXVAL((/(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/))
    1336              :          END DO
    1337              : 
    1338              :          ALLOCATE (max_contraction_k(nsetk))
    1339              :          max_contraction_k = 0.0_dp
    1340              :          DO kset = 1, nsetk
    1341              :             sgfk = first_sgf_k(1, kset)
    1342              :             max_contraction_k(kset) = MAXVAL((/(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/))
    1343              :          END DO
    1344              : 
    1345              :          CALL dbt_blk_sizes(t3c_der_i(jcell, kcell, 1), blk_idx, blk_size)
    1346              : 
    1347              :          ALLOCATE (block_t_i(blk_size(2), blk_size(3), blk_size(1), 3))
    1348              :          ALLOCATE (block_t_j(blk_size(2), blk_size(3), blk_size(1), 3))
    1349              :          ALLOCATE (block_t_k(blk_size(2), blk_size(3), blk_size(1), 3))
    1350              : 
    1351              :          block_t_i = 0.0_dp
    1352              :          block_t_j = 0.0_dp
    1353              :          block_t_k = 0.0_dp
    1354              :          block_j_not_zero = .FALSE.
    1355              :          block_k_not_zero = .FALSE.
    1356              : 
    1357              :          DO iset = 1, nseti
    1358              : 
    1359              :             DO jset = 1, nsetj
    1360              : 
    1361              :                IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE
    1362              : 
    1363              :                DO kset = 1, nsetk
    1364              : 
    1365              :                   IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE
    1366              :                   IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE
    1367              : 
    1368              :                   ncoi = npgfi(iset)*ncoset(lmax_i(iset))
    1369              :                   ncoj = npgfj(jset)*ncoset(lmax_j(jset))
    1370              :                   ncok = npgfk(kset)*ncoset(lmax_k(kset))
    1371              : 
    1372              :                   sgfi = first_sgf_i(1, iset)
    1373              :                   sgfj = first_sgf_j(1, jset)
    1374              :                   sgfk = first_sgf_k(1, kset)
    1375              : 
    1376              :                   IF (ncoj*ncok*ncoi > 0) THEN
    1377              :                      ALLOCATE (dijk_j(ncoj, ncok, ncoi, 3))
    1378              :                      ALLOCATE (dijk_k(ncoj, ncok, ncoi, 3))
    1379              :                      dijk_j(:, :, :, :) = 0.0_dp
    1380              :                      dijk_k(:, :, :, :) = 0.0_dp
    1381              : 
    1382              :                      der_j_zero = .FALSE.
    1383              :                      der_k_zero = .FALSE.
    1384              : 
    1385              :                      !need positions for libint. Only relative positions are needed => set ri to 0.0
    1386              :                      ri = 0.0_dp
    1387              :                      rj = rij ! ri + rij
    1388              :                      rk = rik ! ri + rik
    1389              : 
    1390              :                      IF (op_pos_prv == 1) THEN
    1391              :                         CALL eri_3center_derivs(dijk_j, dijk_k, &
    1392              :                                                 lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
    1393              :                                                 lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
    1394              :                                                 lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
    1395              :                                                 djk, dij, dik, lib, potential_parameter, &
    1396              :                                                 der_abc_1_ext=der_ext_j, der_abc_2_ext=der_ext_k)
    1397              :                      ELSE
    1398              :                         ALLOCATE (tmp_ijk_i(ncoi, ncoj, ncok, 3))
    1399              :                         ALLOCATE (tmp_ijk_j(ncoi, ncoj, ncok, 3))
    1400              :                         tmp_ijk_i(:, :, :, :) = 0.0_dp
    1401              :                         tmp_ijk_j(:, :, :, :) = 0.0_dp
    1402              : 
    1403              :                         CALL eri_3center_derivs(tmp_ijk_i, tmp_ijk_j, &
    1404              :                                                 lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
    1405              :                                                 lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
    1406              :                                                 lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
    1407              :                                                 dij, dik, djk, lib, potential_parameter, &
    1408              :                                                 der_abc_1_ext=der_ext_i, der_abc_2_ext=der_ext_j)
    1409              : 
    1410              :                         !TODO: is that inefficient?
    1411              :                         der_ext_k = 0
    1412              :                         DO i_xyz = 1, 3
    1413              :                            DO i = 1, ncoi
    1414              :                               dijk_j(:, :, i, i_xyz) = tmp_ijk_j(i, :, :, i_xyz)
    1415              :                               dijk_k(:, :, i, i_xyz) = -(dijk_j(:, :, i, i_xyz) + tmp_ijk_i(i, :, :, i_xyz))
    1416              :                               der_ext_k(i_xyz) = MAX(der_ext_k(i_xyz), MAXVAL(ABS(dijk_k(:, :, i, i_xyz))))
    1417              :                            END DO
    1418              :                         END DO
    1419              :                         DEALLOCATE (tmp_ijk_i, tmp_ijk_j)
    1420              : 
    1421              :                      END IF
    1422              : 
    1423              :                      IF (PRESENT(der_eps)) THEN
    1424              :                         DO i_xyz = 1, 3
    1425              :                            IF (der_eps > der_ext_j(i_xyz)*(max_contraction_i(iset)* &
    1426              :                                                            max_contraction_j(jset)* &
    1427              :                                                            max_contraction_k(kset))) THEN
    1428              :                               der_j_zero(i_xyz) = .TRUE.
    1429              :                            END IF
    1430              :                         END DO
    1431              : 
    1432              :                         DO i_xyz = 1, 3
    1433              :                            IF (der_eps > der_ext_k(i_xyz)*(max_contraction_i(iset)* &
    1434              :                                                            max_contraction_j(jset)* &
    1435              :                                                            max_contraction_k(kset))) THEN
    1436              :                               der_k_zero(i_xyz) = .TRUE.
    1437              :                            END IF
    1438              :                         END DO
    1439              :                         IF (ALL(der_j_zero) .AND. ALL(der_k_zero)) THEN
    1440              :                            DEALLOCATE (dijk_j, dijk_k)
    1441              :                            CYCLE
    1442              :                         END IF
    1443              :                      END IF
    1444              : 
    1445              :                      ALLOCATE (dijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
    1446              : 
    1447              :                      block_start_j = sgfj
    1448              :                      block_end_j = sgfj + nsgfj(jset) - 1
    1449              :                      block_start_k = sgfk
    1450              :                      block_end_k = sgfk + nsgfk(kset) - 1
    1451              :                      block_start_i = sgfi
    1452              :                      block_end_i = sgfi + nsgfi(iset) - 1
    1453              : 
    1454              :                      DO i_xyz = 1, 3
    1455              :                         IF (der_j_zero(i_xyz)) CYCLE
    1456              : 
    1457              :                         block_j_not_zero(i_xyz) = .TRUE.
    1458              :                         CALL abc_contract_xsmm(dijk_contr, dijk_j(:, :, :, i_xyz), tspj(jset, jkind)%array, &
    1459              :                                                spk(kset, kkind)%array, spi(iset, ikind)%array, &
    1460              :                                                ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
    1461              :                                                nsgfi(iset), cpp_buffer, ccp_buffer, prefac)
    1462              : 
    1463              :                         block_t_j(block_start_j:block_end_j, &
    1464              :                                   block_start_k:block_end_k, &
    1465              :                                   block_start_i:block_end_i, i_xyz) = &
    1466              :                            block_t_j(block_start_j:block_end_j, &
    1467              :                                      block_start_k:block_end_k, &
    1468              :                                      block_start_i:block_end_i, i_xyz) + &
    1469              :                            dijk_contr(:, :, :)
    1470              : 
    1471              :                      END DO
    1472              : 
    1473              :                      DO i_xyz = 1, 3
    1474              :                         IF (der_k_zero(i_xyz)) CYCLE
    1475              : 
    1476              :                         block_k_not_zero(i_xyz) = .TRUE.
    1477              :                         CALL abc_contract_xsmm(dijk_contr, dijk_k(:, :, :, i_xyz), tspj(jset, jkind)%array, &
    1478              :                                                spk(kset, kkind)%array, spi(iset, ikind)%array, &
    1479              :                                                ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
    1480              :                                                nsgfi(iset), cpp_buffer, ccp_buffer, prefac)
    1481              : 
    1482              :                         block_t_k(block_start_j:block_end_j, &
    1483              :                                   block_start_k:block_end_k, &
    1484              :                                   block_start_i:block_end_i, i_xyz) = &
    1485              :                            block_t_k(block_start_j:block_end_j, &
    1486              :                                      block_start_k:block_end_k, &
    1487              :                                      block_start_i:block_end_i, i_xyz) + &
    1488              :                            dijk_contr(:, :, :)
    1489              : 
    1490              :                      END DO
    1491              : 
    1492              :                      DEALLOCATE (dijk_j, dijk_k, dijk_contr)
    1493              :                   END IF ! number of triples > 0
    1494              :                END DO
    1495              :             END DO
    1496              :          END DO
    1497              : 
    1498              :          CALL timeset(routineN//"_put_dbcsr", handle2)
    1499              : !$OMP CRITICAL
    1500              :          sp = SHAPE(block_t_i(:, :, :, 1))
    1501              :          sp([2, 3, 1]) = sp
    1502              : 
    1503              :          DO i_xyz = 1, 3
    1504              :             !Derivatives wrt to center i can be obtained by translational invariance
    1505              :             IF ((.NOT. block_j_not_zero(i_xyz)) .AND. (.NOT. block_k_not_zero(i_xyz))) CYCLE
    1506              :             block_t_i(:, :, :, i_xyz) = -(block_t_j(:, :, :, i_xyz) + block_t_k(:, :, :, i_xyz))
    1507              : 
    1508              :             CALL dbt_put_block(t3c_der_i(jcell, kcell, i_xyz), blk_idx, sp, &
    1509              :                                RESHAPE(block_t_i(:, :, :, i_xyz), SHAPE=sp, ORDER=[2, 3, 1]), &
    1510              :                                summation=.TRUE.)
    1511              :          END DO
    1512              : 
    1513              :          sp = SHAPE(block_t_k(:, :, :, 1))
    1514              :          sp([2, 3, 1]) = sp
    1515              : 
    1516              :          DO i_xyz = 1, 3
    1517              :             IF (.NOT. block_k_not_zero(i_xyz)) CYCLE
    1518              :             CALL dbt_put_block(t3c_der_k(jcell, kcell, i_xyz), blk_idx, sp, &
    1519              :                                RESHAPE(block_t_k(:, :, :, i_xyz), SHAPE=sp, ORDER=[2, 3, 1]), &
    1520              :                                summation=.TRUE.)
    1521              :          END DO
    1522              : !$OMP END CRITICAL
    1523              : 
    1524              :          IF (nl_3c%sym == symmetric_jk) THEN
    1525              :             sp = SHAPE(block_t_j(:, :, :, 1))
    1526              :             sp([2, 3, 1]) = sp
    1527              : !$OMP CRITICAL
    1528              :             DO i_xyz = 1, 3
    1529              :                IF (.NOT. block_j_not_zero(i_xyz)) CYCLE
    1530              :                CALL dbt_put_block(t3c_der_j(jcell, kcell, i_xyz), blk_idx, sp, &
    1531              :                                   RESHAPE(block_t_j(:, :, :, i_xyz), SHAPE=sp, ORDER=[2, 3, 1]), &
    1532              :                                   summation=.TRUE.)
    1533              :             END DO
    1534              : !$OMP END CRITICAL
    1535              :          END IF
    1536              : 
    1537              :          CALL timestop(handle2)
    1538              : 
    1539              :          DEALLOCATE (block_t_i, block_t_j, block_t_k)
    1540              :          DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k)
    1541              :       END DO
    1542              : 
    1543              :       IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
    1544              :       IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
    1545              : 
    1546              :       CALL cp_libint_cleanup_3eri1(lib)
    1547              :       CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
    1548              : !$OMP END PARALLEL
    1549              : 
    1550          574 :       IF (do_kpoints_prv .AND. .NOT. do_hfx_kpoints_prv) THEN
    1551            0 :          DO i_xyz = 1, 3
    1552            0 :             DO kcell = 1, nimg
    1553            0 :                DO jcell = 1, nimg
    1554              :                   ! need half of filter eps because afterwards we add transposed tensor
    1555            0 :                   CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps/2)
    1556            0 :                   CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps/2)
    1557              :                END DO
    1558              :             END DO
    1559              :          END DO
    1560              : 
    1561          574 :       ELSEIF (nl_3c%sym == symmetric_jk) THEN
    1562              :          !Add the transpose of t3c_der_j to t3c_der_k to get the fully populated tensor
    1563          492 :          CALL dbt_create(t3c_der_k(1, 1, 1), t3c_tmp)
    1564         1968 :          DO i_xyz = 1, 3
    1565         3444 :             DO kcell = 1, nimg
    1566         4428 :                DO jcell = 1, nimg
    1567              :                   CALL dbt_copy(t3c_der_j(jcell, kcell, i_xyz), t3c_der_k(jcell, kcell, i_xyz), &
    1568         1476 :                                 order=[1, 3, 2], move_data=.TRUE., summation=.TRUE.)
    1569         1476 :                   CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
    1570              : 
    1571         1476 :                   CALL dbt_copy(t3c_der_i(jcell, kcell, i_xyz), t3c_tmp)
    1572              :                   CALL dbt_copy(t3c_tmp, t3c_der_i(jcell, kcell, i_xyz), &
    1573         1476 :                                 order=[1, 3, 2], move_data=.TRUE., summation=.TRUE.)
    1574         2952 :                   CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
    1575              :                END DO
    1576              :             END DO
    1577              :          END DO
    1578          492 :          CALL dbt_destroy(t3c_tmp)
    1579              : 
    1580           82 :       ELSEIF (nl_3c%sym == symmetric_none) THEN
    1581          328 :          DO i_xyz = 1, 3
    1582          574 :             DO kcell = 1, ncell_RI
    1583         6804 :                DO jcell = 1, nimg
    1584         6312 :                   CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
    1585         6558 :                   CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
    1586              :                END DO
    1587              :             END DO
    1588              :          END DO
    1589              :       ELSE
    1590            0 :          CPABORT("requested symmetric case not implemented")
    1591              :       END IF
    1592              : 
    1593          574 :       IF (nl_3c%sym == symmetric_jk) THEN
    1594         1968 :          DO i_xyz = 1, 3
    1595         3444 :             DO j_img = 1, nimg
    1596         4428 :                DO i_img = 1, nimg
    1597         2952 :                   CALL dbt_destroy(t3c_der_j(i_img, j_img, i_xyz))
    1598              :                END DO
    1599              :             END DO
    1600              :          END DO
    1601              :       END IF
    1602              : 
    1603         3934 :       DO iset = 1, max_nset
    1604        10258 :          DO ibasis = 1, nbasis
    1605         6324 :             IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
    1606         6324 :             IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
    1607         9684 :             IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
    1608              :          END DO
    1609              :       END DO
    1610              : 
    1611          574 :       DEALLOCATE (spi, tspj, spk)
    1612              : 
    1613          574 :       CALL timestop(handle)
    1614         6642 :    END SUBROUTINE build_3c_derivatives
    1615              : 
    1616              : ! **************************************************************************************************
    1617              : !> \brief Calculates the 3c virial contributions on the fly
    1618              : !> \param work_virial ...
    1619              : !> \param t3c_trace the tensor with which the trace should be taken
    1620              : !> \param pref ...
    1621              : !> \param qs_env ...
    1622              : !> \param nl_3c 3-center neighborlist, with a distribution matching that of t3c_trace
    1623              : !> \param basis_i ...
    1624              : !> \param basis_j ...
    1625              : !> \param basis_k ...
    1626              : !> \param potential_parameter ...
    1627              : !> \param der_eps neglect integrals smaller than der_eps
    1628              : !> \param op_pos operator position.
    1629              : !>        1: calculate (i|jk) integrals,
    1630              : !>        2: calculate (ij|k) integrals
    1631              : !> this routine requires that libint has been static initialised somewhere else
    1632              : ! **************************************************************************************************
    1633           16 :    SUBROUTINE calc_3c_virial(work_virial, t3c_trace, pref, qs_env, &
    1634           16 :                              nl_3c, basis_i, basis_j, basis_k, &
    1635              :                              potential_parameter, der_eps, op_pos)
    1636              : 
    1637              :       REAL(dp), DIMENSION(3, 3), INTENT(INOUT)           :: work_virial
    1638              :       TYPE(dbt_type), INTENT(INOUT)                      :: t3c_trace
    1639              :       REAL(KIND=dp), INTENT(IN)                          :: pref
    1640              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1641              :       TYPE(neighbor_list_3c_type), INTENT(INOUT)         :: nl_3c
    1642              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j, basis_k
    1643              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
    1644              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: der_eps
    1645              :       INTEGER, INTENT(IN), OPTIONAL                      :: op_pos
    1646              : 
    1647              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'calc_3c_virial'
    1648              : 
    1649              :       INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
    1650              :          block_start_k, egfi, handle, i, i_xyz, iatom, ibasis, ikind, ilist, imax, iset, j_xyz, &
    1651              :          jatom, jkind, jset, katom, kkind, kset, m_max, max_ncoj, max_nset, max_nsgfi, maxli, &
    1652              :          maxlj, maxlk, mepos, natom, nbasis, ncoi, ncoj, ncok, nseti, nsetj, nsetk, nthread, &
    1653              :          op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id
    1654              :       INTEGER, DIMENSION(2)                              :: bo
    1655              :       INTEGER, DIMENSION(3)                              :: blk_size, sp
    1656           16 :       INTEGER, DIMENSION(:), POINTER                     :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
    1657           32 :                                                             lmin_k, npgfi, npgfj, npgfk, nsgfi, &
    1658           16 :                                                             nsgfj, nsgfk
    1659           16 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf_i, first_sgf_j, first_sgf_k
    1660              :       LOGICAL                                            :: found, skip
    1661              :       LOGICAL, DIMENSION(3)                              :: block_j_not_zero, block_k_not_zero, &
    1662              :                                                             der_j_zero, der_k_zero
    1663              :       REAL(dp)                                           :: force
    1664              :       REAL(dp), DIMENSION(3)                             :: der_ext_i, der_ext_j, der_ext_k
    1665              :       REAL(KIND=dp)                                      :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
    1666              :                                                             kind_radius_i, kind_radius_j, &
    1667              :                                                             kind_radius_k
    1668           16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ccp_buffer, cpp_buffer, &
    1669           16 :                                                             max_contraction_i, max_contraction_j, &
    1670           16 :                                                             max_contraction_k
    1671           32 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: ablock, dijk_contr, tmp_block
    1672           16 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: block_t_i, block_t_j, block_t_k, dijk_j, &
    1673           16 :                                                             dijk_k
    1674              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rij, rik, rj, rjk, rk, scoord
    1675           16 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_i, set_radius_j, set_radius_k
    1676           16 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
    1677           16 :                                                             sphi_k, zeti, zetj, zetk
    1678           16 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1679              :       TYPE(cell_type), POINTER                           :: cell
    1680           16 :       TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER     :: spi, spk, tspj
    1681              :       TYPE(cp_libint_t)                                  :: lib
    1682              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1683              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
    1684              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1685              :       TYPE(neighbor_list_3c_iterator_type)               :: nl_3c_iter
    1686           16 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1687           16 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1688              : 
    1689           16 :       CALL timeset(routineN, handle)
    1690              : 
    1691           16 :       op_ij = do_potential_id; op_jk = do_potential_id
    1692              : 
    1693           16 :       IF (PRESENT(op_pos)) THEN
    1694           16 :          op_pos_prv = op_pos
    1695              :       ELSE
    1696              :          op_pos_prv = 1
    1697              :       END IF
    1698           16 :       CPASSERT(op_pos == 1)
    1699           16 :       CPASSERT(.NOT. nl_3c%sym == symmetric_jk)
    1700              : 
    1701           16 :       SELECT CASE (op_pos_prv)
    1702              :       CASE (1)
    1703           16 :          op_ij = potential_parameter%potential_type
    1704              :       CASE (2)
    1705           16 :          op_jk = potential_parameter%potential_type
    1706              :       END SELECT
    1707              : 
    1708           16 :       dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
    1709              : 
    1710              :       IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
    1711           16 :           .OR. op_ij == do_potential_mix_cl_trunc) THEN
    1712            6 :          dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
    1713            6 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
    1714           10 :       ELSEIF (op_ij == do_potential_coulomb) THEN
    1715            0 :          dr_ij = 1000000.0_dp
    1716            0 :          dr_ik = 1000000.0_dp
    1717              :       END IF
    1718              : 
    1719              :       IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
    1720           16 :           .OR. op_jk == do_potential_mix_cl_trunc) THEN
    1721            0 :          dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
    1722            0 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
    1723           16 :       ELSEIF (op_jk == do_potential_coulomb) THEN
    1724            0 :          dr_jk = 1000000.0_dp
    1725            0 :          dr_ik = 1000000.0_dp
    1726              :       END IF
    1727              : 
    1728           16 :       NULLIFY (qs_kind_set, atomic_kind_set)
    1729              : 
    1730              :       ! get stuff
    1731              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
    1732              :                       natom=natom, dft_control=dft_control, para_env=para_env, &
    1733           16 :                       particle_set=particle_set, cell=cell)
    1734              : 
    1735              :       !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
    1736           16 :       nbasis = SIZE(basis_i)
    1737           16 :       max_nsgfi = 0
    1738           16 :       max_nset = 0
    1739           16 :       maxli = 0
    1740           48 :       DO ibasis = 1, nbasis
    1741              :          CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
    1742           32 :                                 nset=iset, nsgf_set=nsgfi, npgf=npgfi)
    1743           32 :          maxli = MAX(maxli, imax)
    1744           32 :          max_nset = MAX(max_nset, iset)
    1745          182 :          max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
    1746              :       END DO
    1747              :       max_ncoj = 0
    1748              :       maxlj = 0
    1749           48 :       DO ibasis = 1, nbasis
    1750              :          CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
    1751           32 :                                 nset=jset, nsgf_set=nsgfj, npgf=npgfj)
    1752           32 :          maxlj = MAX(maxlj, imax)
    1753           32 :          max_nset = MAX(max_nset, jset)
    1754          138 :          max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
    1755              :       END DO
    1756              :       maxlk = 0
    1757           48 :       DO ibasis = 1, nbasis
    1758              :          CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
    1759           32 :                                 nset=kset, nsgf_set=nsgfk, npgf=npgfk)
    1760           32 :          maxlk = MAX(maxlk, imax)
    1761           48 :          max_nset = MAX(max_nset, kset)
    1762              :       END DO
    1763           16 :       m_max = maxli + maxlj + maxlk + 1
    1764              : 
    1765              :       !To minimize expensive memory opsand generally optimize contraction, pre-allocate
    1766              :       !contiguous sphi arrays (and transposed in the cas of sphi_i)
    1767              : 
    1768           16 :       NULLIFY (tspj, spi, spk)
    1769          920 :       ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
    1770              : 
    1771           48 :       DO ibasis = 1, nbasis
    1772          280 :          DO iset = 1, max_nset
    1773          232 :             NULLIFY (spi(iset, ibasis)%array)
    1774          232 :             NULLIFY (tspj(iset, ibasis)%array)
    1775              : 
    1776          264 :             NULLIFY (spk(iset, ibasis)%array)
    1777              :          END DO
    1778              :       END DO
    1779              : 
    1780           64 :       DO ilist = 1, 3
    1781          160 :          DO ibasis = 1, nbasis
    1782           96 :             IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
    1783           96 :             IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
    1784           96 :             IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
    1785              : 
    1786          458 :             DO iset = 1, basis_set%nset
    1787              : 
    1788          314 :                ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
    1789          314 :                sgfi = basis_set%first_sgf(1, iset)
    1790          314 :                egfi = sgfi + basis_set%nsgf_set(iset) - 1
    1791              : 
    1792          410 :                IF (ilist == 1) THEN
    1793          536 :                   ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
    1794       177990 :                   spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
    1795              : 
    1796          180 :                ELSE IF (ilist == 2) THEN
    1797          360 :                   ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
    1798         5194 :                   tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
    1799              : 
    1800              :                ELSE
    1801          360 :                   ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
    1802         4442 :                   spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
    1803              :                END IF
    1804              : 
    1805              :             END DO !iset
    1806              :          END DO !ibasis
    1807              :       END DO !ilist
    1808              : 
    1809              :       !Init the truncated Coulomb operator
    1810              :       IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated .OR. &
    1811           16 :           op_ij == do_potential_mix_cl_trunc .OR. op_jk == do_potential_mix_cl_trunc) THEN
    1812              : 
    1813            6 :          IF (m_max > get_lmax_init()) THEN
    1814            0 :             IF (para_env%mepos == 0) THEN
    1815            0 :                CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
    1816              :             END IF
    1817            0 :             CALL init(m_max, unit_id, para_env%mepos, para_env)
    1818            0 :             IF (para_env%mepos == 0) THEN
    1819            0 :                CALL close_file(unit_id)
    1820              :             END IF
    1821              :          END IF
    1822              :       END IF
    1823              : 
    1824           16 :       CALL init_md_ftable(nmax=m_max)
    1825              : 
    1826           16 :       nthread = 1
    1827           16 : !$    nthread = omp_get_max_threads()
    1828              : 
    1829              : !$OMP PARALLEL DEFAULT(NONE) &
    1830              : !$OMP SHARED (nthread,maxli,maxlk,maxlj,work_virial,pref,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,&
    1831              : !$OMP         ncoset,potential_parameter,der_eps,tspj,spi,spk,natom,nl_3c,t3c_trace,cell,particle_set) &
    1832              : !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,i_xyz,j_xyz,first_sgf_i,&
    1833              : !$OMP          lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,sphi_i,zeti,kind_radius_i,first_sgf_j,&
    1834              : !$OMP          lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,&
    1835              : !$OMP          lmax_k,lmin_k,npgfk,nsetk,nsgfk,rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,&
    1836              : !$OMP          ncoi,ncoj,ncok,sgfi,sgfj,sgfk,dijk_j,dijk_k,ri,rj,rk,max_contraction_i,max_contraction_j,&
    1837              : !$OMP          tmp_block,max_contraction_k,iset,jset,kset,block_t_i,blk_size,dijk_contr,found,sp,mepos,&
    1838              : !$OMP          block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,block_t_k,&
    1839              : !$OMP          bo,der_ext_i,der_ext_j,der_ext_k,ablock,force,scoord,skip,cpp_buffer,ccp_buffer,&
    1840           16 : !$OMP          block_k_not_zero,der_k_zero,der_j_zero,block_t_j,block_j_not_zero,i)
    1841              : 
    1842              :       mepos = 0
    1843              : !$    mepos = omp_get_thread_num()
    1844              : 
    1845              :       CALL cp_libint_init_3eri1(lib, MAX(maxli, maxlj, maxlk))
    1846              :       CALL cp_libint_set_contrdepth(lib, 1)
    1847              :       CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
    1848              : 
    1849              :       !We split the provided bounds among the threads such that each threads works on a different set of atoms
    1850              : 
    1851              :       bo = get_limit(natom, nthread, mepos)
    1852              :       CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i=bo)
    1853              : 
    1854              :       skip = .FALSE.
    1855              :       IF (bo(1) > bo(2)) skip = .TRUE.
    1856              : 
    1857              :       DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
    1858              :          CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
    1859              :                                    iatom=iatom, jatom=jatom, katom=katom, &
    1860              :                                    rij=rij, rjk=rjk, rik=rik)
    1861              :          IF (skip) EXIT
    1862              : 
    1863              :          CALL dbt_get_block(t3c_trace, [iatom, jatom, katom], tmp_block, found)
    1864              :          IF (.NOT. found) CYCLE
    1865              : 
    1866              :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
    1867              :                                 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
    1868              :                                 sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
    1869              : 
    1870              :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
    1871              :                                 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
    1872              :                                 sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
    1873              : 
    1874              :          CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
    1875              :                                 npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
    1876              :                                 sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
    1877              : 
    1878              :          djk = NORM2(rjk)
    1879              :          dij = NORM2(rij)
    1880              :          dik = NORM2(rik)
    1881              : 
    1882              :          IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
    1883              :          IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
    1884              :          IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
    1885              : 
    1886              :          ALLOCATE (max_contraction_i(nseti))
    1887              :          max_contraction_i = 0.0_dp
    1888              :          DO iset = 1, nseti
    1889              :             sgfi = first_sgf_i(1, iset)
    1890              :             max_contraction_i(iset) = MAXVAL((/(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/))
    1891              :          END DO
    1892              : 
    1893              :          ALLOCATE (max_contraction_j(nsetj))
    1894              :          max_contraction_j = 0.0_dp
    1895              :          DO jset = 1, nsetj
    1896              :             sgfj = first_sgf_j(1, jset)
    1897              :             max_contraction_j(jset) = MAXVAL((/(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/))
    1898              :          END DO
    1899              : 
    1900              :          ALLOCATE (max_contraction_k(nsetk))
    1901              :          max_contraction_k = 0.0_dp
    1902              :          DO kset = 1, nsetk
    1903              :             sgfk = first_sgf_k(1, kset)
    1904              :             max_contraction_k(kset) = MAXVAL((/(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/))
    1905              :          END DO
    1906              : 
    1907              :          CALL dbt_blk_sizes(t3c_trace, [iatom, jatom, katom], blk_size)
    1908              : 
    1909              :          ALLOCATE (block_t_i(blk_size(2), blk_size(3), blk_size(1), 3))
    1910              :          ALLOCATE (block_t_j(blk_size(2), blk_size(3), blk_size(1), 3))
    1911              :          ALLOCATE (block_t_k(blk_size(2), blk_size(3), blk_size(1), 3))
    1912              : 
    1913              :          ALLOCATE (ablock(blk_size(2), blk_size(3), blk_size(1)))
    1914              :          DO i = 1, blk_size(1)
    1915              :             ablock(:, :, i) = tmp_block(i, :, :)
    1916              :          END DO
    1917              :          DEALLOCATE (tmp_block)
    1918              : 
    1919              :          block_t_i = 0.0_dp
    1920              :          block_t_j = 0.0_dp
    1921              :          block_t_k = 0.0_dp
    1922              :          block_j_not_zero = .FALSE.
    1923              :          block_k_not_zero = .FALSE.
    1924              : 
    1925              :          DO iset = 1, nseti
    1926              : 
    1927              :             DO jset = 1, nsetj
    1928              : 
    1929              :                IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE
    1930              : 
    1931              :                DO kset = 1, nsetk
    1932              : 
    1933              :                   IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE
    1934              :                   IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE
    1935              : 
    1936              :                   ncoi = npgfi(iset)*ncoset(lmax_i(iset))
    1937              :                   ncoj = npgfj(jset)*ncoset(lmax_j(jset))
    1938              :                   ncok = npgfk(kset)*ncoset(lmax_k(kset))
    1939              : 
    1940              :                   sgfi = first_sgf_i(1, iset)
    1941              :                   sgfj = first_sgf_j(1, jset)
    1942              :                   sgfk = first_sgf_k(1, kset)
    1943              : 
    1944              :                   IF (ncoj*ncok*ncoi > 0) THEN
    1945              :                      ALLOCATE (dijk_j(ncoj, ncok, ncoi, 3))
    1946              :                      ALLOCATE (dijk_k(ncoj, ncok, ncoi, 3))
    1947              :                      dijk_j(:, :, :, :) = 0.0_dp
    1948              :                      dijk_k(:, :, :, :) = 0.0_dp
    1949              : 
    1950              :                      der_j_zero = .FALSE.
    1951              :                      der_k_zero = .FALSE.
    1952              : 
    1953              :                      !need positions for libint. Only relative positions are needed => set ri to 0.0
    1954              :                      ri = 0.0_dp
    1955              :                      rj = rij ! ri + rij
    1956              :                      rk = rik ! ri + rik
    1957              : 
    1958              :                      CALL eri_3center_derivs(dijk_j, dijk_k, &
    1959              :                                              lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
    1960              :                                              lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
    1961              :                                              lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
    1962              :                                              djk, dij, dik, lib, potential_parameter, &
    1963              :                                              der_abc_1_ext=der_ext_j, der_abc_2_ext=der_ext_k)
    1964              : 
    1965              :                      IF (PRESENT(der_eps)) THEN
    1966              :                         DO i_xyz = 1, 3
    1967              :                            IF (der_eps > der_ext_j(i_xyz)*(max_contraction_i(iset)* &
    1968              :                                                            max_contraction_j(jset)* &
    1969              :                                                            max_contraction_k(kset))) THEN
    1970              :                               der_j_zero(i_xyz) = .TRUE.
    1971              :                            END IF
    1972              :                         END DO
    1973              : 
    1974              :                         DO i_xyz = 1, 3
    1975              :                            IF (der_eps > der_ext_k(i_xyz)*(max_contraction_i(iset)* &
    1976              :                                                            max_contraction_j(jset)* &
    1977              :                                                            max_contraction_k(kset))) THEN
    1978              :                               der_k_zero(i_xyz) = .TRUE.
    1979              :                            END IF
    1980              :                         END DO
    1981              :                         IF (ALL(der_j_zero) .AND. ALL(der_k_zero)) THEN
    1982              :                            DEALLOCATE (dijk_j, dijk_k)
    1983              :                            CYCLE
    1984              :                         END IF
    1985              :                      END IF
    1986              : 
    1987              :                      ALLOCATE (dijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
    1988              : 
    1989              :                      block_start_j = sgfj
    1990              :                      block_end_j = sgfj + nsgfj(jset) - 1
    1991              :                      block_start_k = sgfk
    1992              :                      block_end_k = sgfk + nsgfk(kset) - 1
    1993              :                      block_start_i = sgfi
    1994              :                      block_end_i = sgfi + nsgfi(iset) - 1
    1995              : 
    1996              :                      DO i_xyz = 1, 3
    1997              :                         IF (der_j_zero(i_xyz)) CYCLE
    1998              : 
    1999              :                         block_j_not_zero(i_xyz) = .TRUE.
    2000              :                         CALL abc_contract_xsmm(dijk_contr, dijk_j(:, :, :, i_xyz), tspj(jset, jkind)%array, &
    2001              :                                                spk(kset, kkind)%array, spi(iset, ikind)%array, &
    2002              :                                                ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
    2003              :                                                nsgfi(iset), cpp_buffer, ccp_buffer)
    2004              : 
    2005              :                         block_t_j(block_start_j:block_end_j, &
    2006              :                                   block_start_k:block_end_k, &
    2007              :                                   block_start_i:block_end_i, i_xyz) = &
    2008              :                            block_t_j(block_start_j:block_end_j, &
    2009              :                                      block_start_k:block_end_k, &
    2010              :                                      block_start_i:block_end_i, i_xyz) + &
    2011              :                            dijk_contr(:, :, :)
    2012              : 
    2013              :                      END DO
    2014              : 
    2015              :                      DO i_xyz = 1, 3
    2016              :                         IF (der_k_zero(i_xyz)) CYCLE
    2017              : 
    2018              :                         block_k_not_zero(i_xyz) = .TRUE.
    2019              :                         CALL abc_contract_xsmm(dijk_contr, dijk_k(:, :, :, i_xyz), tspj(jset, jkind)%array, &
    2020              :                                                spk(kset, kkind)%array, spi(iset, ikind)%array, &
    2021              :                                                ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
    2022              :                                                nsgfi(iset), cpp_buffer, ccp_buffer)
    2023              : 
    2024              :                         block_t_k(block_start_j:block_end_j, &
    2025              :                                   block_start_k:block_end_k, &
    2026              :                                   block_start_i:block_end_i, i_xyz) = &
    2027              :                            block_t_k(block_start_j:block_end_j, &
    2028              :                                      block_start_k:block_end_k, &
    2029              :                                      block_start_i:block_end_i, i_xyz) + &
    2030              :                            dijk_contr(:, :, :)
    2031              : 
    2032              :                      END DO
    2033              : 
    2034              :                      DEALLOCATE (dijk_j, dijk_k, dijk_contr)
    2035              :                   END IF ! number of triples > 0
    2036              :                END DO
    2037              :             END DO
    2038              :          END DO
    2039              : 
    2040              :          !We obtain the derivative wrt to first center using translational invariance
    2041              :          DO i_xyz = 1, 3
    2042              :             block_t_i(:, :, :, i_xyz) = -block_t_j(:, :, :, i_xyz) - block_t_k(:, :, :, i_xyz)
    2043              :          END DO
    2044              : 
    2045              :          !virial contribution coming from deriv wrt to first center
    2046              :          DO i_xyz = 1, 3
    2047              :             force = pref*SUM(ablock(:, :, :)*block_t_i(:, :, :, i_xyz))
    2048              :             CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
    2049              :             DO j_xyz = 1, 3
    2050              : !$OMP ATOMIC
    2051              :                work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
    2052              :             END DO
    2053              :          END DO
    2054              : 
    2055              :          !second center
    2056              :          DO i_xyz = 1, 3
    2057              :             force = pref*SUM(ablock(:, :, :)*block_t_j(:, :, :, i_xyz))
    2058              :             CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
    2059              :             DO j_xyz = 1, 3
    2060              : !$OMP ATOMIC
    2061              :                work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
    2062              :             END DO
    2063              :          END DO
    2064              : 
    2065              :          !third center
    2066              :          DO i_xyz = 1, 3
    2067              :             force = pref*SUM(ablock(:, :, :)*block_t_k(:, :, :, i_xyz))
    2068              :             CALL real_to_scaled(scoord, particle_set(iatom)%r + rik, cell)
    2069              :             DO j_xyz = 1, 3
    2070              : !$OMP ATOMIC
    2071              :                work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
    2072              :             END DO
    2073              :          END DO
    2074              : 
    2075              :          DEALLOCATE (block_t_i, block_t_j, block_t_k)
    2076              :          DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k, ablock)
    2077              :       END DO
    2078              : 
    2079              :       IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
    2080              :       IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
    2081              : 
    2082              :       CALL cp_libint_cleanup_3eri1(lib)
    2083              :       CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
    2084              : !$OMP END PARALLEL
    2085              : 
    2086          132 :       DO iset = 1, max_nset
    2087          364 :          DO ibasis = 1, nbasis
    2088          232 :             IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
    2089          232 :             IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
    2090          348 :             IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
    2091              :          END DO
    2092              :       END DO
    2093              : 
    2094           16 :       DEALLOCATE (spi, tspj, spk)
    2095              : 
    2096           16 :       CALL timestop(handle)
    2097          128 :    END SUBROUTINE calc_3c_virial
    2098              : 
    2099              : ! **************************************************************************************************
    2100              : !> \brief Build 3-center integral tensor
    2101              : !> \param t3c empty DBCSR tensor
    2102              : !>            Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages)
    2103              : !>            if k-points are used
    2104              : !> \param filter_eps Filter threshold for tensor blocks
    2105              : !> \param qs_env ...
    2106              : !> \param nl_3c 3-center neighborlist
    2107              : !> \param basis_i ...
    2108              : !> \param basis_j ...
    2109              : !> \param basis_k ...
    2110              : !> \param potential_parameter ...
    2111              : !> \param int_eps neglect integrals smaller than int_eps
    2112              : !> \param op_pos operator position.
    2113              : !>        1: calculate (i|jk) integrals,
    2114              : !>        2: calculate (ij|k) integrals
    2115              : !> \param do_kpoints ...
    2116              : !> this routine requires that libint has been static initialised somewhere else
    2117              : !> \param do_hfx_kpoints ...
    2118              : !> \param desymmetrize ...
    2119              : !> \param cell_sym ...
    2120              : !> \param bounds_i ...
    2121              : !> \param bounds_j ...
    2122              : !> \param bounds_k ...
    2123              : !> \param RI_range ...
    2124              : !> \param img_to_RI_cell ...
    2125              : !> \param cell_to_index_ext ...
    2126              : ! **************************************************************************************************
    2127         1978 :    SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
    2128         1978 :                                  nl_3c, basis_i, basis_j, basis_k, &
    2129              :                                  potential_parameter, int_eps, &
    2130              :                                  op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, cell_sym, &
    2131              :                                  bounds_i, bounds_j, bounds_k, &
    2132         1978 :                                  RI_range, img_to_RI_cell, cell_to_index_ext)
    2133              : 
    2134              :       TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT)     :: t3c
    2135              :       REAL(KIND=dp), INTENT(IN)                          :: filter_eps
    2136              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2137              :       TYPE(neighbor_list_3c_type), INTENT(INOUT)         :: nl_3c
    2138              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j, basis_k
    2139              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
    2140              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: int_eps
    2141              :       INTEGER, INTENT(IN), OPTIONAL                      :: op_pos
    2142              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_kpoints, do_hfx_kpoints, &
    2143              :                                                             desymmetrize, cell_sym
    2144              :       INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL        :: bounds_i, bounds_j, bounds_k
    2145              :       REAL(dp), INTENT(IN), OPTIONAL                     :: RI_range
    2146              :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: img_to_RI_cell
    2147              :       INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER     :: cell_to_index_ext
    2148              : 
    2149              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integrals'
    2150              : 
    2151              :       INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
    2152              :          block_start_k, egfi, handle, handle2, i, iatom, ibasis, ikind, ilist, imax, iset, jatom, &
    2153              :          jcell, jkind, jset, katom, kcell, kkind, kset, m_max, max_ncoj, max_nset, max_nsgfi, &
    2154              :          maxli, maxlj, maxlk, mepos, natom, nbasis, ncell_RI, ncoi, ncoj, ncok, nimg, nseti, &
    2155              :          nsetj, nsetk, nthread, op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id
    2156         1978 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: img_to_RI_cell_prv
    2157              :       INTEGER, DIMENSION(2)                              :: bo
    2158              :       INTEGER, DIMENSION(3)                              :: blk_idx, blk_size, cell_j, cell_k, &
    2159              :                                                             kp_index_lbounds, kp_index_ubounds, sp
    2160         1978 :       INTEGER, DIMENSION(:), POINTER                     :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
    2161         3956 :                                                             lmin_k, npgfi, npgfj, npgfk, nsgfi, &
    2162         1978 :                                                             nsgfj, nsgfk
    2163         1978 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf_i, first_sgf_j, first_sgf_k
    2164         1978 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    2165              :       LOGICAL                                            :: block_not_zero, cell_sym_prv, debug, &
    2166              :                                                             desymmetrize_prv, do_hfx_kpoints_prv, &
    2167              :                                                             do_kpoints_prv, found, skip
    2168              :       REAL(KIND=dp)                                      :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
    2169              :                                                             kind_radius_i, kind_radius_j, &
    2170              :                                                             kind_radius_k, max_contraction_i, &
    2171              :                                                             prefac, sijk_ext
    2172         1978 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ccp_buffer, cpp_buffer, &
    2173         1978 :                                                             max_contraction_j, max_contraction_k
    2174         3956 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: block_t, dummy_block_t, sijk, &
    2175         3956 :                                                             sijk_contr, tmp_ijk
    2176              :       REAL(KIND=dp), DIMENSION(1, 1, 1)                  :: counter
    2177              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rij, rik, rj, rjk, rk
    2178         1978 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_i, set_radius_j, set_radius_k
    2179         1978 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
    2180         3956 :                                                             sphi_k, zeti, zetj, zetk
    2181         1978 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2182              :       TYPE(cell_type), POINTER                           :: cell
    2183         3956 :       TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER     :: spi, spk, tspj
    2184              :       TYPE(cp_libint_t)                                  :: lib
    2185        13846 :       TYPE(dbt_type)                                     :: t_3c_tmp
    2186              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2187              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
    2188              :       TYPE(kpoint_type), POINTER                         :: kpoints
    2189              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2190              :       TYPE(neighbor_list_3c_iterator_type)               :: nl_3c_iter
    2191         1978 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2192              : 
    2193         1978 :       CALL timeset(routineN, handle)
    2194              : 
    2195         1978 :       debug = .FALSE.
    2196              : 
    2197         1978 :       IF (PRESENT(do_kpoints)) THEN
    2198          382 :          do_kpoints_prv = do_kpoints
    2199              :       ELSE
    2200              :          do_kpoints_prv = .FALSE.
    2201              :       END IF
    2202              : 
    2203         1978 :       IF (PRESENT(do_hfx_kpoints)) THEN
    2204          140 :          do_hfx_kpoints_prv = do_hfx_kpoints
    2205              :       ELSE
    2206              :          do_hfx_kpoints_prv = .FALSE.
    2207              :       END IF
    2208          140 :       IF (do_hfx_kpoints_prv) THEN
    2209          140 :          CPASSERT(do_kpoints_prv)
    2210          140 :          CPASSERT(PRESENT(RI_range))
    2211          140 :          CPASSERT(PRESENT(img_to_RI_cell))
    2212              :       END IF
    2213              : 
    2214         1838 :       IF (PRESENT(img_to_RI_cell)) THEN
    2215          420 :          ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
    2216         6598 :          img_to_RI_cell_prv(:) = img_to_RI_cell
    2217              :       END IF
    2218              : 
    2219         1978 :       IF (PRESENT(desymmetrize)) THEN
    2220         1978 :          desymmetrize_prv = desymmetrize
    2221              :       ELSE
    2222              :          desymmetrize_prv = .TRUE.
    2223              :       END IF
    2224              : 
    2225         1978 :       IF (PRESENT(cell_sym)) THEN
    2226            6 :          cell_sym_prv = cell_sym
    2227              :       ELSE
    2228         1972 :          cell_sym_prv = .FALSE.
    2229              :       END IF
    2230              : 
    2231         1978 :       op_ij = do_potential_id; op_jk = do_potential_id
    2232              : 
    2233         1978 :       IF (PRESENT(op_pos)) THEN
    2234          474 :          op_pos_prv = op_pos
    2235              :       ELSE
    2236         1504 :          op_pos_prv = 1
    2237              :       END IF
    2238              : 
    2239         1838 :       SELECT CASE (op_pos_prv)
    2240              :       CASE (1)
    2241         1838 :          op_ij = potential_parameter%potential_type
    2242              :       CASE (2)
    2243         1978 :          op_jk = potential_parameter%potential_type
    2244              :       END SELECT
    2245              : 
    2246         1978 :       dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
    2247              : 
    2248              :       IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
    2249         1978 :           .OR. op_ij == do_potential_mix_cl_trunc) THEN
    2250         1438 :          dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
    2251         1438 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
    2252          540 :       ELSEIF (op_ij == do_potential_coulomb) THEN
    2253          192 :          dr_ij = 1000000.0_dp
    2254          192 :          dr_ik = 1000000.0_dp
    2255              :       END IF
    2256              : 
    2257              :       IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
    2258         1978 :           .OR. op_jk == do_potential_mix_cl_trunc) THEN
    2259           12 :          dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
    2260           12 :          dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
    2261         1966 :       ELSEIF (op_jk == do_potential_coulomb) THEN
    2262            0 :          dr_jk = 1000000.0_dp
    2263            0 :          dr_ik = 1000000.0_dp
    2264              :       END IF
    2265              : 
    2266         1978 :       NULLIFY (qs_kind_set, atomic_kind_set)
    2267              : 
    2268              :       ! get stuff
    2269              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
    2270         1978 :                       natom=natom, dft_control=dft_control, para_env=para_env)
    2271              : 
    2272         1978 :       IF (do_kpoints_prv) THEN
    2273          152 :          IF (PRESENT(cell_to_index_ext)) THEN
    2274            6 :             cell_to_index => cell_to_index_ext
    2275          288 :             nimg = MAXVAL(cell_to_index)
    2276              :          ELSE
    2277          146 :             CALL get_qs_env(qs_env, kpoints=kpoints)
    2278          146 :             CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
    2279          146 :             nimg = dft_control%nimages
    2280              :          END IF
    2281          152 :          ncell_RI = nimg
    2282          152 :          IF (do_hfx_kpoints_prv) THEN
    2283          140 :             nimg = SIZE(t3c, 1)
    2284          140 :             ncell_RI = SIZE(t3c, 2)
    2285              :          END IF
    2286              :       ELSE
    2287              :          nimg = 1
    2288              :          ncell_RI = 1
    2289              :       END IF
    2290              : 
    2291              :       CALL alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, &
    2292              :                           op_pos=op_pos_prv, do_kpoints=do_kpoints, do_hfx_kpoints=do_hfx_kpoints, &
    2293              :                           bounds_i=bounds_i, bounds_j=bounds_j, bounds_k=bounds_k, &
    2294              :                           RI_range=RI_range, img_to_RI_cell=img_to_RI_cell, cell_sym=cell_sym_prv, &
    2295         3816 :                           cell_to_index=cell_to_index)
    2296              : 
    2297         1978 :       IF (do_hfx_kpoints_prv) THEN
    2298          140 :          CPASSERT(op_pos_prv == 2)
    2299          140 :          CPASSERT(.NOT. desymmetrize_prv)
    2300         1838 :       ELSE IF (do_kpoints_prv) THEN
    2301           36 :          CPASSERT(ALL(SHAPE(t3c) == [nimg, ncell_RI]))
    2302              :       END IF
    2303              : 
    2304              :       !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
    2305         1978 :       nbasis = SIZE(basis_i)
    2306         1978 :       max_nsgfi = 0
    2307         1978 :       max_nset = 0
    2308         1978 :       maxli = 0
    2309         4932 :       DO ibasis = 1, nbasis
    2310              :          CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
    2311         2954 :                                 nset=iset, nsgf_set=nsgfi, npgf=npgfi)
    2312         2954 :          maxli = MAX(maxli, imax)
    2313         2954 :          max_nset = MAX(max_nset, iset)
    2314        17516 :          max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
    2315              :       END DO
    2316              :       max_ncoj = 0
    2317              :       maxlj = 0
    2318         4932 :       DO ibasis = 1, nbasis
    2319              :          CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
    2320         2954 :                                 nset=jset, nsgf_set=nsgfj, npgf=npgfj)
    2321         2954 :          maxlj = MAX(maxlj, imax)
    2322         2954 :          max_nset = MAX(max_nset, jset)
    2323        11314 :          max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
    2324              :       END DO
    2325              :       maxlk = 0
    2326         4932 :       DO ibasis = 1, nbasis
    2327              :          CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
    2328         2954 :                                 nset=kset, nsgf_set=nsgfk, npgf=npgfk)
    2329         2954 :          maxlk = MAX(maxlk, imax)
    2330         4932 :          max_nset = MAX(max_nset, kset)
    2331              :       END DO
    2332         1978 :       m_max = maxli + maxlj + maxlk
    2333              : 
    2334              :       !To minimize expensive memory ops and generally optimize contraction, pre-allocate
    2335              :       !contiguous sphi arrays (and transposed in the case of sphi_i)
    2336              : 
    2337         1978 :       NULLIFY (tspj, spi, spk)
    2338        75842 :       ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
    2339              : 
    2340         4932 :       DO ibasis = 1, nbasis
    2341        21984 :          DO iset = 1, max_nset
    2342        17052 :             NULLIFY (spi(iset, ibasis)%array)
    2343        17052 :             NULLIFY (tspj(iset, ibasis)%array)
    2344              : 
    2345        20006 :             NULLIFY (spk(iset, ibasis)%array)
    2346              :          END DO
    2347              :       END DO
    2348              : 
    2349         7912 :       DO ilist = 1, 3
    2350        16774 :          DO ibasis = 1, nbasis
    2351         8862 :             IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
    2352         8862 :             IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
    2353         8862 :             IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
    2354              : 
    2355        39876 :             DO iset = 1, basis_set%nset
    2356              : 
    2357        25080 :                ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
    2358        25080 :                sgfi = basis_set%first_sgf(1, iset)
    2359        25080 :                egfi = sgfi + basis_set%nsgf_set(iset) - 1
    2360              : 
    2361        33942 :                IF (ilist == 1) THEN
    2362        50336 :                   ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
    2363      3349236 :                   spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
    2364              : 
    2365        12496 :                ELSE IF (ilist == 2) THEN
    2366        25528 :                   ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
    2367       358742 :                   tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
    2368              : 
    2369              :                ELSE
    2370        24456 :                   ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
    2371      1495290 :                   spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
    2372              :                END IF
    2373              : 
    2374              :             END DO !iset
    2375              :          END DO !ibasis
    2376              :       END DO !ilist
    2377              : 
    2378              :       !Init the truncated Coulomb operator
    2379              :       IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated &
    2380         1978 :           .OR. op_ij == do_potential_mix_cl_trunc .OR. op_jk == do_potential_mix_cl_trunc) THEN
    2381              : 
    2382         1444 :          IF (m_max > get_lmax_init()) THEN
    2383           84 :             IF (para_env%mepos == 0) THEN
    2384           42 :                CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
    2385              :             END IF
    2386           84 :             CALL init(m_max, unit_id, para_env%mepos, para_env)
    2387           84 :             IF (para_env%mepos == 0) THEN
    2388           42 :                CALL close_file(unit_id)
    2389              :             END IF
    2390              :          END IF
    2391              :       END IF
    2392              : 
    2393         1978 :       CALL init_md_ftable(nmax=m_max)
    2394              : 
    2395         1978 :       IF (do_kpoints_prv) THEN
    2396          608 :          kp_index_lbounds = LBOUND(cell_to_index)
    2397          608 :          kp_index_ubounds = UBOUND(cell_to_index)
    2398              :       END IF
    2399              : 
    2400         7912 :       counter = 1.0_dp
    2401              : 
    2402         1978 :       nthread = 1
    2403         1978 : !$    nthread = omp_get_max_threads()
    2404              : 
    2405              : !$OMP PARALLEL DEFAULT(NONE) &
    2406              : !$OMP SHARED (nthread,do_kpoints_prv,kp_index_lbounds,kp_index_ubounds,maxli,maxlk,maxlj,bounds_i,&
    2407              : !$OMP         bounds_j,bounds_k,nimg,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,&
    2408              : !$OMP         potential_parameter,int_eps,t3c,tspj,spi,spk,debug,cell_to_index,&
    2409              : !$OMP         natom,nl_3c,cell,op_pos_prv,do_hfx_kpoints_prv,RI_range,ncell_RI, &
    2410              : !$OMP         img_to_RI_cell_prv, cell_sym_prv) &
    2411              : !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,cell_j,cell_k,&
    2412              : !$OMP          prefac,jcell,kcell,first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
    2413              : !$OMP          sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
    2414              : !$OMP          set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
    2415              : !$OMP          rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
    2416              : !$OMP          sgfk,sijk,ri,rj,rk,sijk_ext,block_not_zero,max_contraction_i,max_contraction_j,&
    2417              : !$OMP          max_contraction_k,iset,jset,kset,block_t,blk_size,sijk_contr,cpp_buffer,ccp_buffer,&
    2418              : !$OMP          block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
    2419         1978 : !$OMP          dummy_block_t,sp,handle2,mepos,bo,skip,tmp_ijk,i,blk_idx)
    2420              : 
    2421              :       mepos = 0
    2422              : !$    mepos = omp_get_thread_num()
    2423              : 
    2424              :       CALL cp_libint_init_3eri(lib, MAX(maxli, maxlj, maxlk))
    2425              :       CALL cp_libint_set_contrdepth(lib, 1)
    2426              :       CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
    2427              : 
    2428              :       !We split the provided bounds among the threads such that each threads works on a different set of atoms
    2429              :       IF (PRESENT(bounds_i)) THEN
    2430              :          bo = get_limit(bounds_i(2) - bounds_i(1) + 1, nthread, mepos)
    2431              :          bo(:) = bo(:) + bounds_i(1) - 1
    2432              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
    2433              :       ELSE IF (PRESENT(bounds_j)) THEN
    2434              : 
    2435              :          bo = get_limit(bounds_j(2) - bounds_j(1) + 1, nthread, mepos)
    2436              :          bo(:) = bo(:) + bounds_j(1) - 1
    2437              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bo, bounds_k)
    2438              :       ELSE IF (PRESENT(bounds_k)) THEN
    2439              :          bo = get_limit(bounds_k(2) - bounds_k(1) + 1, nthread, mepos)
    2440              :          bo(:) = bo(:) + bounds_k(1) - 1
    2441              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bo)
    2442              :       ELSE
    2443              :          bo = get_limit(natom, nthread, mepos)
    2444              :          CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
    2445              :       END IF
    2446              : 
    2447              :       skip = .FALSE.
    2448              :       IF (bo(1) > bo(2)) skip = .TRUE.
    2449              : 
    2450              :       DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
    2451              :          CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
    2452              :                                    iatom=iatom, jatom=jatom, katom=katom, &
    2453              :                                    rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
    2454              :          IF (skip) EXIT
    2455              : 
    2456              :          djk = NORM2(rjk)
    2457              :          dij = NORM2(rij)
    2458              :          dik = NORM2(rik)
    2459              : 
    2460              :          IF (nl_3c%sym == symmetric_jk) THEN
    2461              :             IF (jatom == katom) THEN
    2462              :                ! factor 0.5 due to double-counting of diagonal blocks
    2463              :                ! (we desymmetrize by adding transpose)
    2464              :                prefac = 0.5_dp
    2465              :             ELSE
    2466              :                prefac = 1.0_dp
    2467              :             END IF
    2468              :          ELSEIF (nl_3c%sym == symmetric_ij) THEN
    2469              :             IF (iatom == jatom) THEN
    2470              :                ! factor 0.5 due to double-counting of diagonal blocks
    2471              :                ! (we desymmetrize by adding transpose)
    2472              :                prefac = 0.5_dp
    2473              :             ELSE
    2474              :                prefac = 1.0_dp
    2475              :             END IF
    2476              :          ELSE
    2477              :             prefac = 1.0_dp
    2478              :          END IF
    2479              :          IF (do_kpoints_prv) prefac = 1.0_dp
    2480              : 
    2481              :          IF (do_kpoints_prv) THEN
    2482              : 
    2483              :             IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
    2484              :                 ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
    2485              : 
    2486              :             jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
    2487              :             IF (jcell > nimg .OR. jcell < 1) CYCLE
    2488              : 
    2489              :             IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
    2490              :                 ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
    2491              : 
    2492              :             kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
    2493              :             IF (kcell > nimg .OR. kcell < 1) CYCLE
    2494              :             !In case of HFX k-points, we only consider P^k if d_ik <= RI_range
    2495              :             IF (do_hfx_kpoints_prv) THEN
    2496              :                IF (dik > RI_range) CYCLE
    2497              :                kcell = img_to_RI_cell_prv(kcell)
    2498              :             END IF
    2499              :          ELSE
    2500              :             jcell = 1; kcell = 1
    2501              :          END IF
    2502              : 
    2503              :          IF (cell_sym_prv .AND. jcell < kcell) CYCLE
    2504              : 
    2505              :          blk_idx = [iatom, jatom, katom]
    2506              :          IF (do_hfx_kpoints_prv) THEN
    2507              :             blk_idx(3) = (kcell - 1)*natom + katom
    2508              :             kcell = 1
    2509              :          END IF
    2510              : 
    2511              :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
    2512              :                                 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
    2513              :                                 sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
    2514              : 
    2515              :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
    2516              :                                 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
    2517              :                                 sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
    2518              : 
    2519              :          CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
    2520              :                                 npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
    2521              :                                 sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
    2522              : 
    2523              :          IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
    2524              :          IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
    2525              :          IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
    2526              : 
    2527              :          ALLOCATE (max_contraction_j(nsetj))
    2528              :          DO jset = 1, nsetj
    2529              :             sgfj = first_sgf_j(1, jset)
    2530              :             max_contraction_j(jset) = MAXVAL((/(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/))
    2531              :          END DO
    2532              : 
    2533              :          ALLOCATE (max_contraction_k(nsetk))
    2534              :          DO kset = 1, nsetk
    2535              :             sgfk = first_sgf_k(1, kset)
    2536              :             max_contraction_k(kset) = MAXVAL((/(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/))
    2537              :          END DO
    2538              : 
    2539              :          CALL dbt_blk_sizes(t3c(jcell, kcell), blk_idx, blk_size)
    2540              : 
    2541              :          ALLOCATE (block_t(blk_size(2), blk_size(3), blk_size(1)))
    2542              : 
    2543              :          block_t = 0.0_dp
    2544              :          block_not_zero = .FALSE.
    2545              :          DO iset = 1, nseti
    2546              : 
    2547              :             sgfi = first_sgf_i(1, iset)
    2548              :             max_contraction_i = MAXVAL((/(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/))
    2549              : 
    2550              :             DO jset = 1, nsetj
    2551              : 
    2552              :                IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE
    2553              : 
    2554              :                DO kset = 1, nsetk
    2555              : 
    2556              :                   IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE
    2557              :                   IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE
    2558              : 
    2559              :                   ncoi = npgfi(iset)*ncoset(lmax_i(iset))
    2560              :                   ncoj = npgfj(jset)*ncoset(lmax_j(jset))
    2561              :                   ncok = npgfk(kset)*ncoset(lmax_k(kset))
    2562              : 
    2563              :                   !ensure non-zero number of triples below
    2564              :                   IF (ncoj*ncok*ncoi == 0) CYCLE
    2565              : 
    2566              :                   !need positions for libint. Only relative positions are needed => set ri to 0.0
    2567              :                   ri = 0.0_dp
    2568              :                   rj = rij ! ri + rij
    2569              :                   rk = rik ! ri + rik
    2570              : 
    2571              :                   ALLOCATE (sijk(ncoj, ncok, ncoi))
    2572              :                   IF (op_pos_prv == 1) THEN
    2573              :                      sijk(:, :, :) = 0.0_dp
    2574              :                      CALL eri_3center(sijk, &
    2575              :                                       lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
    2576              :                                       lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
    2577              :                                       lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
    2578              :                                       djk, dij, dik, lib, potential_parameter, int_abc_ext=sijk_ext)
    2579              :                   ELSE
    2580              :                      ALLOCATE (tmp_ijk(ncoi, ncoj, ncok))
    2581              :                      tmp_ijk(:, :, :) = 0.0_dp
    2582              :                      CALL eri_3center(tmp_ijk, &
    2583              :                                       lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
    2584              :                                       lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
    2585              :                                       lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
    2586              :                                       dij, dik, djk, lib, potential_parameter, int_abc_ext=sijk_ext)
    2587              : 
    2588              :                      !F08: sijk = RESHAPE(tmp_ijk, [ncoj, ncok, ncoi], ORDER=[2, 3, 1]) with sijk not allocated
    2589              :                      DO i = 1, ncoi !TODO: revise/check for efficiency
    2590              :                         sijk(:, :, i) = tmp_ijk(i, :, :)
    2591              :                      END DO
    2592              :                      DEALLOCATE (tmp_ijk)
    2593              :                   END IF
    2594              : 
    2595              :                   IF (PRESENT(int_eps)) THEN
    2596              :                      IF (int_eps > sijk_ext*(max_contraction_i* &
    2597              :                                              max_contraction_j(jset)* &
    2598              :                                              max_contraction_k(kset))) THEN
    2599              :                         DEALLOCATE (sijk)
    2600              :                         CYCLE
    2601              :                      END IF
    2602              :                   END IF
    2603              : 
    2604              :                   block_not_zero = .TRUE.
    2605              :                   ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
    2606              :                   CALL abc_contract_xsmm(sijk_contr, sijk, tspj(jset, jkind)%array, &
    2607              :                                          spk(kset, kkind)%array, spi(iset, ikind)%array, &
    2608              :                                          ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
    2609              :                                          nsgfi(iset), cpp_buffer, ccp_buffer, prefac)
    2610              :                   DEALLOCATE (sijk)
    2611              : 
    2612              :                   sgfj = first_sgf_j(1, jset)
    2613              :                   sgfk = first_sgf_k(1, kset)
    2614              : 
    2615              :                   block_start_j = sgfj
    2616              :                   block_end_j = sgfj + nsgfj(jset) - 1
    2617              :                   block_start_k = sgfk
    2618              :                   block_end_k = sgfk + nsgfk(kset) - 1
    2619              :                   block_start_i = sgfi
    2620              :                   block_end_i = sgfi + nsgfi(iset) - 1
    2621              : 
    2622              :                   block_t(block_start_j:block_end_j, &
    2623              :                           block_start_k:block_end_k, &
    2624              :                           block_start_i:block_end_i) = &
    2625              :                      block_t(block_start_j:block_end_j, &
    2626              :                              block_start_k:block_end_k, &
    2627              :                              block_start_i:block_end_i) + &
    2628              :                      sijk_contr(:, :, :)
    2629              :                   DEALLOCATE (sijk_contr)
    2630              : 
    2631              :                END DO
    2632              : 
    2633              :             END DO
    2634              : 
    2635              :          END DO
    2636              : 
    2637              :          IF (block_not_zero) THEN
    2638              : !$OMP CRITICAL
    2639              :             CALL timeset(routineN//"_put_dbcsr", handle2)
    2640              :             IF (debug) THEN
    2641              :                CALL dbt_get_block(t3c(jcell, kcell), blk_idx, dummy_block_t, found=found)
    2642              :                CPASSERT(found)
    2643              :             END IF
    2644              : 
    2645              :             sp = SHAPE(block_t)
    2646              :             sp([2, 3, 1]) = sp ! sp = sp([2, 3, 1]) performs worse
    2647              : 
    2648              :             CALL dbt_put_block(t3c(jcell, kcell), blk_idx, sp, &
    2649              :                                RESHAPE(block_t, SHAPE=sp, ORDER=[2, 3, 1]), summation=.TRUE.)
    2650              : 
    2651              :             CALL timestop(handle2)
    2652              : !$OMP END CRITICAL
    2653              :          END IF
    2654              : 
    2655              :          DEALLOCATE (block_t)
    2656              :          DEALLOCATE (max_contraction_j, max_contraction_k)
    2657              :       END DO
    2658              : 
    2659              :       IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
    2660              :       IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
    2661              : 
    2662              :       CALL cp_libint_cleanup_3eri(lib)
    2663              :       CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
    2664              : !$OMP END PARALLEL
    2665              : 
    2666              :       !TODO: deal with hfx_kpoints, because should not filter by 1/2
    2667         1978 :       IF (nl_3c%sym == symmetric_jk .OR. do_kpoints_prv) THEN
    2668              : 
    2669          716 :          IF (.NOT. do_hfx_kpoints_prv) THEN
    2670         1228 :             DO kcell = 1, nimg
    2671         2568 :                DO jcell = 1, nimg
    2672              :                   ! need half of filter eps because afterwards we add transposed tensor
    2673         1992 :                   CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
    2674              :                END DO
    2675              :             END DO
    2676              :          ELSE
    2677          280 :             DO kcell = 1, ncell_RI
    2678         4058 :                DO jcell = 1, nimg
    2679         3918 :                   CALL dbt_filter(t3c(jcell, kcell), filter_eps)
    2680              :                END DO
    2681              :             END DO
    2682              :          END IF
    2683              : 
    2684          716 :          IF (desymmetrize_prv) THEN
    2685              :             ! add transposed of overlap integrals
    2686            0 :             CALL dbt_create(t3c(1, 1), t_3c_tmp)
    2687            0 :             DO kcell = 1, jcell
    2688            0 :                DO jcell = 1, nimg
    2689            0 :                   CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
    2690            0 :                   CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.)
    2691            0 :                   CALL dbt_filter(t3c(kcell, jcell), filter_eps)
    2692              :                END DO
    2693              :             END DO
    2694            0 :             DO kcell = jcell + 1, nimg
    2695            0 :                DO jcell = 1, nimg
    2696            0 :                   CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
    2697            0 :                   CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.)
    2698            0 :                   CALL dbt_filter(t3c(kcell, jcell), filter_eps)
    2699              :                END DO
    2700              :             END DO
    2701            0 :             CALL dbt_destroy(t_3c_tmp)
    2702              :          END IF
    2703         1262 :       ELSEIF (nl_3c%sym == symmetric_ij) THEN
    2704            0 :          DO kcell = 1, nimg
    2705            0 :             DO jcell = 1, nimg
    2706            0 :                CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
    2707              :             END DO
    2708              :          END DO
    2709         1262 :       ELSEIF (nl_3c%sym == symmetric_none) THEN
    2710         2524 :          DO kcell = 1, nimg
    2711         3786 :             DO jcell = 1, nimg
    2712         2524 :                CALL dbt_filter(t3c(jcell, kcell), filter_eps)
    2713              :             END DO
    2714              :          END DO
    2715              :       ELSE
    2716            0 :          CPABORT("requested symmetric case not implemented")
    2717              :       END IF
    2718              : 
    2719        12958 :       DO iset = 1, max_nset
    2720        30010 :          DO ibasis = 1, nbasis
    2721        17052 :             IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
    2722        17052 :             IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
    2723        28032 :             IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
    2724              :          END DO
    2725              :       END DO
    2726         1978 :       DEALLOCATE (spi, tspj, spk)
    2727              : 
    2728         1978 :       CALL timestop(handle)
    2729        17802 :    END SUBROUTINE build_3c_integrals
    2730              : 
    2731              : ! **************************************************************************************************
    2732              : !> \brief Calculates the derivatives of 2-center integrals, wrt to the first center
    2733              : !> \param t2c_der ...
    2734              : !> this routine requires that libint has been static initialised somewhere else
    2735              : !> \param filter_eps Filter threshold for matrix blocks
    2736              : !> \param qs_env ...
    2737              : !> \param nl_2c 2-center neighborlist
    2738              : !> \param basis_i ...
    2739              : !> \param basis_j ...
    2740              : !> \param potential_parameter ...
    2741              : !> \param do_kpoints ...
    2742              : ! **************************************************************************************************
    2743          384 :    SUBROUTINE build_2c_derivatives(t2c_der, filter_eps, qs_env, &
    2744          384 :                                    nl_2c, basis_i, basis_j, &
    2745              :                                    potential_parameter, do_kpoints)
    2746              : 
    2747              :       TYPE(dbcsr_type), DIMENSION(:, :), INTENT(INOUT)   :: t2c_der
    2748              :       REAL(KIND=dp), INTENT(IN)                          :: filter_eps
    2749              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2750              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2751              :          POINTER                                         :: nl_2c
    2752              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j
    2753              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
    2754              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_kpoints
    2755              : 
    2756              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_2c_derivatives'
    2757              : 
    2758              :       INTEGER :: handle, i_xyz, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, &
    2759              :          jset, m_max, maxli, maxlj, natom, ncoi, ncoj, nimg, nseti, nsetj, op_prv, sgfi, sgfj, &
    2760              :          unit_id
    2761              :       INTEGER, DIMENSION(3)                              :: cell_j, kp_index_lbounds, &
    2762              :                                                             kp_index_ubounds
    2763          384 :       INTEGER, DIMENSION(:), POINTER                     :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
    2764          384 :                                                             npgfj, nsgfi, nsgfj
    2765          384 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf_i, first_sgf_j
    2766          384 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    2767              :       LOGICAL                                            :: do_kpoints_prv, do_symmetric, found, &
    2768              :                                                             trans
    2769              :       REAL(KIND=dp)                                      :: dab
    2770          384 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dij_contr, dij_rs
    2771          384 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dij
    2772              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rij, rj
    2773          384 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_i, set_radius_j
    2774          384 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
    2775          384 :                                                             zetj
    2776          384 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2777         1536 :       TYPE(block_p_type), DIMENSION(3)                   :: block_t
    2778              :       TYPE(cp_libint_t)                                  :: lib
    2779              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2780              :       TYPE(kpoint_type), POINTER                         :: kpoints
    2781              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2782              :       TYPE(neighbor_list_iterator_p_type), &
    2783          384 :          DIMENSION(:), POINTER                           :: nl_iterator
    2784          384 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2785              : 
    2786          384 :       CALL timeset(routineN, handle)
    2787              : 
    2788          384 :       IF (PRESENT(do_kpoints)) THEN
    2789           92 :          do_kpoints_prv = do_kpoints
    2790              :       ELSE
    2791              :          do_kpoints_prv = .FALSE.
    2792              :       END IF
    2793              : 
    2794          384 :       op_prv = potential_parameter%potential_type
    2795              : 
    2796          384 :       NULLIFY (qs_kind_set, atomic_kind_set, block_t(1)%block, block_t(2)%block, block_t(3)%block, cell_to_index)
    2797              : 
    2798              :       ! get stuff
    2799              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
    2800          384 :                       natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
    2801              : 
    2802          384 :       IF (do_kpoints_prv) THEN
    2803           92 :          nimg = SIZE(t2c_der, 1)
    2804           92 :          CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
    2805          368 :          kp_index_lbounds = LBOUND(cell_to_index)
    2806          368 :          kp_index_ubounds = UBOUND(cell_to_index)
    2807              :       ELSE
    2808              :          nimg = 1
    2809              :       END IF
    2810              : 
    2811              :       ! check for symmetry
    2812          384 :       CPASSERT(SIZE(nl_2c) > 0)
    2813          384 :       CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
    2814              : 
    2815          384 :       IF (do_symmetric) THEN
    2816          584 :          DO img = 1, nimg
    2817              :             !Derivtive matrix is assymetric
    2818         1460 :             DO i_xyz = 1, 3
    2819         1168 :                CPASSERT(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_antisymmetric)
    2820              :             END DO
    2821              :          END DO
    2822              :       ELSE
    2823         2296 :          DO img = 1, nimg
    2824         8908 :             DO i_xyz = 1, 3
    2825         8816 :                CPASSERT(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_no_symmetry)
    2826              :             END DO
    2827              :          END DO
    2828              :       END IF
    2829              : 
    2830         2880 :       DO img = 1, nimg
    2831        10368 :          DO i_xyz = 1, 3
    2832         9984 :             CALL cp_dbcsr_alloc_block_from_nbl(t2c_der(img, i_xyz), nl_2c)
    2833              :          END DO
    2834              :       END DO
    2835              : 
    2836          384 :       maxli = 0
    2837         1086 :       DO ibasis = 1, SIZE(basis_i)
    2838          702 :          CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
    2839         1086 :          maxli = MAX(maxli, imax)
    2840              :       END DO
    2841          384 :       maxlj = 0
    2842         1086 :       DO ibasis = 1, SIZE(basis_j)
    2843          702 :          CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
    2844         1086 :          maxlj = MAX(maxlj, imax)
    2845              :       END DO
    2846              : 
    2847          384 :       m_max = maxli + maxlj + 1
    2848              : 
    2849              :       !Init the truncated Coulomb operator
    2850          384 :       IF (op_prv == do_potential_truncated .OR. op_prv == do_potential_mix_cl_trunc) THEN
    2851              : 
    2852           60 :          IF (m_max > get_lmax_init()) THEN
    2853           22 :             IF (para_env%mepos == 0) THEN
    2854           11 :                CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
    2855              :             END IF
    2856           22 :             CALL init(m_max, unit_id, para_env%mepos, para_env)
    2857           22 :             IF (para_env%mepos == 0) THEN
    2858           11 :                CALL close_file(unit_id)
    2859              :             END IF
    2860              :          END IF
    2861              :       END IF
    2862              : 
    2863          384 :       CALL init_md_ftable(nmax=m_max)
    2864              : 
    2865          384 :       CALL cp_libint_init_2eri1(lib, MAX(maxli, maxlj))
    2866          384 :       CALL cp_libint_set_contrdepth(lib, 1)
    2867              : 
    2868          384 :       CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
    2869        14908 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    2870              : 
    2871              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
    2872        14524 :                                 iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
    2873        14524 :          IF (do_kpoints_prv) THEN
    2874        11326 :             IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
    2875              :                 ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
    2876         1618 :             img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
    2877         1618 :             IF (img > nimg .OR. img < 1) CYCLE
    2878              :          ELSE
    2879              :             img = 1
    2880              :          END IF
    2881              : 
    2882              :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
    2883              :                                 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
    2884        14496 :                                 sphi=sphi_i, zet=zeti)
    2885              : 
    2886              :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
    2887              :                                 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
    2888        14496 :                                 sphi=sphi_j, zet=zetj)
    2889              : 
    2890        14496 :          IF (do_symmetric) THEN
    2891        12906 :             IF (iatom <= jatom) THEN
    2892         7900 :                irow = iatom
    2893         7900 :                icol = jatom
    2894              :             ELSE
    2895         5006 :                irow = jatom
    2896         5006 :                icol = iatom
    2897              :             END IF
    2898              :          ELSE
    2899         1590 :             irow = iatom
    2900         1590 :             icol = jatom
    2901              :          END IF
    2902              : 
    2903        57984 :          dab = NORM2(rij)
    2904        14496 :          trans = do_symmetric .AND. (iatom > jatom)
    2905              : 
    2906        57984 :          DO i_xyz = 1, 3
    2907              :             CALL dbcsr_get_block_p(matrix=t2c_der(img, i_xyz), &
    2908        43488 :                                    row=irow, col=icol, BLOCK=block_t(i_xyz)%block, found=found)
    2909        57984 :             CPASSERT(found)
    2910              :          END DO
    2911              : 
    2912        49556 :          DO iset = 1, nseti
    2913              : 
    2914        34676 :             ncoi = npgfi(iset)*ncoset(lmax_i(iset))
    2915        34676 :             sgfi = first_sgf_i(1, iset)
    2916              : 
    2917       172652 :             DO jset = 1, nsetj
    2918              : 
    2919       123452 :                ncoj = npgfj(jset)*ncoset(lmax_j(jset))
    2920       123452 :                sgfj = first_sgf_j(1, jset)
    2921              : 
    2922       158128 :                IF (ncoi*ncoj > 0) THEN
    2923       493808 :                   ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
    2924       617260 :                   ALLOCATE (dij(ncoi, ncoj, 3))
    2925    313443755 :                   dij(:, :, :) = 0.0_dp
    2926              : 
    2927       123452 :                   ri = 0.0_dp
    2928       123452 :                   rj = rij
    2929              : 
    2930              :                   CALL eri_2center_derivs(dij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
    2931              :                                           rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
    2932       123452 :                                           rpgf_j(:, jset), rj, dab, lib, potential_parameter)
    2933              : 
    2934       493808 :                   DO i_xyz = 1, 3
    2935              : 
    2936    123412836 :                      dij_contr(:, :) = 0.0_dp
    2937              :                      CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
    2938              :                                       sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
    2939       370356 :                                       ncoi, ncoj, nsgfi(iset), nsgfj(jset))
    2940              : 
    2941       370356 :                      IF (trans) THEN
    2942              :                         !if transpose, then -1 factor for antisymmetry
    2943       449508 :                         ALLOCATE (dij_rs(nsgfj(jset), nsgfi(iset)))
    2944     44948634 :                         dij_rs(:, :) = -1.0_dp*TRANSPOSE(dij_contr)
    2945              :                      ELSE
    2946      1031916 :                         ALLOCATE (dij_rs(nsgfi(iset), nsgfj(jset)))
    2947     78355215 :                         dij_rs(:, :) = dij_contr
    2948              :                      END IF
    2949              : 
    2950              :                      CALL block_add("IN", dij_rs, &
    2951              :                                     nsgfi(iset), nsgfj(jset), block_t(i_xyz)%block, &
    2952       370356 :                                     sgfi, sgfj, trans=trans)
    2953       493808 :                      DEALLOCATE (dij_rs)
    2954              :                   END DO
    2955              : 
    2956       123452 :                   DEALLOCATE (dij, dij_contr)
    2957              :                END IF
    2958              :             END DO
    2959              :          END DO
    2960              :       END DO
    2961              : 
    2962          384 :       CALL cp_libint_cleanup_2eri1(lib)
    2963          384 :       CALL neighbor_list_iterator_release(nl_iterator)
    2964              : 
    2965         2880 :       DO img = 1, nimg
    2966        10368 :          DO i_xyz = 1, 3
    2967         7488 :             CALL dbcsr_finalize(t2c_der(img, i_xyz))
    2968         9984 :             CALL dbcsr_filter(t2c_der(img, i_xyz), filter_eps)
    2969              :          END DO
    2970              :       END DO
    2971              : 
    2972          384 :       CALL timestop(handle)
    2973              : 
    2974         1152 :    END SUBROUTINE build_2c_derivatives
    2975              : 
    2976              : ! **************************************************************************************************
    2977              : !> \brief Calculates the virial coming from 2c derivatives on the fly
    2978              : !> \param work_virial ...
    2979              : !> \param t2c_trace the 2c tensor that we should trace with the derivatives
    2980              : !> \param pref ...
    2981              : !> \param qs_env ...
    2982              : !> \param nl_2c 2-center neighborlist. Assumed to have compatible distribution with t2c_trace,
    2983              : !>              and to be non-symmetric
    2984              : !> \param basis_i ...
    2985              : !> \param basis_j ...
    2986              : !> \param potential_parameter ...
    2987              : ! **************************************************************************************************
    2988           12 :    SUBROUTINE calc_2c_virial(work_virial, t2c_trace, pref, qs_env, nl_2c, basis_i, basis_j, potential_parameter)
    2989              :       REAL(dp), DIMENSION(3, 3), INTENT(INOUT)           :: work_virial
    2990              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: t2c_trace
    2991              :       REAL(KIND=dp), INTENT(IN)                          :: pref
    2992              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2993              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2994              :          POINTER                                         :: nl_2c
    2995              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j
    2996              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
    2997              : 
    2998              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calc_2c_virial'
    2999              : 
    3000              :       INTEGER :: handle, i_xyz, iatom, ibasis, ikind, imax, iset, j_xyz, jatom, jkind, jset, &
    3001              :          m_max, maxli, maxlj, natom, ncoi, ncoj, nseti, nsetj, op_prv, sgfi, sgfj, unit_id
    3002           12 :       INTEGER, DIMENSION(:), POINTER                     :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
    3003           12 :                                                             npgfj, nsgfi, nsgfj
    3004           12 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf_i, first_sgf_j
    3005              :       LOGICAL                                            :: do_symmetric, found
    3006              :       REAL(dp)                                           :: force
    3007           12 :       REAL(dp), DIMENSION(:, :), POINTER                 :: pblock
    3008              :       REAL(KIND=dp)                                      :: dab
    3009           12 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dij_contr
    3010           12 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dij
    3011              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rij, rj, scoord
    3012           12 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_i, set_radius_j
    3013           12 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
    3014           12 :                                                             zetj
    3015           12 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    3016              :       TYPE(cell_type), POINTER                           :: cell
    3017              :       TYPE(cp_libint_t)                                  :: lib
    3018              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3019              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    3020              :       TYPE(neighbor_list_iterator_p_type), &
    3021           12 :          DIMENSION(:), POINTER                           :: nl_iterator
    3022           12 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    3023           12 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    3024              : 
    3025           12 :       CALL timeset(routineN, handle)
    3026              : 
    3027           12 :       op_prv = potential_parameter%potential_type
    3028              : 
    3029           12 :       NULLIFY (qs_kind_set, atomic_kind_set, pblock, particle_set, cell)
    3030              : 
    3031              :       ! get stuff
    3032              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
    3033              :                       natom=natom, dft_control=dft_control, para_env=para_env, &
    3034           12 :                       particle_set=particle_set, cell=cell)
    3035              : 
    3036              :       ! check for symmetry
    3037           12 :       CPASSERT(SIZE(nl_2c) > 0)
    3038           12 :       CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
    3039           12 :       CPASSERT(.NOT. do_symmetric)
    3040              : 
    3041           12 :       maxli = 0
    3042           36 :       DO ibasis = 1, SIZE(basis_i)
    3043           24 :          CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
    3044           36 :          maxli = MAX(maxli, imax)
    3045              :       END DO
    3046           12 :       maxlj = 0
    3047           36 :       DO ibasis = 1, SIZE(basis_j)
    3048           24 :          CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
    3049           36 :          maxlj = MAX(maxlj, imax)
    3050              :       END DO
    3051              : 
    3052           12 :       m_max = maxli + maxlj + 1
    3053              : 
    3054              :       !Init the truncated Coulomb operator
    3055           12 :       IF (op_prv == do_potential_truncated .OR. op_prv == do_potential_mix_cl_trunc) THEN
    3056              : 
    3057            2 :          IF (m_max > get_lmax_init()) THEN
    3058            0 :             IF (para_env%mepos == 0) THEN
    3059            0 :                CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
    3060              :             END IF
    3061            0 :             CALL init(m_max, unit_id, para_env%mepos, para_env)
    3062            0 :             IF (para_env%mepos == 0) THEN
    3063            0 :                CALL close_file(unit_id)
    3064              :             END IF
    3065              :          END IF
    3066              :       END IF
    3067              : 
    3068           12 :       CALL init_md_ftable(nmax=m_max)
    3069              : 
    3070           12 :       CALL cp_libint_init_2eri1(lib, MAX(maxli, maxlj))
    3071           12 :       CALL cp_libint_set_contrdepth(lib, 1)
    3072              : 
    3073           12 :       CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
    3074         1644 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    3075              : 
    3076              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
    3077         1632 :                                 iatom=iatom, jatom=jatom, r=rij)
    3078              : 
    3079              :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
    3080              :                                 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
    3081         1632 :                                 sphi=sphi_i, zet=zeti)
    3082              : 
    3083              :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
    3084              :                                 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
    3085         1632 :                                 sphi=sphi_j, zet=zetj)
    3086              : 
    3087         6528 :          dab = NORM2(rij)
    3088              : 
    3089         1632 :          CALL dbcsr_get_block_p(t2c_trace, iatom, jatom, pblock, found)
    3090         1632 :          IF (.NOT. found) CYCLE
    3091              : 
    3092         6113 :          DO iset = 1, nseti
    3093              : 
    3094         4469 :             ncoi = npgfi(iset)*ncoset(lmax_i(iset))
    3095         4469 :             sgfi = first_sgf_i(1, iset)
    3096              : 
    3097        25808 :             DO jset = 1, nsetj
    3098              : 
    3099        19707 :                ncoj = npgfj(jset)*ncoset(lmax_j(jset))
    3100        19707 :                sgfj = first_sgf_j(1, jset)
    3101              : 
    3102        24176 :                IF (ncoi*ncoj > 0) THEN
    3103        78828 :                   ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
    3104        98535 :                   ALLOCATE (dij(ncoi, ncoj, 3))
    3105      5649186 :                   dij(:, :, :) = 0.0_dp
    3106              : 
    3107        19707 :                   ri = 0.0_dp
    3108        19707 :                   rj = rij
    3109              : 
    3110              :                   CALL eri_2center_derivs(dij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
    3111              :                                           rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
    3112        19707 :                                           rpgf_j(:, jset), rj, dab, lib, potential_parameter)
    3113              : 
    3114        78828 :                   DO i_xyz = 1, 3
    3115              : 
    3116      2522610 :                      dij_contr(:, :) = 0.0_dp
    3117              :                      CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
    3118              :                                       sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
    3119        59121 :                                       ncoi, ncoj, nsgfi(iset), nsgfj(jset))
    3120              : 
    3121      2522610 :                      force = SUM(pblock(sgfi:sgfi + nsgfi(iset) - 1, sgfj:sgfj + nsgfj(jset) - 1)*dij_contr(:, :))
    3122        59121 :                      force = pref*force
    3123              : 
    3124              :                      !iatom virial
    3125        59121 :                      CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
    3126       236484 :                      DO j_xyz = 1, 3
    3127       236484 :                         work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
    3128              :                      END DO
    3129              : 
    3130              :                      !jatom virial
    3131       236484 :                      CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
    3132       256191 :                      DO j_xyz = 1, 3
    3133       236484 :                         work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - force*scoord(j_xyz)
    3134              :                      END DO
    3135              :                   END DO
    3136              : 
    3137        19707 :                   DEALLOCATE (dij, dij_contr)
    3138              :                END IF
    3139              :             END DO
    3140              :          END DO
    3141              :       END DO
    3142              : 
    3143           12 :       CALL neighbor_list_iterator_release(nl_iterator)
    3144           12 :       CALL cp_libint_cleanup_2eri1(lib)
    3145              : 
    3146           12 :       CALL timestop(handle)
    3147              : 
    3148           24 :    END SUBROUTINE calc_2c_virial
    3149              : 
    3150              : ! **************************************************************************************************
    3151              : !> \brief ...
    3152              : !> \param t2c empty DBCSR matrix
    3153              : !> \param filter_eps Filter threshold for matrix blocks
    3154              : !> \param qs_env ...
    3155              : !> \param nl_2c 2-center neighborlist
    3156              : !> \param basis_i ...
    3157              : !> \param basis_j ...
    3158              : !> \param potential_parameter ...
    3159              : !> \param do_kpoints ...
    3160              : !> \param do_hfx_kpoints ...
    3161              : !> this routine requires that libint has been static initialised somewhere else
    3162              : !> \param ext_kpoints ...
    3163              : !> \param regularization_RI ...
    3164              : ! **************************************************************************************************
    3165          882 :    SUBROUTINE build_2c_integrals(t2c, filter_eps, qs_env, &
    3166          882 :                                  nl_2c, basis_i, basis_j, &
    3167              :                                  potential_parameter, do_kpoints, &
    3168              :                                  do_hfx_kpoints, ext_kpoints, regularization_RI)
    3169              : 
    3170              :       TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT)      :: t2c
    3171              :       REAL(KIND=dp), INTENT(IN)                          :: filter_eps
    3172              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    3173              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    3174              :          POINTER                                         :: nl_2c
    3175              :       TYPE(gto_basis_set_p_type), DIMENSION(:)           :: basis_i, basis_j
    3176              :       TYPE(libint_potential_type), INTENT(IN)            :: potential_parameter
    3177              :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_kpoints, do_hfx_kpoints
    3178              :       TYPE(kpoint_type), OPTIONAL, POINTER               :: ext_kpoints
    3179              :       REAL(KIND=dp), OPTIONAL                            :: regularization_RI
    3180              : 
    3181              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_2c_integrals'
    3182              : 
    3183              :       INTEGER :: handle, i_diag, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, &
    3184              :          jset, m_max, maxli, maxlj, natom, ncoi, ncoj, nimg, nseti, nsetj, op_prv, sgfi, sgfj, &
    3185              :          unit_id
    3186              :       INTEGER, DIMENSION(3)                              :: cell_j, kp_index_lbounds, &
    3187              :                                                             kp_index_ubounds
    3188          882 :       INTEGER, DIMENSION(:), POINTER                     :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
    3189          882 :                                                             npgfj, nsgfi, nsgfj
    3190          882 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgf_i, first_sgf_j
    3191          882 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    3192              :       LOGICAL                                            :: do_hfx_kpoints_prv, do_kpoints_prv, &
    3193              :                                                             do_symmetric, found, trans
    3194              :       REAL(KIND=dp)                                      :: dab, min_zet
    3195          882 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: sij, sij_contr, sij_rs
    3196              :       REAL(KIND=dp), DIMENSION(3)                        :: ri, rij, rj
    3197          882 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_i, set_radius_j
    3198          882 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
    3199          882 :                                                             zetj
    3200          882 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    3201              :       TYPE(block_p_type)                                 :: block_t
    3202              :       TYPE(cell_type), POINTER                           :: cell
    3203              :       TYPE(cp_libint_t)                                  :: lib
    3204              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3205              :       TYPE(kpoint_type), POINTER                         :: kpoints
    3206              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    3207              :       TYPE(neighbor_list_iterator_p_type), &
    3208          882 :          DIMENSION(:), POINTER                           :: nl_iterator
    3209          882 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    3210              : 
    3211          882 :       CALL timeset(routineN, handle)
    3212              : 
    3213          882 :       IF (PRESENT(do_kpoints)) THEN
    3214          832 :          do_kpoints_prv = do_kpoints
    3215              :       ELSE
    3216              :          do_kpoints_prv = .FALSE.
    3217              :       END IF
    3218              : 
    3219          882 :       IF (PRESENT(do_hfx_kpoints)) THEN
    3220          302 :          do_hfx_kpoints_prv = do_hfx_kpoints
    3221              :       ELSE
    3222              :          do_hfx_kpoints_prv = .FALSE.
    3223              :       END IF
    3224          302 :       IF (do_hfx_kpoints_prv) THEN
    3225          160 :          CPASSERT(do_kpoints_prv)
    3226              :       END IF
    3227              : 
    3228          882 :       op_prv = potential_parameter%potential_type
    3229              : 
    3230          882 :       NULLIFY (qs_kind_set, atomic_kind_set, block_t%block, cell_to_index)
    3231              : 
    3232              :       ! get stuff
    3233              :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
    3234          882 :                       natom=natom, dft_control=dft_control, para_env=para_env, kpoints=kpoints)
    3235              : 
    3236          882 :       IF (PRESENT(ext_kpoints)) kpoints => ext_kpoints
    3237              : 
    3238          882 :       IF (do_kpoints_prv) THEN
    3239          464 :          nimg = SIZE(t2c)
    3240          464 :          CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
    3241         1856 :          kp_index_lbounds = LBOUND(cell_to_index)
    3242         1856 :          kp_index_ubounds = UBOUND(cell_to_index)
    3243              :       ELSE
    3244              :          nimg = 1
    3245              :       END IF
    3246              : 
    3247              :       ! check for symmetry
    3248          882 :       CPASSERT(SIZE(nl_2c) > 0)
    3249          882 :       CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
    3250              : 
    3251          882 :       IF (do_symmetric) THEN
    3252         2154 :          DO img = 1, nimg
    3253         2154 :             CPASSERT(dbcsr_has_symmetry(t2c(img)))
    3254              :          END DO
    3255              :       ELSE
    3256         4884 :          DO img = 1, nimg
    3257         4884 :             CPASSERT(.NOT. dbcsr_has_symmetry(t2c(img)))
    3258              :          END DO
    3259              :       END IF
    3260              : 
    3261         7038 :       DO img = 1, nimg
    3262         7038 :          CALL cp_dbcsr_alloc_block_from_nbl(t2c(img), nl_2c)
    3263              :       END DO
    3264              : 
    3265          882 :       maxli = 0
    3266         2470 :       DO ibasis = 1, SIZE(basis_i)
    3267         1588 :          CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
    3268         2470 :          maxli = MAX(maxli, imax)
    3269              :       END DO
    3270          882 :       maxlj = 0
    3271         2470 :       DO ibasis = 1, SIZE(basis_j)
    3272         1588 :          CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
    3273         2470 :          maxlj = MAX(maxlj, imax)
    3274              :       END DO
    3275              : 
    3276          882 :       m_max = maxli + maxlj
    3277              : 
    3278              :       !Init the truncated Coulomb operator
    3279          882 :       IF (op_prv == do_potential_truncated .OR. op_prv == do_potential_mix_cl_trunc) THEN
    3280              : 
    3281          544 :          IF (m_max > get_lmax_init()) THEN
    3282          106 :             IF (para_env%mepos == 0) THEN
    3283           53 :                CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
    3284              :             END IF
    3285          106 :             CALL init(m_max, unit_id, para_env%mepos, para_env)
    3286          106 :             IF (para_env%mepos == 0) THEN
    3287           53 :                CALL close_file(unit_id)
    3288              :             END IF
    3289              :          END IF
    3290              :       END IF
    3291              : 
    3292          882 :       CALL init_md_ftable(nmax=m_max)
    3293              : 
    3294          882 :       CALL cp_libint_init_2eri(lib, MAX(maxli, maxlj))
    3295          882 :       CALL cp_libint_set_contrdepth(lib, 1)
    3296              : 
    3297          882 :       CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
    3298        29975 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    3299              : 
    3300              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
    3301        29093 :                                 iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
    3302        29093 :          IF (do_kpoints_prv) THEN
    3303        37009 :             IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
    3304              :                 ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
    3305         5247 :             img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
    3306         5247 :             IF (img > nimg .OR. img < 1) CYCLE
    3307              :          ELSE
    3308              :             img = 1
    3309              :          END IF
    3310              : 
    3311              :          CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
    3312              :                                 npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
    3313        28985 :                                 sphi=sphi_i, zet=zeti)
    3314              : 
    3315              :          CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
    3316              :                                 npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
    3317        28985 :                                 sphi=sphi_j, zet=zetj)
    3318              : 
    3319        28985 :          IF (do_symmetric) THEN
    3320        25105 :             IF (iatom <= jatom) THEN
    3321        14962 :                irow = iatom
    3322        14962 :                icol = jatom
    3323              :             ELSE
    3324        10143 :                irow = jatom
    3325        10143 :                icol = iatom
    3326              :             END IF
    3327              :          ELSE
    3328         3880 :             irow = iatom
    3329         3880 :             icol = jatom
    3330              :          END IF
    3331              : 
    3332       115940 :          dab = NORM2(rij)
    3333              : 
    3334              :          CALL dbcsr_get_block_p(matrix=t2c(img), &
    3335        28985 :                                 row=irow, col=icol, BLOCK=block_t%block, found=found)
    3336        28985 :          CPASSERT(found)
    3337        28985 :          trans = do_symmetric .AND. (iatom > jatom)
    3338              : 
    3339       110025 :          DO iset = 1, nseti
    3340              : 
    3341        80158 :             ncoi = npgfi(iset)*ncoset(lmax_i(iset))
    3342        80158 :             sgfi = first_sgf_i(1, iset)
    3343              : 
    3344       507394 :             DO jset = 1, nsetj
    3345              : 
    3346       398195 :                ncoj = npgfj(jset)*ncoset(lmax_j(jset))
    3347       398195 :                sgfj = first_sgf_j(1, jset)
    3348              : 
    3349       478353 :                IF (ncoi*ncoj > 0) THEN
    3350      1592780 :                   ALLOCATE (sij_contr(nsgfi(iset), nsgfj(jset)))
    3351     79646473 :                   sij_contr(:, :) = 0.0_dp
    3352              : 
    3353      1592780 :                   ALLOCATE (sij(ncoi, ncoj))
    3354    205995414 :                   sij(:, :) = 0.0_dp
    3355              : 
    3356       398195 :                   ri = 0.0_dp
    3357       398195 :                   rj = rij
    3358              : 
    3359              :                   CALL eri_2center(sij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
    3360              :                                    rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
    3361       398195 :                                    rpgf_j(:, jset), rj, dab, lib, potential_parameter)
    3362              : 
    3363              :                   CALL ab_contract(sij_contr, sij, &
    3364              :                                    sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
    3365       398195 :                                    ncoi, ncoj, nsgfi(iset), nsgfj(jset))
    3366              : 
    3367       398195 :                   DEALLOCATE (sij)
    3368       398195 :                   IF (trans) THEN
    3369       492840 :                      ALLOCATE (sij_rs(nsgfj(jset), nsgfi(iset)))
    3370     30206681 :                      sij_rs(:, :) = TRANSPOSE(sij_contr)
    3371              :                   ELSE
    3372      1099940 :                      ALLOCATE (sij_rs(nsgfi(iset), nsgfj(jset)))
    3373     49354022 :                      sij_rs(:, :) = sij_contr
    3374              :                   END IF
    3375              : 
    3376       398195 :                   DEALLOCATE (sij_contr)
    3377              : 
    3378              :                   ! RI regularization
    3379              :                   IF (.NOT. do_hfx_kpoints_prv .AND. PRESENT(regularization_RI) .AND. &
    3380              :                       iatom == jatom .AND. iset == jset .AND. &
    3381       398195 :                       cell_j(1) == 0 .AND. cell_j(2) == 0 .AND. cell_j(3) == 0) THEN
    3382        11990 :                      DO i_diag = 1, nsgfi(iset)
    3383        26016 :                         min_zet = MINVAL(zeti(:, iset))
    3384         8672 :                         CPASSERT(min_zet > 1.0E-10_dp)
    3385              :                         sij_rs(i_diag, i_diag) = sij_rs(i_diag, i_diag) + &
    3386        11990 :                                                  regularization_RI*MAX(1.0_dp, 1.0_dp/min_zet)
    3387              :                      END DO
    3388              :                   END IF
    3389              : 
    3390              :                   CALL block_add("IN", sij_rs, &
    3391              :                                  nsgfi(iset), nsgfj(jset), block_t%block, &
    3392       398195 :                                  sgfi, sgfj, trans=trans)
    3393       398195 :                   DEALLOCATE (sij_rs)
    3394              : 
    3395              :                END IF
    3396              :             END DO
    3397              :          END DO
    3398              :       END DO
    3399              : 
    3400          882 :       CALL cp_libint_cleanup_2eri(lib)
    3401              : 
    3402          882 :       CALL neighbor_list_iterator_release(nl_iterator)
    3403         7038 :       DO img = 1, nimg
    3404         6156 :          CALL dbcsr_finalize(t2c(img))
    3405         7038 :          CALL dbcsr_filter(t2c(img), filter_eps)
    3406              :       END DO
    3407              : 
    3408          882 :       CALL timestop(handle)
    3409              : 
    3410         1764 :    END SUBROUTINE build_2c_integrals
    3411              : 
    3412              : ! **************************************************************************************************
    3413              : !> \brief ...
    3414              : !> \param tensor tensor with data. Data is cleared after compression.
    3415              : !> \param blk_indices ...
    3416              : !> \param compressed compressed tensor data
    3417              : !> \param eps all entries < eps are discarded
    3418              : !> \param memory ...
    3419              : ! **************************************************************************************************
    3420        22202 :    SUBROUTINE compress_tensor(tensor, blk_indices, compressed, eps, memory)
    3421              :       TYPE(dbt_type), INTENT(INOUT)                      :: tensor
    3422              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), &
    3423              :          INTENT(INOUT)                                   :: blk_indices
    3424              :       TYPE(hfx_compression_type), INTENT(INOUT)          :: compressed
    3425              :       REAL(dp), INTENT(IN)                               :: eps
    3426              :       REAL(dp), INTENT(INOUT)                            :: memory
    3427              : 
    3428              :       INTEGER                                            :: buffer_left, buffer_size, buffer_start, &
    3429              :                                                             i, iblk, memory_usage, nbits, nblk, &
    3430              :                                                             nints, offset, shared_offset
    3431              :       INTEGER(int_8)                                     :: estimate_to_store_int, &
    3432              :                                                             storage_counter_integrals
    3433              :       INTEGER, DIMENSION(3)                              :: ind
    3434              :       LOGICAL                                            :: found
    3435              :       REAL(dp)                                           :: spherical_estimate
    3436        22202 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET  :: blk_data
    3437        22202 :       REAL(dp), DIMENSION(:), POINTER                    :: blk_data_1d
    3438              :       TYPE(dbt_iterator_type)                            :: iter
    3439        22202 :       TYPE(hfx_cache_type), DIMENSION(:), POINTER        :: integral_caches
    3440              :       TYPE(hfx_cache_type), POINTER                      :: maxval_cache
    3441        22202 :       TYPE(hfx_container_type), DIMENSION(:), POINTER    :: integral_containers
    3442              :       TYPE(hfx_container_type), POINTER                  :: maxval_container
    3443              : 
    3444        22202 :       CALL dealloc_containers(compressed, memory_usage)
    3445        22202 :       CALL alloc_containers(compressed, 1)
    3446              : 
    3447        22202 :       maxval_container => compressed%maxval_container(1)
    3448        22202 :       integral_containers => compressed%integral_containers(:, 1)
    3449              : 
    3450        22202 :       CALL hfx_init_container(maxval_container, memory_usage, .FALSE.)
    3451      1443130 :       DO i = 1, 64
    3452      1443130 :          CALL hfx_init_container(integral_containers(i), memory_usage, .FALSE.)
    3453              :       END DO
    3454              : 
    3455        22202 :       maxval_cache => compressed%maxval_cache(1)
    3456        22202 :       integral_caches => compressed%integral_caches(:, 1)
    3457              : 
    3458        22202 :       IF (ALLOCATED(blk_indices)) DEALLOCATE (blk_indices)
    3459        59114 :       ALLOCATE (blk_indices(dbt_get_num_blocks(tensor), 3))
    3460        22202 :       shared_offset = 0
    3461              : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices,shared_offset) &
    3462        22202 : !$OMP PRIVATE(iter,ind,offset,nblk,iblk)
    3463              :       CALL dbt_iterator_start(iter, tensor)
    3464              :       nblk = dbt_iterator_num_blocks(iter)
    3465              : !$OMP CRITICAL
    3466              :       offset = shared_offset
    3467              :       shared_offset = shared_offset + nblk
    3468              : !$OMP END CRITICAL
    3469              :       DO iblk = 1, nblk
    3470              :          CALL dbt_iterator_next_block(iter, ind)
    3471              :          blk_indices(offset + iblk, :) = ind(:)
    3472              :       END DO
    3473              :       CALL dbt_iterator_stop(iter)
    3474              : !$OMP END PARALLEL
    3475              : 
    3476              :       ! Can not use the tensor iterator here because the order of the blocks is not guaranteed.
    3477       444593 :       DO i = 1, SIZE(blk_indices, 1)
    3478      1689564 :          ind = blk_indices(i, :)
    3479       422391 :          CALL dbt_get_block(tensor, ind, blk_data, found)
    3480       422391 :          CPASSERT(found)
    3481      1689564 :          nints = SIZE(blk_data)
    3482       422391 :          blk_data_1d(1:nints) => blk_data
    3483    220585911 :          spherical_estimate = MAXVAL(ABS(blk_data_1d))
    3484       422391 :          IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate)
    3485       422391 :          estimate_to_store_int = EXPONENT(spherical_estimate)
    3486       422391 :          estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
    3487              : 
    3488              :          CALL hfx_add_single_cache_element(estimate_to_store_int, 6, &
    3489              :                                            maxval_cache, maxval_container, memory_usage, &
    3490       422391 :                                            .FALSE.)
    3491              : 
    3492       422391 :          spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
    3493              : 
    3494       422391 :          nbits = EXPONENT(ANINT(spherical_estimate/eps)) + 1
    3495       422391 :          IF (nbits > 64) THEN
    3496              :             CALL cp_abort(__LOCATION__, &
    3497            0 :                           "Overflow during tensor compression. Please use a larger EPS_FILTER or EPS_STORAGE_SCALING")
    3498              :          END IF
    3499              : 
    3500              :          buffer_left = nints
    3501              :          buffer_start = 1
    3502              : 
    3503       940678 :          DO WHILE (buffer_left > 0)
    3504       518287 :             buffer_size = MIN(buffer_left, cache_size)
    3505              :             CALL hfx_add_mult_cache_elements(blk_data_1d(buffer_start:), &
    3506              :                                              buffer_size, nbits, &
    3507              :                                              integral_caches(nbits), &
    3508              :                                              integral_containers(nbits), &
    3509              :                                              eps, 1.0_dp, &
    3510              :                                              memory_usage, &
    3511       518287 :                                              .FALSE.)
    3512       518287 :             buffer_left = buffer_left - buffer_size
    3513       518287 :             buffer_start = buffer_start + buffer_size
    3514              :          END DO
    3515              : 
    3516       866984 :          NULLIFY (blk_data_1d); DEALLOCATE (blk_data)
    3517              :       END DO
    3518              : 
    3519        22202 :       CALL dbt_clear(tensor)
    3520              : 
    3521        22202 :       storage_counter_integrals = memory_usage*cache_size
    3522        22202 :       memory = memory + REAL(storage_counter_integrals, dp)/1024/128
    3523              :       !WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
    3524              :       !   "HFX_MEM_INFO| Total memory consumption ERI's RAM [MiB]:            ", memory
    3525              : 
    3526              :       CALL hfx_flush_last_cache(6, maxval_cache, maxval_container, memory_usage, &
    3527        22202 :                                 .FALSE.)
    3528      1443130 :       DO i = 1, 64
    3529              :          CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
    3530      1443130 :                                    memory_usage, .FALSE.)
    3531              :       END DO
    3532              : 
    3533        22202 :       CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .FALSE.)
    3534      1443130 :       DO i = 1, 64
    3535              :          CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
    3536      1443130 :                                             memory_usage, .FALSE.)
    3537              :       END DO
    3538              : 
    3539        44404 :    END SUBROUTINE
    3540              : 
    3541              : ! **************************************************************************************************
    3542              : !> \brief ...
    3543              : !> \param tensor empty tensor which is filled by decompressed data
    3544              : !> \param blk_indices indices of blocks to be reserved
    3545              : !> \param compressed compressed data
    3546              : !> \param eps all entries < eps are discarded
    3547              : ! **************************************************************************************************
    3548        63866 :    SUBROUTINE decompress_tensor(tensor, blk_indices, compressed, eps)
    3549              : 
    3550              :       TYPE(dbt_type), INTENT(INOUT)                      :: tensor
    3551              :       INTEGER, DIMENSION(:, :)                           :: blk_indices
    3552              :       TYPE(hfx_compression_type), INTENT(INOUT)          :: compressed
    3553              :       REAL(dp), INTENT(IN)                               :: eps
    3554              : 
    3555              :       INTEGER                                            :: A, b, buffer_left, buffer_size, &
    3556              :                                                             buffer_start, i, memory_usage, nbits, &
    3557              :                                                             nblk_per_thread, nints
    3558              :       INTEGER(int_8)                                     :: estimate_to_store_int
    3559              :       INTEGER, DIMENSION(3)                              :: blk_size, ind
    3560              :       REAL(dp)                                           :: spherical_estimate
    3561        63866 :       REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET        :: blk_data
    3562        63866 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: blk_data_3d
    3563        63866 :       TYPE(hfx_cache_type), DIMENSION(:), POINTER        :: integral_caches
    3564              :       TYPE(hfx_cache_type), POINTER                      :: maxval_cache
    3565        63866 :       TYPE(hfx_container_type), DIMENSION(:), POINTER    :: integral_containers
    3566              :       TYPE(hfx_container_type), POINTER                  :: maxval_container
    3567              : 
    3568        63866 :       maxval_cache => compressed%maxval_cache(1)
    3569        63866 :       maxval_container => compressed%maxval_container(1)
    3570        63866 :       integral_caches => compressed%integral_caches(:, 1)
    3571        63866 :       integral_containers => compressed%integral_containers(:, 1)
    3572              : 
    3573        63866 :       memory_usage = 0
    3574              : 
    3575        63866 :       CALL hfx_decompress_first_cache(6, maxval_cache, maxval_container, memory_usage, .FALSE.)
    3576              : 
    3577      4151290 :       DO i = 1, 64
    3578              :          CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), &
    3579      4151290 :                                          memory_usage, .FALSE.)
    3580              :       END DO
    3581              : 
    3582              : !TODO: Parallelize creation of block list.
    3583        63866 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices) PRIVATE(nblk_per_thread,A,b)
    3584              :       nblk_per_thread = SIZE(blk_indices, 1)/omp_get_num_threads() + 1
    3585              :       a = omp_get_thread_num()*nblk_per_thread + 1
    3586              :       b = MIN(a + nblk_per_thread, SIZE(blk_indices, 1))
    3587              :       CALL dbt_reserve_blocks(tensor, blk_indices(a:b, :))
    3588              : !$OMP END PARALLEL
    3589              : 
    3590              :       ! Can not use the tensor iterator here because the order of the blocks is not guaranteed.
    3591      1517272 :       DO i = 1, SIZE(blk_indices, 1)
    3592      5813624 :          ind = blk_indices(i, :)
    3593      1453406 :          CALL dbt_blk_sizes(tensor, ind, blk_size)
    3594      5813624 :          nints = PRODUCT(blk_size)
    3595              :          CALL hfx_get_single_cache_element( &
    3596              :             estimate_to_store_int, 6, &
    3597              :             maxval_cache, maxval_container, memory_usage, &
    3598      1453406 :             .FALSE.)
    3599              : 
    3600      1453406 :          spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
    3601              : 
    3602      1453406 :          nbits = EXPONENT(ANINT(spherical_estimate/eps)) + 1
    3603              : 
    3604      1453406 :          buffer_left = nints
    3605      1453406 :          buffer_start = 1
    3606              : 
    3607      4360218 :          ALLOCATE (blk_data(nints))
    3608      3095835 :          DO WHILE (buffer_left > 0)
    3609      1642429 :             buffer_size = MIN(buffer_left, cache_size)
    3610              :             CALL hfx_get_mult_cache_elements(blk_data(buffer_start), &
    3611              :                                              buffer_size, nbits, &
    3612              :                                              integral_caches(nbits), &
    3613              :                                              integral_containers(nbits), &
    3614              :                                              eps, 1.0_dp, &
    3615              :                                              memory_usage, &
    3616      1642429 :                                              .FALSE.)
    3617      1642429 :             buffer_left = buffer_left - buffer_size
    3618      1642429 :             buffer_start = buffer_start + buffer_size
    3619              :          END DO
    3620              : 
    3621      1453406 :          blk_data_3d(1:blk_size(1), 1:blk_size(2), 1:blk_size(3)) => blk_data
    3622      1453406 :          CALL dbt_put_block(tensor, ind, blk_size, blk_data_3d)
    3623      1517272 :          NULLIFY (blk_data_3d); DEALLOCATE (blk_data)
    3624              :       END DO
    3625              : 
    3626        63866 :       CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .FALSE.)
    3627      4151290 :       DO i = 1, 64
    3628              :          CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
    3629      4151290 :                                             memory_usage, .FALSE.)
    3630              :       END DO
    3631       127732 :    END SUBROUTINE
    3632              : 
    3633              : ! **************************************************************************************************
    3634              : !> \brief ...
    3635              : !> \param tensor ...
    3636              : !> \param nze ...
    3637              : !> \param occ ...
    3638              : ! **************************************************************************************************
    3639       139344 :    SUBROUTINE get_tensor_occupancy(tensor, nze, occ)
    3640              :       TYPE(dbt_type), INTENT(IN)                         :: tensor
    3641              :       INTEGER(int_8), INTENT(OUT)                        :: nze
    3642              :       REAL(dp), INTENT(OUT)                              :: occ
    3643              : 
    3644       278688 :       INTEGER, DIMENSION(dbt_ndims(tensor))              :: dims
    3645              : 
    3646       139344 :       nze = dbt_get_nze_total(tensor)
    3647       139344 :       CALL dbt_get_info(tensor, nfull_total=dims)
    3648       529578 :       occ = REAL(nze, dp)/PRODUCT(REAL(dims, dp))
    3649              : 
    3650       139344 :    END SUBROUTINE
    3651              : 
    3652            0 : END MODULE
        

Generated by: LCOV version 2.0-1