rpa_util.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Utility functions for RPA calculations
!> \par History
!>      06.2019 Moved from rpa_ri_gpw.F [Frederick Stein]
! **************************************************************************************************
MODULE rpa_util
 
   USE cell_types,                      ONLY: cell_type,&
                                              get_cell
   USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
                                              cp_blacs_env_release,&
                                              cp_blacs_env_type
   USE cp_cfm_types,                    ONLY: cp_cfm_create,&
                                              cp_cfm_release,&
                                              cp_cfm_set_all,&
                                              cp_cfm_type
   USE cp_dbcsr_api,                    ONLY: &
        dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, &
        dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
                                              copy_fm_to_dbcsr,&
                                              dbcsr_allocate_matrix_set,&
                                              dbcsr_deallocate_matrix_set
   USE cp_fm_basic_linalg,              ONLY: cp_fm_syrk,&
                                              cp_fm_transpose
   USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose
   USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
                                              cp_fm_struct_release,&
                                              cp_fm_struct_type
   USE cp_fm_types,                     ONLY: &
        cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_release, cp_fm_set_all, &
        cp_fm_set_element, cp_fm_to_fm, cp_fm_to_fm_submat_general, cp_fm_type
   USE dbt_api,                         ONLY: dbt_destroy,&
                                              dbt_type
   USE dgemm_counter_types,             ONLY: dgemm_counter_start,&
                                              dgemm_counter_stop,&
                                              dgemm_counter_type
   USE hfx_types,                       ONLY: block_ind_type,&
                                              dealloc_containers,&
                                              hfx_compression_type
   USE input_constants,                 ONLY: wfc_mm_style_gemm,&
                                              wfc_mm_style_syrk
   USE kinds,                           ONLY: dp
   USE kpoint_types,                    ONLY: get_kpoint_info,&
                                              kpoint_release,&
                                              kpoint_type
   USE mathconstants,                   ONLY: z_zero
   USE message_passing,                 ONLY: mp_para_env_release,&
                                              mp_para_env_type
   USE mp2_laplace,                     ONLY: calc_fm_mat_s_laplace
   USE parallel_gemm_api,               ONLY: parallel_gemm
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE rpa_gw_kpoints_util,             ONLY: compute_wkp_w
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util'
 
   PUBLIC :: compute_erpa_by_freq_int, alloc_im_time, calc_mat_q, q_trace_and_add_unit_matrix, &
             dealloc_im_time, contract_p_omega_with_mat_l, calc_fm_mat_s_rpa, remove_scaling_factor_rpa
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param para_env ...
!> \param dimen_RI ...
!> \param dimen_RI_red ...
!> \param num_integ_points ...
!> \param nspins ...
!> \param fm_mat_Q ...
!> \param fm_mo_coeff_occ ...
!> \param fm_mo_coeff_virt ...
!> \param fm_matrix_Minv_L_kpoints ...
!> \param fm_matrix_L_kpoints ...
!> \param mat_P_global ...
!> \param t_3c_O ...
!> \param matrix_s ...
!> \param kpoints ...
!> \param eps_filter_im_time ...
!> \param cut_memory ...
!> \param nkp ...
!> \param num_cells_dm ...
!> \param num_3c_repl ...
!> \param size_P ...
!> \param ikp_local ...
!> \param index_to_cell_3c ...
!> \param cell_to_index_3c ...
!> \param col_blk_size ...
!> \param do_ic_model ...
!> \param do_kpoints_cubic_RPA ...
!> \param do_kpoints_from_Gamma ...
!> \param do_ri_Sigma_x ...
!> \param my_open_shell ...
!> \param has_mat_P_blocks ...
!> \param wkp_W ...
!> \param cfm_mat_Q ...
!> \param fm_mat_Minv_L_kpoints ...
!> \param fm_mat_L_kpoints ...
!> \param fm_mat_RI_global_work ...
!> \param fm_mat_work ...
!> \param fm_mo_coeff_occ_scaled ...
!> \param fm_mo_coeff_virt_scaled ...
!> \param mat_dm ...
!> \param mat_L ...
!> \param mat_M_P_munu_occ ...
!> \param mat_M_P_munu_virt ...
!> \param mat_MinvVMinv ...
!> \param mat_P_omega ...
!> \param mat_P_omega_kp ...
!> \param mat_work ...
!> \param mo_coeff ...
!> \param fm_scaled_dm_occ_tau ...
!> \param fm_scaled_dm_virt_tau ...
!> \param homo ...
!> \param nmo ...
! **************************************************************************************************
   SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_points, nspins, &
                            fm_mat_Q, fm_mo_coeff_occ, fm_mo_coeff_virt, &
                            fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, mat_P_global, &
                            t_3c_O, matrix_s, kpoints, eps_filter_im_time, &
                            cut_memory, nkp, num_cells_dm, num_3c_repl, &
                            size_P, ikp_local, &
                            index_to_cell_3c, &
                            cell_to_index_3c, &
                            col_blk_size, &
                            do_ic_model, do_kpoints_cubic_RPA, &
                            do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, &
                            has_mat_P_blocks, wkp_W, &
                            cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
                            fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, &
                            fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, &
                            mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
                            mat_work, mo_coeff, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
 
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(mp_para_env_type), POINTER                    :: para_env
      INTEGER, INTENT(IN)                                :: dimen_ri, dimen_ri_red, &
                                                            num_integ_points, nspins
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_q
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_mo_coeff_occ, fm_mo_coeff_virt
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_matrix_minv_l_kpoints, &
                                                            fm_matrix_l_kpoints
      TYPE(dbcsr_p_type), INTENT(IN)                     :: mat_p_global
      TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
         INTENT(INOUT)                                   :: t_3c_o
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      TYPE(kpoint_type), POINTER                         :: kpoints
      REAL(kind=dp), INTENT(IN)                          :: eps_filter_im_time
      INTEGER, INTENT(IN)                                :: cut_memory
      INTEGER, INTENT(OUT)                               :: nkp, num_cells_dm, num_3c_repl, size_p, &
                                                            ikp_local
      INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell_3c
      INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
         INTENT(OUT)                                     :: cell_to_index_3c
      INTEGER, DIMENSION(:), POINTER                     :: col_blk_size
      LOGICAL, INTENT(IN)                                :: do_ic_model, do_kpoints_cubic_rpa, &
                                                            do_kpoints_from_gamma, do_ri_sigma_x, &
                                                            my_open_shell
      LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
         INTENT(OUT)                                     :: has_mat_p_blocks
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
         INTENT(OUT)                                     :: wkp_w
      TYPE(cp_cfm_type), INTENT(OUT)                     :: cfm_mat_q
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_mat_minv_l_kpoints, fm_mat_l_kpoints
      TYPE(cp_fm_type), INTENT(OUT)                      :: fm_mat_ri_global_work, fm_mat_work, &
                                                            fm_mo_coeff_occ_scaled, &
                                                            fm_mo_coeff_virt_scaled
      TYPE(dbcsr_p_type), INTENT(OUT)                    :: mat_dm, mat_l, mat_m_p_munu_occ, &
                                                            mat_m_p_munu_virt, mat_minvvminv
      TYPE(dbcsr_p_type), ALLOCATABLE, &
         DIMENSION(:, :, :)                              :: mat_p_omega
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_p_omega_kp
      TYPE(dbcsr_type), POINTER                          :: mat_work
      TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff
      TYPE(cp_fm_type), INTENT(OUT)                      :: fm_scaled_dm_occ_tau, &
                                                            fm_scaled_dm_virt_tau
      INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
      INTEGER, INTENT(IN)                                :: nmo
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'alloc_im_time'
 
      INTEGER                                            :: cell_grid_dm(3), first_ikp_local, &
                                                            handle, i_dim, i_kp, ispin, jquad, &
                                                            nspins_p_omega, periodic(3)
      INTEGER, DIMENSION(:), POINTER                     :: row_blk_size
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: wkp_v
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_sub_kp
 
      CALL timeset(routinen, handle)
 
      ALLOCATE (fm_mo_coeff_occ(nspins), fm_mo_coeff_virt(nspins))
 
      CALL cp_fm_create(fm_scaled_dm_occ_tau, mo_coeff(1)%matrix_struct)
      CALL cp_fm_set_all(fm_scaled_dm_occ_tau, 0.0_dp)
 
      CALL cp_fm_create(fm_scaled_dm_virt_tau, mo_coeff(1)%matrix_struct)
      CALL cp_fm_set_all(fm_scaled_dm_virt_tau, 0.0_dp)
 
      DO ispin = 1, SIZE(mo_coeff)
         CALL create_occ_virt_mo_coeffs(fm_mo_coeff_occ(ispin), fm_mo_coeff_virt(ispin), mo_coeff(ispin), &
                                        nmo, homo(ispin))
      END DO
 
      num_3c_repl = SIZE(t_3c_o, 2)
 
      IF (do_kpoints_cubic_rpa) THEN
         ! we always use an odd number of image cells
         ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
         DO i_dim = 1, 3
            cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
         END DO
         num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
         ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2)))
         cpassert(SIZE(kpoints%index_to_cell, 1) == 3)
         index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :)
         ALLOCATE (cell_to_index_3c(lbound(kpoints%cell_to_index, 1):ubound(kpoints%cell_to_index, 1), &
                                    lbound(kpoints%cell_to_index, 2):ubound(kpoints%cell_to_index, 2), &
                                    lbound(kpoints%cell_to_index, 3):ubound(kpoints%cell_to_index, 3)))
         cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :)
 
      ELSE
         ALLOCATE (index_to_cell_3c(3, 1))
         index_to_cell_3c(:, 1) = 0
         ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0))
         cell_to_index_3c(0, 0, 0) = 1
         num_cells_dm = 1
      END IF
 
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
 
         CALL get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_ri, ikp_local, first_ikp_local)
 
         CALL cp_cfm_create(cfm_mat_q, fm_struct_sub_kp)
         CALL cp_cfm_set_all(cfm_mat_q, z_zero)
      ELSE
         first_ikp_local = 1
      END IF
 
      ! if we do kpoints, mat_P has a kpoint and mat_P_omega has the inted
      ! mat_P(tau, kpoint)
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
 
         NULLIFY (cell)
         CALL get_qs_env(qs_env, cell=cell)
         CALL get_cell(cell=cell, periodic=periodic)
 
         CALL get_kpoint_info(kpoints, nkp=nkp)
         ! compute k-point weights such that functions 1/k^2, 1/k and const function are
         ! integrated correctly
         CALL compute_wkp_w(qs_env, wkp_w, wkp_v, kpoints, cell%h_inv, periodic)
         DEALLOCATE (wkp_v)
 
      ELSE
         nkp = 1
      END IF
 
      IF (do_kpoints_cubic_rpa) THEN
         size_p = max(num_cells_dm/2 + 1, nkp)
      ELSE IF (do_kpoints_from_gamma) THEN
         size_p = max(3**(periodic(1) + periodic(2) + periodic(3)), nkp)
      ELSE
         size_p = 1
      END IF
 
      nspins_p_omega = 1
      IF (my_open_shell) nspins_p_omega = 2
 
      ALLOCATE (mat_p_omega(num_integ_points, size_p, nspins_p_omega))
      DO ispin = 1, nspins_p_omega
         DO i_kp = 1, size_p
            DO jquad = 1, num_integ_points
               NULLIFY (mat_p_omega(jquad, i_kp, ispin)%matrix)
               ALLOCATE (mat_p_omega(jquad, i_kp, ispin)%matrix)
               CALL dbcsr_create(matrix=mat_p_omega(jquad, i_kp, ispin)%matrix, &
                                 template=mat_p_global%matrix)
               CALL dbcsr_set(mat_p_omega(jquad, i_kp, ispin)%matrix, 0.0_dp)
            END DO
         END DO
      END DO
 
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
         CALL alloc_mat_p_omega(mat_p_omega_kp, 2, size_p, mat_p_global%matrix)
      END IF
 
      CALL cp_fm_create(fm_mo_coeff_occ_scaled, fm_mo_coeff_occ(1)%matrix_struct)
      CALL cp_fm_to_fm(fm_mo_coeff_occ(1), fm_mo_coeff_occ_scaled)
      CALL cp_fm_set_all(matrix=fm_mo_coeff_occ_scaled, alpha=0.0_dp)
 
      CALL cp_fm_create(fm_mo_coeff_virt_scaled, fm_mo_coeff_virt(1)%matrix_struct)
      CALL cp_fm_to_fm(fm_mo_coeff_virt(1), fm_mo_coeff_virt_scaled)
      CALL cp_fm_set_all(matrix=fm_mo_coeff_virt_scaled, alpha=0.0_dp)
 
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
         CALL cp_fm_create(fm_mat_ri_global_work, fm_matrix_minv_l_kpoints(1, 1)%matrix_struct)
         CALL cp_fm_to_fm(fm_matrix_minv_l_kpoints(1, 1), fm_mat_ri_global_work)
         CALL cp_fm_set_all(fm_mat_ri_global_work, 0.0_dp)
      END IF
 
      ALLOCATE (has_mat_p_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl))
      has_mat_p_blocks = .true.
 
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
         CALL reorder_mat_l(fm_mat_minv_l_kpoints, fm_matrix_minv_l_kpoints, fm_mat_q%matrix_struct, para_env, mat_l, &
                            mat_p_global%matrix, dimen_ri, dimen_ri_red, first_ikp_local, ikp_local, fm_struct_sub_kp, &
                            allocate_mat_l=.false.)
 
         CALL reorder_mat_l(fm_mat_l_kpoints, fm_matrix_l_kpoints, fm_mat_q%matrix_struct, para_env, mat_l, &
                            mat_p_global%matrix, dimen_ri, dimen_ri_red, first_ikp_local, ikp_local, fm_struct_sub_kp)
 
         CALL cp_fm_struct_release(fm_struct_sub_kp)
 
      ELSE
         CALL reorder_mat_l(fm_mat_minv_l_kpoints, fm_matrix_minv_l_kpoints, fm_mat_q%matrix_struct, para_env, mat_l, &
                            mat_p_global%matrix, dimen_ri, dimen_ri_red, first_ikp_local)
      END IF
 
      ! Create Scalapack working matrix for the contraction with the metric
      IF (dimen_ri == dimen_ri_red) THEN
         CALL cp_fm_create(fm_mat_work, fm_mat_q%matrix_struct)
         CALL cp_fm_to_fm(fm_mat_q, fm_mat_work)
         CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
 
      ELSE
         NULLIFY (fm_struct)
         CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_mat_q%matrix_struct, &
                                  ncol_global=dimen_ri_red, nrow_global=dimen_ri)
 
         CALL cp_fm_create(fm_mat_work, fm_struct)
         CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
 
         CALL cp_fm_struct_release(fm_struct)
 
      END IF
 
      ! Then its DBCSR counter part
      IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
         CALL dbcsr_get_info(mat_l%matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
 
         ! Create mat_work having the shape of the transposed of mat_L (compare with contract_P_omega_with_mat_L)
         NULLIFY (mat_work)
         ALLOCATE (mat_work)
         CALL dbcsr_create(mat_work, template=mat_l%matrix, row_blk_size=col_blk_size, col_blk_size=row_blk_size)
      END IF
 
      IF (do_ri_sigma_x .OR. do_ic_model) THEN
 
         NULLIFY (mat_minvvminv%matrix)
         ALLOCATE (mat_minvvminv%matrix)
         CALL dbcsr_create(mat_minvvminv%matrix, template=mat_p_global%matrix)
         CALL dbcsr_set(mat_minvvminv%matrix, 0.0_dp)
 
         ! for kpoints we compute SinvVSinv later with kpoints
         IF (.NOT. do_kpoints_from_gamma) THEN
 
            !  get the Coulomb matrix for Sigma_x = G*V
            CALL dbcsr_multiply("T", "N", 1.0_dp, mat_l%matrix, mat_l%matrix, &
                                0.0_dp, mat_minvvminv%matrix, filter_eps=eps_filter_im_time)
 
         END IF
 
      END IF
 
      IF (do_ri_sigma_x) THEN
 
         NULLIFY (mat_dm%matrix)
         ALLOCATE (mat_dm%matrix)
         CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix)
 
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE alloc_im_time
 
! **************************************************************************************************
!> \brief ...
!> \param fm_mo_coeff_occ ...
!> \param fm_mo_coeff_virt ...
!> \param mo_coeff ...
!> \param nmo ...
!> \param homo ...
! **************************************************************************************************
   SUBROUTINE create_occ_virt_mo_coeffs(fm_mo_coeff_occ, fm_mo_coeff_virt, mo_coeff, &
                                        nmo, homo)
 
      TYPE(cp_fm_type), INTENT(OUT)                      :: fm_mo_coeff_occ, fm_mo_coeff_virt
      TYPE(cp_fm_type), INTENT(IN)                       :: mo_coeff
      INTEGER, INTENT(IN)                                :: nmo, homo
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'create_occ_virt_mo_coeffs'
 
      INTEGER                                            :: handle, icol_global, irow_global
 
      CALL timeset(routinen, handle)
 
      CALL cp_fm_create(fm_mo_coeff_occ, mo_coeff%matrix_struct)
      CALL cp_fm_set_all(fm_mo_coeff_occ, 0.0_dp)
      CALL cp_fm_to_fm(mo_coeff, fm_mo_coeff_occ)
 
      ! set all virtual MO coeffs to zero
      DO irow_global = 1, nmo
         DO icol_global = homo + 1, nmo
            CALL cp_fm_set_element(fm_mo_coeff_occ, irow_global, icol_global, 0.0_dp)
         END DO
      END DO
 
      CALL cp_fm_create(fm_mo_coeff_virt, mo_coeff%matrix_struct)
      CALL cp_fm_set_all(fm_mo_coeff_virt, 0.0_dp)
      CALL cp_fm_to_fm(mo_coeff, fm_mo_coeff_virt)
 
      ! set all occupied MO coeffs to zero
      DO irow_global = 1, nmo
         DO icol_global = 1, homo
            CALL cp_fm_set_element(fm_mo_coeff_virt, irow_global, icol_global, 0.0_dp)
         END DO
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE create_occ_virt_mo_coeffs
 
! **************************************************************************************************
!> \brief ...
!> \param mat_P_omega ...
!> \param num_integ_points ...
!> \param size_P ...
!> \param template ...
! **************************************************************************************************
   SUBROUTINE alloc_mat_p_omega(mat_P_omega, num_integ_points, size_P, template)
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_p_omega
      INTEGER, INTENT(IN)                                :: num_integ_points, size_p
      TYPE(dbcsr_type), POINTER                          :: template
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'alloc_mat_P_omega'
 
      INTEGER                                            :: handle, i_kp, jquad
 
      CALL timeset(routinen, handle)
 
      NULLIFY (mat_p_omega)
      CALL dbcsr_allocate_matrix_set(mat_p_omega, num_integ_points, size_p)
      DO i_kp = 1, size_p
         DO jquad = 1, num_integ_points
            ALLOCATE (mat_p_omega(jquad, i_kp)%matrix)
            CALL dbcsr_create(matrix=mat_p_omega(jquad, i_kp)%matrix, &
                              template=template)
            CALL dbcsr_set(mat_p_omega(jquad, i_kp)%matrix, 0.0_dp)
         END DO
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE alloc_mat_p_omega
 
! **************************************************************************************************
!> \brief ...
!> \param fm_mat_L ...
!> \param fm_matrix_Minv_L_kpoints ...
!> \param fm_struct_template ...
!> \param para_env ...
!> \param mat_L ...
!> \param mat_template ...
!> \param dimen_RI ...
!> \param dimen_RI_red ...
!> \param first_ikp_local ...
!> \param ikp_local ...
!> \param fm_struct_sub_kp ...
!> \param allocate_mat_L ...
! **************************************************************************************************
   SUBROUTINE reorder_mat_l(fm_mat_L, fm_matrix_Minv_L_kpoints, fm_struct_template, para_env, mat_L, mat_template, &
                            dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp, allocate_mat_L)
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_mat_l, fm_matrix_minv_l_kpoints
      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_template
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(dbcsr_p_type), INTENT(OUT)                    :: mat_l
      TYPE(dbcsr_type), INTENT(IN)                       :: mat_template
      INTEGER, INTENT(IN)                                :: dimen_ri, dimen_ri_red, first_ikp_local
      INTEGER, OPTIONAL                                  :: ikp_local
      TYPE(cp_fm_struct_type), OPTIONAL, POINTER         :: fm_struct_sub_kp
      LOGICAL, INTENT(IN), OPTIONAL                      :: allocate_mat_l
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'reorder_mat_L'
 
      INTEGER                                            :: handle, ikp, j_size, nblk
      INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, row_blk_size
      LOGICAL                                            :: do_kpoints, my_allocate_mat_l
      TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
      TYPE(cp_fm_type)                                   :: fm_mat_l_transposed, fmdummy
 
      CALL timeset(routinen, handle)
 
      do_kpoints = .false.
      IF (PRESENT(ikp_local) .AND. PRESENT(fm_struct_sub_kp)) THEN
         do_kpoints = .true.
      END IF
 
      ! Get the fm_struct for fm_mat_L
      NULLIFY (fm_struct)
      IF (dimen_ri == dimen_ri_red) THEN
         fm_struct => fm_struct_template
      ELSE
         ! The template is assumed to be square such that we need a new fm_struct if dimensions are not equal
         CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_ri_red, ncol_global=dimen_ri, template_fmstruct=fm_struct_template)
      END IF
 
      ! Start to allocate the new full matrix
      ALLOCATE (fm_mat_l(SIZE(fm_matrix_minv_l_kpoints, 1), SIZE(fm_matrix_minv_l_kpoints, 2)))
      DO ikp = 1, SIZE(fm_matrix_minv_l_kpoints, 1)
         DO j_size = 1, SIZE(fm_matrix_minv_l_kpoints, 2)
            IF (do_kpoints) THEN
               IF (ikp == first_ikp_local .OR. ikp_local == -1) THEN
                  CALL cp_fm_create(fm_mat_l(ikp, j_size), fm_struct_sub_kp)
                  CALL cp_fm_set_all(fm_mat_l(ikp, j_size), 0.0_dp)
               END IF
            ELSE
               CALL cp_fm_create(fm_mat_l(ikp, j_size), fm_struct)
               CALL cp_fm_set_all(fm_mat_l(ikp, j_size), 0.0_dp)
            END IF
         END DO
      END DO
 
      ! For the transposed matric we need a different fm_struct
      IF (dimen_ri == dimen_ri_red) THEN
         fm_struct => fm_mat_l(first_ikp_local, 1)%matrix_struct
      ELSE
         CALL cp_fm_struct_release(fm_struct)
 
         ! Create a fm_struct with transposed sizes
         NULLIFY (fm_struct)
         CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_ri, ncol_global=dimen_ri_red, &
                                  template_fmstruct=fm_mat_l(first_ikp_local, 1)%matrix_struct) !, force_block=.TRUE.)
      END IF
 
      ! Allocate buffer matrix
      CALL cp_fm_create(fm_mat_l_transposed, fm_struct)
      CALL cp_fm_set_all(matrix=fm_mat_l_transposed, alpha=0.0_dp)
 
      IF (dimen_ri /= dimen_ri_red) CALL cp_fm_struct_release(fm_struct)
 
      CALL cp_fm_get_info(fm_mat_l_transposed, context=blacs_env)
 
      ! For k-points copy matrices of your group
      ! Without kpoints, transpose matrix
      ! without kpoints, the size of fm_mat_L is 1x1. with kpoints, the size is N_kpoints x 2 (2 for real/complex)
      DO ikp = 1, SIZE(fm_matrix_minv_l_kpoints, 1)
      DO j_size = 1, SIZE(fm_matrix_minv_l_kpoints, 2)
         IF (do_kpoints) THEN
            IF (ikp_local == ikp .OR. ikp_local == -1) THEN
               CALL cp_fm_copy_general(fm_matrix_minv_l_kpoints(ikp, j_size), fm_mat_l_transposed, para_env)
               CALL cp_fm_to_fm(fm_mat_l_transposed, fm_mat_l(ikp, j_size))
            ELSE
               CALL cp_fm_copy_general(fm_matrix_minv_l_kpoints(ikp, j_size), fmdummy, para_env)
            END IF
         ELSE
            CALL cp_fm_copy_general(fm_matrix_minv_l_kpoints(ikp, j_size), fm_mat_l_transposed, blacs_env%para_env)
            CALL cp_fm_transpose(fm_mat_l_transposed, fm_mat_l(ikp, j_size))
         END IF
      END DO
      END DO
 
      ! Release old matrix
      CALL cp_fm_release(fm_matrix_minv_l_kpoints)
      ! Release buffer
      CALL cp_fm_release(fm_mat_l_transposed)
 
      my_allocate_mat_l = .true.
      IF (PRESENT(allocate_mat_l)) my_allocate_mat_l = allocate_mat_l
 
      IF (my_allocate_mat_l) THEN
         ! Create sparse variant of L
         NULLIFY (mat_l%matrix)
         ALLOCATE (mat_l%matrix)
         IF (dimen_ri == dimen_ri_red) THEN
            CALL dbcsr_create(mat_l%matrix, template=mat_template)
         ELSE
            CALL dbcsr_get_info(mat_template, nblkrows_total=nblk, col_blk_size=col_blk_size)
 
            CALL calculate_equal_blk_size(row_blk_size, dimen_ri_red, nblk)
 
            CALL dbcsr_create(mat_l%matrix, template=mat_template, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
 
            DEALLOCATE (row_blk_size)
         END IF
 
         IF (.NOT. (do_kpoints)) THEN
            CALL copy_fm_to_dbcsr(fm_mat_l(1, 1), mat_l%matrix)
         END IF
 
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE reorder_mat_l
 
! **************************************************************************************************
!> \brief ...
!> \param blk_size_new ...
!> \param dimen_RI_red ...
!> \param nblk ...
! **************************************************************************************************
   SUBROUTINE calculate_equal_blk_size(blk_size_new, dimen_RI_red, nblk)
      INTEGER, DIMENSION(:), POINTER                     :: blk_size_new
      INTEGER, INTENT(IN)                                :: dimen_ri_red, nblk
 
      INTEGER                                            :: col_per_blk, remainder
 
      NULLIFY (blk_size_new)
      ALLOCATE (blk_size_new(nblk))
 
      remainder = mod(dimen_ri_red, nblk)
      col_per_blk = dimen_ri_red/nblk
 
      ! Determine a new distribution for the columns (corresponding to the number of columns)
      IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1
      blk_size_new(remainder + 1:nblk) = col_per_blk
 
   END SUBROUTINE calculate_equal_blk_size
 
! **************************************************************************************************
!> \brief ...
!> \param fm_mat_S ...
!> \param do_ri_sos_laplace_mp2 ...
!> \param first_cycle ...
!> \param virtual ...
!> \param Eigenval ...
!> \param homo ...
!> \param omega ...
!> \param omega_old ...
!> \param jquad ...
!> \param mm_style ...
!> \param dimen_RI ...
!> \param dimen_ia ...
!> \param alpha ...
!> \param fm_mat_Q ...
!> \param fm_mat_Q_gemm ...
!> \param do_bse ...
!> \param fm_mat_Q_static_bse_gemm ...
!> \param dgemm_counter ...
!> \param num_integ_points ...
!> \param count_ev_sc_GW ...
! **************************************************************************************************
   SUBROUTINE calc_mat_q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, virtual, &
                         Eigenval, homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, &
                         do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
                         num_integ_points, count_ev_sc_GW)
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_s
      LOGICAL, INTENT(IN)                                :: do_ri_sos_laplace_mp2, first_cycle
      INTEGER, INTENT(IN)                                :: virtual
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: eigenval
      INTEGER, INTENT(IN)                                :: homo
      REAL(kind=dp), INTENT(IN)                          :: omega, omega_old
      INTEGER, INTENT(IN)                                :: jquad, mm_style, dimen_ri, dimen_ia
      REAL(kind=dp), INTENT(IN)                          :: alpha
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_q, fm_mat_q_gemm
      LOGICAL, INTENT(IN)                                :: do_bse
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_q_static_bse_gemm
      TYPE(dgemm_counter_type), INTENT(INOUT)            :: dgemm_counter
      INTEGER, INTENT(IN)                                :: num_integ_points, count_ev_sc_gw
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'calc_mat_Q'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
 
      IF (do_ri_sos_laplace_mp2) THEN
         ! the first index of tau_tj starts with 0 (see mp2_weights)
         CALL calc_fm_mat_s_laplace(fm_mat_s, homo, virtual, eigenval, omega - omega_old)
      ELSE
         CALL calc_fm_mat_s_rpa(fm_mat_s, first_cycle, virtual, eigenval, &
                                homo, omega, omega_old)
      END IF
 
      CALL contract_s_to_q(mm_style, dimen_ri, dimen_ia, alpha, fm_mat_s, fm_mat_q_gemm, &
                           fm_mat_q, dgemm_counter)
      ! fm_mat_Q_static_bse_gemm does not enter W_ijab (A matrix in TDA), but only full ABBA
      ! (since only B_ij_bar enters W_ijab)
      ! Changing jquad, since omega=0 is at last idx
      ! We enforce W0 for BSE in case of evGW
      IF (do_bse .AND. jquad == num_integ_points .AND. count_ev_sc_gw == 1) THEN
         CALL cp_fm_to_fm(fm_mat_q_gemm, fm_mat_q_static_bse_gemm)
      END IF
      CALL timestop(handle)
 
   END SUBROUTINE calc_mat_q
 
! **************************************************************************************************
!> \brief ...
!> \param fm_mat_S ...
!> \param virtual ...
!> \param Eigenval_last ...
!> \param homo ...
!> \param omega_old ...
! **************************************************************************************************
   SUBROUTINE remove_scaling_factor_rpa(fm_mat_S, virtual, Eigenval_last, homo, omega_old)
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_s
      INTEGER, INTENT(IN)                                :: virtual
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: eigenval_last
      INTEGER, INTENT(IN)                                :: homo
      REAL(kind=dp), INTENT(IN)                          :: omega_old
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'remove_scaling_factor_rpa'
 
      INTEGER                                            :: avirt, handle, i_global, iib, iocc, &
                                                            ncol_local
      INTEGER, DIMENSION(:), POINTER                     :: col_indices
      REAL(kind=dp)                                      :: eigen_diff
 
      CALL timeset(routinen, handle)
 
      ! get info of fm_mat_S
      CALL cp_fm_get_info(matrix=fm_mat_s, &
                          ncol_local=ncol_local, &
                          col_indices=col_indices)
 
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
!$OMP             SHARED(ncol_local,col_indices,Eigenval_last,fm_mat_S,virtual,homo,omega_old)
      DO iib = 1, ncol_local
         i_global = col_indices(iib)
 
         iocc = max(1, i_global - 1)/virtual + 1
         avirt = i_global - (iocc - 1)*virtual
         eigen_diff = eigenval_last(avirt + homo) - eigenval_last(iocc)
 
         fm_mat_s%local_data(:, iib) = fm_mat_s%local_data(:, iib)/ &
                                       sqrt(eigen_diff/(eigen_diff**2 + omega_old**2))
 
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE remove_scaling_factor_rpa
 
! **************************************************************************************************
!> \brief ...
!> \param fm_mat_S ...
!> \param first_cycle ...
!> \param virtual ...
!> \param Eigenval ...
!> \param homo ...
!> \param omega ...
!> \param omega_old ...
! **************************************************************************************************
   SUBROUTINE calc_fm_mat_s_rpa(fm_mat_S, first_cycle, virtual, Eigenval, homo, &
                                omega, omega_old)
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_s
      LOGICAL, INTENT(IN)                                :: first_cycle
      INTEGER, INTENT(IN)                                :: virtual
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: eigenval
      INTEGER, INTENT(IN)                                :: homo
      REAL(kind=dp), INTENT(IN)                          :: omega, omega_old
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'calc_fm_mat_S_rpa'
 
      INTEGER                                            :: avirt, handle, i_global, iib, iocc, &
                                                            ncol_local
      INTEGER, DIMENSION(:), POINTER                     :: col_indices
      REAL(kind=dp)                                      :: eigen_diff
 
      CALL timeset(routinen, handle)
 
      ! get info of fm_mat_S
      CALL cp_fm_get_info(matrix=fm_mat_s, &
                          ncol_local=ncol_local, &
                          col_indices=col_indices)
 
      ! update G matrix with the new value of omega
      IF (first_cycle) THEN
         ! In this case just update the matrix (symmetric form) with
         ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2))
         !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
         !$OMP             SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega)
         DO iib = 1, ncol_local
            i_global = col_indices(iib)
 
            iocc = max(1, i_global - 1)/virtual + 1
            avirt = i_global - (iocc - 1)*virtual
            eigen_diff = eigenval(avirt + homo) - eigenval(iocc)
 
            fm_mat_s%local_data(:, iib) = fm_mat_s%local_data(:, iib)* &
                                          sqrt(eigen_diff/(eigen_diff**2 + omega**2))
 
         END DO
      ELSE
         ! In this case the update has to remove the old omega component thus
         ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2))
         !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
         !$OMP             SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega,omega_old)
         DO iib = 1, ncol_local
            i_global = col_indices(iib)
 
            iocc = max(1, i_global - 1)/virtual + 1
            avirt = i_global - (iocc - 1)*virtual
            eigen_diff = eigenval(avirt + homo) - eigenval(iocc)
 
            fm_mat_s%local_data(:, iib) = fm_mat_s%local_data(:, iib)* &
                                          sqrt((eigen_diff**2 + omega_old**2)/(eigen_diff**2 + omega**2))
 
         END DO
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE calc_fm_mat_s_rpa
 
! **************************************************************************************************
!> \brief ...
!> \param mm_style ...
!> \param dimen_RI ...
!> \param dimen_ia ...
!> \param alpha ...
!> \param fm_mat_S ...
!> \param fm_mat_Q_gemm ...
!> \param fm_mat_Q ...
!> \param dgemm_counter ...
! **************************************************************************************************
   SUBROUTINE contract_s_to_q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
                              fm_mat_Q, dgemm_counter)
 
      INTEGER, INTENT(IN)                                :: mm_style, dimen_ri, dimen_ia
      REAL(kind=dp), INTENT(IN)                          :: alpha
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_s, fm_mat_q_gemm, fm_mat_q
      TYPE(dgemm_counter_type), INTENT(INOUT)            :: dgemm_counter
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'contract_S_to_Q'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
 
      CALL dgemm_counter_start(dgemm_counter)
      SELECT CASE (mm_style)
      CASE (wfc_mm_style_gemm)
         ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm !!!
         CALL parallel_gemm(transa="N", transb="T", m=dimen_ri, n=dimen_ri, k=dimen_ia, alpha=alpha, &
                            matrix_a=fm_mat_s, matrix_b=fm_mat_s, beta=0.0_dp, &
                            matrix_c=fm_mat_q_gemm)
      CASE (wfc_mm_style_syrk)
         ! will only compute the upper half of the matrix, which is fine, since we only use it for cholesky later
         CALL cp_fm_syrk(uplo='U', trans='N', k=dimen_ia, alpha=alpha, matrix_a=fm_mat_s, &
                         ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_q_gemm)
      CASE DEFAULT
         cpabort("")
      END SELECT
      CALL dgemm_counter_stop(dgemm_counter, dimen_ri, dimen_ri, dimen_ia)
 
      ! copy/redistribute fm_mat_Q_gemm to fm_mat_Q
      CALL cp_fm_set_all(matrix=fm_mat_q, alpha=0.0_dp)
      CALL cp_fm_to_fm_submat_general(fm_mat_q_gemm, fm_mat_q, dimen_ri, dimen_ri, 1, 1, 1, 1, &
                                      fm_mat_q_gemm%matrix_struct%context)
 
      CALL timestop(handle)
 
   END SUBROUTINE contract_s_to_q
 
! **************************************************************************************************
!> \brief ...
!> \param dimen_RI ...
!> \param trace_Qomega ...
!> \param fm_mat_Q ...
! **************************************************************************************************
   SUBROUTINE q_trace_and_add_unit_matrix(dimen_RI, trace_Qomega, fm_mat_Q)
 
      INTEGER, INTENT(IN)                                :: dimen_ri
      REAL(kind=dp), DIMENSION(dimen_RI), INTENT(OUT)    :: trace_qomega
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_q
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'Q_trace_and_add_unit_matrix'
 
      INTEGER                                            :: handle, i_global, iib, j_global, jjb, &
                                                            ncol_local, nrow_local
      INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
      TYPE(mp_para_env_type), POINTER                    :: para_env
 
      CALL timeset(routinen, handle)
 
      CALL cp_fm_get_info(matrix=fm_mat_q, &
                          nrow_local=nrow_local, &
                          ncol_local=ncol_local, &
                          row_indices=row_indices, &
                          col_indices=col_indices, &
                          para_env=para_env)
 
      ! calculate the trace of Q and add 1 on the diagonal
      trace_qomega = 0.0_dp
!$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
!$OMP                       SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,fm_mat_Q,dimen_RI)
      DO jjb = 1, ncol_local
         j_global = col_indices(jjb)
         DO iib = 1, nrow_local
            i_global = row_indices(iib)
            IF (j_global == i_global .AND. i_global <= dimen_ri) THEN
               trace_qomega(i_global) = fm_mat_q%local_data(iib, jjb)
               fm_mat_q%local_data(iib, jjb) = fm_mat_q%local_data(iib, jjb) + 1.0_dp
            END IF
         END DO
      END DO
      CALL para_env%sum(trace_qomega)
 
      CALL timestop(handle)
 
   END SUBROUTINE q_trace_and_add_unit_matrix
 
! **************************************************************************************************
!> \brief ...
!> \param dimen_RI ...
!> \param trace_Qomega ...
!> \param fm_mat_Q ...
!> \param para_env_RPA ...
!> \param Erpa ...
!> \param wjquad ...
! **************************************************************************************************
   SUBROUTINE compute_erpa_by_freq_int(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA, Erpa, wjquad)
 
      INTEGER, INTENT(IN)                                :: dimen_ri
      REAL(kind=dp), DIMENSION(dimen_RI), INTENT(IN)     :: trace_qomega
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_q
      TYPE(mp_para_env_type), INTENT(IN)                 :: para_env_rpa
      REAL(kind=dp), INTENT(INOUT)                       :: erpa
      REAL(kind=dp), INTENT(IN)                          :: wjquad
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_Erpa_by_freq_int'
 
      INTEGER                                            :: handle, i_global, iib, info_chol, &
                                                            j_global, jjb, ncol_local, nrow_local
      INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
      REAL(kind=dp)                                      :: fcomega
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: q_log
 
      CALL timeset(routinen, handle)
 
      CALL cp_fm_get_info(matrix=fm_mat_q, &
                          nrow_local=nrow_local, &
                          ncol_local=ncol_local, &
                          row_indices=row_indices, &
                          col_indices=col_indices)
 
      ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition
      CALL cp_fm_cholesky_decompose(matrix=fm_mat_q, n=dimen_ri, info_out=info_chol)
      IF (info_chol .NE. 0) THEN
         CALL cp_warn(__location__, &
                      "The Cholesky decomposition before inverting the RPA matrix / dielectric "// &
                      "function failed. "// &
                      "In case of low-scaling RPA/GW, decreasing EPS_FILTER in the &LOW_SCALING "// &
                      "section might "// &
                      "increase the overall accuracy making the matrix positive definite. "// &
                      "Code will abort.")
      END IF
 
      cpassert(info_chol == 0)
 
      ALLOCATE (q_log(dimen_ri))
      q_log = 0.0_dp
!$OMP             PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
!$OMP                         SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,fm_mat_Q,dimen_RI)
      DO jjb = 1, ncol_local
         j_global = col_indices(jjb)
         DO iib = 1, nrow_local
            i_global = row_indices(iib)
            IF (j_global == i_global .AND. i_global <= dimen_ri) THEN
               q_log(i_global) = 2.0_dp*log(fm_mat_q%local_data(iib, jjb))
            END IF
         END DO
      END DO
      CALL para_env_rpa%sum(q_log)
 
      ! the following frequency integration is Eq. (27) in M. Del Ben et al., JCTC 9, 2654 (2013)
      ! (https://doi.org/10.1021/ct4002202)
      fcomega = 0.0_dp
      DO iib = 1, dimen_ri
         IF (modulo(iib, para_env_rpa%num_pe) /= para_env_rpa%mepos) cycle
         fcomega = fcomega + (q_log(iib) - trace_qomega(iib))/2.0_dp
      END DO
      erpa = erpa + fcomega*wjquad
 
      DEALLOCATE (q_log)
 
      CALL timestop(handle)
 
   END SUBROUTINE compute_erpa_by_freq_int
 
! **************************************************************************************************
!> \brief ...
!> \param fm_struct_sub_kp ...
!> \param para_env ...
!> \param dimen_RI ...
!> \param ikp_local ...
!> \param first_ikp_local ...
! **************************************************************************************************
   SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, &
                              ikp_local, first_ikp_local)
      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_sub_kp
      TYPE(mp_para_env_type), POINTER                    :: para_env
      INTEGER, INTENT(IN)                                :: dimen_ri
      INTEGER, INTENT(OUT)                               :: ikp_local, first_ikp_local
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'get_sub_para_kp'
 
      INTEGER                                            :: color_sub_kp, handle, num_proc_per_kp
      TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_sub_kp
      TYPE(mp_para_env_type), POINTER                    :: para_env_sub_kp
 
      CALL timeset(routinen, handle)
 
      ! we use all processors for every k-point, subgroups for cp_cfm_heevd only seems to work for
      ! very small subgroups with 1, 2, or 3 MPI ranks. For more MPI-ranks, eigenvalues and
      ! eigenvectors coming out of cp_cfm_heevd are totally wrong unfortunately.
      num_proc_per_kp = para_env%num_pe
 
      ! IF(nkp > para_env%num_pe) THEN
      !   num_proc_per_kp = para_env%num_pe
      ! ELSE
      !   num_proc_per_kp = para_env%num_pe/nkp
      ! END IF
 
      color_sub_kp = para_env%mepos/num_proc_per_kp
      ALLOCATE (para_env_sub_kp)
      CALL para_env_sub_kp%from_split(para_env, color_sub_kp)
 
      ! grid_2d(1) = 1
      ! grid_2d(2) = para_env_sub_kp%num_pe
 
      NULLIFY (blacs_env_sub_kp)
      ! CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp, grid_2d=grid_2d)
      CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp)
 
      NULLIFY (fm_struct_sub_kp)
      CALL cp_fm_struct_create(fm_struct_sub_kp, context=blacs_env_sub_kp, nrow_global=dimen_ri, &
                               ncol_global=dimen_ri, para_env=para_env_sub_kp)
 
      CALL cp_blacs_env_release(blacs_env_sub_kp)
 
      ! IF(nkp > para_env%num_pe) THEN
      ! every processor has all ikp's
      ikp_local = -1
      first_ikp_local = 1
      ! ELSE
      !    ikp_local = 0
      !    first_ikp_local = 1
      !    DO ikp = 1, nkp
      !      IF(MOD(ikp-1, para_env%num_pe/num_proc_per_kp) == color_sub_kp) THEN
      !        ikp_local = ikp
      !        first_ikp_local = ikp
      !      END IF
      !    END DO
      ! END IF
 
      CALL mp_para_env_release(para_env_sub_kp)
 
      CALL timestop(handle)
 
   END SUBROUTINE get_sub_para_kp
 
! **************************************************************************************************
!> \brief ...
!> \param fm_mo_coeff_occ ...
!> \param fm_mo_coeff_virt ...
!> \param fm_scaled_dm_occ_tau ...
!> \param fm_scaled_dm_virt_tau ...
!> \param index_to_cell_3c ...
!> \param cell_to_index_3c ...
!> \param do_ic_model ...
!> \param do_kpoints_cubic_RPA ...
!> \param do_kpoints_from_Gamma ...
!> \param do_ri_Sigma_x ...
!> \param has_mat_P_blocks ...
!> \param wkp_W ...
!> \param cfm_mat_Q ...
!> \param fm_mat_Minv_L_kpoints ...
!> \param fm_mat_L_kpoints ...
!> \param fm_matrix_Minv ...
!> \param fm_matrix_Minv_Vtrunc_Minv ...
!> \param fm_mat_RI_global_work ...
!> \param fm_mat_work ...
!> \param fm_mo_coeff_occ_scaled ...
!> \param fm_mo_coeff_virt_scaled ...
!> \param mat_dm ...
!> \param mat_L ...
!> \param mat_MinvVMinv ...
!> \param mat_P_omega ...
!> \param mat_P_omega_kp ...
!> \param t_3c_M ...
!> \param t_3c_O ...
!> \param t_3c_O_compressed ...
!> \param t_3c_O_ind ...
!> \param mat_work ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, &
                              fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, &
                              cell_to_index_3c, do_ic_model, &
                              do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, &
                              has_mat_P_blocks, &
                              wkp_W, cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
                              fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
                              fm_mat_RI_global_work, fm_mat_work, &
                              fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, &
                              mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
                              t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
                              mat_work, qs_env)
 
      TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: fm_mo_coeff_occ, fm_mo_coeff_virt
      TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_scaled_dm_occ_tau, &
                                                            fm_scaled_dm_virt_tau
      INTEGER, ALLOCATABLE, DIMENSION(:, :), &
         INTENT(INOUT)                                   :: index_to_cell_3c
      INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
         INTENT(INOUT)                                   :: cell_to_index_3c
      LOGICAL, INTENT(IN)                                :: do_ic_model, do_kpoints_cubic_rpa, &
                                                            do_kpoints_from_gamma, do_ri_sigma_x
      LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
         INTENT(INOUT)                                   :: has_mat_p_blocks
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
         INTENT(INOUT)                                   :: wkp_w
      TYPE(cp_cfm_type), INTENT(INOUT)                   :: cfm_mat_q
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_mat_minv_l_kpoints, fm_mat_l_kpoints, &
                                                            fm_matrix_minv, &
                                                            fm_matrix_minv_vtrunc_minv
      TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_mat_ri_global_work, fm_mat_work, &
                                                            fm_mo_coeff_occ_scaled, &
                                                            fm_mo_coeff_virt_scaled
      TYPE(dbcsr_p_type), INTENT(INOUT)                  :: mat_dm, mat_l, mat_minvvminv
      TYPE(dbcsr_p_type), ALLOCATABLE, &
         DIMENSION(:, :, :), INTENT(INOUT)               :: mat_p_omega
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_p_omega_kp
      TYPE(dbt_type)                                     :: t_3c_m
      TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_o
      TYPE(hfx_compression_type), ALLOCATABLE, &
         DIMENSION(:, :, :), INTENT(INOUT)               :: t_3c_o_compressed
      TYPE(block_ind_type), ALLOCATABLE, &
         DIMENSION(:, :, :), INTENT(INOUT)               :: t_3c_o_ind
      TYPE(dbcsr_type), POINTER                          :: mat_work
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'dealloc_im_time'
 
      INTEGER                                            :: cut_memory, handle, i_kp, i_mem, i_size, &
                                                            ispin, j_size, jquad, nspins, unused
      LOGICAL                                            :: my_open_shell
 
      CALL timeset(routinen, handle)
 
      nspins = SIZE(fm_mo_coeff_occ)
      my_open_shell = (nspins == 2)
 
      CALL cp_fm_release(fm_scaled_dm_occ_tau)
      CALL cp_fm_release(fm_scaled_dm_virt_tau)
      DO ispin = 1, SIZE(fm_mo_coeff_occ)
         CALL cp_fm_release(fm_mo_coeff_occ(ispin))
         CALL cp_fm_release(fm_mo_coeff_virt(ispin))
      END DO
      CALL cp_fm_release(fm_mo_coeff_occ_scaled)
      CALL cp_fm_release(fm_mo_coeff_virt_scaled)
 
      CALL cp_fm_release(fm_mat_minv_l_kpoints)
      CALL cp_fm_release(fm_mat_l_kpoints)
 
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
         CALL cp_fm_release(fm_matrix_minv_vtrunc_minv)
         CALL cp_fm_release(fm_matrix_minv)
      END IF
 
      CALL cp_fm_release(fm_mat_work)
 
      IF (.NOT. (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma)) THEN
         CALL dbcsr_release(mat_work)
         DEALLOCATE (mat_work)
      END IF
 
      CALL dbcsr_release(mat_l%matrix)
      DEALLOCATE (mat_l%matrix)
 
      IF (do_ri_sigma_x .OR. do_ic_model) THEN
         CALL dbcsr_release(mat_minvvminv%matrix)
         DEALLOCATE (mat_minvvminv%matrix)
      END IF
      IF (do_ri_sigma_x) THEN
         CALL dbcsr_release(mat_dm%matrix)
         DEALLOCATE (mat_dm%matrix)
      END IF
 
      DEALLOCATE (index_to_cell_3c, cell_to_index_3c)
 
      IF (ALLOCATED(mat_p_omega)) THEN
         DO ispin = 1, SIZE(mat_p_omega, 3)
            DO i_kp = 1, SIZE(mat_p_omega, 2)
               DO jquad = 1, SIZE(mat_p_omega, 1)
                  CALL dbcsr_deallocate_matrix(mat_p_omega(jquad, i_kp, ispin)%matrix)
               END DO
            END DO
         END DO
         DEALLOCATE (mat_p_omega)
      END IF
 
      DO i_size = 1, SIZE(t_3c_o, 1)
         DO j_size = 1, SIZE(t_3c_o, 2)
            CALL dbt_destroy(t_3c_o(i_size, j_size))
         END DO
      END DO
 
      DEALLOCATE (t_3c_o)
      CALL dbt_destroy(t_3c_m)
 
      DEALLOCATE (has_mat_p_blocks)
 
      IF (do_kpoints_cubic_rpa .OR. do_kpoints_from_gamma) THEN
         CALL cp_cfm_release(cfm_mat_q)
         CALL cp_fm_release(fm_mat_ri_global_work)
         CALL dbcsr_deallocate_matrix_set(mat_p_omega_kp)
         DEALLOCATE (wkp_w)
      END IF
 
      cut_memory = SIZE(t_3c_o_compressed, 3)
 
      DEALLOCATE (t_3c_o_ind)
      DO i_size = 1, SIZE(t_3c_o_compressed, 1)
         DO j_size = 1, SIZE(t_3c_o_compressed, 2)
            DO i_mem = 1, cut_memory
               CALL dealloc_containers(t_3c_o_compressed(i_size, j_size, i_mem), unused)
            END DO
         END DO
      END DO
      DEALLOCATE (t_3c_o_compressed)
 
      IF (do_kpoints_from_gamma) THEN
         CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_G)
         IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
            CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma)
            CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc)
         END IF
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE dealloc_im_time
 
! **************************************************************************************************
!> \brief ...
!> \param mat_P_omega ...
!> \param mat_L ...
!> \param mat_work ...
!> \param eps_filter_im_time ...
!> \param fm_mat_work ...
!> \param dimen_RI ...
!> \param dimen_RI_red ...
!> \param fm_mat_L ...
!> \param fm_mat_Q ...
! **************************************************************************************************
   SUBROUTINE contract_p_omega_with_mat_l(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, &
                                          dimen_RI_red, fm_mat_L, fm_mat_Q)
 
      TYPE(dbcsr_type), INTENT(IN)                       :: mat_p_omega, mat_l
      TYPE(dbcsr_type), INTENT(INOUT)                    :: mat_work
      REAL(kind=dp), INTENT(IN)                          :: eps_filter_im_time
      TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_mat_work
      INTEGER, INTENT(IN)                                :: dimen_ri, dimen_ri_red
      TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_l, fm_mat_q
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'contract_P_omega_with_mat_L'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
 
      ! multiplication with RI metric/Coulomb operator
      CALL dbcsr_multiply("N", "T", 1.0_dp, mat_p_omega, mat_l, &
                          0.0_dp, mat_work, filter_eps=eps_filter_im_time)
 
      CALL copy_dbcsr_to_fm(mat_work, fm_mat_work)
 
      CALL parallel_gemm('N', 'N', dimen_ri_red, dimen_ri_red, dimen_ri, 1.0_dp, fm_mat_l, fm_mat_work, &
                         0.0_dp, fm_mat_q)
 
      ! Reset mat_work to save memory
      CALL dbcsr_set(mat_work, 0.0_dp)
      CALL dbcsr_filter(mat_work, 1.0_dp)
 
      CALL timestop(handle)
 
   END SUBROUTINE contract_p_omega_with_mat_l
 
END MODULE rpa_util