LCOV - code coverage report
Current view: top level - src - mp2_integrals.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 520 544 95.6 %
Date: 2025-02-18 08:24:35 Functions: 7 9 77.8 %

          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 Routines to calculate and distribute 2c- and 3c- integrals for RI
      10             : !> \par History
      11             : !>      06.2012 created [Mauro Del Ben]
      12             : !>      03.2019 separated from mp2_ri_gpw [Frederick Stein]
      13             : ! **************************************************************************************************
      14             : MODULE mp2_integrals
      15             :    USE OMP_LIB,                         ONLY: omp_get_num_threads,&
      16             :                                               omp_get_thread_num
      17             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      18             :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      19             :                                               gto_basis_set_type
      20             :    USE bibliography,                    ONLY: DelBen2013,&
      21             :                                               cite_reference
      22             :    USE cell_types,                      ONLY: cell_type,&
      23             :                                               get_cell
      24             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      25             :    USE cp_control_types,                ONLY: dft_control_type
      26             :    USE cp_dbcsr_api,                    ONLY: &
      27             :         dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
      28             :         dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      29             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      30             :                                               cp_dbcsr_m_by_n_from_template
      31             :    USE cp_eri_mme_interface,            ONLY: cp_eri_mme_param,&
      32             :                                               cp_eri_mme_set_params
      33             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      34             :                                               cp_fm_struct_release,&
      35             :                                               cp_fm_struct_type
      36             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      37             :                                               cp_fm_get_info,&
      38             :                                               cp_fm_release,&
      39             :                                               cp_fm_type
      40             :    USE cp_log_handling,                 ONLY: cp_to_string
      41             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      42             :    USE dbt_api,                         ONLY: &
      43             :         dbt_clear, dbt_contract, dbt_copy, dbt_create, dbt_destroy, dbt_distribution_destroy, &
      44             :         dbt_distribution_new, dbt_distribution_type, dbt_filter, dbt_get_block, dbt_get_info, &
      45             :         dbt_get_stored_coordinates, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
      46             :         dbt_pgrid_type, dbt_put_block, dbt_reserve_blocks, dbt_scale, dbt_split_blocks, dbt_type
      47             :    USE group_dist_types,                ONLY: create_group_dist,&
      48             :                                               get_group_dist,&
      49             :                                               group_dist_d1_type
      50             :    USE hfx_types,                       ONLY: alloc_containers,&
      51             :                                               block_ind_type,&
      52             :                                               hfx_compression_type
      53             :    USE input_constants,                 ONLY: &
      54             :         do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, &
      55             :         do_potential_long, do_potential_short, do_potential_truncated, kp_weights_W_auto, &
      56             :         kp_weights_W_tailored, kp_weights_W_uniform
      57             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      58             :                                               section_vals_type,&
      59             :                                               section_vals_val_get
      60             :    USE kinds,                           ONLY: default_string_length,&
      61             :                                               dp,&
      62             :                                               int_8
      63             :    USE kpoint_methods,                  ONLY: kpoint_init_cell_index
      64             :    USE kpoint_types,                    ONLY: kpoint_type
      65             :    USE libint_2c_3c,                    ONLY: compare_potential_types,&
      66             :                                               libint_potential_type
      67             :    USE machine,                         ONLY: m_flush
      68             :    USE message_passing,                 ONLY: mp_cart_type,&
      69             :                                               mp_comm_type,&
      70             :                                               mp_para_env_type
      71             :    USE mp2_eri,                         ONLY: mp2_eri_3c_integrate
      72             :    USE mp2_eri_gpw,                     ONLY: cleanup_gpw,&
      73             :                                               mp2_eri_3c_integrate_gpw,&
      74             :                                               prepare_gpw
      75             :    USE mp2_ri_2c,                       ONLY: get_2c_integrals
      76             :    USE mp2_types,                       ONLY: three_dim_real_array
      77             :    USE particle_methods,                ONLY: get_particle_set
      78             :    USE particle_types,                  ONLY: particle_type
      79             :    USE pw_env_types,                    ONLY: pw_env_type
      80             :    USE pw_poisson_types,                ONLY: pw_poisson_type
      81             :    USE pw_pool_types,                   ONLY: pw_pool_type
      82             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      83             :                                               pw_r3d_rs_type
      84             :    USE qs_environment_types,            ONLY: get_qs_env,&
      85             :                                               qs_environment_type,&
      86             :                                               set_qs_env
      87             :    USE qs_integral_utils,               ONLY: basis_set_list_setup
      88             :    USE qs_interactions,                 ONLY: init_interaction_radii_orb_basis
      89             :    USE qs_kind_types,                   ONLY: qs_kind_type
      90             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      91             :    USE qs_tensors,                      ONLY: build_3c_integrals,&
      92             :                                               build_3c_neighbor_lists,&
      93             :                                               compress_tensor,&
      94             :                                               get_tensor_occupancy,&
      95             :                                               neighbor_list_3c_destroy
      96             :    USE qs_tensors_types,                ONLY: create_3c_tensor,&
      97             :                                               create_tensor_batches,&
      98             :                                               distribution_3d_create,&
      99             :                                               distribution_3d_type,&
     100             :                                               neighbor_list_3c_type,&
     101             :                                               pgf_block_sizes
     102             :    USE task_list_types,                 ONLY: task_list_type
     103             :    USE util,                            ONLY: get_limit
     104             : #include "./base/base_uses.f90"
     105             : 
     106             :    IMPLICIT NONE
     107             : 
     108             :    PRIVATE
     109             : 
     110             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_integrals'
     111             : 
     112             :    PUBLIC :: mp2_ri_gpw_compute_in, compute_kpoints
     113             : 
     114             :    TYPE intermediate_matrix_type
     115             :       TYPE(dbcsr_type) :: matrix_ia_jnu, matrix_ia_jb
     116             :       INTEGER :: max_row_col_local = 0
     117             :       INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info
     118             :       TYPE(cp_fm_type) :: fm_BIb_jb = cp_fm_type()
     119             :       CHARACTER(LEN=default_string_length) :: descr = ""
     120             :    END TYPE intermediate_matrix_type
     121             : 
     122             : CONTAINS
     123             : 
     124             : ! **************************************************************************************************
     125             : !> \brief with ri mp2 gpw
     126             : !> \param BIb_C ...
     127             : !> \param BIb_C_gw ...
     128             : !> \param BIb_C_bse_ij ...
     129             : !> \param BIb_C_bse_ab ...
     130             : !> \param gd_array ...
     131             : !> \param gd_B_virtual ...
     132             : !> \param dimen_RI ...
     133             : !> \param dimen_RI_red ...
     134             : !> \param qs_env ...
     135             : !> \param para_env ...
     136             : !> \param para_env_sub ...
     137             : !> \param color_sub ...
     138             : !> \param cell ...
     139             : !> \param particle_set ...
     140             : !> \param atomic_kind_set ...
     141             : !> \param qs_kind_set ...
     142             : !> \param fm_matrix_PQ ...
     143             : !> \param fm_matrix_L_kpoints ...
     144             : !> \param fm_matrix_Minv_L_kpoints ...
     145             : !> \param fm_matrix_Minv ...
     146             : !> \param fm_matrix_Minv_Vtrunc_Minv ...
     147             : !> \param nmo ...
     148             : !> \param homo ...
     149             : !> \param mat_munu ...
     150             : !> \param sab_orb_sub ...
     151             : !> \param mo_coeff_o ...
     152             : !> \param mo_coeff_v ...
     153             : !> \param mo_coeff_all ...
     154             : !> \param mo_coeff_gw ...
     155             : !> \param mo_coeff_o_bse ...
     156             : !> \param mo_coeff_v_bse ...
     157             : !> \param eps_filter ...
     158             : !> \param unit_nr ...
     159             : !> \param mp2_memory ...
     160             : !> \param calc_PQ_cond_num ...
     161             : !> \param calc_forces ...
     162             : !> \param blacs_env_sub ...
     163             : !> \param my_do_gw ...
     164             : !> \param do_bse ...
     165             : !> \param gd_B_all ...
     166             : !> \param starts_array_mc ...
     167             : !> \param ends_array_mc ...
     168             : !> \param starts_array_mc_block ...
     169             : !> \param ends_array_mc_block ...
     170             : !> \param gw_corr_lev_occ ...
     171             : !> \param gw_corr_lev_virt ...
     172             : !> \param bse_lev_virt ...
     173             : !> \param do_im_time ...
     174             : !> \param do_kpoints_cubic_RPA ...
     175             : !> \param kpoints ...
     176             : !> \param t_3c_M ...
     177             : !> \param t_3c_O ...
     178             : !> \param t_3c_O_compressed ...
     179             : !> \param t_3c_O_ind ...
     180             : !> \param ri_metric ...
     181             : !> \param gd_B_occ_bse ...
     182             : !> \param gd_B_virt_bse ...
     183             : !> \author Mauro Del Ben
     184             : ! **************************************************************************************************
     185        5264 :    SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, &
     186             :                                     dimen_RI, dimen_RI_red, qs_env, para_env, para_env_sub, color_sub, &
     187             :                                     cell, particle_set, atomic_kind_set, qs_kind_set, &
     188             :                                     fm_matrix_PQ, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
     189             :                                     fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
     190         658 :                                     nmo, homo, mat_munu, &
     191        1316 :                                     sab_orb_sub, mo_coeff_o, mo_coeff_v, mo_coeff_all, &
     192         658 :                                     mo_coeff_gw, mo_coeff_o_bse, mo_coeff_v_bse, eps_filter, unit_nr, &
     193             :                                     mp2_memory, calc_PQ_cond_num, calc_forces, blacs_env_sub, my_do_gw, do_bse, &
     194             :                                     gd_B_all, starts_array_mc, ends_array_mc, &
     195             :                                     starts_array_mc_block, ends_array_mc_block, &
     196             :                                     gw_corr_lev_occ, gw_corr_lev_virt, &
     197             :                                     bse_lev_virt, &
     198             :                                     do_im_time, do_kpoints_cubic_RPA, kpoints, &
     199             :                                     t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
     200             :                                     ri_metric, gd_B_occ_bse, gd_B_virt_bse)
     201             : 
     202             :       TYPE(three_dim_real_array), ALLOCATABLE, &
     203             :          DIMENSION(:), INTENT(OUT)                       :: BIb_C, BIb_C_gw
     204             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
     205             :          INTENT(OUT)                                     :: BIb_C_bse_ij, BIb_C_bse_ab
     206             :       TYPE(group_dist_d1_type), INTENT(OUT)              :: gd_array
     207             :       TYPE(group_dist_d1_type), ALLOCATABLE, &
     208             :          DIMENSION(:), INTENT(OUT)                       :: gd_B_virtual
     209             :       INTEGER, INTENT(OUT)                               :: dimen_RI, dimen_RI_red
     210             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     211             :       TYPE(mp_para_env_type), POINTER                    :: para_env, para_env_sub
     212             :       INTEGER, INTENT(IN)                                :: color_sub
     213             :       TYPE(cell_type), POINTER                           :: cell
     214             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     215             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     216             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     217             :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_matrix_PQ
     218             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_matrix_L_kpoints, &
     219             :                                                             fm_matrix_Minv_L_kpoints, &
     220             :                                                             fm_matrix_Minv, &
     221             :                                                             fm_matrix_Minv_Vtrunc_Minv
     222             :       INTEGER, INTENT(IN)                                :: nmo
     223             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
     224             :       TYPE(dbcsr_p_type), INTENT(INOUT)                  :: mat_munu
     225             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     226             :          INTENT(IN), POINTER                             :: sab_orb_sub
     227             :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN)       :: mo_coeff_o, mo_coeff_v, mo_coeff_all, &
     228             :                                                             mo_coeff_gw, mo_coeff_o_bse, &
     229             :                                                             mo_coeff_v_bse
     230             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
     231             :       INTEGER, INTENT(IN)                                :: unit_nr
     232             :       REAL(KIND=dp), INTENT(IN)                          :: mp2_memory
     233             :       LOGICAL, INTENT(IN)                                :: calc_PQ_cond_num, calc_forces
     234             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_sub
     235             :       LOGICAL, INTENT(IN)                                :: my_do_gw, do_bse
     236             :       TYPE(group_dist_d1_type), INTENT(OUT)              :: gd_B_all
     237             :       INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT)    :: starts_array_mc, ends_array_mc, &
     238             :                                                             starts_array_mc_block, &
     239             :                                                             ends_array_mc_block
     240             :       INTEGER, INTENT(IN)                                :: gw_corr_lev_occ, gw_corr_lev_virt, &
     241             :                                                             bse_lev_virt
     242             :       LOGICAL, INTENT(IN)                                :: do_im_time, do_kpoints_cubic_RPA
     243             :       TYPE(kpoint_type), POINTER                         :: kpoints
     244             :       TYPE(dbt_type), INTENT(OUT)                        :: t_3c_M
     245             :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
     246             :          INTENT(OUT)                                     :: t_3c_O
     247             :       TYPE(hfx_compression_type), ALLOCATABLE, &
     248             :          DIMENSION(:, :, :), INTENT(INOUT)               :: t_3c_O_compressed
     249             :       TYPE(block_ind_type), ALLOCATABLE, &
     250             :          DIMENSION(:, :, :)                              :: t_3c_O_ind
     251             :       TYPE(libint_potential_type), INTENT(IN)            :: ri_metric
     252             :       TYPE(group_dist_d1_type), INTENT(OUT)              :: gd_B_occ_bse, gd_B_virt_bse
     253             : 
     254             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_ri_gpw_compute_in'
     255             : 
     256             :       INTEGER :: cm, cut_memory, cut_memory_int, eri_method, gw_corr_lev_total, handle, handle2, &
     257             :          handle4, i, i_counter, i_mem, ibasis, ispin, itmp(2), j, jcell, kcell, LLL, min_bsize, &
     258             :          my_B_all_end, my_B_all_size, my_B_all_start, my_B_occ_bse_end, my_B_occ_bse_size, &
     259             :          my_B_occ_bse_start, my_B_virt_bse_end, my_B_virt_bse_size, my_B_virt_bse_start, &
     260             :          my_group_L_end, my_group_L_size, my_group_L_start, n_rep, natom, ngroup, nimg, nkind, &
     261             :          nspins, potential_type, ri_metric_type
     262             :       INTEGER(int_8)                                     :: nze
     263         658 :       INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI, &
     264         658 :          ends_array_mc_block_int, ends_array_mc_int, my_B_size, my_B_virtual_end, &
     265        1316 :          my_B_virtual_start, sizes_AO, sizes_AO_split, sizes_RI, sizes_RI_split, &
     266        1316 :          starts_array_mc_block_int, starts_array_mc_int, virtual
     267             :       INTEGER, DIMENSION(2, 3)                           :: bounds
     268             :       INTEGER, DIMENSION(3)                              :: bounds_3c, pcoord, pdims, pdims_t3c, &
     269             :                                                             periodic
     270             :       LOGICAL                                            :: do_gpw, do_kpoints_from_Gamma, do_svd, &
     271             :                                                             memory_info
     272             :       REAL(KIND=dp) :: compression_factor, cutoff_old, eps_pgf_orb, eps_pgf_orb_old, eps_svd, &
     273             :          mem_for_abK, mem_for_iaK, mem_for_ijK, memory_3c, occ, omega_pot, rc_ang, &
     274             :          relative_cutoff_old
     275         658 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: e_cutoff_old
     276         658 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: my_Lrows, my_Vrows
     277             :       TYPE(cp_eri_mme_param), POINTER                    :: eri_param
     278         658 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_munu_local_L
     279        3290 :       TYPE(dbt_pgrid_type)                               :: pgrid_t3c_M, pgrid_t3c_overl
     280        8554 :       TYPE(dbt_type)                                     :: t_3c_overl_int_template, t_3c_tmp
     281         658 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_overl_int
     282             :       TYPE(dft_control_type), POINTER                    :: dft_control
     283             :       TYPE(distribution_3d_type)                         :: dist_3d
     284         658 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_ao, basis_set_ri_aux
     285             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis, ri_basis
     286         658 :       TYPE(intermediate_matrix_type)                     :: intermed_mat_bse_ab, intermed_mat_bse_ij
     287             :       TYPE(intermediate_matrix_type), ALLOCATABLE, &
     288         658 :          DIMENSION(:)                                    :: intermed_mat, intermed_mat_gw
     289         658 :       TYPE(mp_cart_type)                                 :: mp_comm_t3c_2
     290             :       TYPE(neighbor_list_3c_type)                        :: nl_3c
     291             :       TYPE(pw_c1d_gs_type)                               :: pot_g, rho_g
     292             :       TYPE(pw_env_type), POINTER                         :: pw_env_sub
     293             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     294             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     295             :       TYPE(pw_r3d_rs_type)                               :: psi_L, rho_r
     296             :       TYPE(section_vals_type), POINTER                   :: qs_section
     297             :       TYPE(task_list_type), POINTER                      :: task_list_sub
     298             : 
     299         658 :       CALL timeset(routineN, handle)
     300             : 
     301         658 :       CALL cite_reference(DelBen2013)
     302             : 
     303         658 :       nspins = SIZE(homo)
     304             : 
     305        1974 :       ALLOCATE (virtual(nspins))
     306        1462 :       virtual(:) = nmo - homo(:)
     307         658 :       gw_corr_lev_total = gw_corr_lev_virt + gw_corr_lev_occ
     308             : 
     309         658 :       eri_method = qs_env%mp2_env%eri_method
     310         658 :       eri_param => qs_env%mp2_env%eri_mme_param
     311         658 :       do_svd = qs_env%mp2_env%do_svd
     312         658 :       eps_svd = qs_env%mp2_env%eps_svd
     313         658 :       potential_type = qs_env%mp2_env%potential_parameter%potential_type
     314         658 :       ri_metric_type = ri_metric%potential_type
     315         658 :       omega_pot = qs_env%mp2_env%potential_parameter%omega
     316             : 
     317             :       ! whether we need gpw integrals (plus pw stuff)
     318             :       do_gpw = (eri_method == do_eri_gpw) .OR. &
     319             :                ((potential_type == do_potential_long .OR. ri_metric_type == do_potential_long) &
     320             :                 .AND. qs_env%mp2_env%eri_method == do_eri_os) &
     321         658 :                .OR. (ri_metric_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)
     322             : 
     323         658 :       IF (do_svd .AND. calc_forces) THEN
     324           0 :          CPABORT("SVD not implemented for forces.!")
     325             :       END IF
     326             : 
     327         658 :       do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
     328         658 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
     329             :          CALL get_qs_env(qs_env=qs_env, &
     330          22 :                          kpoints=kpoints)
     331             :       END IF
     332          22 :       IF (do_kpoints_from_Gamma) THEN
     333          18 :          CALL compute_kpoints(qs_env, kpoints, unit_nr)
     334             :       END IF
     335             : 
     336         658 :       IF (do_bse) THEN
     337          42 :          IF (.NOT. my_do_gw) THEN
     338           0 :             CALL cp_abort(__LOCATION__, "BSE calculations require prior GW calculations.")
     339             :          END IF
     340          42 :          IF (do_im_time) THEN
     341           0 :             CALL cp_abort(__LOCATION__, "BSE calculations are not implemented for low-scaling GW.")
     342             :          END IF
     343             :          ! GPW integrals have to be implemented later
     344          42 :          IF (eri_method == do_eri_gpw) THEN
     345             :             CALL cp_abort(__LOCATION__, &
     346             :                           "BSE calculations are not implemented for GPW integrals. "// &
     347             :                           "This is probably caused by invoking a periodic calculation. "// &
     348           0 :                           "Use PERIODIC NONE for BSE calculations.")
     349             :          END IF
     350             :       END IF
     351             : 
     352         658 :       ngroup = para_env%num_pe/para_env_sub%num_pe
     353             : 
     354             :       ! Preparations for MME method to compute ERIs
     355         658 :       IF (qs_env%mp2_env%eri_method .EQ. do_eri_mme) THEN
     356             :          ! cell might have changed, so we need to reset parameters
     357         126 :          CALL cp_eri_mme_set_params(eri_param, cell, qs_kind_set, basis_type_1="ORB", basis_type_2="RI_AUX", para_env=para_env)
     358             :       END IF
     359             : 
     360         658 :       CALL get_cell(cell=cell, periodic=periodic)
     361             :       ! for minimax Ewald summation, full periodicity is required
     362         658 :       IF (eri_method == do_eri_mme) THEN
     363         126 :          CPASSERT(periodic(1) == 1 .AND. periodic(2) == 1 .AND. periodic(3) == 1)
     364             :       END IF
     365             : 
     366         658 :       IF (do_svd .AND. (do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA)) THEN
     367           0 :          CPABORT("SVD with kpoints not implemented yet!")
     368             :       END IF
     369             : 
     370             :       CALL get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, &
     371             :                             my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, &
     372             :                             kpoints, my_group_L_size, my_group_L_start, my_group_L_end, &
     373             :                             gd_array, calc_PQ_cond_num .AND. .NOT. do_svd, do_svd, eps_svd, &
     374             :                             qs_env%mp2_env%potential_parameter, ri_metric, &
     375             :                             fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
     376             :                             do_im_time, do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA, qs_env%mp2_env%mp2_gpw%eps_pgf_orb_S, &
     377        1952 :                             qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
     378             : 
     379         658 :       IF (unit_nr > 0) THEN
     380             :          ASSOCIATE (ri_metric => qs_env%mp2_env%ri_metric)
     381         578 :             SELECT CASE (ri_metric%potential_type)
     382             :             CASE (do_potential_coulomb)
     383             :                WRITE (unit_nr, FMT="(/T3,A,T74,A)") &
     384         249 :                   "RI_INFO| RI metric: ", "COULOMB"
     385             :             CASE (do_potential_short)
     386             :                WRITE (unit_nr, FMT="(T3,A,T71,A)") &
     387           0 :                   "RI_INFO| RI metric: ", "SHORTRANGE"
     388             :                WRITE (unit_nr, '(T3,A,T61,F20.10)') &
     389           0 :                   "RI_INFO| Omega:     ", ri_metric%omega
     390           0 :                rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
     391             :                WRITE (unit_nr, '(T3,A,T61,F20.10)') &
     392           0 :                   "RI_INFO| Cutoff Radius [angstrom]:     ", rc_ang
     393             :             CASE (do_potential_long)
     394             :                WRITE (unit_nr, FMT="(T3,A,T72,A)") &
     395           8 :                   "RI_INFO| RI metric: ", "LONGRANGE"
     396             :                WRITE (unit_nr, '(T3,A,T61,F20.10)') &
     397           8 :                   "RI_INFO| Omega:     ", ri_metric%omega
     398             :             CASE (do_potential_id)
     399             :                WRITE (unit_nr, FMT="(T3,A,T74,A)") &
     400          40 :                   "RI_INFO| RI metric: ", "OVERLAP"
     401             :             CASE (do_potential_truncated)
     402             :                WRITE (unit_nr, FMT="(T3,A,T64,A)") &
     403          32 :                   "RI_INFO| RI metric: ", "TRUNCATED COULOMB"
     404          32 :                rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
     405             :                WRITE (unit_nr, '(T3,A,T61,F20.2)') &
     406         361 :                   "RI_INFO| Cutoff Radius [angstrom]:     ", rc_ang
     407             :             END SELECT
     408             :          END ASSOCIATE
     409             :       END IF
     410             : 
     411         658 :       IF (calc_forces .AND. .NOT. do_im_time) THEN
     412             :          ! we need (P|Q)^(-1/2) for future use, just save it
     413             :          ! in a fully (home made) distributed way
     414         264 :          itmp = get_limit(dimen_RI, para_env_sub%num_pe, para_env_sub%mepos)
     415         264 :          lll = itmp(2) - itmp(1) + 1
     416        1056 :          ALLOCATE (qs_env%mp2_env%ri_grad%PQ_half(lll, my_group_L_size))
     417      944212 :          qs_env%mp2_env%ri_grad%PQ_half(:, :) = my_Lrows(itmp(1):itmp(2), 1:my_group_L_size)
     418         264 :          IF (.NOT. compare_potential_types(qs_env%mp2_env%ri_metric, qs_env%mp2_env%potential_parameter)) THEN
     419          36 :             ALLOCATE (qs_env%mp2_env%ri_grad%operator_half(lll, my_group_L_size))
     420       41844 :             qs_env%mp2_env%ri_grad%operator_half(:, :) = my_Vrows(itmp(1):itmp(2), 1:my_group_L_size)
     421          12 :             DEALLOCATE (my_Vrows)
     422             :          END IF
     423             :       END IF
     424             : 
     425         658 :       IF (unit_nr > 0) THEN
     426             :          WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
     427         329 :             "RI_INFO| Number of auxiliary basis functions:", dimen_RI, &
     428         329 :             "GENERAL_INFO| Number of basis functions:", nmo, &
     429         329 :             "GENERAL_INFO| Number of occupied orbitals:", homo(1), &
     430         658 :             "GENERAL_INFO| Number of virtual orbitals:", virtual(1)
     431         329 :          IF (do_svd) THEN
     432             :             WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
     433          22 :                "RI_INFO| Reduced auxiliary basis set size:", dimen_RI_red
     434             :          END IF
     435             : 
     436         731 :          mem_for_iaK = dimen_RI*REAL(SUM(homo*virtual), KIND=dp)*8.0_dp/(1024_dp**2)
     437         658 :          mem_for_ijK = dimen_RI*REAL(SUM([homo(1)]**2), KIND=dp)*8.0_dp/(1024_dp**2)
     438         658 :          mem_for_abK = dimen_RI*REAL(SUM([bse_lev_virt]**2), KIND=dp)*8.0_dp/(1024_dp**2)
     439             : 
     440         329 :          IF (.NOT. do_im_time) THEN
     441         261 :             WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ia|K) integrals:', &
     442         522 :                mem_for_iaK, ' MiB'
     443         261 :             IF (my_do_gw .AND. .NOT. do_im_time) THEN
     444          35 :                mem_for_iaK = dimen_RI*REAL(nmo, KIND=dp)*gw_corr_lev_total*8.0_dp/(1024_dp**2)
     445             : 
     446          35 :                WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for G0W0-(nm|K) integrals:', &
     447          70 :                   mem_for_iaK, ' MiB'
     448             :             END IF
     449             :          END IF
     450         329 :          IF (do_bse) THEN
     451          21 :             WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ij|K) integrals:', &
     452          42 :                mem_for_ijK, ' MiB'
     453          21 :             WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ab|K) integrals:', &
     454          42 :                mem_for_abK, ' MiB'
     455             :          END IF
     456         329 :          CALL m_flush(unit_nr)
     457             :       END IF
     458             : 
     459         658 :       CALL para_env%sync() ! sync to see memory output
     460             : 
     461             :       ! in case we do imaginary time, we need the overlap tensor (alpha beta P) or trunc. Coulomb tensor
     462         658 :       IF (.NOT. do_im_time) THEN
     463             : 
     464        3886 :          ALLOCATE (gd_B_virtual(nspins), intermed_mat(nspins))
     465        2088 :          ALLOCATE (my_B_virtual_start(nspins), my_B_virtual_end(nspins), my_B_size(nspins))
     466        1160 :          DO ispin = 1, nspins
     467             : 
     468             :             CALL create_intermediate_matrices(intermed_mat(ispin), mo_coeff_o(ispin)%matrix, virtual(ispin), homo(ispin), &
     469         638 :                                               TRIM(ADJUSTL(cp_to_string(ispin))), blacs_env_sub, para_env_sub)
     470             : 
     471         638 :             CALL create_group_dist(gd_B_virtual(ispin), para_env_sub%num_pe, virtual(ispin))
     472             :             CALL get_group_dist(gd_B_virtual(ispin), para_env_sub%mepos, my_B_virtual_start(ispin), my_B_virtual_end(ispin), &
     473        1160 :                                 my_B_size(ispin))
     474             : 
     475             :          END DO
     476             : 
     477             :          ! in the case of G0W0, we need (K|nm), n,m may be occ or virt (m restricted to corrected levels)
     478         522 :          IF (my_do_gw) THEN
     479             : 
     480         214 :             ALLOCATE (intermed_mat_gw(nspins))
     481         144 :             DO ispin = 1, nspins
     482             :                CALL create_intermediate_matrices(intermed_mat_gw(ispin), mo_coeff_gw(ispin)%matrix, &
     483             :                                                  nmo, gw_corr_lev_total, &
     484             :                                                  "gw_"//TRIM(ADJUSTL(cp_to_string(ispin))), &
     485         144 :                                                  blacs_env_sub, para_env_sub)
     486             : 
     487             :             END DO
     488             : 
     489          70 :             CALL create_group_dist(gd_B_all, para_env_sub%num_pe, nmo)
     490          70 :             CALL get_group_dist(gd_B_all, para_env_sub%mepos, my_B_all_start, my_B_all_end, my_B_all_size)
     491             : 
     492          70 :             IF (do_bse) THEN
     493             :                ! virt x virt matrices
     494             :                CALL create_intermediate_matrices(intermed_mat_bse_ab, mo_coeff_v_bse(1)%matrix, bse_lev_virt, bse_lev_virt, &
     495          42 :                                                  "bse_ab", blacs_env_sub, para_env_sub)
     496             : 
     497          42 :                CALL create_group_dist(gd_B_virt_bse, para_env_sub%num_pe, bse_lev_virt)
     498          42 :                CALL get_group_dist(gd_B_virt_bse, para_env_sub%mepos, my_B_virt_bse_start, my_B_virt_bse_end, my_B_virt_bse_size)
     499             : 
     500             :                ! occ x occ matrices
     501             :                ! We do not implement bse_lev_occ here, because the small number of occupied levels
     502             :                ! does not critically influence the memory
     503             :                CALL create_intermediate_matrices(intermed_mat_bse_ij, mo_coeff_o_bse(1)%matrix, homo(1), homo(1), &
     504          42 :                                                  "bse_ij", blacs_env_sub, para_env_sub)
     505             : 
     506          42 :                CALL create_group_dist(gd_B_occ_bse, para_env_sub%num_pe, homo(1))
     507          42 :                CALL get_group_dist(gd_B_occ_bse, para_env_sub%mepos, my_B_occ_bse_start, my_B_occ_bse_end, my_B_occ_bse_size)
     508             : 
     509             :             END IF
     510             :          END IF
     511             : 
     512             :          ! array that will store the (ia|K) integrals
     513        2204 :          ALLOCATE (BIb_C(nspins))
     514        1160 :          DO ispin = 1, nspins
     515        3184 :             ALLOCATE (BIb_C(ispin)%array(my_group_L_size, my_B_size(ispin), homo(ispin)))
     516     1502432 :             BIb_C(ispin)%array = 0.0_dp
     517             :          END DO
     518             : 
     519             :          ! in the case of GW, we also need (nm|K)
     520         522 :          IF (my_do_gw) THEN
     521             : 
     522         214 :             ALLOCATE (BIb_C_gw(nspins))
     523         144 :             DO ispin = 1, nspins
     524         370 :                ALLOCATE (BIb_C_gw(ispin)%array(my_group_L_size, my_B_all_size, gw_corr_lev_total))
     525     1079572 :                BIb_C_gw(ispin)%array = 0.0_dp
     526             :             END DO
     527             : 
     528             :          END IF
     529             : 
     530         522 :          IF (do_bse) THEN
     531             : 
     532         210 :             ALLOCATE (BIb_C_bse_ij(my_group_L_size, my_B_occ_bse_size, homo(1)))
     533       28770 :             BIb_C_bse_ij = 0.0_dp
     534             : 
     535         210 :             ALLOCATE (BIb_C_bse_ab(my_group_L_size, my_B_virt_bse_size, bse_lev_virt))
     536      257586 :             BIb_C_bse_ab = 0.0_dp
     537             : 
     538             :          END IF
     539             : 
     540         522 :          CALL timeset(routineN//"_loop", handle2)
     541             : 
     542             :          IF (eri_method == do_eri_mme .AND. &
     543         522 :              (ri_metric%potential_type == do_potential_coulomb .OR. ri_metric%potential_type == do_potential_long) .OR. &
     544             :              eri_method == do_eri_os .AND. ri_metric%potential_type == do_potential_coulomb) THEN
     545             : 
     546             :             ! Add a warning for automatically generated RI_AUX basis sets
     547             :             ! Tend to be not sufficiently converged
     548         182 :             IF (qs_env%mp2_env%ri_aux_auto_generated) THEN
     549             :                CALL cp_warn(__LOCATION__, &
     550             :                             "At least one RI_AUX basis set was not explicitly invoked in &KIND-section. "// &
     551             :                             "Automatically RI-basis sets and ERI_METHOD OS tend to be not converged. "// &
     552           0 :                             "Consider specifying BASIS_SET RI_AUX explicitly with a sufficiently large basis.")
     553             :             END IF
     554             : 
     555         182 :             NULLIFY (mat_munu_local_L)
     556        6476 :             ALLOCATE (mat_munu_local_L(my_group_L_size))
     557        6112 :             DO LLL = 1, my_group_L_size
     558        5930 :                NULLIFY (mat_munu_local_L(LLL)%matrix)
     559        5930 :                ALLOCATE (mat_munu_local_L(LLL)%matrix)
     560        5930 :                CALL dbcsr_copy(mat_munu_local_L(LLL)%matrix, mat_munu%matrix)
     561        6112 :                CALL dbcsr_set(mat_munu_local_L(LLL)%matrix, 0.0_dp)
     562             :             END DO
     563             :             CALL mp2_eri_3c_integrate(eri_param, ri_metric, para_env_sub, qs_env, &
     564             :                                       first_c=my_group_L_start, last_c=my_group_L_end, &
     565             :                                       mat_ab=mat_munu_local_L, &
     566             :                                       basis_type_a="ORB", basis_type_b="ORB", &
     567             :                                       basis_type_c="RI_AUX", &
     568         182 :                                       sab_nl=sab_orb_sub, eri_method=eri_method)
     569             : 
     570         370 :             DO ispin = 1, nspins
     571        6476 :                DO LLL = 1, my_group_L_size
     572             :                   CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat(ispin), &
     573             :                                             BIb_C(ispin)%array(LLL, :, :), &
     574             :                                             mo_coeff_o(ispin)%matrix, mo_coeff_v(ispin)%matrix, &
     575             :                                             eps_filter, &
     576        6476 :                                             my_B_virtual_end(ispin), my_B_virtual_start(ispin))
     577             :                END DO
     578             :                CALL contract_B_L(BIb_C(ispin)%array, my_Lrows, gd_B_virtual(ispin)%sizes, &
     579             :                                  gd_array%sizes, qs_env%mp2_env%eri_blksize, &
     580         370 :                                  ngroup, color_sub, para_env, para_env_sub)
     581             :             END DO
     582             : 
     583         182 :             IF (my_do_gw) THEN
     584             : 
     585         144 :                DO ispin = 1, nspins
     586        3089 :                   DO LLL = 1, my_group_L_size
     587             :                      CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat_gw(ispin), &
     588             :                                                BIb_C_gw(ispin)%array(LLL, :, :), &
     589             :                                                mo_coeff_gw(ispin)%matrix, mo_coeff_all(ispin)%matrix, eps_filter, &
     590        3089 :                                                my_B_all_end, my_B_all_start)
     591             :                   END DO
     592             :                   CALL contract_B_L(BIb_C_gw(ispin)%array, my_Lrows, gd_B_all%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, &
     593         144 :                                     ngroup, color_sub, para_env, para_env_sub)
     594             :                END DO
     595             :             END IF
     596             : 
     597         182 :             IF (do_bse) THEN
     598             : 
     599             :                ! B^ab_P matrix elements for BSE
     600        1785 :                DO LLL = 1, my_group_L_size
     601             :                   CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat_bse_ab, &
     602             :                                             BIb_C_bse_ab(LLL, :, :), &
     603             :                                             mo_coeff_v_bse(1)%matrix, mo_coeff_v_bse(1)%matrix, eps_filter, &
     604        1785 :                                             my_B_all_end, my_B_all_start)
     605             :                END DO
     606             :                CALL contract_B_L(BIb_C_bse_ab, my_Lrows, gd_B_virt_bse%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, &
     607          42 :                                  ngroup, color_sub, para_env, para_env_sub)
     608             : 
     609             :                ! B^ij_P matrix elements for BSE
     610        1785 :                DO LLL = 1, my_group_L_size
     611             :                   CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat_bse_ij, &
     612             :                                             BIb_C_bse_ij(LLL, :, :), &
     613             :                                             mo_coeff_o(1)%matrix, mo_coeff_o(1)%matrix, eps_filter, &
     614        1785 :                                             my_B_occ_bse_end, my_B_occ_bse_start)
     615             :                END DO
     616             :                CALL contract_B_L(BIb_C_bse_ij, my_Lrows, gd_B_occ_bse%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, &
     617          42 :                                  ngroup, color_sub, para_env, para_env_sub)
     618             : 
     619             :             END IF
     620             : 
     621        6112 :             DO LLL = 1, my_group_L_size
     622        6112 :                CALL dbcsr_release_p(mat_munu_local_L(LLL)%matrix)
     623             :             END DO
     624         182 :             DEALLOCATE (mat_munu_local_L)
     625             : 
     626         340 :          ELSE IF (do_gpw) THEN
     627             : 
     628             :             CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
     629         340 :                              auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
     630             : 
     631       14739 :             DO i_counter = 1, my_group_L_size
     632             : 
     633             :                CALL mp2_eri_3c_integrate_gpw(psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, &
     634             :                                              particle_set, pw_env_sub, my_Lrows(:, i_counter), poisson_env, rho_r, pot_g, &
     635       14399 :                                              ri_metric, mat_munu, qs_env, task_list_sub)
     636             : 
     637       33524 :                DO ispin = 1, nspins
     638             :                   CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, intermed_mat(ispin), &
     639             :                                             BIb_C(ispin)%array(i_counter, :, :), &
     640             :                                             mo_coeff_o(ispin)%matrix, mo_coeff_v(ispin)%matrix, eps_filter, &
     641       33524 :                                             my_B_virtual_end(ispin), my_B_virtual_start(ispin))
     642             : 
     643             :                END DO
     644             : 
     645       14739 :                IF (my_do_gw) THEN
     646             :                   ! transform (K|mu nu) to (K|nm), n corresponds to corrected GW levels, m is in nmo
     647           0 :                   DO ispin = 1, nspins
     648             :                      CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, intermed_mat_gw(ispin), &
     649             :                                                BIb_C_gw(ispin)%array(i_counter, :, :), &
     650             :                                                mo_coeff_gw(ispin)%matrix, mo_coeff_all(ispin)%matrix, eps_filter, &
     651           0 :                                                my_B_all_end, my_B_all_start)
     652             : 
     653             :                   END DO
     654             :                END IF
     655             : 
     656             :             END DO
     657             : 
     658             :             CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
     659         340 :                              task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
     660             :          ELSE
     661           0 :             CPABORT("Integration method not implemented!")
     662             :          END IF
     663             : 
     664         522 :          CALL timestop(handle2)
     665             : 
     666         522 :          DEALLOCATE (my_Lrows)
     667             : 
     668        1160 :          DO ispin = 1, nspins
     669        1160 :             CALL release_intermediate_matrices(intermed_mat(ispin))
     670             :          END DO
     671        1160 :          DEALLOCATE (intermed_mat)
     672             : 
     673         522 :          IF (my_do_gw) THEN
     674         144 :             DO ispin = 1, nspins
     675         144 :                CALL release_intermediate_matrices(intermed_mat_gw(ispin))
     676             :             END DO
     677         144 :             DEALLOCATE (intermed_mat_gw)
     678             :          END IF
     679             : 
     680        1044 :          IF (do_bse) THEN
     681          42 :             CALL release_intermediate_matrices(intermed_mat_bse_ab)
     682          42 :             CALL release_intermediate_matrices(intermed_mat_bse_ij)
     683             :          END IF
     684             : 
     685             :          ! imag. time = low-scaling SOS-MP2, RPA, GW
     686             :       ELSE
     687             : 
     688         136 :          memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info
     689             : 
     690             :          ! we need 3 tensors:
     691             :          ! 1) t_3c_overl_int: 3c overlap integrals, optimized for easy access to integral blocks
     692             :          !                   (atomic blocks)
     693             :          ! 2) t_3c_O: 3c overlap integrals, optimized for contraction (split blocks)
     694             :          ! 3) t_3c_M: tensor M, optimized for contraction
     695             : 
     696         136 :          CALL get_qs_env(qs_env, natom=natom, nkind=nkind, dft_control=dft_control)
     697             : 
     698         136 :          pdims_t3c = 0
     699         136 :          CALL dbt_pgrid_create(para_env, pdims_t3c, pgrid_t3c_overl)
     700             : 
     701             :          ! set up basis
     702         544 :          ALLOCATE (sizes_RI(natom), sizes_AO(natom))
     703        1080 :          ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
     704         136 :          CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
     705         136 :          CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux)
     706         136 :          CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
     707         136 :          CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao)
     708             : 
     709             :          ! make sure we use the QS%EPS_PGF_ORB
     710         136 :          qs_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS")
     711         136 :          CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
     712         136 :          IF (n_rep /= 0) THEN
     713          82 :             CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=eps_pgf_orb)
     714             :          ELSE
     715          54 :             CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=eps_pgf_orb)
     716          54 :             eps_pgf_orb = SQRT(eps_pgf_orb)
     717             :          END IF
     718         136 :          eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
     719             : 
     720         404 :          DO ibasis = 1, SIZE(basis_set_ao)
     721         268 :             orb_basis => basis_set_ao(ibasis)%gto_basis_set
     722         268 :             CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb)
     723         268 :             ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
     724         404 :             CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb)
     725             :          END DO
     726             : 
     727         136 :          cut_memory_int = qs_env%mp2_env%ri_rpa_im_time%cut_memory
     728             :          CALL create_tensor_batches(sizes_RI, cut_memory_int, starts_array_mc_int, ends_array_mc_int, &
     729         136 :                                     starts_array_mc_block_int, ends_array_mc_block_int)
     730             : 
     731         136 :          DEALLOCATE (starts_array_mc_int, ends_array_mc_int)
     732             : 
     733             :          CALL create_3c_tensor(t_3c_overl_int_template, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, &
     734             :                                sizes_RI, sizes_AO, sizes_AO, map1=[1, 2], map2=[3], &
     735         136 :                                name="O (RI AO | AO)")
     736             : 
     737         136 :          CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
     738         136 :          CALL dbt_mp_environ_pgrid(pgrid_t3c_overl, pdims, pcoord)
     739         136 :          CALL mp_comm_t3c_2%create(pgrid_t3c_overl%mp_comm_2d, 3, pdims)
     740             :          CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, &
     741         136 :                                      nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.)
     742         136 :          DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
     743             : 
     744             :          CALL build_3c_neighbor_lists(nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, &
     745             :                                       dist_3d, ri_metric, "RPA_3c_nl", qs_env, &
     746         136 :                                       sym_jk=.NOT. do_kpoints_cubic_RPA, own_dist=.TRUE.)
     747             : 
     748             :          ! init k points
     749         136 :          IF (do_kpoints_cubic_RPA) THEN
     750             :             ! set up new kpoint type with periodic images according to eps_grid from MP2 section
     751             :             ! instead of eps_pgf_orb from QS section
     752           4 :             CALL kpoint_init_cell_index(kpoints, nl_3c%jk_list, para_env, dft_control)
     753           4 :             IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
     754           2 :                "3C_OVERLAP_INTEGRALS_INFO| Number of periodic images considered:", dft_control%nimages
     755             : 
     756           4 :             nimg = dft_control%nimages
     757             :          ELSE
     758             :             nimg = 1
     759             :          END IF
     760             : 
     761        1744 :          ALLOCATE (t_3c_overl_int(nimg, nimg))
     762             : 
     763         288 :          DO i = 1, SIZE(t_3c_overl_int, 1)
     764         520 :             DO j = 1, SIZE(t_3c_overl_int, 2)
     765         384 :                CALL dbt_create(t_3c_overl_int_template, t_3c_overl_int(i, j))
     766             :             END DO
     767             :          END DO
     768             : 
     769         136 :          CALL dbt_destroy(t_3c_overl_int_template)
     770             : 
     771             :          ! split blocks to improve load balancing for tensor contraction
     772         136 :          min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize
     773             : 
     774         136 :          CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_AO_split)
     775         136 :          CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_RI_split)
     776             : 
     777         136 :          pdims_t3c = 0
     778         136 :          CALL dbt_pgrid_create(para_env, pdims_t3c, pgrid_t3c_M)
     779             : 
     780             :          ASSOCIATE (cut_memory => qs_env%mp2_env%ri_rpa_im_time%cut_memory)
     781             :             CALL create_tensor_batches(sizes_AO_split, cut_memory, starts_array_mc, ends_array_mc, &
     782         136 :                                        starts_array_mc_block, ends_array_mc_block)
     783             :             CALL create_tensor_batches(sizes_RI_split, cut_memory, &
     784             :                                        qs_env%mp2_env%ri_rpa_im_time%starts_array_mc_RI, &
     785             :                                        qs_env%mp2_env%ri_rpa_im_time%ends_array_mc_RI, &
     786             :                                        qs_env%mp2_env%ri_rpa_im_time%starts_array_mc_block_RI, &
     787         272 :                                        qs_env%mp2_env%ri_rpa_im_time%ends_array_mc_block_RI)
     788             : 
     789             :          END ASSOCIATE
     790         136 :          cut_memory = qs_env%mp2_env%ri_rpa_im_time%cut_memory
     791             : 
     792             :          CALL create_3c_tensor(t_3c_M, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_M, &
     793             :                                sizes_RI_split, sizes_AO_split, sizes_AO_split, &
     794             :                                map1=[1], map2=[2, 3], &
     795         136 :                                name="M (RI | AO AO)")
     796         136 :          DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
     797         136 :          CALL dbt_pgrid_destroy(pgrid_t3c_M)
     798             : 
     799        1744 :          ALLOCATE (t_3c_O(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2)))
     800      288902 :          ALLOCATE (t_3c_O_compressed(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2), cut_memory))
     801        1670 :          ALLOCATE (t_3c_O_ind(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2), cut_memory))
     802             :          CALL create_3c_tensor(t_3c_O(1, 1), dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, &
     803             :                                sizes_RI_split, sizes_AO_split, sizes_AO_split, &
     804             :                                map1=[1, 2], map2=[3], &
     805         136 :                                name="O (RI AO | AO)")
     806         136 :          DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
     807         136 :          CALL dbt_pgrid_destroy(pgrid_t3c_overl)
     808             : 
     809         288 :          DO i = 1, SIZE(t_3c_O, 1)
     810         520 :             DO j = 1, SIZE(t_3c_O, 2)
     811         384 :                IF (i > 1 .OR. j > 1) CALL dbt_create(t_3c_O(1, 1), t_3c_O(i, j))
     812             :             END DO
     813             :          END DO
     814             : 
     815             :          ! build integrals in batches and copy to optimized format
     816             :          ! note: integrals are stored in terms of atomic blocks. To avoid a memory bottleneck,
     817             :          ! integrals are calculated in batches and copied to optimized format with subatomic blocks
     818             : 
     819         378 :          DO cm = 1, cut_memory_int
     820             :             CALL build_3c_integrals(t_3c_overl_int, &
     821             :                                     qs_env%mp2_env%ri_rpa_im_time%eps_filter/2, &
     822             :                                     qs_env, &
     823             :                                     nl_3c, &
     824             :                                     int_eps=qs_env%mp2_env%ri_rpa_im_time%eps_filter/2, &
     825             :                                     basis_i=basis_set_ri_aux, &
     826             :                                     basis_j=basis_set_ao, basis_k=basis_set_ao, &
     827             :                                     potential_parameter=ri_metric, &
     828             :                                     do_kpoints=do_kpoints_cubic_RPA, &
     829         726 :                                     bounds_i=[starts_array_mc_block_int(cm), ends_array_mc_block_int(cm)], desymmetrize=.FALSE.)
     830         242 :             CALL timeset(routineN//"_copy_3c", handle4)
     831             :             ! copy integral tensor t_3c_overl_int to t_3c_O tensor optimized for contraction
     832         508 :             DO i = 1, SIZE(t_3c_overl_int, 1)
     833         894 :                DO j = 1, SIZE(t_3c_overl_int, 2)
     834             : 
     835             :                   CALL dbt_copy(t_3c_overl_int(i, j), t_3c_O(i, j), order=[1, 3, 2], &
     836         386 :                                 summation=.TRUE., move_data=.TRUE.)
     837         386 :                   CALL dbt_clear(t_3c_overl_int(i, j))
     838         386 :                   CALL dbt_filter(t_3c_O(i, j), qs_env%mp2_env%ri_rpa_im_time%eps_filter/2)
     839             :                   ! rescaling, probably because of neighbor list
     840         652 :                   IF (do_kpoints_cubic_RPA .AND. cm == cut_memory_int) THEN
     841         100 :                      CALL dbt_scale(t_3c_O(i, j), 0.5_dp)
     842             :                   END IF
     843             :                END DO
     844             :             END DO
     845         620 :             CALL timestop(handle4)
     846             :          END DO
     847             : 
     848         288 :          DO i = 1, SIZE(t_3c_overl_int, 1)
     849         520 :             DO j = 1, SIZE(t_3c_overl_int, 2)
     850         384 :                CALL dbt_destroy(t_3c_overl_int(i, j))
     851             :             END DO
     852             :          END DO
     853         368 :          DEALLOCATE (t_3c_overl_int)
     854             : 
     855         136 :          CALL timeset(routineN//"_copy_3c", handle4)
     856             :          ! desymmetrize
     857         136 :          CALL dbt_create(t_3c_O(1, 1), t_3c_tmp)
     858         288 :          DO jcell = 1, nimg
     859         480 :             DO kcell = 1, jcell
     860         192 :                CALL dbt_copy(t_3c_O(jcell, kcell), t_3c_tmp)
     861         192 :                CALL dbt_copy(t_3c_tmp, t_3c_O(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.)
     862         344 :                CALL dbt_filter(t_3c_O(kcell, jcell), qs_env%mp2_env%ri_rpa_im_time%eps_filter)
     863             :             END DO
     864             :          END DO
     865         288 :          DO jcell = 1, nimg
     866         328 :             DO kcell = jcell + 1, nimg
     867          40 :                CALL dbt_copy(t_3c_O(jcell, kcell), t_3c_tmp)
     868          40 :                CALL dbt_copy(t_3c_tmp, t_3c_O(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.)
     869         192 :                CALL dbt_filter(t_3c_O(kcell, jcell), qs_env%mp2_env%ri_rpa_im_time%eps_filter)
     870             :             END DO
     871             :          END DO
     872             : 
     873         136 :          CALL dbt_get_info(t_3c_O(1, 1), nfull_total=bounds_3c)
     874         136 :          CALL get_tensor_occupancy(t_3c_O(1, 1), nze, occ)
     875         136 :          memory_3c = 0.0_dp
     876             : 
     877         408 :          bounds(:, 1) = [1, bounds_3c(1)]
     878         408 :          bounds(:, 3) = [1, bounds_3c(3)]
     879         288 :          DO i = 1, SIZE(t_3c_O, 1)
     880         520 :             DO j = 1, SIZE(t_3c_O, 2)
     881         650 :                DO i_mem = 1, cut_memory
     882        1254 :                   bounds(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
     883         418 :                   CALL dbt_copy(t_3c_O(i, j), t_3c_tmp, bounds=bounds)
     884             : 
     885         418 :                   CALL alloc_containers(t_3c_O_compressed(i, j, i_mem), 1)
     886             :                   CALL compress_tensor(t_3c_tmp, t_3c_O_ind(i, j, i_mem)%ind, &
     887             :                                        t_3c_O_compressed(i, j, i_mem), &
     888         650 :                                        qs_env%mp2_env%ri_rpa_im_time%eps_compress, memory_3c)
     889             :                END DO
     890         384 :                CALL dbt_clear(t_3c_O(i, j))
     891             :             END DO
     892             :          END DO
     893             : 
     894         136 :          CALL para_env%sum(memory_3c)
     895             : 
     896         136 :          compression_factor = REAL(nze, dp)*1.0E-06*8.0_dp/memory_3c
     897             : 
     898         136 :          IF (unit_nr > 0) THEN
     899             :             WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") &
     900          68 :                "MEMORY_INFO| Memory for 3-center integrals (compressed):", memory_3c, ' MiB'
     901             : 
     902             :             WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") &
     903          68 :                "MEMORY_INFO| Compression factor:                  ", compression_factor
     904             :          END IF
     905             : 
     906         136 :          CALL dbt_destroy(t_3c_tmp)
     907             : 
     908         136 :          CALL timestop(handle4)
     909             : 
     910         404 :          DO ibasis = 1, SIZE(basis_set_ao)
     911         268 :             orb_basis => basis_set_ao(ibasis)%gto_basis_set
     912         268 :             CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb_old)
     913         268 :             ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
     914         404 :             CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb_old)
     915             :          END DO
     916             : 
     917         136 :          DEALLOCATE (basis_set_ri_aux, basis_set_ao)
     918             : 
     919        1224 :          CALL neighbor_list_3c_destroy(nl_3c)
     920             : 
     921             :       END IF
     922             : 
     923         658 :       CALL timestop(handle)
     924             : 
     925        2632 :    END SUBROUTINE mp2_ri_gpw_compute_in
     926             : 
     927             : ! **************************************************************************************************
     928             : !> \brief Contract (P|ai) = (R|P) x (R|ai)
     929             : !> \param BIb_C (R|ai)
     930             : !> \param my_Lrows (R|P)
     931             : !> \param sizes_B number of a (virtual) indices per subgroup process
     932             : !> \param sizes_L number of P / R (auxiliary) indices per subgroup
     933             : !> \param blk_size ...
     934             : !> \param ngroup how many subgroups (NG)
     935             : !> \param igroup subgroup color
     936             : !> \param mp_comm communicator
     937             : !> \param para_env_sub ...
     938             : ! **************************************************************************************************
     939         346 :    SUBROUTINE contract_B_L(BIb_C, my_Lrows, sizes_B, sizes_L, blk_size, ngroup, igroup, mp_comm, para_env_sub)
     940             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: BIb_C
     941             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: my_Lrows
     942             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: sizes_B, sizes_L
     943             :       INTEGER, DIMENSION(2), INTENT(IN)                  :: blk_size
     944             :       INTEGER, INTENT(IN)                                :: ngroup, igroup
     945             : 
     946             :       CLASS(mp_comm_type), INTENT(IN)                    :: mp_comm
     947             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env_sub
     948             : 
     949             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'contract_B_L'
     950             :       LOGICAL, PARAMETER                                 :: debug = .FALSE.
     951             : 
     952             :       INTEGER                                            :: check_proc, handle, i, iend, ii, ioff, &
     953             :                                                             istart, loc_a, loc_P, nblk_per_thread
     954         346 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: block_ind_L_P, block_ind_L_R
     955             :       INTEGER, DIMENSION(1)                              :: dist_B_i, map_B_1, map_L_1, map_L_2, &
     956             :                                                             sizes_i
     957             :       INTEGER, DIMENSION(2)                              :: map_B_2, pdims_L
     958             :       INTEGER, DIMENSION(3)                              :: pdims_B
     959             :       LOGICAL                                            :: found
     960         692 :       INTEGER, DIMENSION(ngroup)                         :: dist_L_P, dist_L_R
     961         692 :       INTEGER, DIMENSION(para_env_sub%num_pe)            :: dist_B_a
     962        5882 :       TYPE(dbt_distribution_type)                        :: dist_B, dist_L
     963        1730 :       TYPE(dbt_pgrid_type)                               :: mp_comm_B, mp_comm_L
     964        8650 :       TYPE(dbt_type)                                     :: tB_in, tB_in_split, tB_out, &
     965        8650 :                                                             tB_out_split, tL, tL_split
     966             : 
     967         346 :       CALL timeset(routineN, handle)
     968             : 
     969         346 :       sizes_i(1) = SIZE(BIb_C, 3)
     970             : 
     971             :       ASSOCIATE (nproc => para_env_sub%num_pe, iproc => para_env_sub%mepos, iproc_glob => mp_comm%mepos)
     972             : 
     973             :          ! local block index for R/P and a
     974         346 :          loc_P = igroup + 1; loc_a = iproc + 1
     975             : 
     976           0 :          CPASSERT(SIZE(sizes_L) .EQ. ngroup)
     977         346 :          CPASSERT(SIZE(sizes_B) .EQ. nproc)
     978         346 :          CPASSERT(sizes_L(loc_P) .EQ. SIZE(BIb_C, 1))
     979         346 :          CPASSERT(sizes_L(loc_P) .EQ. SIZE(my_Lrows, 2))
     980         346 :          CPASSERT(sizes_B(loc_a) .EQ. SIZE(BIb_C, 2))
     981             : 
     982             :          ! Tensor distributions as follows:
     983             :          ! Process grid NG x Nw
     984             :          ! Each process has coordinates (np, nw)
     985             :          ! tB_in: (R|ai): R distributed (np), a distributed (nw)
     986             :          ! tB_out: (P|ai): P distributed (np), a distributed (nw)
     987             :          ! tL: (R|P): R distributed (nw), P distributed (np)
     988             : 
     989             :          ! define mappings between tensor index and matrix index:
     990             :          ! (R|ai) and (P|ai):
     991         346 :          map_B_1 = [1] ! index 1 (R or P) maps to 1st matrix index (np distributed)
     992         346 :          map_B_2 = [2, 3] ! indices 2, 3 (a, i) map to 2nd matrix index (nw distributed)
     993             :          ! (R|P):
     994         346 :          map_L_1 = [2] ! index 2 (P) maps to 1st matrix index (np distributed)
     995         346 :          map_L_2 = [1] ! index 1 (R) maps to 2nd matrix index (nw distributed)
     996             : 
     997             :          ! derive nd process grid that is compatible with distributions and 2d process grid
     998             :          ! (R|ai) / (P|ai) on process grid NG x Nw x 1
     999             :          ! (R|P) on process grid NG x Nw
    1000        1384 :          pdims_B = [ngroup, nproc, 1]
    1001        1038 :          pdims_L = [nproc, ngroup]
    1002             : 
    1003         346 :          CALL dbt_pgrid_create(mp_comm, pdims_B, mp_comm_B)
    1004         346 :          CALL dbt_pgrid_create(mp_comm, pdims_L, mp_comm_L)
    1005             : 
    1006             :          ! setup distribution vectors such that distribution matches parallel data layout of BIb_C and my_Lrows
    1007         346 :          dist_B_i = [0]
    1008        1396 :          dist_B_a = (/(i, i=0, nproc - 1)/)
    1009        2064 :          dist_L_R = (/(MODULO(i, nproc), i=0, ngroup - 1)/) ! R index is replicated in my_Lrows, we impose a cyclic distribution
    1010        2064 :          dist_L_P = (/(i, i=0, ngroup - 1)/)
    1011             : 
    1012             :          ! create distributions and tensors
    1013         346 :          CALL dbt_distribution_new(dist_B, mp_comm_B, dist_L_P, dist_B_a, dist_B_i)
    1014         346 :          CALL dbt_distribution_new(dist_L, mp_comm_L, dist_L_R, dist_L_P)
    1015             : 
    1016         346 :          CALL dbt_create(tB_in, "(R|ai)", dist_B, map_B_1, map_B_2, sizes_L, sizes_B, sizes_i)
    1017         346 :          CALL dbt_create(tB_out, "(P|ai)", dist_B, map_B_1, map_B_2, sizes_L, sizes_B, sizes_i)
    1018         346 :          CALL dbt_create(tL, "(R|P)", dist_L, map_L_1, map_L_2, sizes_L, sizes_L)
    1019             : 
    1020             :          IF (debug) THEN
    1021             :             ! check that tensor distribution is correct
    1022             :             CALL dbt_get_stored_coordinates(tB_in, [loc_P, loc_a, 1], check_proc)
    1023             :             CPASSERT(check_proc .EQ. iproc_glob)
    1024             :          END IF
    1025             : 
    1026             :          ! reserve (R|ai) block
    1027         346 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tB_in,loc_P,loc_a)
    1028             :          CALL dbt_reserve_blocks(tB_in, [loc_P], [loc_a], [1])
    1029             : !$OMP END PARALLEL
    1030             : 
    1031             :          ! reserve (R|P) blocks
    1032             :          ! in my_Lrows, R index is replicated. For (R|P), we distribute quadratic blocks cyclically over
    1033             :          ! the processes in a subgroup.
    1034             :          ! There are NG blocks, so each process holds at most NG/Nw+1 blocks.
    1035        1038 :          ALLOCATE (block_ind_L_R(ngroup/nproc + 1))
    1036         692 :          ALLOCATE (block_ind_L_P(ngroup/nproc + 1))
    1037        2744 :          block_ind_L_R(:) = 0; block_ind_L_P(:) = 0
    1038             :          ii = 0
    1039        1032 :          DO i = 1, ngroup
    1040        2058 :             CALL dbt_get_stored_coordinates(tL, [i, loc_P], check_proc)
    1041        1032 :             IF (check_proc == iproc_glob) THEN
    1042         683 :                ii = ii + 1
    1043         683 :                block_ind_L_R(ii) = i
    1044         683 :                block_ind_L_P(ii) = loc_P
    1045             :             END IF
    1046             :          END DO
    1047             : 
    1048             : !TODO: Parallelize creation of block list.
    1049             : !$OMP PARALLEL DEFAULT(NONE) SHARED(tL,block_ind_L_R,block_ind_L_P,ii) &
    1050         346 : !$OMP PRIVATE(nblk_per_thread,istart,iend)
    1051             :          nblk_per_thread = ii/omp_get_num_threads() + 1
    1052             :          istart = omp_get_thread_num()*nblk_per_thread + 1
    1053             :          iend = MIN(istart + nblk_per_thread, ii)
    1054             :          CALL dbt_reserve_blocks(tL, block_ind_L_R(istart:iend), block_ind_L_P(istart:iend))
    1055             : !$OMP END PARALLEL
    1056             : 
    1057             :          ! insert (R|ai) block
    1058        2422 :          CALL dbt_put_block(tB_in, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C)
    1059             : 
    1060             :          ! insert (R|P) blocks
    1061         346 :          ioff = 0
    1062        1378 :          DO i = 1, ngroup
    1063         686 :             istart = ioff + 1; iend = ioff + sizes_L(i)
    1064         686 :             ioff = ioff + sizes_L(i)
    1065        2058 :             CALL dbt_get_stored_coordinates(tL, [i, loc_P], check_proc)
    1066        1032 :             IF (check_proc == iproc_glob) THEN
    1067      980422 :                CALL dbt_put_block(tL, [i, loc_P], [sizes_L(i), sizes_L(loc_P)], my_Lrows(istart:iend, :))
    1068             :             END IF
    1069             :          END DO
    1070             :       END ASSOCIATE
    1071             : 
    1072        1384 :       CALL dbt_split_blocks(tB_in, tB_in_split, [blk_size(2), blk_size(1), blk_size(1)])
    1073        1038 :       CALL dbt_split_blocks(tL, tL_split, [blk_size(2), blk_size(2)])
    1074        1384 :       CALL dbt_split_blocks(tB_out, tB_out_split, [blk_size(2), blk_size(1), blk_size(1)])
    1075             : 
    1076             :       ! contract
    1077             :       CALL dbt_contract(alpha=1.0_dp, tensor_1=tB_in_split, tensor_2=tL_split, &
    1078             :                         beta=0.0_dp, tensor_3=tB_out_split, &
    1079             :                         contract_1=[1], notcontract_1=[2, 3], &
    1080             :                         contract_2=[1], notcontract_2=[2], &
    1081         346 :                         map_1=[2, 3], map_2=[1], optimize_dist=.TRUE.)
    1082             : 
    1083             :       ! retrieve local block of contraction result (P|ai)
    1084         346 :       CALL dbt_copy(tB_out_split, tB_out)
    1085             : 
    1086        2422 :       CALL dbt_get_block(tB_out, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C, found)
    1087         346 :       CPASSERT(found)
    1088             : 
    1089             :       ! cleanup
    1090         346 :       CALL dbt_destroy(tB_in)
    1091         346 :       CALL dbt_destroy(tB_in_split)
    1092         346 :       CALL dbt_destroy(tB_out)
    1093         346 :       CALL dbt_destroy(tB_out_split)
    1094         346 :       CALL dbt_destroy(tL)
    1095         346 :       CALL dbt_destroy(tL_split)
    1096             : 
    1097         346 :       CALL dbt_distribution_destroy(dist_B)
    1098         346 :       CALL dbt_distribution_destroy(dist_L)
    1099             : 
    1100         346 :       CALL dbt_pgrid_destroy(mp_comm_B)
    1101         346 :       CALL dbt_pgrid_destroy(mp_comm_L)
    1102             : 
    1103         346 :       CALL timestop(handle)
    1104             : 
    1105         692 :    END SUBROUTINE contract_B_L
    1106             : 
    1107             : ! **************************************************************************************************
    1108             : !> \brief Encapsulate building of intermediate matrices matrix_ia_jnu(_beta
    1109             : !>         matrix_ia_jb(_beta),fm_BIb_jb(_beta),matrix_in_jnu(for G0W0) and
    1110             : !>         fm_BIb_all(for G0W0)
    1111             : !> \param intermed_mat ...
    1112             : !> \param mo_coeff_templ ...
    1113             : !> \param size_1 ...
    1114             : !> \param size_2 ...
    1115             : !> \param matrix_name_2 ...
    1116             : !> \param blacs_env_sub ...
    1117             : !> \param para_env_sub ...
    1118             : !> \author Jan Wilhelm
    1119             : ! **************************************************************************************************
    1120           0 :    SUBROUTINE create_intermediate_matrices(intermed_mat, mo_coeff_templ, size_1, size_2, &
    1121             :                                            matrix_name_2, blacs_env_sub, para_env_sub)
    1122             : 
    1123             :       TYPE(intermediate_matrix_type), INTENT(OUT)        :: intermed_mat
    1124             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: mo_coeff_templ
    1125             :       INTEGER, INTENT(IN)                                :: size_1, size_2
    1126             :       CHARACTER(LEN=*), INTENT(IN)                       :: matrix_name_2
    1127             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_sub
    1128             :       TYPE(mp_para_env_type), POINTER                    :: para_env_sub
    1129             : 
    1130             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_intermediate_matrices'
    1131             : 
    1132             :       INTEGER                                            :: handle, ncol_local, nfullcols_total, &
    1133             :                                                             nfullrows_total, nrow_local
    1134         796 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1135             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    1136             : 
    1137         796 :       CALL timeset(routineN, handle)
    1138             : 
    1139             :       ! initialize and create the matrix (K|jnu)
    1140         796 :       CALL dbcsr_create(intermed_mat%matrix_ia_jnu, template=mo_coeff_templ)
    1141             : 
    1142             :       ! Allocate Sparse matrices: (K|jb)
    1143             :       CALL cp_dbcsr_m_by_n_from_template(intermed_mat%matrix_ia_jb, template=mo_coeff_templ, m=size_2, n=size_1, &
    1144         796 :                                          sym=dbcsr_type_no_symmetry)
    1145             : 
    1146             :       ! set all to zero in such a way that the memory is actually allocated
    1147         796 :       CALL dbcsr_set(intermed_mat%matrix_ia_jnu, 0.0_dp)
    1148         796 :       CALL dbcsr_set(intermed_mat%matrix_ia_jb, 0.0_dp)
    1149             : 
    1150             :       ! create the analogous of matrix_ia_jb in fm type
    1151         796 :       NULLIFY (fm_struct)
    1152         796 :       CALL dbcsr_get_info(intermed_mat%matrix_ia_jb, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
    1153             :       CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, nrow_global=nfullrows_total, &
    1154         796 :                                ncol_global=nfullcols_total, para_env=para_env_sub)
    1155         796 :       CALL cp_fm_create(intermed_mat%fm_BIb_jb, fm_struct, name="fm_BIb_jb_"//matrix_name_2)
    1156             : 
    1157         796 :       CALL copy_dbcsr_to_fm(intermed_mat%matrix_ia_jb, intermed_mat%fm_BIb_jb)
    1158         796 :       CALL cp_fm_struct_release(fm_struct)
    1159             : 
    1160             :       CALL cp_fm_get_info(matrix=intermed_mat%fm_BIb_jb, &
    1161             :                           nrow_local=nrow_local, &
    1162             :                           ncol_local=ncol_local, &
    1163             :                           row_indices=row_indices, &
    1164         796 :                           col_indices=col_indices)
    1165             : 
    1166         796 :       intermed_mat%max_row_col_local = MAX(nrow_local, ncol_local)
    1167         796 :       CALL para_env_sub%max(intermed_mat%max_row_col_local)
    1168             : 
    1169        3184 :       ALLOCATE (intermed_mat%local_col_row_info(0:intermed_mat%max_row_col_local, 2))
    1170       28576 :       intermed_mat%local_col_row_info = 0
    1171             :       ! 0,1 nrows
    1172         796 :       intermed_mat%local_col_row_info(0, 1) = nrow_local
    1173        4637 :       intermed_mat%local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
    1174             :       ! 0,2 ncols
    1175         796 :       intermed_mat%local_col_row_info(0, 2) = ncol_local
    1176       12986 :       intermed_mat%local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
    1177             : 
    1178         796 :       intermed_mat%descr = matrix_name_2
    1179             : 
    1180         796 :       CALL timestop(handle)
    1181             : 
    1182        2388 :    END SUBROUTINE create_intermediate_matrices
    1183             : 
    1184             : ! **************************************************************************************************
    1185             : !> \brief Encapsulate ERI postprocessing: AO to MO transformation and store in B matrix.
    1186             : !> \param para_env ...
    1187             : !> \param mat_munu ...
    1188             : !> \param intermed_mat ...
    1189             : !> \param BIb_jb ...
    1190             : !> \param mo_coeff_o ...
    1191             : !> \param mo_coeff_v ...
    1192             : !> \param eps_filter ...
    1193             : !> \param my_B_end ...
    1194             : !> \param my_B_start ...
    1195             : ! **************************************************************************************************
    1196       31914 :    SUBROUTINE ao_to_mo_and_store_B(para_env, mat_munu, intermed_mat, BIb_jb, &
    1197             :                                    mo_coeff_o, mo_coeff_v, eps_filter, &
    1198             :                                    my_B_end, my_B_start)
    1199             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1200             :       TYPE(dbcsr_p_type), INTENT(IN)                     :: mat_munu
    1201             :       TYPE(intermediate_matrix_type), INTENT(INOUT)      :: intermed_mat
    1202             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: BIb_jb
    1203             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff_o, mo_coeff_v
    1204             :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter
    1205             :       INTEGER, INTENT(IN)                                :: my_B_end, my_B_start
    1206             : 
    1207             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ao_to_mo_and_store_B'
    1208             : 
    1209             :       INTEGER                                            :: handle
    1210             : 
    1211       31914 :       CALL timeset(routineN//"_mult_"//TRIM(intermed_mat%descr), handle)
    1212             : 
    1213             :       CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu%matrix, mo_coeff_o, &
    1214       31914 :                           0.0_dp, intermed_mat%matrix_ia_jnu, filter_eps=eps_filter)
    1215             :       CALL dbcsr_multiply("T", "N", 1.0_dp, intermed_mat%matrix_ia_jnu, mo_coeff_v, &
    1216       31914 :                           0.0_dp, intermed_mat%matrix_ia_jb, filter_eps=eps_filter)
    1217       31914 :       CALL timestop(handle)
    1218             : 
    1219       31914 :       CALL timeset(routineN//"_E_Ex_"//TRIM(intermed_mat%descr), handle)
    1220       31914 :       CALL copy_dbcsr_to_fm(intermed_mat%matrix_ia_jb, intermed_mat%fm_BIb_jb)
    1221             : 
    1222             :       CALL grep_my_integrals(para_env, intermed_mat%fm_BIb_jb, BIb_jb, intermed_mat%max_row_col_local, &
    1223             :                              intermed_mat%local_col_row_info, &
    1224       31914 :                              my_B_end, my_B_start)
    1225             : 
    1226       31914 :       CALL timestop(handle)
    1227       31914 :    END SUBROUTINE ao_to_mo_and_store_B
    1228             : 
    1229             : ! **************************************************************************************************
    1230             : !> \brief ...
    1231             : !> \param intermed_mat ...
    1232             : ! **************************************************************************************************
    1233         796 :    SUBROUTINE release_intermediate_matrices(intermed_mat)
    1234             :       TYPE(intermediate_matrix_type), INTENT(INOUT)      :: intermed_mat
    1235             : 
    1236         796 :       CALL dbcsr_release(intermed_mat%matrix_ia_jnu)
    1237         796 :       CALL dbcsr_release(intermed_mat%matrix_ia_jb)
    1238         796 :       CALL cp_fm_release(intermed_mat%fm_BIb_jb)
    1239         796 :       DEALLOCATE (intermed_mat%local_col_row_info)
    1240             : 
    1241         796 :    END SUBROUTINE
    1242             : 
    1243             : ! **************************************************************************************************
    1244             : !> \brief ...
    1245             : !> \param qs_env ...
    1246             : !> \param kpoints ...
    1247             : !> \param unit_nr ...
    1248             : ! **************************************************************************************************
    1249          18 :    SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr)
    1250             : 
    1251             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1252             :       TYPE(kpoint_type), POINTER                         :: kpoints
    1253             :       INTEGER                                            :: unit_nr
    1254             : 
    1255             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_kpoints'
    1256             : 
    1257             :       INTEGER                                            :: handle, i, i_dim, ix, iy, iz, nkp, &
    1258             :                                                             nkp_extra, nkp_orig
    1259             :       INTEGER, DIMENSION(3)                              :: nkp_grid, nkp_grid_extra, periodic
    1260             :       LOGICAL                                            :: do_extrapolate_kpoints
    1261             :       TYPE(cell_type), POINTER                           :: cell
    1262             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1263             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1264             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1265          18 :          POINTER                                         :: sab_orb
    1266             : 
    1267          18 :       CALL timeset(routineN, handle)
    1268             : 
    1269          18 :       NULLIFY (cell, dft_control, para_env)
    1270          18 :       CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb)
    1271          18 :       CALL get_cell(cell=cell, periodic=periodic)
    1272             : 
    1273             :       ! general because we augment a Monkhorst-Pack mesh by additional points in the BZ
    1274          18 :       kpoints%kp_scheme = "GENERAL"
    1275          18 :       kpoints%symmetry = .FALSE.
    1276          18 :       kpoints%verbose = .FALSE.
    1277          18 :       kpoints%full_grid = .TRUE.
    1278          18 :       kpoints%use_real_wfn = .FALSE.
    1279          18 :       kpoints%eps_geo = 1.e-6_dp
    1280          72 :       nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid(1:3)
    1281          18 :       do_extrapolate_kpoints = qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints
    1282             : 
    1283          72 :       DO i_dim = 1, 3
    1284          54 :          IF (periodic(i_dim) == 1) THEN
    1285          36 :             CPASSERT(MODULO(nkp_grid(i_dim), 2) == 0)
    1286             :          END IF
    1287          72 :          IF (periodic(i_dim) == 0) THEN
    1288          18 :             CPASSERT(nkp_grid(i_dim) == 1)
    1289             :          END IF
    1290             :       END DO
    1291             : 
    1292          18 :       nkp_orig = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2
    1293             : 
    1294          18 :       IF (do_extrapolate_kpoints) THEN
    1295             : 
    1296          18 :          CPASSERT(qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method == kp_weights_W_uniform)
    1297             : 
    1298          72 :          DO i_dim = 1, 3
    1299          54 :             IF (periodic(i_dim) == 1) nkp_grid_extra(i_dim) = nkp_grid(i_dim) + 2
    1300          72 :             IF (periodic(i_dim) == 0) nkp_grid_extra(i_dim) = 1
    1301             :          END DO
    1302             : 
    1303          72 :          qs_env%mp2_env%ri_rpa_im_time%kp_grid_extra(1:3) = nkp_grid_extra(1:3)
    1304             : 
    1305          18 :          nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
    1306             : 
    1307             :       ELSE
    1308             : 
    1309           0 :          nkp_grid_extra(1:3) = 0
    1310           0 :          nkp_extra = 0
    1311             : 
    1312             :       END IF
    1313             : 
    1314          18 :       nkp = nkp_orig + nkp_extra
    1315             : 
    1316          18 :       qs_env%mp2_env%ri_rpa_im_time%nkp_orig = nkp_orig
    1317          18 :       qs_env%mp2_env%ri_rpa_im_time%nkp_extra = nkp_extra
    1318             : 
    1319          90 :       ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
    1320             : 
    1321          72 :       kpoints%nkp_grid(1:3) = nkp_grid(1:3)
    1322          18 :       kpoints%nkp = nkp
    1323             : 
    1324          36 :       ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%wkp_V(nkp))
    1325          18 :       IF (do_extrapolate_kpoints) THEN
    1326             :          kpoints%wkp(1:nkp_orig) = 1.0_dp/REAL(nkp_orig, KIND=dp) &
    1327         162 :                                    /(1.0_dp - SQRT(REAL(nkp_extra, KIND=dp)/REAL(nkp_orig, KIND=dp)))
    1328             :          kpoints%wkp(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp) &
    1329         342 :                                          /(1.0_dp - SQRT(REAL(nkp_orig, KIND=dp)/REAL(nkp_extra, KIND=dp)))
    1330         162 :          qs_env%mp2_env%ri_rpa_im_time%wkp_V(1:nkp_orig) = 0.0_dp
    1331         342 :          qs_env%mp2_env%ri_rpa_im_time%wkp_V(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp)
    1332             :       ELSE
    1333           0 :          kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp)
    1334           0 :          qs_env%mp2_env%ri_rpa_im_time%wkp_V(:) = kpoints%wkp(:)
    1335             :       END IF
    1336             : 
    1337          18 :       i = 0
    1338          48 :       DO ix = 1, nkp_grid(1)
    1339         120 :          DO iy = 1, nkp_grid(2)
    1340         390 :             DO iz = 1, nkp_grid(3)
    1341             : 
    1342         288 :                IF (i == nkp_orig) CYCLE
    1343         144 :                i = i + 1
    1344             : 
    1345         144 :                kpoints%xkp(1, i) = REAL(2*ix - nkp_grid(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(1), KIND=dp))
    1346         144 :                kpoints%xkp(2, i) = REAL(2*iy - nkp_grid(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(2), KIND=dp))
    1347         360 :                kpoints%xkp(3, i) = REAL(2*iz - nkp_grid(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(3), KIND=dp))
    1348             : 
    1349             :             END DO
    1350             :          END DO
    1351             :       END DO
    1352             : 
    1353          56 :       DO ix = 1, nkp_grid_extra(1)
    1354         164 :          DO iy = 1, nkp_grid_extra(2)
    1355         794 :             DO iz = 1, nkp_grid_extra(3)
    1356             : 
    1357         648 :                i = i + 1
    1358         648 :                IF (i > nkp) CYCLE
    1359             : 
    1360         324 :                kpoints%xkp(1, i) = REAL(2*ix - nkp_grid_extra(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid_extra(1), KIND=dp))
    1361         324 :                kpoints%xkp(2, i) = REAL(2*iy - nkp_grid_extra(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid_extra(2), KIND=dp))
    1362         756 :                kpoints%xkp(3, i) = REAL(2*iz - nkp_grid_extra(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid_extra(3), KIND=dp))
    1363             : 
    1364             :             END DO
    1365             :          END DO
    1366             :       END DO
    1367             : 
    1368          18 :       CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control)
    1369             : 
    1370          18 :       CALL set_qs_env(qs_env, kpoints=kpoints)
    1371             : 
    1372          18 :       IF (unit_nr > 0) THEN
    1373             : 
    1374           9 :          IF (do_extrapolate_kpoints) THEN
    1375           9 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh for V (leading to Sigma^x):", nkp_grid(1:3)
    1376           9 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T69)") "KPOINT_INFO| K-point extrapolation for W^c is used (W^c leads to Sigma^c):"
    1377           9 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh 1 for W^c:", nkp_grid(1:3)
    1378           9 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh 2 for W^c:", nkp_grid_extra(1:3)
    1379             :          ELSE
    1380           0 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh for V and W:", nkp_grid(1:3)
    1381           0 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T75,I6)") "KPOINT_INFO| Number of kpoints for V and W:", nkp
    1382             :          END IF
    1383             : 
    1384           9 :          SELECT CASE (qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method)
    1385             :          CASE (kp_weights_W_tailored)
    1386             :             WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") &
    1387           0 :                "KPOINT_INFO| K-point weights for W:                                   TAILORED"
    1388             :          CASE (kp_weights_W_auto)
    1389             :             WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") &
    1390           0 :                "KPOINT_INFO| K-point weights for W:                                       AUTO"
    1391             :          CASE (kp_weights_W_uniform)
    1392             :             WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") &
    1393           9 :                "KPOINT_INFO| K-point weights for W:                                    UNIFORM"
    1394             :          END SELECT
    1395             : 
    1396             :       END IF
    1397             : 
    1398          18 :       CALL timestop(handle)
    1399             : 
    1400          18 :    END SUBROUTINE compute_kpoints
    1401             : 
    1402             : ! **************************************************************************************************
    1403             : !> \brief ...
    1404             : !> \param para_env_sub ...
    1405             : !> \param fm_BIb_jb ...
    1406             : !> \param BIb_jb ...
    1407             : !> \param max_row_col_local ...
    1408             : !> \param local_col_row_info ...
    1409             : !> \param my_B_virtual_end ...
    1410             : !> \param my_B_virtual_start ...
    1411             : ! **************************************************************************************************
    1412       31914 :    SUBROUTINE grep_my_integrals(para_env_sub, fm_BIb_jb, BIb_jb, max_row_col_local, &
    1413             :                                 local_col_row_info, &
    1414             :                                 my_B_virtual_end, my_B_virtual_start)
    1415             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env_sub
    1416             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_BIb_jb
    1417             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: BIb_jb
    1418             :       INTEGER, INTENT(IN)                                :: max_row_col_local
    1419             :       INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN)  :: local_col_row_info
    1420             :       INTEGER, INTENT(IN)                                :: my_B_virtual_end, my_B_virtual_start
    1421             : 
    1422             :       INTEGER                                            :: i_global, iiB, j_global, jjB, ncol_rec, &
    1423             :                                                             nrow_rec, proc_receive, proc_send, &
    1424             :                                                             proc_shift
    1425             :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: rec_col_row_info
    1426       31914 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices_rec, row_indices_rec
    1427       31914 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: local_BI, rec_BI
    1428             : 
    1429      127656 :       ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
    1430             : 
    1431     1201746 :       rec_col_row_info(:, :) = local_col_row_info
    1432             : 
    1433       31914 :       nrow_rec = rec_col_row_info(0, 1)
    1434       31914 :       ncol_rec = rec_col_row_info(0, 2)
    1435             : 
    1436       95700 :       ALLOCATE (row_indices_rec(nrow_rec))
    1437      192995 :       row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
    1438             : 
    1439       95742 :       ALLOCATE (col_indices_rec(ncol_rec))
    1440      548408 :       col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
    1441             : 
    1442             :       ! accumulate data on BIb_jb buffer starting from myself
    1443      548408 :       DO jjB = 1, ncol_rec
    1444      516494 :          j_global = col_indices_rec(jjB)
    1445      548408 :          IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN
    1446     3236970 :             DO iiB = 1, nrow_rec
    1447     2745076 :                i_global = row_indices_rec(iiB)
    1448     3236970 :                BIb_jb(j_global - my_B_virtual_start + 1, i_global) = fm_BIb_jb%local_data(iiB, jjB)
    1449             :             END DO
    1450             :          END IF
    1451             :       END DO
    1452             : 
    1453       31914 :       DEALLOCATE (row_indices_rec)
    1454       31914 :       DEALLOCATE (col_indices_rec)
    1455             : 
    1456       31914 :       IF (para_env_sub%num_pe > 1) THEN
    1457        9816 :          ALLOCATE (local_BI(nrow_rec, ncol_rec))
    1458      156209 :          local_BI(1:nrow_rec, 1:ncol_rec) = fm_BIb_jb%local_data(1:nrow_rec, 1:ncol_rec)
    1459             : 
    1460        4908 :          DO proc_shift = 1, para_env_sub%num_pe - 1
    1461        2454 :             proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
    1462        2454 :             proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
    1463             : 
    1464             :             ! first exchange information on the local data
    1465      110670 :             rec_col_row_info = 0
    1466        2454 :             CALL para_env_sub%sendrecv(local_col_row_info, proc_send, rec_col_row_info, proc_receive)
    1467        2454 :             nrow_rec = rec_col_row_info(0, 1)
    1468        2454 :             ncol_rec = rec_col_row_info(0, 2)
    1469             : 
    1470        7362 :             ALLOCATE (row_indices_rec(nrow_rec))
    1471        7705 :             row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
    1472             : 
    1473        7362 :             ALLOCATE (col_indices_rec(ncol_rec))
    1474       51654 :             col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
    1475             : 
    1476        9816 :             ALLOCATE (rec_BI(nrow_rec, ncol_rec))
    1477      156209 :             rec_BI = 0.0_dp
    1478             : 
    1479             :             ! then send and receive the real data
    1480      309964 :             CALL para_env_sub%sendrecv(local_BI, proc_send, rec_BI, proc_receive)
    1481             : 
    1482             :             ! accumulate the received data on BIb_jb buffer
    1483       51654 :             DO jjB = 1, ncol_rec
    1484       49200 :                j_global = col_indices_rec(jjB)
    1485       51654 :                IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN
    1486       76719 :                   DO iiB = 1, nrow_rec
    1487       52119 :                      i_global = row_indices_rec(iiB)
    1488       76719 :                      BIb_jb(j_global - my_B_virtual_start + 1, i_global) = rec_BI(iiB, jjB)
    1489             :                   END DO
    1490             :                END IF
    1491             :             END DO
    1492             : 
    1493        2454 :             DEALLOCATE (col_indices_rec)
    1494        2454 :             DEALLOCATE (row_indices_rec)
    1495        4908 :             DEALLOCATE (rec_BI)
    1496             :          END DO
    1497             : 
    1498        2454 :          DEALLOCATE (local_BI)
    1499             :       END IF
    1500             : 
    1501       31914 :       DEALLOCATE (rec_col_row_info)
    1502             : 
    1503       31914 :    END SUBROUTINE grep_my_integrals
    1504             : 
    1505           0 : END MODULE mp2_integrals

Generated by: LCOV version 1.15