LCOV - code coverage report
Current view: top level - src - mp2_ri_2c.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 611 642 95.2 %
Date: 2025-02-18 08:24:35 Functions: 16 16 100.0 %

          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 Framework for 2c-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_ri_2c
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16             :                                               get_atomic_kind_set
      17             :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      18             :                                               gto_basis_set_type
      19             :    USE cell_types,                      ONLY: cell_type,&
      20             :                                               get_cell,&
      21             :                                               plane_distance
      22             :    USE constants_operator,              ONLY: operator_coulomb
      23             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
      24             :                                               cp_blacs_env_release,&
      25             :                                               cp_blacs_env_type
      26             :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_scale_and_add_fm
      27             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      28             :                                               cp_cfm_release,&
      29             :                                               cp_cfm_to_fm,&
      30             :                                               cp_cfm_type
      31             :    USE cp_control_types,                ONLY: dft_control_type
      32             :    USE cp_dbcsr_api,                    ONLY: &
      33             :         dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
      34             :         dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, dbcsr_release, &
      35             :         dbcsr_reserve_all_blocks, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
      36             :         dbcsr_type_no_symmetry, dbcsr_type_symmetric
      37             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      38             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      39             :                                               cp_dbcsr_dist2d_to_dist,&
      40             :                                               dbcsr_allocate_matrix_set,&
      41             :                                               dbcsr_deallocate_matrix_set
      42             :    USE cp_eri_mme_interface,            ONLY: cp_eri_mme_param
      43             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      44             :                                               cp_fm_triangular_invert
      45             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose
      46             :    USE cp_fm_diag,                      ONLY: choose_eigv_solver,&
      47             :                                               cp_fm_power,&
      48             :                                               cp_fm_syevx
      49             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      50             :                                               cp_fm_struct_release,&
      51             :                                               cp_fm_struct_type
      52             :    USE cp_fm_types,                     ONLY: cp_fm_copy_general,&
      53             :                                               cp_fm_create,&
      54             :                                               cp_fm_get_info,&
      55             :                                               cp_fm_release,&
      56             :                                               cp_fm_set_all,&
      57             :                                               cp_fm_to_fm,&
      58             :                                               cp_fm_type
      59             :    USE distribution_2d_types,           ONLY: distribution_2d_type
      60             :    USE group_dist_types,                ONLY: create_group_dist,&
      61             :                                               get_group_dist,&
      62             :                                               group_dist_d1_type,&
      63             :                                               release_group_dist
      64             :    USE input_constants,                 ONLY: do_eri_gpw,&
      65             :                                               do_eri_mme,&
      66             :                                               do_eri_os,&
      67             :                                               do_potential_coulomb,&
      68             :                                               do_potential_id,&
      69             :                                               do_potential_long,&
      70             :                                               do_potential_truncated
      71             :    USE kinds,                           ONLY: dp
      72             :    USE kpoint_coulomb_2c,               ONLY: build_2c_coulomb_matrix_kp
      73             :    USE kpoint_methods,                  ONLY: kpoint_init_cell_index,&
      74             :                                               rskp_transform
      75             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      76             :                                               kpoint_type
      77             :    USE libint_2c_3c,                    ONLY: compare_potential_types,&
      78             :                                               libint_potential_type
      79             :    USE machine,                         ONLY: m_flush
      80             :    USE message_passing,                 ONLY: mp_comm_type,&
      81             :                                               mp_para_env_release,&
      82             :                                               mp_para_env_type,&
      83             :                                               mp_request_type
      84             :    USE mp2_eri,                         ONLY: mp2_eri_2c_integrate
      85             :    USE mp2_eri_gpw,                     ONLY: mp2_eri_2c_integrate_gpw
      86             :    USE mp2_types,                       ONLY: integ_mat_buffer_type
      87             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      88             :    USE particle_methods,                ONLY: get_particle_set
      89             :    USE particle_types,                  ONLY: particle_type
      90             :    USE qs_environment_types,            ONLY: get_qs_env,&
      91             :                                               qs_environment_type
      92             :    USE qs_integral_utils,               ONLY: basis_set_list_setup
      93             :    USE qs_interactions,                 ONLY: init_interaction_radii
      94             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      95             :                                               qs_kind_type
      96             :    USE qs_ks_types,                     ONLY: get_ks_env,&
      97             :                                               set_ks_env
      98             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type,&
      99             :                                               release_neighbor_list_sets
     100             :    USE qs_tensors,                      ONLY: build_2c_integrals,&
     101             :                                               build_2c_neighbor_lists
     102             :    USE rpa_communication,               ONLY: communicate_buffer
     103             :    USE rpa_gw_kpoints_util,             ONLY: cp_cfm_power
     104             : 
     105             : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
     106             : #include "./base/base_uses.f90"
     107             : 
     108             :    IMPLICIT NONE
     109             : 
     110             :    PRIVATE
     111             : 
     112             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_2c'
     113             : 
     114             :    PUBLIC :: get_2c_integrals, trunc_coulomb_for_exchange, RI_2c_integral_mat, &
     115             :              inversion_of_M_and_mult_with_chol_dec_of_V
     116             : 
     117             : CONTAINS
     118             : 
     119             : ! **************************************************************************************************
     120             : !> \brief ...
     121             : !> \param qs_env ...
     122             : !> \param eri_method ...
     123             : !> \param eri_param ...
     124             : !> \param para_env ...
     125             : !> \param para_env_sub ...
     126             : !> \param mp2_memory ...
     127             : !> \param my_Lrows ...
     128             : !> \param my_Vrows ...
     129             : !> \param fm_matrix_PQ ...
     130             : !> \param ngroup ...
     131             : !> \param color_sub ...
     132             : !> \param dimen_RI ...
     133             : !> \param dimen_RI_red reduced RI dimension,  needed if we perform SVD instead of Cholesky
     134             : !> \param kpoints ...
     135             : !> \param my_group_L_size ...
     136             : !> \param my_group_L_start ...
     137             : !> \param my_group_L_end ...
     138             : !> \param gd_array ...
     139             : !> \param calc_PQ_cond_num ...
     140             : !> \param do_svd ...
     141             : !> \param eps_svd ...
     142             : !> \param potential ...
     143             : !> \param ri_metric ...
     144             : !> \param fm_matrix_L_kpoints ...
     145             : !> \param fm_matrix_Minv_L_kpoints ...
     146             : !> \param fm_matrix_Minv ...
     147             : !> \param fm_matrix_Minv_Vtrunc_Minv ...
     148             : !> \param do_im_time ...
     149             : !> \param do_kpoints ...
     150             : !> \param mp2_eps_pgf_orb_S ...
     151             : !> \param qs_kind_set ...
     152             : !> \param sab_orb_sub ...
     153             : !> \param calc_forces ...
     154             : !> \param unit_nr ...
     155             : ! **************************************************************************************************
     156         658 :    SUBROUTINE get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, &
     157             :                                my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, &
     158             :                                kpoints, my_group_L_size, my_group_L_start, my_group_L_end, &
     159             :                                gd_array, calc_PQ_cond_num, do_svd, eps_svd, potential, ri_metric, &
     160             :                                fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
     161             :                                fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
     162             :                                do_im_time, do_kpoints, mp2_eps_pgf_orb_S, qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
     163             : 
     164             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     165             :       INTEGER, INTENT(IN)                                :: eri_method
     166             :       TYPE(cp_eri_mme_param), POINTER                    :: eri_param
     167             :       TYPE(mp_para_env_type), POINTER                    :: para_env, para_env_sub
     168             :       REAL(KIND=dp), INTENT(IN)                          :: mp2_memory
     169             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     170             :          INTENT(OUT)                                     :: my_Lrows, my_Vrows
     171             :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_matrix_PQ
     172             :       INTEGER, INTENT(IN)                                :: ngroup, color_sub
     173             :       INTEGER, INTENT(OUT)                               :: dimen_RI, dimen_RI_red
     174             :       TYPE(kpoint_type), POINTER                         :: kpoints
     175             :       INTEGER, INTENT(OUT)                               :: my_group_L_size, my_group_L_start, &
     176             :                                                             my_group_L_end
     177             :       TYPE(group_dist_d1_type), INTENT(OUT)              :: gd_array
     178             :       LOGICAL, INTENT(IN)                                :: calc_PQ_cond_num, do_svd
     179             :       REAL(KIND=dp), INTENT(IN)                          :: eps_svd
     180             :       TYPE(libint_potential_type)                        :: potential, ri_metric
     181             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_matrix_L_kpoints, &
     182             :                                                             fm_matrix_Minv_L_kpoints, &
     183             :                                                             fm_matrix_Minv, &
     184             :                                                             fm_matrix_Minv_Vtrunc_Minv
     185             :       LOGICAL, INTENT(IN)                                :: do_im_time, do_kpoints
     186             :       REAL(KIND=dp), INTENT(IN)                          :: mp2_eps_pgf_orb_S
     187             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     188             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     189             :          POINTER                                         :: sab_orb_sub
     190             :       LOGICAL, INTENT(IN)                                :: calc_forces
     191             :       INTEGER, INTENT(IN)                                :: unit_nr
     192             : 
     193             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_2c_integrals'
     194             : 
     195             :       INTEGER                                            :: handle, num_small_eigen
     196             :       REAL(KIND=dp)                                      :: cond_num, eps_pgf_orb_old
     197             :       TYPE(cp_fm_type)                                   :: fm_matrix_L_work, fm_matrix_M_inv_work, &
     198             :                                                             fm_matrix_V
     199             :       TYPE(dft_control_type), POINTER                    :: dft_control
     200             :       TYPE(libint_potential_type)                        :: trunc_coulomb
     201             :       TYPE(mp_para_env_type), POINTER                    :: para_env_L
     202             : 
     203         658 :       CALL timeset(routineN, handle)
     204             : 
     205             :       ! calculate V and store it in fm_matrix_L_work
     206             :       CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
     207             :                                 fm_matrix_L_work, ngroup, color_sub, dimen_RI, &
     208             :                                 my_group_L_size, my_group_L_start, my_group_L_end, &
     209             :                                 gd_array, calc_PQ_cond_num, cond_num, &
     210         658 :                                 num_small_eigen, potential, sab_orb_sub, do_im_time=do_im_time)
     211             : 
     212         658 :       IF (do_im_time .AND. calc_forces) THEN
     213             :          !need a copy of the (P|Q) integrals
     214          50 :          CALL cp_fm_create(fm_matrix_PQ, fm_matrix_L_work%matrix_struct)
     215          50 :          CALL cp_fm_to_fm(fm_matrix_L_work, fm_matrix_PQ)
     216             :       END IF
     217             : 
     218         658 :       dimen_RI_red = dimen_RI
     219             : 
     220         658 :       IF (compare_potential_types(ri_metric, potential) .AND. .NOT. do_im_time) THEN
     221             :          CALL decomp_mat_L(fm_matrix_L_work, do_svd, eps_svd, num_small_eigen, cond_num, .TRUE., gd_array, ngroup, &
     222         486 :                            dimen_RI, dimen_RI_red, para_env)
     223             :       ELSE
     224             : 
     225             :          ! RI-metric matrix (depending on RI metric: Coulomb, overlap or truncated Coulomb)
     226         172 :          IF (do_im_time) THEN
     227         136 :             CALL get_qs_env(qs_env, dft_control=dft_control)
     228             : 
     229             :             ! re-init the radii to be able to generate pair lists with appropriate screening for overlap matrix
     230         136 :             eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
     231         136 :             dft_control%qs_control%eps_pgf_orb = mp2_eps_pgf_orb_S
     232         136 :             CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
     233             : 
     234             :             CALL RI_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L_work, dimen_RI, ri_metric, &
     235             :                                     do_kpoints, kpoints, put_mat_KS_env=.TRUE., &
     236         136 :                                     regularization_RI=qs_env%mp2_env%ri_rpa_im_time%regularization_RI)
     237             : 
     238             :             ! re-init the radii to previous values
     239         136 :             dft_control%qs_control%eps_pgf_orb = eps_pgf_orb_old
     240         136 :             CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
     241             : 
     242             :             ! GPW overlap matrix
     243             :          ELSE
     244             : 
     245          36 :             CALL mp_para_env_release(para_env_L)
     246          36 :             CALL release_group_dist(gd_array)
     247             : 
     248         108 :             ALLOCATE (fm_matrix_Minv_L_kpoints(1, 1))
     249             : 
     250             :             ! Calculate matrix of RI operator (for overlap metric: S), store it in fm_matrix_Minv_L_kpoints
     251             :             CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
     252             :                                       fm_matrix_Minv_L_kpoints(1, 1), ngroup, color_sub, dimen_RI, &
     253             :                                       my_group_L_size, my_group_L_start, my_group_L_end, &
     254             :                                       gd_array, calc_PQ_cond_num, cond_num, &
     255             :                                       num_small_eigen, ri_metric, sab_orb_sub, &
     256          36 :                                       fm_matrix_L_extern=fm_matrix_L_work)
     257             : 
     258             :          END IF
     259             : 
     260         172 :          IF (do_kpoints) THEN
     261             : 
     262             :             ! k-dependent 1/r Coulomb matrix V_PQ(k)
     263          22 :             CALL compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, kpoints)
     264             : 
     265             :             CALL inversion_of_M_and_mult_with_chol_dec_of_V(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, dimen_RI, &
     266          22 :                                                             kpoints, qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S)
     267             : 
     268          22 :             CALL trunc_coulomb_for_exchange(qs_env, trunc_coulomb)
     269             : 
     270             :             ! Gamma-only truncated Coulomb matrix V^tr with cutoff radius = half the unit cell size; for exchange self-energy
     271             :             CALL RI_2c_integral_mat(qs_env, fm_matrix_Minv_Vtrunc_Minv, fm_matrix_L_work, dimen_RI, trunc_coulomb, &
     272          22 :                                     do_kpoints=.FALSE., kpoints=kpoints, put_mat_KS_env=.FALSE., regularization_RI=0.0_dp)
     273             : 
     274             :             ! Gamma-only RI-metric matrix; for computing Gamma-only/MIC self-energy
     275             :             CALL RI_2c_integral_mat(qs_env, fm_matrix_Minv, fm_matrix_L_work, dimen_RI, ri_metric, &
     276          22 :                                     do_kpoints=.FALSE., kpoints=kpoints, put_mat_KS_env=.FALSE., regularization_RI=0.0_dp)
     277             : 
     278             :             CALL Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
     279          22 :                                                                             fm_matrix_Minv, qs_env)
     280             :          ELSE
     281         150 :             IF (calc_forces .AND. (.NOT. do_im_time)) THEN
     282             :                ! For the gradients, we make a copy of the inverse root of L
     283          12 :                CALL cp_fm_create(fm_matrix_V, fm_matrix_L_work%matrix_struct)
     284          12 :                CALL cp_fm_to_fm(fm_matrix_L_work, fm_matrix_V)
     285             : 
     286             :                CALL decomp_mat_L(fm_matrix_V, do_svd, eps_svd, num_small_eigen, cond_num, .TRUE., gd_array, ngroup, &
     287          12 :                                  dimen_RI, dimen_RI_red, para_env)
     288             :             END IF
     289             : 
     290             :             CALL decomp_mat_L(fm_matrix_L_work, do_svd, eps_svd, num_small_eigen, cond_num, .FALSE., gd_array, ngroup, &
     291         150 :                               dimen_RI, dimen_RI_red, para_env)
     292             : 
     293             :             CALL decomp_mat_L(fm_matrix_Minv_L_kpoints(1, 1), .FALSE., 0.0_dp, num_small_eigen, cond_num, .TRUE., &
     294         150 :                               gd_array, ngroup, dimen_RI, dimen_RI_red, para_env)
     295             : 
     296         150 :             CALL cp_fm_create(fm_matrix_M_inv_work, fm_matrix_Minv_L_kpoints(1, 1)%matrix_struct)
     297         150 :             CALL cp_fm_set_all(fm_matrix_M_inv_work, 0.0_dp)
     298             : 
     299             :             CALL parallel_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_Minv_L_kpoints(1, 1), &
     300         150 :                                fm_matrix_Minv_L_kpoints(1, 1), 0.0_dp, fm_matrix_M_inv_work)
     301             : 
     302         150 :             IF (do_svd) THEN
     303             :                ! We have to reset the size of fm_matrix_Minv_L_kpoints
     304          32 :                CALL reset_size_matrix(fm_matrix_Minv_L_kpoints(1, 1), dimen_RI_red, fm_matrix_L_work%matrix_struct)
     305             : 
     306             :                ! L (=fm_matrix_Minv_L_kpoints) = S^(-1)*chol_dec(V)
     307             :                CALL parallel_gemm('T', 'N', dimen_RI, dimen_RI_red, dimen_RI, 1.0_dp, fm_matrix_M_inv_work, &
     308          32 :                                   fm_matrix_L_work, 0.0_dp, fm_matrix_Minv_L_kpoints(1, 1))
     309             :             ELSE
     310             :                ! L (=fm_matrix_Minv_L_kpoints) = S^(-1)*chol_dec(V)
     311             :                CALL parallel_gemm('T', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_M_inv_work, &
     312         118 :                                   fm_matrix_L_work, 0.0_dp, fm_matrix_Minv_L_kpoints(1, 1))
     313             :             END IF
     314             : 
     315             :             ! clean the S_inv matrix
     316         150 :             CALL cp_fm_release(fm_matrix_M_inv_work)
     317             :          END IF
     318             : 
     319         172 :          IF (.NOT. do_im_time) THEN
     320             : 
     321          36 :             CALL cp_fm_to_fm(fm_matrix_Minv_L_kpoints(1, 1), fm_matrix_L_work)
     322          36 :             CALL cp_fm_release(fm_matrix_Minv_L_kpoints)
     323             : 
     324             :          END IF
     325             : 
     326             :       END IF
     327             : 
     328             :       ! If the group distribution changed because of an SVD, we get the new values here
     329         658 :       CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
     330             : 
     331        1180 :       IF (.NOT. do_im_time) THEN
     332         522 :          IF (unit_nr > 0) THEN
     333         261 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Cholesky decomposition group size:", para_env_L%num_pe
     334         261 :             WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Number of groups for auxiliary basis functions", ngroup
     335         261 :             IF (calc_PQ_cond_num .OR. do_svd) THEN
     336             :                WRITE (UNIT=unit_nr, FMT="(T3,A,T67,ES14.5)") &
     337          16 :                   "RI_INFO| Condition number of the (P|Q):", cond_num
     338             :                WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
     339          16 :                   "RI_INFO| Number of non-positive Eigenvalues of (P|Q):", num_small_eigen
     340             :             END IF
     341         261 :             CALL m_flush(unit_nr)
     342             :          END IF
     343             : 
     344             :          ! replicate the necessary row of the L^{-1} matrix on each proc
     345         522 :          CALL grep_Lcols(fm_matrix_L_work, my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
     346         534 :          IF (calc_forces .AND. .NOT. compare_potential_types(qs_env%mp2_env%ri_metric, &
     347             :                                                              qs_env%mp2_env%potential_parameter)) THEN
     348          12 :             CALL grep_Lcols(fm_matrix_V, my_group_L_start, my_group_L_end, my_group_L_size, my_Vrows)
     349             :          END IF
     350             :       END IF
     351         658 :       CALL cp_fm_release(fm_matrix_L_work)
     352         658 :       IF (calc_forces .AND. .NOT. (do_im_time .OR. compare_potential_types(qs_env%mp2_env%ri_metric, &
     353          12 :                                                                qs_env%mp2_env%potential_parameter))) CALL cp_fm_release(fm_matrix_V)
     354         658 :       CALL mp_para_env_release(para_env_L)
     355             : 
     356         658 :       CALL timestop(handle)
     357             : 
     358         658 :    END SUBROUTINE get_2c_integrals
     359             : 
     360             : ! **************************************************************************************************
     361             : !> \brief ...
     362             : !> \param fm_matrix_L ...
     363             : !> \param do_svd ...
     364             : !> \param eps_svd ...
     365             : !> \param num_small_eigen ...
     366             : !> \param cond_num ...
     367             : !> \param do_inversion ...
     368             : !> \param gd_array ...
     369             : !> \param ngroup ...
     370             : !> \param dimen_RI ...
     371             : !> \param dimen_RI_red ...
     372             : !> \param para_env ...
     373             : ! **************************************************************************************************
     374         798 :    SUBROUTINE decomp_mat_L(fm_matrix_L, do_svd, eps_svd, num_small_eigen, cond_num, do_inversion, gd_array, ngroup, &
     375             :                            dimen_RI, dimen_RI_red, para_env)
     376             : 
     377             :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_matrix_L
     378             :       LOGICAL, INTENT(IN)                                :: do_svd
     379             :       REAL(KIND=dp), INTENT(IN)                          :: eps_svd
     380             :       INTEGER, INTENT(INOUT)                             :: num_small_eigen
     381             :       REAL(KIND=dp), INTENT(INOUT)                       :: cond_num
     382             :       LOGICAL, INTENT(IN)                                :: do_inversion
     383             :       TYPE(group_dist_d1_type), INTENT(INOUT)            :: gd_array
     384             :       INTEGER, INTENT(IN)                                :: ngroup, dimen_RI
     385             :       INTEGER, INTENT(INOUT)                             :: dimen_RI_red
     386             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     387             : 
     388         798 :       IF (do_svd) THEN
     389          44 :          CALL matrix_root_with_svd(fm_matrix_L, num_small_eigen, cond_num, eps_svd, do_inversion, para_env)
     390             : 
     391          44 :          dimen_RI_red = dimen_RI - num_small_eigen
     392             : 
     393             :          ! We changed the size of fm_matrix_L in matrix_root_with_svd.
     394             :          ! So, we have to get new group sizes
     395          44 :          CALL release_group_dist(gd_array)
     396             : 
     397          44 :          CALL create_group_dist(gd_array, ngroup, dimen_RI_red)
     398             :       ELSE
     399             : 
     400             :          ! calculate Cholesky decomposition L (V=LL^T)
     401         754 :          CALL cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion=do_inversion)
     402             : 
     403         754 :          IF (do_inversion) CALL invert_mat(fm_matrix_L)
     404             :       END IF
     405             : 
     406         798 :    END SUBROUTINE decomp_mat_L
     407             : 
     408             : ! **************************************************************************************************
     409             : !> \brief ...
     410             : !> \param qs_env ...
     411             : !> \param fm_matrix_L_kpoints ...
     412             : !> \param fm_matrix_Minv_L_kpoints ...
     413             : !> \param kpoints ...
     414             : ! **************************************************************************************************
     415          22 :    SUBROUTINE compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, kpoints)
     416             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     417             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_matrix_L_kpoints, &
     418             :                                                             fm_matrix_Minv_L_kpoints
     419             :       TYPE(kpoint_type), POINTER                         :: kpoints
     420             : 
     421             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_by_lattice_sum'
     422             : 
     423             :       INTEGER                                            :: handle, i_dim, i_real_imag, ikp, nkp, &
     424             :                                                             nkp_extra, nkp_orig
     425             :       INTEGER, DIMENSION(3)                              :: nkp_grid_orig, periodic
     426          22 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     427             :       TYPE(cell_type), POINTER                           :: cell
     428          22 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_RI_aux_transl, matrix_v_RI_kp
     429          22 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     430          22 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     431             : 
     432          22 :       CALL timeset(routineN, handle)
     433             : 
     434          22 :       NULLIFY (matrix_s_RI_aux_transl, particle_set, cell, qs_kind_set)
     435             : 
     436             :       CALL get_qs_env(qs_env=qs_env, &
     437             :                       matrix_s_RI_aux_kp=matrix_s_RI_aux_transl, &
     438             :                       particle_set=particle_set, &
     439             :                       cell=cell, &
     440             :                       qs_kind_set=qs_kind_set, &
     441          22 :                       atomic_kind_set=atomic_kind_set)
     442             : 
     443             :       ! check that we have a 2n x 2m x 2k mesh to guarantee good convergence
     444          22 :       CALL get_cell(cell=cell, periodic=periodic)
     445          88 :       DO i_dim = 1, 3
     446          88 :          IF (periodic(i_dim) == 1) THEN
     447          40 :             CPASSERT(MODULO(kpoints%nkp_grid(i_dim), 2) == 0)
     448             :          END IF
     449             :       END DO
     450             : 
     451          22 :       nkp = kpoints%nkp
     452             : 
     453        1062 :       ALLOCATE (fm_matrix_L_kpoints(nkp, 2))
     454         498 :       DO ikp = 1, nkp
     455        1450 :          DO i_real_imag = 1, 2
     456             :             CALL cp_fm_create(fm_matrix_L_kpoints(ikp, i_real_imag), &
     457         952 :                               fm_matrix_Minv_L_kpoints(1, i_real_imag)%matrix_struct)
     458        1428 :             CALL cp_fm_set_all(fm_matrix_L_kpoints(ikp, i_real_imag), 0.0_dp)
     459             :          END DO
     460             :       END DO
     461             : 
     462          22 :       CALL allocate_matrix_v_RI_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
     463             : 
     464          22 :       IF (qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints) THEN
     465             : 
     466          18 :          nkp_orig = qs_env%mp2_env%ri_rpa_im_time%nkp_orig
     467          18 :          nkp_extra = qs_env%mp2_env%ri_rpa_im_time%nkp_extra
     468             : 
     469             :          CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", &
     470             :                                          cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
     471             :                                          atomic_kind_set=atomic_kind_set, &
     472             :                                          size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
     473          18 :                                          operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp_orig)
     474             : 
     475          72 :          nkp_grid_orig = kpoints%nkp_grid
     476         144 :          kpoints%nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid_extra(1:3)
     477             : 
     478             :          CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", &
     479             :                                          cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
     480             :                                          atomic_kind_set=atomic_kind_set, &
     481             :                                          size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
     482          18 :                                          operator_type=operator_coulomb, ikp_start=nkp_orig + 1, ikp_end=nkp)
     483             : 
     484          72 :          kpoints%nkp_grid = nkp_grid_orig
     485             : 
     486             :       ELSE
     487             : 
     488             :          CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", &
     489             :                                          cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
     490             :                                          atomic_kind_set=atomic_kind_set, &
     491             :                                          size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
     492           4 :                                          operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp)
     493             : 
     494             :       END IF
     495             : 
     496         498 :       DO ikp = 1, nkp
     497             : 
     498         476 :          CALL copy_dbcsr_to_fm(matrix_v_RI_kp(ikp, 1)%matrix, fm_matrix_L_kpoints(ikp, 1))
     499         498 :          CALL copy_dbcsr_to_fm(matrix_v_RI_kp(ikp, 2)%matrix, fm_matrix_L_kpoints(ikp, 2))
     500             : 
     501             :       END DO
     502             : 
     503          22 :       CALL dbcsr_deallocate_matrix_set(matrix_v_RI_kp)
     504             : 
     505          22 :       CALL timestop(handle)
     506             : 
     507          22 :    END SUBROUTINE compute_V_by_lattice_sum
     508             : 
     509             : ! **************************************************************************************************
     510             : !> \brief ...
     511             : !> \param matrix_v_RI_kp ...
     512             : !> \param matrix_s_RI_aux_transl ...
     513             : !> \param nkp ...
     514             : ! **************************************************************************************************
     515          22 :    SUBROUTINE allocate_matrix_v_RI_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
     516             : 
     517             :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_v_RI_kp, matrix_s_RI_aux_transl
     518             :       INTEGER                                            :: nkp
     519             : 
     520             :       INTEGER                                            :: ikp
     521             : 
     522          22 :       NULLIFY (matrix_v_RI_kp)
     523          22 :       CALL dbcsr_allocate_matrix_set(matrix_v_RI_kp, nkp, 2)
     524             : 
     525         498 :       DO ikp = 1, nkp
     526             : 
     527         476 :          ALLOCATE (matrix_v_RI_kp(ikp, 1)%matrix)
     528             :          CALL dbcsr_create(matrix_v_RI_kp(ikp, 1)%matrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
     529         476 :                            matrix_type=dbcsr_type_no_symmetry)
     530         476 :          CALL dbcsr_reserve_all_blocks(matrix_v_RI_kp(ikp, 1)%matrix)
     531         476 :          CALL dbcsr_set(matrix_v_RI_kp(ikp, 1)%matrix, 0.0_dp)
     532             : 
     533         476 :          ALLOCATE (matrix_v_RI_kp(ikp, 2)%matrix)
     534             :          CALL dbcsr_create(matrix_v_RI_kp(ikp, 2)%matrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
     535         476 :                            matrix_type=dbcsr_type_no_symmetry)
     536         476 :          CALL dbcsr_reserve_all_blocks(matrix_v_RI_kp(ikp, 2)%matrix)
     537         498 :          CALL dbcsr_set(matrix_v_RI_kp(ikp, 2)%matrix, 0.0_dp)
     538             : 
     539             :       END DO
     540             : 
     541          22 :    END SUBROUTINE allocate_matrix_v_RI_kp
     542             : 
     543             : ! **************************************************************************************************
     544             : !> \brief ...
     545             : !> \param qs_env ...
     546             : !> \param fm_matrix_Minv_L_kpoints ...
     547             : !> \param fm_matrix_L ...
     548             : !> \param dimen_RI ...
     549             : !> \param ri_metric ...
     550             : !> \param do_kpoints ...
     551             : !> \param kpoints ...
     552             : !> \param put_mat_KS_env ...
     553             : !> \param regularization_RI ...
     554             : !> \param ikp_ext ...
     555             : !> \param do_build_cell_index ...
     556             : ! **************************************************************************************************
     557         494 :    SUBROUTINE RI_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L, dimen_RI, ri_metric, &
     558             :                                  do_kpoints, kpoints, put_mat_KS_env, regularization_RI, ikp_ext, &
     559             :                                  do_build_cell_index)
     560             : 
     561             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     562             :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_matrix_Minv_L_kpoints
     563             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_matrix_L
     564             :       INTEGER, INTENT(IN)                                :: dimen_RI
     565             :       TYPE(libint_potential_type)                        :: ri_metric
     566             :       LOGICAL, INTENT(IN)                                :: do_kpoints
     567             :       TYPE(kpoint_type), OPTIONAL, POINTER               :: kpoints
     568             :       LOGICAL, OPTIONAL                                  :: put_mat_KS_env
     569             :       REAL(KIND=dp), OPTIONAL                            :: regularization_RI
     570             :       INTEGER, OPTIONAL                                  :: ikp_ext
     571             :       LOGICAL, OPTIONAL                                  :: do_build_cell_index
     572             : 
     573             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'RI_2c_integral_mat'
     574             : 
     575             :       INTEGER                                            :: handle, i_real_imag, i_size, ikp, &
     576             :                                                             ikp_for_xkp, img, n_real_imag, natom, &
     577             :                                                             nimg, nkind, nkp
     578         494 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sizes_RI
     579         494 :       INTEGER, DIMENSION(:), POINTER                     :: col_bsize, row_bsize
     580         494 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     581             :       LOGICAL                                            :: my_do_build_cell_index, my_put_mat_KS_env
     582             :       REAL(KIND=dp)                                      :: my_regularization_RI
     583         494 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
     584             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     585             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     586             :       TYPE(cp_fm_type)                                   :: fm_matrix_S_global
     587             :       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
     588         494 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s_RI_aux_transl
     589         494 :       TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: mat_2c
     590             :       TYPE(dbcsr_type), POINTER                          :: cmatrix, matrix_s_RI_aux_desymm, &
     591             :                                                             rmatrix, tmpmat
     592             :       TYPE(dft_control_type), POINTER                    :: dft_control
     593             :       TYPE(distribution_2d_type), POINTER                :: dist_2d
     594         494 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_RI
     595             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     596             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     597         494 :          POINTER                                         :: sab_RI
     598         494 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     599         494 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     600             : 
     601         494 :       CALL timeset(routineN, handle)
     602             : 
     603         494 :       NULLIFY (sab_RI, matrix_s_RI_aux_transl, dist_2d)
     604             : 
     605         494 :       IF (PRESENT(regularization_RI)) THEN
     606         446 :          my_regularization_RI = regularization_RI
     607             :       ELSE
     608          48 :          my_regularization_RI = 0.0_dp
     609             :       END IF
     610             : 
     611         494 :       IF (PRESENT(put_mat_KS_env)) THEN
     612         180 :          my_put_mat_KS_env = put_mat_KS_env
     613             :       ELSE
     614             :          my_put_mat_KS_env = .FALSE.
     615             :       END IF
     616             : 
     617         494 :       IF (PRESENT(do_build_cell_index)) THEN
     618         266 :          my_do_build_cell_index = do_build_cell_index
     619             :       ELSE
     620             :          my_do_build_cell_index = .FALSE.
     621             :       END IF
     622             : 
     623             :       CALL get_qs_env(qs_env=qs_env, &
     624             :                       para_env=para_env, &
     625             :                       blacs_env=blacs_env, &
     626             :                       nkind=nkind, &
     627             :                       distribution_2d=dist_2d, &
     628             :                       qs_kind_set=qs_kind_set, &
     629             :                       particle_set=particle_set, &
     630             :                       dft_control=dft_control, &
     631         494 :                       natom=natom)
     632             : 
     633        1482 :       ALLOCATE (sizes_RI(natom))
     634        2434 :       ALLOCATE (basis_set_RI(nkind))
     635             : 
     636         494 :       CALL basis_set_list_setup(basis_set_RI, 'RI_AUX', qs_kind_set)
     637         494 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_RI)
     638             : 
     639             :       CALL build_2c_neighbor_lists(sab_RI, basis_set_RI, basis_set_RI, ri_metric, "RPA_2c_nl_RI", qs_env, &
     640         494 :                                    sym_ij=.TRUE., dist_2d=dist_2d)
     641             : 
     642         494 :       CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
     643        1482 :       ALLOCATE (row_bsize(SIZE(sizes_RI)))
     644        1482 :       ALLOCATE (col_bsize(SIZE(sizes_RI)))
     645        1796 :       row_bsize(:) = sizes_RI
     646        1796 :       col_bsize(:) = sizes_RI
     647             : 
     648         494 :       IF (do_kpoints) THEN
     649         288 :          CPASSERT(PRESENT(kpoints))
     650         288 :          IF (my_do_build_cell_index) THEN
     651          16 :             CALL kpoint_init_cell_index(kpoints, sab_RI, para_env, dft_control)
     652             :          END IF
     653             :          CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, &
     654         288 :                               cell_to_index=cell_to_index)
     655         288 :          n_real_imag = 2
     656         288 :          nimg = dft_control%nimages
     657             :       ELSE
     658         206 :          nkp = 1
     659         206 :          n_real_imag = 1
     660         206 :          nimg = 1
     661             :       END IF
     662             : 
     663        3094 :       ALLOCATE (mat_2c(nimg))
     664             :       CALL dbcsr_create(mat_2c(1), "(RI|RI)", dbcsr_dist, dbcsr_type_symmetric, &
     665         494 :                         row_bsize, col_bsize)
     666         494 :       DEALLOCATE (row_bsize, col_bsize)
     667             : 
     668        1612 :       DO img = 2, nimg
     669        1612 :          CALL dbcsr_create(mat_2c(img), template=mat_2c(1))
     670             :       END DO
     671             : 
     672             :       CALL build_2c_integrals(mat_2c, 0.0_dp, qs_env, sab_RI, basis_set_RI, basis_set_RI, &
     673             :                               ri_metric, do_kpoints=do_kpoints, ext_kpoints=kpoints, &
     674         494 :                               regularization_RI=my_regularization_RI)
     675             : 
     676         494 :       CALL dbcsr_distribution_release(dbcsr_dist)
     677         494 :       DEALLOCATE (basis_set_RI)
     678             : 
     679         494 :       IF (my_put_mat_KS_env) THEN
     680         136 :          CALL get_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
     681         136 :          IF (ASSOCIATED(matrix_s_RI_aux_transl)) CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
     682             :       END IF
     683         494 :       CALL dbcsr_allocate_matrix_set(matrix_s_RI_aux_transl, 1, nimg)
     684        2106 :       DO img = 1, nimg
     685        1612 :          ALLOCATE (matrix_s_RI_aux_transl(1, img)%matrix)
     686        1612 :          CALL dbcsr_copy(matrix_s_RI_aux_transl(1, img)%matrix, mat_2c(img))
     687        2106 :          CALL dbcsr_release(mat_2c(img))
     688             :       END DO
     689             : 
     690         494 :       IF (my_put_mat_KS_env) THEN
     691         136 :          CALL set_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
     692             :       END IF
     693             : 
     694         494 :       IF (PRESENT(ikp_ext)) nkp = 1
     695             : 
     696        4448 :       ALLOCATE (fm_matrix_Minv_L_kpoints(nkp, n_real_imag))
     697        1442 :       DO i_size = 1, nkp
     698        3132 :          DO i_real_imag = 1, n_real_imag
     699        1690 :             CALL cp_fm_create(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), fm_matrix_L%matrix_struct)
     700        2638 :             CALL cp_fm_set_all(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), 0.0_dp)
     701             :          END DO
     702             :       END DO
     703             : 
     704         494 :       NULLIFY (fm_struct)
     705             :       CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
     706         494 :                                ncol_global=dimen_RI, para_env=para_env)
     707             : 
     708         494 :       CALL cp_fm_create(fm_matrix_S_global, fm_struct)
     709         494 :       CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
     710             : 
     711         494 :       IF (do_kpoints) THEN
     712             : 
     713         288 :          ALLOCATE (rmatrix, cmatrix, tmpmat)
     714             :          CALL dbcsr_create(rmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
     715         288 :                            matrix_type=dbcsr_type_symmetric)
     716             :          CALL dbcsr_create(cmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
     717         288 :                            matrix_type=dbcsr_type_antisymmetric)
     718             :          CALL dbcsr_create(tmpmat, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
     719         288 :                            matrix_type=dbcsr_type_no_symmetry)
     720         288 :          CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_RI)
     721         288 :          CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_RI)
     722             : 
     723        1030 :          DO ikp = 1, nkp
     724             : 
     725             :             ! s matrix is not spin dependent, double the work
     726         742 :             CALL dbcsr_set(rmatrix, 0.0_dp)
     727         742 :             CALL dbcsr_set(cmatrix, 0.0_dp)
     728             : 
     729         742 :             IF (PRESENT(ikp_ext)) THEN
     730         266 :                ikp_for_xkp = ikp_ext
     731             :             ELSE
     732             :                ikp_for_xkp = ikp
     733             :             END IF
     734             : 
     735             :             CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s_RI_aux_transl, ispin=1, &
     736         742 :                                 xkp=xkp(1:3, ikp_for_xkp), cell_to_index=cell_to_index, sab_nl=sab_RI)
     737             : 
     738         742 :             CALL dbcsr_set(tmpmat, 0.0_dp)
     739         742 :             CALL dbcsr_desymmetrize(rmatrix, tmpmat)
     740             : 
     741         742 :             CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
     742         742 :             CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
     743         742 :             CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 1), para_env)
     744             : 
     745         742 :             CALL dbcsr_set(tmpmat, 0.0_dp)
     746         742 :             CALL dbcsr_desymmetrize(cmatrix, tmpmat)
     747             : 
     748         742 :             CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
     749         742 :             CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
     750        1030 :             CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 2), para_env)
     751             : 
     752             :          END DO
     753             : 
     754         288 :          CALL dbcsr_deallocate_matrix(rmatrix)
     755         288 :          CALL dbcsr_deallocate_matrix(cmatrix)
     756         288 :          CALL dbcsr_deallocate_matrix(tmpmat)
     757             : 
     758             :       ELSE
     759             : 
     760             :          NULLIFY (matrix_s_RI_aux_desymm)
     761         206 :          ALLOCATE (matrix_s_RI_aux_desymm)
     762             :          CALL dbcsr_create(matrix_s_RI_aux_desymm, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
     763         206 :                            name='S_RI non_symm', matrix_type=dbcsr_type_no_symmetry)
     764             : 
     765         206 :          CALL dbcsr_desymmetrize(matrix_s_RI_aux_transl(1, 1)%matrix, matrix_s_RI_aux_desymm)
     766             : 
     767         206 :          CALL copy_dbcsr_to_fm(matrix_s_RI_aux_desymm, fm_matrix_S_global)
     768             : 
     769         206 :          CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(1, 1), para_env)
     770             : 
     771         206 :          CALL dbcsr_deallocate_matrix(matrix_s_RI_aux_desymm)
     772             : 
     773             :       END IF
     774             : 
     775         494 :       CALL release_neighbor_list_sets(sab_RI)
     776             : 
     777         494 :       CALL cp_fm_struct_release(fm_struct)
     778             : 
     779         494 :       CALL cp_fm_release(fm_matrix_S_global)
     780             : 
     781         494 :       IF (.NOT. my_put_mat_KS_env) THEN
     782         358 :          CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
     783             :       END IF
     784             : 
     785         494 :       CALL timestop(handle)
     786             : 
     787        1976 :    END SUBROUTINE RI_2c_integral_mat
     788             : 
     789             : ! **************************************************************************************************
     790             : !> \brief ...
     791             : !> \param qs_env ...
     792             : !> \param eri_method ...
     793             : !> \param eri_param ...
     794             : !> \param para_env ...
     795             : !> \param para_env_sub ...
     796             : !> \param para_env_L ...
     797             : !> \param mp2_memory ...
     798             : !> \param fm_matrix_L ...
     799             : !> \param ngroup ...
     800             : !> \param color_sub ...
     801             : !> \param dimen_RI ...
     802             : !> \param my_group_L_size ...
     803             : !> \param my_group_L_start ...
     804             : !> \param my_group_L_end ...
     805             : !> \param gd_array ...
     806             : !> \param calc_PQ_cond_num ...
     807             : !> \param cond_num ...
     808             : !> \param num_small_eigen ...
     809             : !> \param potential ...
     810             : !> \param sab_orb_sub ...
     811             : !> \param do_im_time ...
     812             : !> \param fm_matrix_L_extern ...
     813             : ! **************************************************************************************************
     814         694 :    SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
     815             :                                    fm_matrix_L, ngroup, color_sub, dimen_RI, &
     816             :                                    my_group_L_size, my_group_L_start, my_group_L_end, &
     817             :                                    gd_array, calc_PQ_cond_num, cond_num, num_small_eigen, potential, &
     818             :                                    sab_orb_sub, do_im_time, fm_matrix_L_extern)
     819             : 
     820             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     821             :       INTEGER, INTENT(IN)                                :: eri_method
     822             :       TYPE(cp_eri_mme_param), POINTER                    :: eri_param
     823             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     824             :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env_sub
     825             :       TYPE(mp_para_env_type), INTENT(OUT), POINTER       :: para_env_L
     826             :       REAL(KIND=dp), INTENT(IN)                          :: mp2_memory
     827             :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_matrix_L
     828             :       INTEGER, INTENT(IN)                                :: ngroup, color_sub
     829             :       INTEGER, INTENT(OUT)                               :: dimen_RI, my_group_L_size, &
     830             :                                                             my_group_L_start, my_group_L_end
     831             :       TYPE(group_dist_d1_type), INTENT(OUT)              :: gd_array
     832             :       LOGICAL, INTENT(IN)                                :: calc_PQ_cond_num
     833             :       REAL(KIND=dp), INTENT(OUT)                         :: cond_num
     834             :       INTEGER, INTENT(OUT)                               :: num_small_eigen
     835             :       TYPE(libint_potential_type), INTENT(IN)            :: potential
     836             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     837             :          POINTER                                         :: sab_orb_sub
     838             :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_im_time
     839             :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: fm_matrix_L_extern
     840             : 
     841             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_2c_integrals'
     842             : 
     843             :       INTEGER :: best_group_size, color_L, group_size, handle, handle2, i_global, iatom, iiB, &
     844             :          ikind, iproc, j_global, jjB, natom, ncol_local, nkind, nrow_local, nsgf, potential_type, &
     845             :          proc_receive, proc_receive_static, proc_send_static, proc_shift, rec_L_end, rec_L_size, &
     846             :          rec_L_start, strat_group_size
     847         694 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     848         694 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     849             :       LOGICAL                                            :: my_do_im_time
     850             :       REAL(KIND=dp)                                      :: min_mem_for_QK
     851         694 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: egen_L
     852         694 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: L_external_col, L_local_col
     853         694 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     854             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_L
     855             :       TYPE(cp_fm_type)                                   :: fm_matrix_L_diag
     856         694 :       TYPE(group_dist_d1_type)                           :: gd_sub_array
     857             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a
     858         694 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     859         694 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     860             : 
     861         694 :       CALL timeset(routineN, handle)
     862             : 
     863         694 :       my_do_im_time = .FALSE.
     864         694 :       IF (PRESENT(do_im_time)) THEN
     865         658 :          my_do_im_time = do_im_time
     866             :       END IF
     867             : 
     868             :       ! get stuff
     869             :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
     870         694 :                       particle_set=particle_set)
     871             : 
     872         694 :       nkind = SIZE(qs_kind_set)
     873         694 :       natom = SIZE(particle_set)
     874             : 
     875         694 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     876             : 
     877        1974 :       DO ikind = 1, nkind
     878        1280 :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
     879        1974 :          CPASSERT(ASSOCIATED(basis_set_a))
     880             :       END DO
     881             : 
     882         694 :       dimen_RI = 0
     883        2734 :       DO iatom = 1, natom
     884        2040 :          ikind = kind_of(iatom)
     885        2040 :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
     886        2734 :          dimen_RI = dimen_RI + nsgf
     887             :       END DO
     888             : 
     889             :       ! check that very small systems are not running on too many processes
     890         694 :       IF (dimen_RI < ngroup) THEN
     891             :          CALL cp_abort(__LOCATION__, "Product of block size and number "// &
     892           0 :                        "of RI functions should not exceed total number of processes")
     893             :       END IF
     894             : 
     895         694 :       CALL create_group_dist(gd_array, ngroup, dimen_RI)
     896             : 
     897         694 :       CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
     898             : 
     899         694 :       CALL timeset(routineN//"_loop_lm", handle2)
     900             : 
     901        2776 :       ALLOCATE (L_local_col(dimen_RI, my_group_L_size))
     902     2385102 :       L_local_col = 0.0_dp
     903             : 
     904         694 :       potential_type = potential%potential_type
     905             : 
     906             :       IF ((eri_method .EQ. do_eri_mme .OR. eri_method .EQ. do_eri_os) &
     907         694 :           .AND. potential_type .EQ. do_potential_coulomb .OR. &
     908             :           (eri_method .EQ. do_eri_mme .AND. potential_type .EQ. do_potential_long)) THEN
     909             : 
     910             :          CALL mp2_eri_2c_integrate(eri_param, potential, para_env_sub, qs_env, &
     911             :                                    basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
     912             :                                    hab=L_local_col, first_b=my_group_L_start, last_b=my_group_L_end, &
     913         306 :                                    eri_method=eri_method)
     914             : 
     915             :       ELSEIF (eri_method .EQ. do_eri_gpw .OR. &
     916             :               (potential_type == do_potential_long .AND. qs_env%mp2_env%eri_method == do_eri_os) &
     917         388 :               .OR. (potential_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)) THEN
     918             : 
     919             :          CALL mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, &
     920         388 :                                        natom, potential, sab_orb_sub, L_local_col, kind_of)
     921             : 
     922             :       ELSE
     923           0 :          CPABORT("unknown ERI method")
     924             :       END IF
     925             : 
     926         694 :       CALL timestop(handle2)
     927             : 
     928             :       ! split the total number of proc in a subgroup of the size of ~1/10 of the
     929             :       ! total num of proc
     930         694 :       best_group_size = para_env%num_pe
     931             : 
     932         694 :       strat_group_size = MAX(1, para_env%num_pe/10)
     933             : 
     934         694 :       min_mem_for_QK = REAL(dimen_RI, KIND=dp)*dimen_RI*3.0_dp*8.0_dp/1024_dp/1024_dp
     935             : 
     936         694 :       group_size = strat_group_size - 1
     937         728 :       DO iproc = strat_group_size, para_env%num_pe
     938         728 :          group_size = group_size + 1
     939             :          ! check that group_size is a multiple of sub_group_size and a divisor of
     940             :          ! the total num of proc
     941         728 :          IF (MOD(para_env%num_pe, group_size) /= 0 .OR. MOD(group_size, para_env_sub%num_pe) /= 0) CYCLE
     942             : 
     943             :          ! check for memory
     944         694 :          IF (REAL(group_size, KIND=dp)*mp2_memory < min_mem_for_QK) CYCLE
     945             : 
     946             :          best_group_size = group_size
     947          34 :          EXIT
     948             :       END DO
     949             : 
     950         694 :       IF (my_do_im_time) THEN
     951             :          ! para_env_L is the para_env for im. time to avoid bug
     952         136 :          best_group_size = para_env%num_pe
     953             :       END IF
     954             : 
     955             :       ! create the L group
     956         694 :       color_L = para_env%mepos/best_group_size
     957         694 :       ALLOCATE (para_env_L)
     958         694 :       CALL para_env_L%from_split(para_env, color_L)
     959             : 
     960             :       ! create the blacs_L
     961         694 :       NULLIFY (blacs_env_L)
     962         694 :       CALL cp_blacs_env_create(blacs_env=blacs_env_L, para_env=para_env_L)
     963             : 
     964             :       ! create the full matrix L defined in the L group
     965         694 :       CALL create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, "fm_matrix_L", fm_matrix_L_extern)
     966             : 
     967         694 :       IF (my_do_im_time .AND. para_env%num_pe > 1) THEN
     968             : 
     969             :          CALL fill_fm_L_from_L_loc_non_blocking(fm_matrix_L, L_local_col, para_env, &
     970             :                                                 my_group_L_start, my_group_L_end, &
     971         136 :                                                 dimen_RI)
     972             : 
     973             :       ELSE
     974         558 :          BLOCK
     975             :             TYPE(mp_comm_type) :: comm_exchange
     976             : 
     977         558 :             CALL comm_exchange%from_split(para_env_L, para_env_sub%mepos)
     978             : 
     979             :             CALL create_group_dist(gd_sub_array, my_group_L_start, &
     980         558 :                                    my_group_L_end, my_group_L_size, comm_exchange)
     981             : 
     982             :             CALL cp_fm_get_info(matrix=fm_matrix_L, &
     983             :                                 nrow_local=nrow_local, &
     984             :                                 ncol_local=ncol_local, &
     985             :                                 row_indices=row_indices, &
     986         558 :                                 col_indices=col_indices)
     987             : 
     988       42336 :             DO jjB = 1, ncol_local
     989       41778 :                j_global = col_indices(jjB)
     990       42336 :                IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
     991     1723644 :                   DO iiB = 1, nrow_local
     992     1701635 :                      i_global = row_indices(iiB)
     993     1723644 :                      fm_matrix_L%local_data(iiB, jjB) = L_local_col(i_global, j_global - my_group_L_start + 1)
     994             :                   END DO
     995             :                END IF
     996             :             END DO
     997             : 
     998         558 :             proc_send_static = MODULO(comm_exchange%mepos + 1, comm_exchange%num_pe)
     999         558 :             proc_receive_static = MODULO(comm_exchange%mepos - 1, comm_exchange%num_pe)
    1000             : 
    1001         558 :             DO proc_shift = 1, comm_exchange%num_pe - 1
    1002           0 :                proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
    1003             : 
    1004           0 :                CALL get_group_dist(gd_sub_array, proc_receive, rec_L_start, rec_L_end, rec_L_size)
    1005             : 
    1006           0 :                ALLOCATE (L_external_col(dimen_RI, rec_L_size))
    1007           0 :                L_external_col = 0.0_dp
    1008             : 
    1009           0 :                CALL comm_exchange%sendrecv(L_local_col, proc_send_static, L_external_col, proc_receive_static)
    1010             : 
    1011           0 :                DO jjB = 1, ncol_local
    1012           0 :                   j_global = col_indices(jjB)
    1013           0 :                   IF (j_global >= rec_L_start .AND. j_global <= rec_L_end) THEN
    1014           0 :                      DO iiB = 1, nrow_local
    1015           0 :                         i_global = row_indices(iiB)
    1016           0 :                         fm_matrix_L%local_data(iiB, jjB) = L_external_col(i_global, j_global - rec_L_start + 1)
    1017             :                      END DO
    1018             :                   END IF
    1019             :                END DO
    1020             : 
    1021         558 :                CALL MOVE_ALLOC(L_external_col, L_local_col)
    1022             :             END DO
    1023             : 
    1024         558 :             CALL release_group_dist(gd_sub_array)
    1025        1116 :             CALL comm_exchange%free()
    1026             :          END BLOCK
    1027             :       END IF
    1028             : 
    1029         694 :       DEALLOCATE (L_local_col)
    1030             : 
    1031             :       ! create the new group for the sum of the local data over all processes
    1032             :       BLOCK
    1033             :          TYPE(mp_comm_type) :: comm_exchange
    1034         694 :          comm_exchange = fm_matrix_L%matrix_struct%context%interconnect(para_env)
    1035         694 :          CALL comm_exchange%sum(fm_matrix_L%local_data)
    1036        1388 :          CALL comm_exchange%free()
    1037             :       END BLOCK
    1038             : 
    1039         694 :       cond_num = 1.0_dp
    1040         694 :       num_small_eigen = 0
    1041         694 :       IF (calc_PQ_cond_num) THEN
    1042             :          ! calculate the condition number of the (P|Q) matrix
    1043             :          ! create a copy of the matrix
    1044             : 
    1045           0 :          CALL create_matrix_L(fm_matrix_L_diag, blacs_env_L, dimen_RI, para_env_L, "fm_matrix_L_diag", fm_matrix_L_extern)
    1046             : 
    1047           0 :          CALL cp_fm_to_fm(source=fm_matrix_L, destination=fm_matrix_L_diag)
    1048             : 
    1049           0 :          ALLOCATE (egen_L(dimen_RI))
    1050             : 
    1051           0 :          egen_L = 0.0_dp
    1052           0 :          CALL cp_fm_syevx(matrix=fm_matrix_L_diag, eigenvalues=egen_L)
    1053             : 
    1054           0 :          num_small_eigen = 0
    1055           0 :          DO iiB = 1, dimen_RI
    1056           0 :             IF (ABS(egen_L(iiB)) < 0.001_dp) num_small_eigen = num_small_eigen + 1
    1057             :          END DO
    1058             : 
    1059           0 :          cond_num = MAXVAL(ABS(egen_L))/MINVAL(ABS(egen_L))
    1060             : 
    1061           0 :          CALL cp_fm_release(fm_matrix_L_diag)
    1062             : 
    1063           0 :          DEALLOCATE (egen_L)
    1064             :       END IF
    1065             : 
    1066             :       ! release blacs_env
    1067         694 :       CALL cp_blacs_env_release(blacs_env_L)
    1068             : 
    1069         694 :       CALL timestop(handle)
    1070             : 
    1071        3470 :    END SUBROUTINE compute_2c_integrals
    1072             : 
    1073             : ! **************************************************************************************************
    1074             : !> \brief ...
    1075             : !> \param matrix ...
    1076             : !> \param num_small_evals ...
    1077             : !> \param cond_num ...
    1078             : !> \param eps_svd ...
    1079             : !> \param do_inversion ...
    1080             : !> \param para_env ...
    1081             : ! **************************************************************************************************
    1082          44 :    SUBROUTINE matrix_root_with_svd(matrix, num_small_evals, cond_num, eps_svd, do_inversion, para_env)
    1083             :       TYPE(cp_fm_type), INTENT(INOUT)                    :: matrix
    1084             :       INTEGER, INTENT(OUT)                               :: num_small_evals
    1085             :       REAL(KIND=dp), INTENT(OUT)                         :: cond_num
    1086             :       REAL(KIND=dp), INTENT(IN)                          :: eps_svd
    1087             :       LOGICAL, INTENT(IN)                                :: do_inversion
    1088             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1089             : 
    1090             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_root_with_svd'
    1091             : 
    1092             :       INTEGER                                            :: group_size_L, handle, ii, needed_evals, &
    1093             :                                                             nrow, pos_max(1)
    1094          44 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: num_eval
    1095             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: evals
    1096             :       TYPE(cp_fm_type)                                   :: evecs
    1097             :       TYPE(mp_comm_type)                                 :: comm_exchange
    1098             : 
    1099          44 :       CALL timeset(routineN, handle)
    1100             : 
    1101             :       ! Create arrays carrying eigenvalues and eigenvectors later
    1102          44 :       CALL cp_fm_get_info(matrix=matrix, nrow_global=nrow)
    1103             : 
    1104         132 :       ALLOCATE (evals(nrow))
    1105        3724 :       evals = 0
    1106             : 
    1107          44 :       CALL cp_fm_create(evecs, matrix%matrix_struct)
    1108             : 
    1109             :       ! Perform the EVD (=SVD of a positive semidefinite matrix)
    1110          44 :       CALL choose_eigv_solver(matrix, evecs, evals)
    1111             : 
    1112             :       ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
    1113          44 :       num_small_evals = 0
    1114         202 :       DO ii = 1, nrow
    1115         202 :          IF (evals(ii) > eps_svd) THEN
    1116          44 :             num_small_evals = ii - 1
    1117          44 :             EXIT
    1118             :          END IF
    1119             :       END DO
    1120          44 :       needed_evals = nrow - num_small_evals
    1121             : 
    1122             :       ! Get the condition number w.r.t. considered eigenvalues
    1123          44 :       cond_num = evals(nrow)/evals(num_small_evals + 1)
    1124             : 
    1125             :       ! Determine the eigenvalues of the request matrix root or its inverse
    1126         202 :       evals(1:num_small_evals) = 0.0_dp
    1127          44 :       IF (do_inversion) THEN
    1128        1008 :          evals(num_small_evals + 1:nrow) = 1.0_dp/SQRT(evals(num_small_evals + 1:nrow))
    1129             :       ELSE
    1130        2558 :          evals(num_small_evals + 1:nrow) = SQRT(evals(num_small_evals + 1:nrow))
    1131             :       END IF
    1132             : 
    1133          44 :       CALL cp_fm_column_scale(evecs, evals)
    1134             : 
    1135             :       ! As it turns out, the results in the subgroups differ.
    1136             :       ! Thus, we have to equilibrate the results.
    1137             :       ! Step 1: Get a communicator connecting processes with same id within subgroups
    1138             :       ! use that group_size_L is divisible by the total number of processes (see above)
    1139          44 :       group_size_L = para_env%num_pe/matrix%matrix_struct%para_env%num_pe
    1140          44 :       comm_exchange = matrix%matrix_struct%context%interconnect(para_env)
    1141             : 
    1142         132 :       ALLOCATE (num_eval(0:group_size_L - 1))
    1143          44 :       CALL comm_exchange%allgather(num_small_evals, num_eval)
    1144             : 
    1145         118 :       num_small_evals = MINVAL(num_eval)
    1146             : 
    1147         118 :       IF (num_small_evals /= MAXVAL(num_eval)) THEN
    1148             :          ! Step 2: Get position of maximum value
    1149           0 :          pos_max = MAXLOC(num_eval)
    1150           0 :          num_small_evals = num_eval(pos_max(1))
    1151           0 :          needed_evals = nrow - num_small_evals
    1152             : 
    1153             :          ! Step 3: Broadcast your local data to all other processes
    1154           0 :          CALL comm_exchange%bcast(evecs%local_data, pos_max(1))
    1155           0 :          CALL comm_exchange%bcast(cond_num, pos_max(1))
    1156             :       END IF
    1157             : 
    1158          44 :       DEALLOCATE (num_eval)
    1159             : 
    1160          44 :       CALL comm_exchange%free()
    1161             : 
    1162          44 :       CALL reset_size_matrix(matrix, needed_evals, matrix%matrix_struct)
    1163             : 
    1164             :       ! Copy the needed eigenvectors
    1165          44 :       CALL cp_fm_to_fm(evecs, matrix, needed_evals, num_small_evals + 1)
    1166             : 
    1167          44 :       CALL cp_fm_release(evecs)
    1168             : 
    1169          44 :       CALL timestop(handle)
    1170             : 
    1171         132 :    END SUBROUTINE matrix_root_with_svd
    1172             : 
    1173             : ! **************************************************************************************************
    1174             : !> \brief ...
    1175             : !> \param matrix ...
    1176             : !> \param new_size ...
    1177             : !> \param fm_struct_template ...
    1178             : ! **************************************************************************************************
    1179         152 :    SUBROUTINE reset_size_matrix(matrix, new_size, fm_struct_template)
    1180             :       TYPE(cp_fm_type), INTENT(INOUT)                    :: matrix
    1181             :       INTEGER, INTENT(IN)                                :: new_size
    1182             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_template
    1183             : 
    1184             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'reset_size_matrix'
    1185             : 
    1186             :       INTEGER                                            :: handle
    1187             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    1188             : 
    1189          76 :       CALL timeset(routineN, handle)
    1190             : 
    1191             :       ! Choose the block sizes as large as in the template for the later steps
    1192          76 :       NULLIFY (fm_struct)
    1193          76 :       CALL cp_fm_struct_create(fm_struct, ncol_global=new_size, template_fmstruct=fm_struct_template, force_block=.TRUE.)
    1194             : 
    1195          76 :       CALL cp_fm_release(matrix)
    1196             : 
    1197          76 :       CALL cp_fm_create(matrix, fm_struct)
    1198          76 :       CALL cp_fm_set_all(matrix, 0.0_dp)
    1199             : 
    1200          76 :       CALL cp_fm_struct_release(fm_struct)
    1201             : 
    1202          76 :       CALL timestop(handle)
    1203             : 
    1204          76 :    END SUBROUTINE reset_size_matrix
    1205             : 
    1206             : ! **************************************************************************************************
    1207             : !> \brief ...
    1208             : !> \param fm_matrix_L ...
    1209             : !> \param L_local_col ...
    1210             : !> \param para_env ...
    1211             : !> \param my_group_L_start ...
    1212             : !> \param my_group_L_end ...
    1213             : !> \param dimen_RI ...
    1214             : ! **************************************************************************************************
    1215         136 :    SUBROUTINE fill_fm_L_from_L_loc_non_blocking(fm_matrix_L, L_local_col, para_env, my_group_L_start, &
    1216             :                                                 my_group_L_end, dimen_RI)
    1217             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_matrix_L
    1218             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
    1219             :          INTENT(IN)                                      :: L_local_col
    1220             :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
    1221             :       INTEGER, INTENT(IN)                                :: my_group_L_start, my_group_L_end, &
    1222             :                                                             dimen_RI
    1223             : 
    1224             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_fm_L_from_L_loc_non_blocking'
    1225             : 
    1226             :       INTEGER :: dummy_proc, handle, handle2, i_entry_rec, i_row, i_row_global, iproc, j_col, &
    1227             :          j_col_global, LLL, MMM, ncol_local, nrow_local, proc_send, send_pcol, send_prow
    1228         136 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: entry_counter, num_entries_rec, &
    1229             :                                                             num_entries_send
    1230         136 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1231             :       TYPE(integ_mat_buffer_type), ALLOCATABLE, &
    1232         136 :          DIMENSION(:)                                    :: buffer_rec, buffer_send
    1233         136 :       TYPE(mp_request_type), DIMENSION(:, :), POINTER    :: req_array
    1234             : 
    1235         136 :       CALL timeset(routineN, handle)
    1236             : 
    1237         136 :       CALL timeset(routineN//"_1", handle2)
    1238             : 
    1239             :       ! get info for the dest
    1240             :       CALL cp_fm_get_info(matrix=fm_matrix_L, &
    1241             :                           nrow_local=nrow_local, &
    1242             :                           ncol_local=ncol_local, &
    1243             :                           row_indices=row_indices, &
    1244         136 :                           col_indices=col_indices)
    1245             : 
    1246         408 :       ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
    1247         272 :       ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
    1248             : 
    1249         408 :       num_entries_rec(:) = 0
    1250         408 :       num_entries_send(:) = 0
    1251             : 
    1252         136 :       dummy_proc = 0
    1253             : 
    1254             :       ! get the process, where the elements from L_local_col have to be sent
    1255       11548 :       DO LLL = 1, dimen_RI
    1256             : 
    1257       11412 :          send_prow = fm_matrix_L%matrix_struct%g2p_row(LLL)
    1258             : 
    1259      568986 :          DO MMM = my_group_L_start, my_group_L_end
    1260             : 
    1261      557438 :             send_pcol = fm_matrix_L%matrix_struct%g2p_col(MMM)
    1262             : 
    1263      557438 :             proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
    1264             : 
    1265      568850 :             num_entries_send(proc_send) = num_entries_send(proc_send) + 1
    1266             : 
    1267             :          END DO
    1268             : 
    1269             :       END DO
    1270             : 
    1271         136 :       CALL timestop(handle2)
    1272             : 
    1273         136 :       CALL timeset(routineN//"_2", handle2)
    1274             : 
    1275         136 :       CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
    1276             : 
    1277         136 :       CALL timestop(handle2)
    1278             : 
    1279         136 :       CALL timeset(routineN//"_3", handle2)
    1280             : 
    1281             :       ! allocate buffers to send the entries and the information of the entries
    1282         680 :       ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
    1283         544 :       ALLOCATE (buffer_send(0:para_env%num_pe - 1))
    1284             : 
    1285             :       ! allocate data message and corresponding indices
    1286         408 :       DO iproc = 0, para_env%num_pe - 1
    1287             : 
    1288         816 :          ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
    1289      557846 :          buffer_rec(iproc)%msg = 0.0_dp
    1290             : 
    1291             :       END DO
    1292             : 
    1293         136 :       CALL timestop(handle2)
    1294             : 
    1295         136 :       CALL timeset(routineN//"_4", handle2)
    1296             : 
    1297         408 :       DO iproc = 0, para_env%num_pe - 1
    1298             : 
    1299         816 :          ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
    1300      557846 :          buffer_send(iproc)%msg = 0.0_dp
    1301             : 
    1302             :       END DO
    1303             : 
    1304         136 :       CALL timestop(handle2)
    1305             : 
    1306         136 :       CALL timeset(routineN//"_5", handle2)
    1307             : 
    1308         408 :       DO iproc = 0, para_env%num_pe - 1
    1309             : 
    1310         816 :          ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
    1311     1115828 :          buffer_rec(iproc)%indx = 0
    1312             : 
    1313             :       END DO
    1314             : 
    1315         136 :       CALL timestop(handle2)
    1316             : 
    1317         136 :       CALL timeset(routineN//"_6", handle2)
    1318             : 
    1319         408 :       DO iproc = 0, para_env%num_pe - 1
    1320             : 
    1321         816 :          ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
    1322     1115828 :          buffer_send(iproc)%indx = 0
    1323             : 
    1324             :       END DO
    1325             : 
    1326         136 :       CALL timestop(handle2)
    1327             : 
    1328         136 :       CALL timeset(routineN//"_7", handle2)
    1329             : 
    1330         272 :       ALLOCATE (entry_counter(0:para_env%num_pe - 1))
    1331         408 :       entry_counter(:) = 0
    1332             : 
    1333             :       ! get the process, where the elements from L_local_col have to be sent and
    1334             :       ! write the data and the info to buffer_send
    1335       11548 :       DO LLL = 1, dimen_RI
    1336             : 
    1337       11412 :          send_prow = fm_matrix_L%matrix_struct%g2p_row(LLL)
    1338             : 
    1339      568986 :          DO MMM = my_group_L_start, my_group_L_end
    1340             : 
    1341      557438 :             send_pcol = fm_matrix_L%matrix_struct%g2p_col(MMM)
    1342             : 
    1343      557438 :             proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
    1344             : 
    1345      557438 :             entry_counter(proc_send) = entry_counter(proc_send) + 1
    1346             : 
    1347             :             buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
    1348      557438 :                L_local_col(LLL, MMM - my_group_L_start + 1)
    1349             : 
    1350      557438 :             buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = LLL
    1351      568850 :             buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = MMM
    1352             : 
    1353             :          END DO
    1354             : 
    1355             :       END DO
    1356             : 
    1357        2040 :       ALLOCATE (req_array(1:para_env%num_pe, 4))
    1358             : 
    1359         136 :       CALL timestop(handle2)
    1360             : 
    1361         136 :       CALL timeset(routineN//"_8", handle2)
    1362             : 
    1363             :       ! communicate the buffer
    1364             :       CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, &
    1365         136 :                               buffer_send, req_array)
    1366             : 
    1367      541430 :       fm_matrix_L%local_data = 0.0_dp
    1368             : 
    1369         136 :       CALL timestop(handle2)
    1370             : 
    1371         136 :       CALL timeset(routineN//"_9", handle2)
    1372             : 
    1373             :       ! fill fm_matrix_L with the entries from buffer_rec
    1374         408 :       DO iproc = 0, para_env%num_pe - 1
    1375             : 
    1376      557846 :          DO i_entry_rec = 1, num_entries_rec(iproc)
    1377             : 
    1378    31532766 :             DO i_row = 1, nrow_local
    1379             : 
    1380    30975056 :                i_row_global = row_indices(i_row)
    1381             : 
    1382  4692723120 :                DO j_col = 1, ncol_local
    1383             : 
    1384  4661190626 :                   j_col_global = col_indices(j_col)
    1385             : 
    1386  4661190626 :                   IF (i_row_global == buffer_rec(iproc)%indx(i_entry_rec, 1) .AND. &
    1387    30975056 :                       j_col_global == buffer_rec(iproc)%indx(i_entry_rec, 2)) THEN
    1388             : 
    1389      557438 :                      fm_matrix_L%local_data(i_row, j_col) = buffer_rec(iproc)%msg(i_entry_rec)
    1390             : 
    1391             :                   END IF
    1392             : 
    1393             :                END DO
    1394             : 
    1395             :             END DO
    1396             : 
    1397             :          END DO
    1398             : 
    1399             :       END DO
    1400             : 
    1401         136 :       CALL timestop(handle2)
    1402             : 
    1403         136 :       CALL timeset(routineN//"_10", handle2)
    1404             : 
    1405         408 :       DO iproc = 0, para_env%num_pe - 1
    1406         272 :          DEALLOCATE (buffer_rec(iproc)%msg)
    1407         272 :          DEALLOCATE (buffer_rec(iproc)%indx)
    1408         272 :          DEALLOCATE (buffer_send(iproc)%msg)
    1409         408 :          DEALLOCATE (buffer_send(iproc)%indx)
    1410             :       END DO
    1411             : 
    1412         680 :       DEALLOCATE (buffer_rec, buffer_send)
    1413         136 :       DEALLOCATE (req_array)
    1414         136 :       DEALLOCATE (entry_counter)
    1415         136 :       DEALLOCATE (num_entries_rec, num_entries_send)
    1416             : 
    1417         136 :       CALL timestop(handle2)
    1418             : 
    1419         136 :       CALL timestop(handle)
    1420             : 
    1421        1360 :    END SUBROUTINE fill_fm_L_from_L_loc_non_blocking
    1422             : 
    1423             : ! **************************************************************************************************
    1424             : !> \brief ...
    1425             : !> \param fm_matrix_L ...
    1426             : !> \param dimen_RI ...
    1427             : !> \param do_inversion ...
    1428             : ! **************************************************************************************************
    1429        1552 :    SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion)
    1430             : 
    1431             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_matrix_L
    1432             :       INTEGER, INTENT(IN)                                :: dimen_RI
    1433             :       LOGICAL, INTENT(IN)                                :: do_inversion
    1434             : 
    1435             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cholesky_decomp'
    1436             : 
    1437             :       INTEGER                                            :: handle, i_global, iiB, info_chol, &
    1438             :                                                             j_global, jjB, ncol_local, nrow_local
    1439         776 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1440             : 
    1441         776 :       CALL timeset(routineN, handle)
    1442             : 
    1443             :       ! do cholesky decomposition
    1444         776 :       CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=dimen_RI, info_out=info_chol)
    1445         776 :       CPASSERT(info_chol == 0)
    1446             : 
    1447         776 :       IF (.NOT. do_inversion) THEN
    1448             :          ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
    1449             :          CALL cp_fm_get_info(matrix=fm_matrix_L, &
    1450             :                              nrow_local=nrow_local, &
    1451             :                              ncol_local=ncol_local, &
    1452             :                              row_indices=row_indices, &
    1453         118 :                              col_indices=col_indices)
    1454        5888 :          DO iiB = 1, nrow_local
    1455        5770 :             i_global = row_indices(iiB)
    1456      545304 :             DO jjB = 1, ncol_local
    1457      539416 :                j_global = col_indices(jjB)
    1458      545186 :                IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
    1459             :             END DO
    1460             :          END DO
    1461             : 
    1462             :       END IF
    1463             : 
    1464         776 :       CALL timestop(handle)
    1465             : 
    1466         776 :    END SUBROUTINE cholesky_decomp
    1467             : 
    1468             : ! **************************************************************************************************
    1469             : !> \brief ...
    1470             : !> \param fm_matrix_Minv_L_kpoints ...
    1471             : !> \param fm_matrix_L_kpoints ...
    1472             : !> \param dimen_RI ...
    1473             : !> \param kpoints ...
    1474             : !> \param eps_eigval_S ...
    1475             : ! **************************************************************************************************
    1476          22 :    SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, &
    1477             :                                                          dimen_RI, kpoints, eps_eigval_S)
    1478             : 
    1479             :       TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN)      :: fm_matrix_Minv_L_kpoints, &
    1480             :                                                             fm_matrix_L_kpoints
    1481             :       INTEGER, INTENT(IN)                                :: dimen_RI
    1482             :       TYPE(kpoint_type), POINTER                         :: kpoints
    1483             :       REAL(KIND=dp), INTENT(IN)                          :: eps_eigval_S
    1484             : 
    1485             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'inversion_of_M_and_mult_with_chol_dec_of_V'
    1486             :       COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
    1487             :          czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp), ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
    1488             : 
    1489             :       INTEGER                                            :: handle, ikp, nkp
    1490             :       TYPE(cp_cfm_type)                                  :: cfm_matrix_K_tmp, cfm_matrix_M_tmp, &
    1491             :                                                             cfm_matrix_V_tmp, cfm_matrix_Vtrunc_tmp
    1492             :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
    1493             : 
    1494          22 :       CALL timeset(routineN, handle)
    1495             : 
    1496          22 :       CALL cp_fm_get_info(fm_matrix_Minv_L_kpoints(1, 1), matrix_struct=matrix_struct)
    1497             : 
    1498          22 :       CALL cp_cfm_create(cfm_matrix_M_tmp, matrix_struct)
    1499          22 :       CALL cp_cfm_create(cfm_matrix_V_tmp, matrix_struct)
    1500          22 :       CALL cp_cfm_create(cfm_matrix_K_tmp, matrix_struct)
    1501          22 :       CALL cp_cfm_create(cfm_matrix_Vtrunc_tmp, matrix_struct)
    1502             : 
    1503          22 :       CALL get_kpoint_info(kpoints, nkp=nkp)
    1504             : 
    1505         498 :       DO ikp = 1, nkp
    1506             : 
    1507         476 :          CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_M_tmp, cone, fm_matrix_Minv_L_kpoints(ikp, 1))
    1508         476 :          CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_M_tmp, ione, fm_matrix_Minv_L_kpoints(ikp, 2))
    1509             : 
    1510         476 :          CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_V_tmp, cone, fm_matrix_L_kpoints(ikp, 1))
    1511         476 :          CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_V_tmp, ione, fm_matrix_L_kpoints(ikp, 2))
    1512             : 
    1513         476 :          CALL cp_cfm_power(cfm_matrix_M_tmp, threshold=eps_eigval_S, exponent=-1.0_dp)
    1514             : 
    1515         476 :          CALL cp_cfm_power(cfm_matrix_V_tmp, threshold=0.0_dp, exponent=0.5_dp)
    1516             : 
    1517             :          ! save L(k) = SQRT(V(k))
    1518         476 :          CALL cp_cfm_to_fm(cfm_matrix_V_tmp, fm_matrix_L_kpoints(ikp, 1), fm_matrix_L_kpoints(ikp, 2))
    1519             : 
    1520             :          ! get K = M^(-1)*L, where L is the Cholesky decomposition of V = L^T*L
    1521             :          CALL parallel_gemm("N", "C", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_M_tmp, cfm_matrix_V_tmp, &
    1522         476 :                             czero, cfm_matrix_K_tmp)
    1523             : 
    1524             :          ! move
    1525         498 :          CALL cp_cfm_to_fm(cfm_matrix_K_tmp, fm_matrix_Minv_L_kpoints(ikp, 1), fm_matrix_Minv_L_kpoints(ikp, 2))
    1526             : 
    1527             :       END DO
    1528             : 
    1529          22 :       CALL cp_cfm_release(cfm_matrix_M_tmp)
    1530          22 :       CALL cp_cfm_release(cfm_matrix_V_tmp)
    1531          22 :       CALL cp_cfm_release(cfm_matrix_K_tmp)
    1532          22 :       CALL cp_cfm_release(cfm_matrix_Vtrunc_tmp)
    1533             : 
    1534          22 :       CALL timestop(handle)
    1535             : 
    1536          22 :    END SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V
    1537             : 
    1538             : ! **************************************************************************************************
    1539             : !> \brief ...
    1540             : !> \param fm_matrix_Minv_Vtrunc_Minv ...
    1541             : !> \param fm_matrix_Minv ...
    1542             : !> \param qs_env ...
    1543             : ! **************************************************************************************************
    1544          88 :    SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
    1545          22 :                                                                          fm_matrix_Minv, qs_env)
    1546             : 
    1547             :       TYPE(cp_fm_type), DIMENSION(:, :)                  :: fm_matrix_Minv_Vtrunc_Minv, &
    1548             :                                                             fm_matrix_Minv
    1549             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1550             : 
    1551             :       CHARACTER(LEN=*), PARAMETER :: &
    1552             :          routineN = 'Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc'
    1553             : 
    1554             :       INTEGER                                            :: dimen_RI, handle, ndep
    1555             :       REAL(KIND=dp)                                      :: eps_eigval_S_Gamma
    1556             :       TYPE(cp_fm_type)                                   :: fm_matrix_RI_metric_inv_work, fm_work
    1557             : 
    1558          22 :       CALL timeset(routineN, handle)
    1559             : 
    1560          22 :       CALL cp_fm_create(fm_work, fm_matrix_Minv(1, 1)%matrix_struct)
    1561          22 :       CALL cp_fm_set_all(fm_work, 0.0_dp)
    1562             : 
    1563          22 :       CALL cp_fm_create(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1)%matrix_struct)
    1564          22 :       CALL cp_fm_set_all(fm_matrix_RI_metric_inv_work, 0.0_dp)
    1565             : 
    1566          22 :       CALL cp_fm_get_info(matrix=fm_matrix_Minv(1, 1), nrow_global=dimen_RI)
    1567             : 
    1568          22 :       eps_eigval_S_Gamma = qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S_Gamma
    1569             : 
    1570          22 :       IF (eps_eigval_S_Gamma > 1.0E-18) THEN
    1571             : 
    1572             :          ! remove small eigenvalues
    1573             :          CALL cp_fm_power(fm_matrix_Minv(1, 1), fm_matrix_RI_metric_inv_work, -0.5_dp, &
    1574           0 :                           eps_eigval_S_Gamma, ndep)
    1575             : 
    1576             :       ELSE
    1577             : 
    1578          22 :          CALL cholesky_decomp(fm_matrix_Minv(1, 1), dimen_RI, do_inversion=.TRUE.)
    1579             : 
    1580          22 :          CALL invert_mat(fm_matrix_Minv(1, 1))
    1581             : 
    1582             :       END IF
    1583             : 
    1584             :       CALL parallel_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_Minv(1, 1), &
    1585          22 :                          fm_matrix_Minv(1, 1), 0.0_dp, fm_matrix_RI_metric_inv_work)
    1586             : 
    1587          22 :       CALL cp_fm_to_fm(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1))
    1588             : 
    1589             :       CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_RI_metric_inv_work, &
    1590          22 :                          fm_matrix_Minv_Vtrunc_Minv(1, 1), 0.0_dp, fm_work)
    1591             : 
    1592             :       CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work, &
    1593          22 :                          fm_matrix_RI_metric_inv_work, 0.0_dp, fm_matrix_Minv_Vtrunc_Minv(1, 1))
    1594             : 
    1595          22 :       CALL cp_fm_release(fm_work)
    1596          22 :       CALL cp_fm_release(fm_matrix_RI_metric_inv_work)
    1597             : 
    1598          22 :       CALL timestop(handle)
    1599             : 
    1600          22 :    END SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc
    1601             : 
    1602             : ! **************************************************************************************************
    1603             : !> \brief ...
    1604             : !> \param qs_env ...
    1605             : !> \param trunc_coulomb ...
    1606             : !> \param rel_cutoff_trunc_coulomb_ri_x ...
    1607             : !> \param cell_grid ...
    1608             : !> \param do_BvK_cell ...
    1609             : ! **************************************************************************************************
    1610          58 :    SUBROUTINE trunc_coulomb_for_exchange(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x, &
    1611             :                                          cell_grid, do_BvK_cell)
    1612             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1613             :       TYPE(libint_potential_type), OPTIONAL              :: trunc_coulomb
    1614             :       REAL(KIND=dp), OPTIONAL                            :: rel_cutoff_trunc_coulomb_ri_x
    1615             :       INTEGER, DIMENSION(3), OPTIONAL                    :: cell_grid
    1616             :       LOGICAL, OPTIONAL                                  :: do_BvK_cell
    1617             : 
    1618             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'trunc_coulomb_for_exchange'
    1619             : 
    1620             :       INTEGER                                            :: handle, i_dim
    1621             :       INTEGER, DIMENSION(3)                              :: periodic
    1622             :       LOGICAL                                            :: my_do_BvK_cell
    1623             :       REAL(KIND=dp) :: kp_fac, kp_fac_idim, my_rel_cutoff_trunc_coulomb_ri_x, &
    1624             :          shortest_dist_cell_planes
    1625             :       TYPE(cell_type), POINTER                           :: cell
    1626             :       TYPE(kpoint_type), POINTER                         :: kpoints_scf
    1627             : 
    1628          58 :       CALL timeset(routineN, handle)
    1629             : 
    1630          58 :       NULLIFY (cell)
    1631          58 :       CALL get_qs_env(qs_env, cell=cell, kpoints=kpoints_scf)
    1632          58 :       CALL get_cell(cell=cell, periodic=periodic)
    1633             : 
    1634          58 :       my_do_BvK_cell = .FALSE.
    1635          58 :       IF (PRESENT(do_BvK_cell)) my_do_BvK_cell = do_BvK_cell
    1636          36 :       IF (my_do_BvK_cell) THEN
    1637             :          kp_fac = 1.0E10_dp
    1638          32 :          DO i_dim = 1, 3
    1639             :             ! look for smallest k-point mesh in periodic direction
    1640          32 :             IF (periodic(i_dim) == 1) THEN
    1641          16 :                IF (PRESENT(cell_grid)) THEN
    1642          16 :                   kp_fac_idim = REAL(cell_grid(i_dim), KIND=dp)
    1643             :                ELSE
    1644           0 :                   kp_fac_idim = REAL(kpoints_scf%nkp_grid(i_dim), KIND=dp)
    1645             :                END IF
    1646          16 :                IF (kp_fac > kp_fac_idim) kp_fac = kp_fac_idim
    1647             :             END IF
    1648             :          END DO
    1649             :       ELSE
    1650             :          kp_fac = 1.0_dp
    1651             :       END IF
    1652             : 
    1653          58 :       shortest_dist_cell_planes = 1.0E4_dp
    1654          58 :       IF (periodic(1) == 1) THEN
    1655          16 :          IF (shortest_dist_cell_planes > plane_distance(1, 0, 0, cell)) THEN
    1656          16 :             shortest_dist_cell_planes = plane_distance(1, 0, 0, cell)
    1657             :          END IF
    1658             :       END IF
    1659          58 :       IF (periodic(2) == 1) THEN
    1660          38 :          IF (shortest_dist_cell_planes > plane_distance(0, 1, 0, cell)) THEN
    1661          26 :             shortest_dist_cell_planes = plane_distance(0, 1, 0, cell)
    1662             :          END IF
    1663             :       END IF
    1664          58 :       IF (periodic(3) == 1) THEN
    1665          34 :          IF (shortest_dist_cell_planes > plane_distance(0, 0, 1, cell)) THEN
    1666           8 :             shortest_dist_cell_planes = plane_distance(0, 0, 1, cell)
    1667             :          END IF
    1668             :       END IF
    1669             : 
    1670          58 :       IF (PRESENT(rel_cutoff_trunc_coulomb_ri_x)) THEN
    1671          36 :          my_rel_cutoff_trunc_coulomb_ri_x = rel_cutoff_trunc_coulomb_ri_x
    1672             :       ELSE
    1673          22 :          my_rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
    1674             :       END IF
    1675             : 
    1676          58 :       IF (PRESENT(trunc_coulomb)) THEN
    1677          58 :          trunc_coulomb%potential_type = do_potential_truncated
    1678             :          trunc_coulomb%cutoff_radius = shortest_dist_cell_planes* &
    1679             :                                        my_rel_cutoff_trunc_coulomb_ri_x* &
    1680          58 :                                        kp_fac
    1681          58 :          trunc_coulomb%filename = "t_c_g.dat"
    1682             :          ! dummy
    1683          58 :          trunc_coulomb%omega = 0.0_dp
    1684             :       END IF
    1685             : 
    1686          58 :       CALL timestop(handle)
    1687             : 
    1688          58 :    END SUBROUTINE trunc_coulomb_for_exchange
    1689             : 
    1690             : ! **************************************************************************************************
    1691             : !> \brief ...
    1692             : !> \param fm_matrix_L ...
    1693             : ! **************************************************************************************************
    1694        1316 :    SUBROUTINE invert_mat(fm_matrix_L)
    1695             : 
    1696             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_matrix_L
    1697             : 
    1698             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'invert_mat'
    1699             : 
    1700             :       INTEGER                                            :: handle, i_global, iiB, j_global, jjB, &
    1701             :                                                             ncol_local, nrow_local
    1702         658 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1703             : 
    1704         658 :       CALL timeset(routineN, handle)
    1705             : 
    1706         658 :       CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U')
    1707             : 
    1708             :       ! clear the lower part of the L^{-1} matrix (just to have no surprises afterwards)
    1709             :       CALL cp_fm_get_info(matrix=fm_matrix_L, &
    1710             :                           nrow_local=nrow_local, &
    1711             :                           ncol_local=ncol_local, &
    1712             :                           row_indices=row_indices, &
    1713         658 :                           col_indices=col_indices)
    1714             :       ! assume upper triangular format
    1715       50832 :       DO jjB = 1, ncol_local
    1716       50174 :          j_global = col_indices(jjB)
    1717     3639757 :          DO iiB = 1, nrow_local
    1718     3588925 :             i_global = row_indices(iiB)
    1719     3639099 :             IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
    1720             :          END DO
    1721             :       END DO
    1722             : 
    1723         658 :       CALL timestop(handle)
    1724             : 
    1725         658 :    END SUBROUTINE invert_mat
    1726             : 
    1727             : ! **************************************************************************************************
    1728             : !> \brief ...
    1729             : !> \param fm_matrix_L ...
    1730             : !> \param blacs_env_L ...
    1731             : !> \param dimen_RI ...
    1732             : !> \param para_env_L ...
    1733             : !> \param name ...
    1734             : !> \param fm_matrix_L_extern ...
    1735             : ! **************************************************************************************************
    1736         694 :    SUBROUTINE create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, fm_matrix_L_extern)
    1737             :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_matrix_L
    1738             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_L
    1739             :       INTEGER, INTENT(IN)                                :: dimen_RI
    1740             :       TYPE(mp_para_env_type), POINTER                    :: para_env_L
    1741             :       CHARACTER(LEN=*), INTENT(IN)                       :: name
    1742             :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: fm_matrix_L_extern
    1743             : 
    1744             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'create_matrix_L'
    1745             : 
    1746             :       INTEGER                                            :: handle
    1747             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    1748             : 
    1749         694 :       CALL timeset(routineN, handle)
    1750             : 
    1751             :       ! create the full matrix L defined in the L group
    1752         694 :       IF (PRESENT(fm_matrix_L_extern)) THEN
    1753          36 :          CALL cp_fm_create(fm_matrix_L, fm_matrix_L_extern%matrix_struct, name=name)
    1754             :       ELSE
    1755         658 :          NULLIFY (fm_struct)
    1756             :          CALL cp_fm_struct_create(fm_struct, context=blacs_env_L, nrow_global=dimen_RI, &
    1757         658 :                                   ncol_global=dimen_RI, para_env=para_env_L)
    1758             : 
    1759         658 :          CALL cp_fm_create(fm_matrix_L, fm_struct, name=name)
    1760             : 
    1761         658 :          CALL cp_fm_struct_release(fm_struct)
    1762             :       END IF
    1763             : 
    1764         694 :       CALL cp_fm_set_all(matrix=fm_matrix_L, alpha=0.0_dp)
    1765             : 
    1766         694 :       CALL timestop(handle)
    1767             : 
    1768         694 :    END SUBROUTINE create_matrix_L
    1769             : 
    1770             : ! **************************************************************************************************
    1771             : !> \brief ...
    1772             : !> \param fm_matrix_L ...
    1773             : !> \param my_group_L_start ...
    1774             : !> \param my_group_L_end ...
    1775             : !> \param my_group_L_size ...
    1776             : !> \param my_Lrows ...
    1777             : ! **************************************************************************************************
    1778         534 :    SUBROUTINE grep_Lcols(fm_matrix_L, &
    1779             :                          my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
    1780             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_matrix_L
    1781             :       INTEGER, INTENT(IN)                                :: my_group_L_start, my_group_L_end, &
    1782             :                                                             my_group_L_size
    1783             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
    1784             :          INTENT(OUT)                                     :: my_Lrows
    1785             : 
    1786             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'grep_Lcols'
    1787             : 
    1788             :       INTEGER :: dimen_RI, handle, handle2, i_global, iiB, j_global, jjB, max_row_col_local, &
    1789             :          ncol_local, ncol_rec, nrow_local, nrow_rec, proc_receive_static, proc_send_static, &
    1790             :          proc_shift
    1791         534 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: local_col_row_info, rec_col_row_info
    1792         534 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, col_indices_rec, &
    1793         534 :                                                             row_indices, row_indices_rec
    1794         534 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: local_L, rec_L
    1795             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
    1796         534 :          POINTER                                         :: local_L_internal
    1797             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1798             : 
    1799         534 :       CALL timeset(routineN, handle)
    1800             : 
    1801             :       CALL cp_fm_get_info(matrix=fm_matrix_L, &
    1802             :                           nrow_local=nrow_local, &
    1803             :                           ncol_local=ncol_local, &
    1804             :                           row_indices=row_indices, &
    1805             :                           col_indices=col_indices, &
    1806             :                           nrow_global=dimen_RI, &
    1807             :                           local_data=local_L_internal, &
    1808         534 :                           para_env=para_env)
    1809             : 
    1810        2136 :       ALLOCATE (my_Lrows(dimen_RI, my_group_L_size))
    1811     1720484 :       my_Lrows = 0.0_dp
    1812             : 
    1813        2136 :       ALLOCATE (local_L(nrow_local, ncol_local))
    1814     3171430 :       local_L(:, :) = local_L_internal(1:nrow_local, 1:ncol_local)
    1815             : 
    1816         534 :       max_row_col_local = MAX(nrow_local, ncol_local)
    1817         534 :       CALL para_env%max(max_row_col_local)
    1818             : 
    1819        2136 :       ALLOCATE (local_col_row_info(0:max_row_col_local, 2))
    1820       82118 :       local_col_row_info = 0
    1821             :       ! 0,1 nrows
    1822         534 :       local_col_row_info(0, 1) = nrow_local
    1823       39255 :       local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
    1824             :       ! 0,2 ncols
    1825         534 :       local_col_row_info(0, 2) = ncol_local
    1826       40148 :       local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
    1827             : 
    1828        1068 :       ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
    1829             : 
    1830             :       ! accumulate data on my_Lrows starting from myself
    1831       40148 :       DO jjB = 1, ncol_local
    1832       39614 :          j_global = col_indices(jjB)
    1833       40148 :          IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
    1834     1630962 :             DO iiB = 1, nrow_local
    1835     1610135 :                i_global = row_indices(iiB)
    1836     1630962 :                my_Lrows(i_global, j_global - my_group_L_start + 1) = local_L(iiB, jjB)
    1837             :             END DO
    1838             :          END IF
    1839             :       END DO
    1840             : 
    1841         534 :       proc_send_static = MODULO(para_env%mepos + 1, para_env%num_pe)
    1842         534 :       proc_receive_static = MODULO(para_env%mepos - 1, para_env%num_pe)
    1843             : 
    1844         534 :       CALL timeset(routineN//"_comm", handle2)
    1845             : 
    1846         558 :       DO proc_shift = 1, para_env%num_pe - 1
    1847             :          ! first exchange information on the local data
    1848        4200 :          rec_col_row_info = 0
    1849          24 :          CALL para_env%sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static)
    1850          24 :          nrow_rec = rec_col_row_info(0, 1)
    1851          24 :          ncol_rec = rec_col_row_info(0, 2)
    1852             : 
    1853          72 :          ALLOCATE (row_indices_rec(nrow_rec))
    1854        1061 :          row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
    1855             : 
    1856          72 :          ALLOCATE (col_indices_rec(ncol_rec))
    1857        2064 :          col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
    1858             : 
    1859          96 :          ALLOCATE (rec_L(nrow_rec, ncol_rec))
    1860       91052 :          rec_L = 0.0_dp
    1861             : 
    1862             :          ! then send and receive the real data
    1863          24 :          CALL para_env%sendrecv(local_L, proc_send_static, rec_L, proc_receive_static)
    1864             : 
    1865             :          ! accumulate the received data on my_Lrows
    1866        2064 :          DO jjB = 1, ncol_rec
    1867        2040 :             j_global = col_indices_rec(jjB)
    1868        2064 :             IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
    1869       91028 :                DO iiB = 1, nrow_rec
    1870       88988 :                   i_global = row_indices_rec(iiB)
    1871       91028 :                   my_Lrows(i_global, j_global - my_group_L_start + 1) = rec_L(iiB, jjB)
    1872             :                END DO
    1873             :             END IF
    1874             :          END DO
    1875             : 
    1876        4200 :          local_col_row_info(:, :) = rec_col_row_info
    1877          24 :          CALL MOVE_ALLOC(rec_L, local_L)
    1878             : 
    1879          24 :          DEALLOCATE (col_indices_rec)
    1880         558 :          DEALLOCATE (row_indices_rec)
    1881             :       END DO
    1882         534 :       CALL timestop(handle2)
    1883             : 
    1884         534 :       DEALLOCATE (local_col_row_info)
    1885         534 :       DEALLOCATE (rec_col_row_info)
    1886         534 :       DEALLOCATE (local_L)
    1887             : 
    1888         534 :       CALL timestop(handle)
    1889             : 
    1890        1602 :    END SUBROUTINE grep_Lcols
    1891             : 
    1892             : END MODULE mp2_ri_2c

Generated by: LCOV version 1.15