LCOV - code coverage report
Current view: top level - src - rpa_util.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.1 % 341 338
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 14 14

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Utility functions for RPA calculations
      10              : !> \par History
      11              : !>      06.2019 Moved from rpa_ri_gpw.F [Frederick Stein]
      12              : ! **************************************************************************************************
      13              : MODULE rpa_util
      14              : 
      15              :    USE cell_types,                      ONLY: cell_type,&
      16              :                                               get_cell
      17              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
      18              :                                               cp_blacs_env_release,&
      19              :                                               cp_blacs_env_type
      20              :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      21              :                                               cp_cfm_release,&
      22              :                                               cp_cfm_set_all,&
      23              :                                               cp_cfm_type
      24              :    USE cp_dbcsr_api,                    ONLY: &
      25              :         dbcsr_create, dbcsr_deallocate_matrix, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, &
      26              :         dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
      27              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      28              :                                               copy_fm_to_dbcsr,&
      29              :                                               dbcsr_allocate_matrix_set,&
      30              :                                               dbcsr_deallocate_matrix_set
      31              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_syrk,&
      32              :                                               cp_fm_transpose
      33              :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose
      34              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      35              :                                               cp_fm_struct_release,&
      36              :                                               cp_fm_struct_type
      37              :    USE cp_fm_types,                     ONLY: &
      38              :         cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_release, cp_fm_set_all, &
      39              :         cp_fm_set_element, cp_fm_to_fm, cp_fm_to_fm_submat_general, cp_fm_type
      40              :    USE dbt_api,                         ONLY: dbt_destroy,&
      41              :                                               dbt_type
      42              :    USE dgemm_counter_types,             ONLY: dgemm_counter_start,&
      43              :                                               dgemm_counter_stop,&
      44              :                                               dgemm_counter_type
      45              :    USE hfx_types,                       ONLY: block_ind_type,&
      46              :                                               dealloc_containers,&
      47              :                                               hfx_compression_type
      48              :    USE input_constants,                 ONLY: wfc_mm_style_gemm,&
      49              :                                               wfc_mm_style_syrk
      50              :    USE kinds,                           ONLY: dp
      51              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      52              :                                               kpoint_release,&
      53              :                                               kpoint_type
      54              :    USE mathconstants,                   ONLY: z_zero
      55              :    USE message_passing,                 ONLY: mp_para_env_release,&
      56              :                                               mp_para_env_type
      57              :    USE mp2_laplace,                     ONLY: calc_fm_mat_S_laplace
      58              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      59              :    USE qs_environment_types,            ONLY: get_qs_env,&
      60              :                                               qs_environment_type
      61              :    USE rpa_gw_kpoints_util,             ONLY: compute_wkp_W
      62              : #include "./base/base_uses.f90"
      63              : 
      64              :    IMPLICIT NONE
      65              : 
      66              :    PRIVATE
      67              : 
      68              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_util'
      69              : 
      70              :    PUBLIC :: compute_Erpa_by_freq_int, alloc_im_time, calc_mat_Q, Q_trace_and_add_unit_matrix, &
      71              :              dealloc_im_time, contract_P_omega_with_mat_L, calc_fm_mat_S_rpa, remove_scaling_factor_rpa
      72              : 
      73              : CONTAINS
      74              : 
      75              : ! **************************************************************************************************
      76              : !> \brief ...
      77              : !> \param qs_env ...
      78              : !> \param para_env ...
      79              : !> \param dimen_RI ...
      80              : !> \param dimen_RI_red ...
      81              : !> \param num_integ_points ...
      82              : !> \param nspins ...
      83              : !> \param fm_mat_Q ...
      84              : !> \param fm_mo_coeff_occ ...
      85              : !> \param fm_mo_coeff_virt ...
      86              : !> \param fm_matrix_Minv_L_kpoints ...
      87              : !> \param fm_matrix_L_kpoints ...
      88              : !> \param mat_P_global ...
      89              : !> \param t_3c_O ...
      90              : !> \param matrix_s ...
      91              : !> \param kpoints ...
      92              : !> \param eps_filter_im_time ...
      93              : !> \param cut_memory ...
      94              : !> \param nkp ...
      95              : !> \param num_cells_dm ...
      96              : !> \param num_3c_repl ...
      97              : !> \param size_P ...
      98              : !> \param ikp_local ...
      99              : !> \param index_to_cell_3c ...
     100              : !> \param cell_to_index_3c ...
     101              : !> \param col_blk_size ...
     102              : !> \param do_ic_model ...
     103              : !> \param do_kpoints_cubic_RPA ...
     104              : !> \param do_kpoints_from_Gamma ...
     105              : !> \param do_ri_Sigma_x ...
     106              : !> \param my_open_shell ...
     107              : !> \param has_mat_P_blocks ...
     108              : !> \param wkp_W ...
     109              : !> \param cfm_mat_Q ...
     110              : !> \param fm_mat_Minv_L_kpoints ...
     111              : !> \param fm_mat_L_kpoints ...
     112              : !> \param fm_mat_RI_global_work ...
     113              : !> \param fm_mat_work ...
     114              : !> \param fm_mo_coeff_occ_scaled ...
     115              : !> \param fm_mo_coeff_virt_scaled ...
     116              : !> \param mat_dm ...
     117              : !> \param mat_L ...
     118              : !> \param mat_M_P_munu_occ ...
     119              : !> \param mat_M_P_munu_virt ...
     120              : !> \param mat_MinvVMinv ...
     121              : !> \param mat_P_omega ...
     122              : !> \param mat_P_omega_kp ...
     123              : !> \param mat_work ...
     124              : !> \param mo_coeff ...
     125              : !> \param fm_scaled_dm_occ_tau ...
     126              : !> \param fm_scaled_dm_virt_tau ...
     127              : !> \param homo ...
     128              : !> \param nmo ...
     129              : ! **************************************************************************************************
     130          268 :    SUBROUTINE alloc_im_time(qs_env, para_env, dimen_RI, dimen_RI_red, num_integ_points, nspins, &
     131              :                             fm_mat_Q, fm_mo_coeff_occ, fm_mo_coeff_virt, &
     132              :                             fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, mat_P_global, &
     133              :                             t_3c_O, matrix_s, kpoints, eps_filter_im_time, &
     134              :                             cut_memory, nkp, num_cells_dm, num_3c_repl, &
     135              :                             size_P, ikp_local, &
     136              :                             index_to_cell_3c, &
     137              :                             cell_to_index_3c, &
     138              :                             col_blk_size, &
     139              :                             do_ic_model, do_kpoints_cubic_RPA, &
     140              :                             do_kpoints_from_Gamma, do_ri_Sigma_x, my_open_shell, &
     141              :                             has_mat_P_blocks, wkp_W, &
     142              :                             cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
     143              :                             fm_mat_RI_global_work, fm_mat_work, fm_mo_coeff_occ_scaled, &
     144              :                             fm_mo_coeff_virt_scaled, mat_dm, mat_L, mat_M_P_munu_occ, mat_M_P_munu_virt, &
     145              :                             mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
     146          134 :                             mat_work, mo_coeff, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, homo, nmo)
     147              : 
     148              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     149              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     150              :       INTEGER, INTENT(IN)                                :: dimen_RI, dimen_RI_red, &
     151              :                                                             num_integ_points, nspins
     152              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_Q
     153              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_mo_coeff_occ, fm_mo_coeff_virt
     154              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_matrix_Minv_L_kpoints, &
     155              :                                                             fm_matrix_L_kpoints
     156              :       TYPE(dbcsr_p_type), INTENT(IN)                     :: mat_P_global
     157              :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
     158              :          INTENT(INOUT)                                   :: t_3c_O
     159              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     160              :       TYPE(kpoint_type), POINTER                         :: kpoints
     161              :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter_im_time
     162              :       INTEGER, INTENT(IN)                                :: cut_memory
     163              :       INTEGER, INTENT(OUT)                               :: nkp, num_cells_dm, num_3c_repl, size_P, &
     164              :                                                             ikp_local
     165              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(OUT) :: index_to_cell_3c
     166              :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
     167              :          INTENT(OUT)                                     :: cell_to_index_3c
     168              :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size
     169              :       LOGICAL, INTENT(IN)                                :: do_ic_model, do_kpoints_cubic_RPA, &
     170              :                                                             do_kpoints_from_Gamma, do_ri_Sigma_x, &
     171              :                                                             my_open_shell
     172              :       LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
     173              :          INTENT(OUT)                                     :: has_mat_P_blocks
     174              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     175              :          INTENT(OUT)                                     :: wkp_W
     176              :       TYPE(cp_cfm_type), INTENT(OUT)                     :: cfm_mat_Q
     177              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_mat_Minv_L_kpoints, fm_mat_L_kpoints
     178              :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_mat_RI_global_work, fm_mat_work, &
     179              :                                                             fm_mo_coeff_occ_scaled, &
     180              :                                                             fm_mo_coeff_virt_scaled
     181              :       TYPE(dbcsr_p_type), INTENT(OUT)                    :: mat_dm, mat_L, mat_M_P_munu_occ, &
     182              :                                                             mat_M_P_munu_virt, mat_MinvVMinv
     183              :       TYPE(dbcsr_p_type), ALLOCATABLE, &
     184              :          DIMENSION(:, :, :)                              :: mat_P_omega
     185              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_P_omega_kp
     186              :       TYPE(dbcsr_type), POINTER                          :: mat_work
     187              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff
     188              :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_scaled_dm_occ_tau, &
     189              :                                                             fm_scaled_dm_virt_tau
     190              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
     191              :       INTEGER, INTENT(IN)                                :: nmo
     192              : 
     193              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'alloc_im_time'
     194              : 
     195              :       INTEGER                                            :: cell_grid_dm(3), first_ikp_local, &
     196              :                                                             handle, i_dim, i_kp, ispin, jquad, &
     197              :                                                             nspins_P_omega, periodic(3)
     198          134 :       INTEGER, DIMENSION(:), POINTER                     :: row_blk_size
     199          134 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: wkp_V
     200              :       TYPE(cell_type), POINTER                           :: cell
     201              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_sub_kp
     202              : 
     203          134 :       CALL timeset(routineN, handle)
     204              : 
     205          860 :       ALLOCATE (fm_mo_coeff_occ(nspins), fm_mo_coeff_virt(nspins))
     206              : 
     207          134 :       CALL cp_fm_create(fm_scaled_dm_occ_tau, mo_coeff(1)%matrix_struct)
     208          134 :       CALL cp_fm_set_all(fm_scaled_dm_occ_tau, 0.0_dp)
     209              : 
     210          134 :       CALL cp_fm_create(fm_scaled_dm_virt_tau, mo_coeff(1)%matrix_struct)
     211          134 :       CALL cp_fm_set_all(fm_scaled_dm_virt_tau, 0.0_dp)
     212              : 
     213          296 :       DO ispin = 1, SIZE(mo_coeff)
     214              :          CALL create_occ_virt_mo_coeffs(fm_mo_coeff_occ(ispin), fm_mo_coeff_virt(ispin), mo_coeff(ispin), &
     215          296 :                                         nmo, homo(ispin))
     216              :       END DO
     217              : 
     218          134 :       num_3c_repl = SIZE(t_3c_O, 2)
     219              : 
     220          134 :       IF (do_kpoints_cubic_RPA) THEN
     221              :          ! we always use an odd number of image cells
     222              :          ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
     223           16 :          DO i_dim = 1, 3
     224           16 :             cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
     225              :          END DO
     226            4 :          num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
     227           12 :          ALLOCATE (index_to_cell_3c(3, SIZE(kpoints%index_to_cell, 2)))
     228            4 :          CPASSERT(SIZE(kpoints%index_to_cell, 1) == 3)
     229           84 :          index_to_cell_3c(:, :) = kpoints%index_to_cell(:, :)
     230            0 :          ALLOCATE (cell_to_index_3c(LBOUND(kpoints%cell_to_index, 1):UBOUND(kpoints%cell_to_index, 1), &
     231              :                                     LBOUND(kpoints%cell_to_index, 2):UBOUND(kpoints%cell_to_index, 2), &
     232           20 :                                     LBOUND(kpoints%cell_to_index, 3):UBOUND(kpoints%cell_to_index, 3)))
     233           64 :          cell_to_index_3c(:, :, :) = kpoints%cell_to_index(:, :, :)
     234              : 
     235              :       ELSE
     236          130 :          ALLOCATE (index_to_cell_3c(3, 1))
     237          520 :          index_to_cell_3c(:, 1) = 0
     238          130 :          ALLOCATE (cell_to_index_3c(0:0, 0:0, 0:0))
     239          130 :          cell_to_index_3c(0, 0, 0) = 1
     240          130 :          num_cells_dm = 1
     241              :       END IF
     242              : 
     243          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
     244              : 
     245           20 :          CALL get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, ikp_local, first_ikp_local)
     246              : 
     247           20 :          CALL cp_cfm_create(cfm_mat_Q, fm_struct_sub_kp)
     248           20 :          CALL cp_cfm_set_all(cfm_mat_Q, z_zero)
     249              :       ELSE
     250          114 :          first_ikp_local = 1
     251              :       END IF
     252              : 
     253              :       ! if we do kpoints, mat_P has a kpoint and mat_P_omega has the inted
     254              :       ! mat_P(tau, kpoint)
     255          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
     256              : 
     257           20 :          NULLIFY (cell)
     258           20 :          CALL get_qs_env(qs_env, cell=cell)
     259           20 :          CALL get_cell(cell=cell, periodic=periodic)
     260              : 
     261           20 :          CALL get_kpoint_info(kpoints, nkp=nkp)
     262              :          ! compute k-point weights such that functions 1/k^2, 1/k and const function are
     263              :          ! integrated correctly
     264           20 :          CALL compute_wkp_W(qs_env, wkp_W, wkp_V, kpoints, cell%h_inv, periodic)
     265           20 :          DEALLOCATE (wkp_V)
     266              : 
     267              :       ELSE
     268          114 :          nkp = 1
     269              :       END IF
     270              : 
     271          134 :       IF (do_kpoints_cubic_RPA) THEN
     272            4 :          size_P = MAX(num_cells_dm/2 + 1, nkp)
     273          130 :       ELSE IF (do_kpoints_from_Gamma) THEN
     274           16 :          size_P = MAX(3**(periodic(1) + periodic(2) + periodic(3)), nkp)
     275              :       ELSE
     276          114 :          size_P = 1
     277              :       END IF
     278              : 
     279          134 :       nspins_P_omega = 1
     280          134 :       IF (my_open_shell) nspins_P_omega = 2
     281              : 
     282         5352 :       ALLOCATE (mat_P_omega(num_integ_points, size_P, nspins_P_omega))
     283          296 :       DO ispin = 1, nspins_P_omega
     284          912 :          DO i_kp = 1, size_P
     285         4682 :             DO jquad = 1, num_integ_points
     286         3904 :                NULLIFY (mat_P_omega(jquad, i_kp, ispin)%matrix)
     287         3904 :                ALLOCATE (mat_P_omega(jquad, i_kp, ispin)%matrix)
     288              :                CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp, ispin)%matrix, &
     289         3904 :                                  template=mat_P_global%matrix)
     290         4520 :                CALL dbcsr_set(mat_P_omega(jquad, i_kp, ispin)%matrix, 0.0_dp)
     291              :             END DO
     292              :          END DO
     293              :       END DO
     294              : 
     295          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
     296           20 :          CALL alloc_mat_P_omega(mat_P_omega_kp, 2, size_P, mat_P_global%matrix)
     297              :       END IF
     298              : 
     299          134 :       CALL cp_fm_create(fm_mo_coeff_occ_scaled, fm_mo_coeff_occ(1)%matrix_struct)
     300          134 :       CALL cp_fm_to_fm(fm_mo_coeff_occ(1), fm_mo_coeff_occ_scaled)
     301          134 :       CALL cp_fm_set_all(matrix=fm_mo_coeff_occ_scaled, alpha=0.0_dp)
     302              : 
     303          134 :       CALL cp_fm_create(fm_mo_coeff_virt_scaled, fm_mo_coeff_virt(1)%matrix_struct)
     304          134 :       CALL cp_fm_to_fm(fm_mo_coeff_virt(1), fm_mo_coeff_virt_scaled)
     305          134 :       CALL cp_fm_set_all(matrix=fm_mo_coeff_virt_scaled, alpha=0.0_dp)
     306              : 
     307          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
     308           20 :          CALL cp_fm_create(fm_mat_RI_global_work, fm_matrix_Minv_L_kpoints(1, 1)%matrix_struct)
     309           20 :          CALL cp_fm_to_fm(fm_matrix_Minv_L_kpoints(1, 1), fm_mat_RI_global_work)
     310           20 :          CALL cp_fm_set_all(fm_mat_RI_global_work, 0.0_dp)
     311              :       END IF
     312              : 
     313          938 :       ALLOCATE (has_mat_P_blocks(num_cells_dm/2 + 1, cut_memory, cut_memory, num_3c_repl, num_3c_repl))
     314         3000 :       has_mat_P_blocks = .TRUE.
     315              : 
     316          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
     317              :          CALL reorder_mat_L(fm_mat_Minv_L_kpoints, fm_matrix_Minv_L_kpoints, fm_mat_Q%matrix_struct, para_env, mat_L, &
     318              :                             mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp, &
     319           20 :                             allocate_mat_L=.FALSE.)
     320              : 
     321              :          CALL reorder_mat_L(fm_mat_L_kpoints, fm_matrix_L_kpoints, fm_mat_Q%matrix_struct, para_env, mat_L, &
     322           20 :                             mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp)
     323              : 
     324           20 :          CALL cp_fm_struct_release(fm_struct_sub_kp)
     325              : 
     326              :       ELSE
     327              :          CALL reorder_mat_L(fm_mat_Minv_L_kpoints, fm_matrix_Minv_L_kpoints, fm_mat_Q%matrix_struct, para_env, mat_L, &
     328          114 :                             mat_P_global%matrix, dimen_RI, dimen_RI_red, first_ikp_local)
     329              :       END IF
     330              : 
     331              :       ! Create Scalapack working matrix for the contraction with the metric
     332          134 :       IF (dimen_RI == dimen_RI_red) THEN
     333          130 :          CALL cp_fm_create(fm_mat_work, fm_mat_Q%matrix_struct)
     334          130 :          CALL cp_fm_to_fm(fm_mat_Q, fm_mat_work)
     335          130 :          CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
     336              : 
     337              :       ELSE
     338            4 :          NULLIFY (fm_struct)
     339              :          CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_mat_Q%matrix_struct, &
     340            4 :                                   ncol_global=dimen_RI_red, nrow_global=dimen_RI)
     341              : 
     342            4 :          CALL cp_fm_create(fm_mat_work, fm_struct)
     343            4 :          CALL cp_fm_set_all(matrix=fm_mat_work, alpha=0.0_dp)
     344              : 
     345            4 :          CALL cp_fm_struct_release(fm_struct)
     346              : 
     347              :       END IF
     348              : 
     349              :       ! Then its DBCSR counter part
     350          134 :       IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
     351          114 :          CALL dbcsr_get_info(mat_L%matrix, col_blk_size=col_blk_size, row_blk_size=row_blk_size)
     352              : 
     353              :          ! Create mat_work having the shape of the transposed of mat_L (compare with contract_P_omega_with_mat_L)
     354              :          NULLIFY (mat_work)
     355          114 :          ALLOCATE (mat_work)
     356          114 :          CALL dbcsr_create(mat_work, template=mat_L%matrix, row_blk_size=col_blk_size, col_blk_size=row_blk_size)
     357              :       END IF
     358              : 
     359          134 :       IF (do_ri_Sigma_x .OR. do_ic_model) THEN
     360              : 
     361              :          NULLIFY (mat_MinvVMinv%matrix)
     362          104 :          ALLOCATE (mat_MinvVMinv%matrix)
     363          104 :          CALL dbcsr_create(mat_MinvVMinv%matrix, template=mat_P_global%matrix)
     364          104 :          CALL dbcsr_set(mat_MinvVMinv%matrix, 0.0_dp)
     365              : 
     366              :          ! for kpoints we compute SinvVSinv later with kpoints
     367          104 :          IF (.NOT. do_kpoints_from_Gamma) THEN
     368              : 
     369              :             !  get the Coulomb matrix for Sigma_x = G*V
     370              :             CALL dbcsr_multiply("T", "N", 1.0_dp, mat_L%matrix, mat_L%matrix, &
     371           92 :                                 0.0_dp, mat_MinvVMinv%matrix, filter_eps=eps_filter_im_time)
     372              : 
     373              :          END IF
     374              : 
     375              :       END IF
     376              : 
     377          134 :       IF (do_ri_Sigma_x) THEN
     378              : 
     379              :          NULLIFY (mat_dm%matrix)
     380          104 :          ALLOCATE (mat_dm%matrix)
     381          104 :          CALL dbcsr_create(mat_dm%matrix, template=matrix_s(1)%matrix)
     382              : 
     383              :       END IF
     384              : 
     385          134 :       CALL timestop(handle)
     386              : 
     387          804 :    END SUBROUTINE alloc_im_time
     388              : 
     389              : ! **************************************************************************************************
     390              : !> \brief ...
     391              : !> \param fm_mo_coeff_occ ...
     392              : !> \param fm_mo_coeff_virt ...
     393              : !> \param mo_coeff ...
     394              : !> \param nmo ...
     395              : !> \param homo ...
     396              : ! **************************************************************************************************
     397          486 :    SUBROUTINE create_occ_virt_mo_coeffs(fm_mo_coeff_occ, fm_mo_coeff_virt, mo_coeff, &
     398              :                                         nmo, homo)
     399              : 
     400              :       TYPE(cp_fm_type), INTENT(OUT)                      :: fm_mo_coeff_occ, fm_mo_coeff_virt
     401              :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_coeff
     402              :       INTEGER, INTENT(IN)                                :: nmo, homo
     403              : 
     404              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_occ_virt_mo_coeffs'
     405              : 
     406              :       INTEGER                                            :: handle, icol_global, irow_global
     407              : 
     408          162 :       CALL timeset(routineN, handle)
     409              : 
     410          162 :       CALL cp_fm_create(fm_mo_coeff_occ, mo_coeff%matrix_struct)
     411          162 :       CALL cp_fm_set_all(fm_mo_coeff_occ, 0.0_dp)
     412          162 :       CALL cp_fm_to_fm(mo_coeff, fm_mo_coeff_occ)
     413              : 
     414              :       ! set all virtual MO coeffs to zero
     415         3552 :       DO icol_global = homo + 1, nmo
     416        98204 :          DO irow_global = 1, nmo
     417        98042 :             CALL cp_fm_set_element(fm_mo_coeff_occ, irow_global, icol_global, 0.0_dp)
     418              :          END DO
     419              :       END DO
     420              : 
     421          162 :       CALL cp_fm_create(fm_mo_coeff_virt, mo_coeff%matrix_struct)
     422          162 :       CALL cp_fm_set_all(fm_mo_coeff_virt, 0.0_dp)
     423          162 :       CALL cp_fm_to_fm(mo_coeff, fm_mo_coeff_virt)
     424              : 
     425              :       ! set all occupied MO coeffs to zero
     426          766 :       DO icol_global = 1, homo
     427        15928 :          DO irow_global = 1, nmo
     428        15766 :             CALL cp_fm_set_element(fm_mo_coeff_virt, irow_global, icol_global, 0.0_dp)
     429              :          END DO
     430              :       END DO
     431              : 
     432          162 :       CALL timestop(handle)
     433              : 
     434          162 :    END SUBROUTINE create_occ_virt_mo_coeffs
     435              : 
     436              : ! **************************************************************************************************
     437              : !> \brief ...
     438              : !> \param mat_P_omega ...
     439              : !> \param num_integ_points ...
     440              : !> \param size_P ...
     441              : !> \param template ...
     442              : ! **************************************************************************************************
     443           20 :    SUBROUTINE alloc_mat_P_omega(mat_P_omega, num_integ_points, size_P, template)
     444              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_P_omega
     445              :       INTEGER, INTENT(IN)                                :: num_integ_points, size_P
     446              :       TYPE(dbcsr_type), POINTER                          :: template
     447              : 
     448              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'alloc_mat_P_omega'
     449              : 
     450              :       INTEGER                                            :: handle, i_kp, jquad
     451              : 
     452           20 :       CALL timeset(routineN, handle)
     453              : 
     454           20 :       NULLIFY (mat_P_omega)
     455           20 :       CALL dbcsr_allocate_matrix_set(mat_P_omega, num_integ_points, size_P)
     456          444 :       DO i_kp = 1, size_P
     457         1292 :          DO jquad = 1, num_integ_points
     458          848 :             ALLOCATE (mat_P_omega(jquad, i_kp)%matrix)
     459              :             CALL dbcsr_create(matrix=mat_P_omega(jquad, i_kp)%matrix, &
     460          848 :                               template=template)
     461         1272 :             CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp)
     462              :          END DO
     463              :       END DO
     464              : 
     465           20 :       CALL timestop(handle)
     466              : 
     467           20 :    END SUBROUTINE alloc_mat_P_omega
     468              : 
     469              : ! **************************************************************************************************
     470              : !> \brief ...
     471              : !> \param fm_mat_L ...
     472              : !> \param fm_matrix_Minv_L_kpoints ...
     473              : !> \param fm_struct_template ...
     474              : !> \param para_env ...
     475              : !> \param mat_L ...
     476              : !> \param mat_template ...
     477              : !> \param dimen_RI ...
     478              : !> \param dimen_RI_red ...
     479              : !> \param first_ikp_local ...
     480              : !> \param ikp_local ...
     481              : !> \param fm_struct_sub_kp ...
     482              : !> \param allocate_mat_L ...
     483              : ! **************************************************************************************************
     484          154 :    SUBROUTINE reorder_mat_L(fm_mat_L, fm_matrix_Minv_L_kpoints, fm_struct_template, para_env, mat_L, mat_template, &
     485              :                             dimen_RI, dimen_RI_red, first_ikp_local, ikp_local, fm_struct_sub_kp, allocate_mat_L)
     486              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_mat_L, fm_matrix_Minv_L_kpoints
     487              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_template
     488              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     489              :       TYPE(dbcsr_p_type), INTENT(OUT)                    :: mat_L
     490              :       TYPE(dbcsr_type), INTENT(IN)                       :: mat_template
     491              :       INTEGER, INTENT(IN)                                :: dimen_RI, dimen_RI_red, first_ikp_local
     492              :       INTEGER, OPTIONAL                                  :: ikp_local
     493              :       TYPE(cp_fm_struct_type), OPTIONAL, POINTER         :: fm_struct_sub_kp
     494              :       LOGICAL, INTENT(IN), OPTIONAL                      :: allocate_mat_L
     495              : 
     496              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'reorder_mat_L'
     497              : 
     498              :       INTEGER                                            :: handle, ikp, j_size, nblk
     499          154 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_size, row_blk_size
     500              :       LOGICAL                                            :: do_kpoints, my_allocate_mat_L
     501              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     502              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
     503              :       TYPE(cp_fm_type)                                   :: fm_mat_L_transposed, fmdummy
     504              : 
     505          154 :       CALL timeset(routineN, handle)
     506              : 
     507          154 :       do_kpoints = .FALSE.
     508          154 :       IF (PRESENT(ikp_local) .AND. PRESENT(fm_struct_sub_kp)) THEN
     509           40 :          do_kpoints = .TRUE.
     510              :       END IF
     511              : 
     512              :       ! Get the fm_struct for fm_mat_L
     513          154 :       NULLIFY (fm_struct)
     514          154 :       IF (dimen_RI == dimen_RI_red) THEN
     515          150 :          fm_struct => fm_struct_template
     516              :       ELSE
     517              :          ! The template is assumed to be square such that we need a new fm_struct if dimensions are not equal
     518            4 :          CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI_red, ncol_global=dimen_RI, template_fmstruct=fm_struct_template)
     519              :       END IF
     520              : 
     521              :       ! Start to allocate the new full matrix
     522         2620 :       ALLOCATE (fm_mat_L(SIZE(fm_matrix_Minv_L_kpoints, 1), SIZE(fm_matrix_Minv_L_kpoints, 2)))
     523          348 :       DO j_size = 1, SIZE(fm_matrix_Minv_L_kpoints, 2)
     524         2158 :          DO ikp = 1, SIZE(fm_matrix_Minv_L_kpoints, 1)
     525         2004 :             IF (do_kpoints) THEN
     526         1696 :                IF (ikp == first_ikp_local .OR. ikp_local == -1) THEN
     527         1696 :                   CALL cp_fm_create(fm_mat_L(ikp, j_size), fm_struct_sub_kp)
     528         1696 :                   CALL cp_fm_set_all(fm_mat_L(ikp, j_size), 0.0_dp)
     529              :                END IF
     530              :             ELSE
     531          114 :                CALL cp_fm_create(fm_mat_L(ikp, j_size), fm_struct)
     532          114 :                CALL cp_fm_set_all(fm_mat_L(ikp, j_size), 0.0_dp)
     533              :             END IF
     534              :          END DO
     535              :       END DO
     536              : 
     537              :       ! For the transposed matric we need a different fm_struct
     538          154 :       IF (dimen_RI == dimen_RI_red) THEN
     539          150 :          fm_struct => fm_mat_L(first_ikp_local, 1)%matrix_struct
     540              :       ELSE
     541            4 :          CALL cp_fm_struct_release(fm_struct)
     542              : 
     543              :          ! Create a fm_struct with transposed sizes
     544            4 :          NULLIFY (fm_struct)
     545              :          CALL cp_fm_struct_create(fm_struct, nrow_global=dimen_RI, ncol_global=dimen_RI_red, &
     546            4 :                                   template_fmstruct=fm_mat_L(first_ikp_local, 1)%matrix_struct) !, force_block=.TRUE.)
     547              :       END IF
     548              : 
     549              :       ! Allocate buffer matrix
     550          154 :       CALL cp_fm_create(fm_mat_L_transposed, fm_struct)
     551          154 :       CALL cp_fm_set_all(matrix=fm_mat_L_transposed, alpha=0.0_dp)
     552              : 
     553          154 :       IF (dimen_RI /= dimen_RI_red) CALL cp_fm_struct_release(fm_struct)
     554              : 
     555          154 :       CALL cp_fm_get_info(fm_mat_L_transposed, context=blacs_env)
     556              : 
     557              :       ! For k-points copy matrices of your group
     558              :       ! Without kpoints, transpose matrix
     559              :       ! without kpoints, the size of fm_mat_L is 1x1. with kpoints, the size is N_kpoints x 2 (2 for real/complex)
     560          348 :       DO j_size = 1, SIZE(fm_matrix_Minv_L_kpoints, 2)
     561         2158 :       DO ikp = 1, SIZE(fm_matrix_Minv_L_kpoints, 1)
     562         2004 :          IF (do_kpoints) THEN
     563         1696 :             IF (ikp_local == ikp .OR. ikp_local == -1) THEN
     564         1696 :                CALL cp_fm_copy_general(fm_matrix_Minv_L_kpoints(ikp, j_size), fm_mat_L_transposed, para_env)
     565         1696 :                CALL cp_fm_to_fm(fm_mat_L_transposed, fm_mat_L(ikp, j_size))
     566              :             ELSE
     567            0 :                CALL cp_fm_copy_general(fm_matrix_Minv_L_kpoints(ikp, j_size), fmdummy, para_env)
     568              :             END IF
     569              :          ELSE
     570          114 :             CALL cp_fm_copy_general(fm_matrix_Minv_L_kpoints(ikp, j_size), fm_mat_L_transposed, blacs_env%para_env)
     571          114 :             CALL cp_fm_transpose(fm_mat_L_transposed, fm_mat_L(ikp, j_size))
     572              :          END IF
     573              :       END DO
     574              :       END DO
     575              : 
     576              :       ! Release old matrix
     577          154 :       CALL cp_fm_release(fm_matrix_Minv_L_kpoints)
     578              :       ! Release buffer
     579          154 :       CALL cp_fm_release(fm_mat_L_transposed)
     580              : 
     581          154 :       my_allocate_mat_L = .TRUE.
     582          154 :       IF (PRESENT(allocate_mat_L)) my_allocate_mat_L = allocate_mat_L
     583              : 
     584           20 :       IF (my_allocate_mat_L) THEN
     585              :          ! Create sparse variant of L
     586              :          NULLIFY (mat_L%matrix)
     587          134 :          ALLOCATE (mat_L%matrix)
     588          134 :          IF (dimen_RI == dimen_RI_red) THEN
     589          130 :             CALL dbcsr_create(mat_L%matrix, template=mat_template)
     590              :          ELSE
     591            4 :             CALL dbcsr_get_info(mat_template, nblkrows_total=nblk, col_blk_size=col_blk_size)
     592              : 
     593            4 :             CALL calculate_equal_blk_size(row_blk_size, dimen_RI_red, nblk)
     594              : 
     595            4 :             CALL dbcsr_create(mat_L%matrix, template=mat_template, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
     596              : 
     597            4 :             DEALLOCATE (row_blk_size)
     598              :          END IF
     599              : 
     600          134 :          IF (.NOT. (do_kpoints)) THEN
     601          114 :             CALL copy_fm_to_dbcsr(fm_mat_L(1, 1), mat_L%matrix)
     602              :          END IF
     603              : 
     604              :       END IF
     605              : 
     606          154 :       CALL timestop(handle)
     607              : 
     608          308 :    END SUBROUTINE reorder_mat_L
     609              : 
     610              : ! **************************************************************************************************
     611              : !> \brief ...
     612              : !> \param blk_size_new ...
     613              : !> \param dimen_RI_red ...
     614              : !> \param nblk ...
     615              : ! **************************************************************************************************
     616            4 :    SUBROUTINE calculate_equal_blk_size(blk_size_new, dimen_RI_red, nblk)
     617              :       INTEGER, DIMENSION(:), POINTER                     :: blk_size_new
     618              :       INTEGER, INTENT(IN)                                :: dimen_RI_red, nblk
     619              : 
     620              :       INTEGER                                            :: col_per_blk, remainder
     621              : 
     622            4 :       NULLIFY (blk_size_new)
     623           12 :       ALLOCATE (blk_size_new(nblk))
     624              : 
     625            4 :       remainder = MOD(dimen_RI_red, nblk)
     626            4 :       col_per_blk = dimen_RI_red/nblk
     627              : 
     628              :       ! Determine a new distribution for the columns (corresponding to the number of columns)
     629           10 :       IF (remainder > 0) blk_size_new(1:remainder) = col_per_blk + 1
     630           10 :       blk_size_new(remainder + 1:nblk) = col_per_blk
     631              : 
     632            4 :    END SUBROUTINE calculate_equal_blk_size
     633              : 
     634              : ! **************************************************************************************************
     635              : !> \brief ...
     636              : !> \param fm_mat_S ...
     637              : !> \param do_ri_sos_laplace_mp2 ...
     638              : !> \param first_cycle ...
     639              : !> \param virtual ...
     640              : !> \param Eigenval ...
     641              : !> \param homo ...
     642              : !> \param omega ...
     643              : !> \param omega_old ...
     644              : !> \param jquad ...
     645              : !> \param mm_style ...
     646              : !> \param dimen_RI ...
     647              : !> \param dimen_ia ...
     648              : !> \param alpha ...
     649              : !> \param fm_mat_Q ...
     650              : !> \param fm_mat_Q_gemm ...
     651              : !> \param do_bse ...
     652              : !> \param fm_mat_Q_static_bse_gemm ...
     653              : !> \param dgemm_counter ...
     654              : !> \param num_integ_points ...
     655              : !> \param count_ev_sc_GW ...
     656              : ! **************************************************************************************************
     657        21132 :    SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, virtual, &
     658        10566 :                          Eigenval, homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, &
     659              :                          do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
     660              :                          num_integ_points, count_ev_sc_GW)
     661              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_S
     662              :       LOGICAL, INTENT(IN)                                :: do_ri_sos_laplace_mp2, first_cycle
     663              :       INTEGER, INTENT(IN)                                :: virtual
     664              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval
     665              :       INTEGER, INTENT(IN)                                :: homo
     666              :       REAL(KIND=dp), INTENT(IN)                          :: omega, omega_old
     667              :       INTEGER, INTENT(IN)                                :: jquad, mm_style, dimen_RI, dimen_ia
     668              :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     669              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_Q, fm_mat_Q_gemm
     670              :       LOGICAL, INTENT(IN)                                :: do_bse
     671              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_Q_static_bse_gemm
     672              :       TYPE(dgemm_counter_type), INTENT(INOUT)            :: dgemm_counter
     673              :       INTEGER, INTENT(IN)                                :: num_integ_points, count_ev_sc_GW
     674              : 
     675              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'calc_mat_Q'
     676              : 
     677              :       INTEGER                                            :: handle
     678              : 
     679        10566 :       CALL timeset(routineN, handle)
     680              : 
     681        10566 :       IF (do_ri_sos_laplace_mp2) THEN
     682              :          ! the first index of tau_tj starts with 0 (see mp2_weights)
     683          128 :          CALL calc_fm_mat_S_laplace(fm_mat_S, homo, virtual, Eigenval, omega - omega_old)
     684              :       ELSE
     685              :          CALL calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, &
     686        10438 :                                 homo, omega, omega_old)
     687              :       END IF
     688              : 
     689              :       CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
     690        10566 :                            fm_mat_Q, dgemm_counter)
     691              :       ! fm_mat_Q_static_bse_gemm does not enter W_ijab (A matrix in TDA), but only full ABBA
     692              :       ! (since only B_ij_bar enters W_ijab)
     693              :       ! Changing jquad, since omega=0 is at last idx
     694              :       ! We enforce W0 for BSE in case of evGW
     695        10566 :       IF (do_bse .AND. jquad == num_integ_points .AND. count_ev_sc_GW == 1) THEN
     696           30 :          CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm)
     697              :       END IF
     698        10566 :       CALL timestop(handle)
     699              : 
     700        10566 :    END SUBROUTINE calc_mat_Q
     701              : 
     702              : ! **************************************************************************************************
     703              : !> \brief ...
     704              : !> \param fm_mat_S ...
     705              : !> \param virtual ...
     706              : !> \param Eigenval_last ...
     707              : !> \param homo ...
     708              : !> \param omega_old ...
     709              : ! **************************************************************************************************
     710          244 :    SUBROUTINE remove_scaling_factor_rpa(fm_mat_S, virtual, Eigenval_last, homo, omega_old)
     711              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_S
     712              :       INTEGER, INTENT(IN)                                :: virtual
     713              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval_last
     714              :       INTEGER, INTENT(IN)                                :: homo
     715              :       REAL(KIND=dp), INTENT(IN)                          :: omega_old
     716              : 
     717              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'remove_scaling_factor_rpa'
     718              : 
     719              :       INTEGER                                            :: avirt, handle, i_global, iiB, iocc, &
     720              :                                                             ncol_local
     721          244 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices
     722              :       REAL(KIND=dp)                                      :: eigen_diff
     723              : 
     724          244 :       CALL timeset(routineN, handle)
     725              : 
     726              :       ! get info of fm_mat_S
     727              :       CALL cp_fm_get_info(matrix=fm_mat_S, &
     728              :                           ncol_local=ncol_local, &
     729          244 :                           col_indices=col_indices)
     730              : 
     731              : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
     732          244 : !$OMP             SHARED(ncol_local,col_indices,Eigenval_last,fm_mat_S,virtual,homo,omega_old)
     733              :       DO iiB = 1, ncol_local
     734              :          i_global = col_indices(iiB)
     735              : 
     736              :          iocc = MAX(1, i_global - 1)/virtual + 1
     737              :          avirt = i_global - (iocc - 1)*virtual
     738              :          eigen_diff = Eigenval_last(avirt + homo) - Eigenval_last(iocc)
     739              : 
     740              :          fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)/ &
     741              :                                        SQRT(eigen_diff/(eigen_diff**2 + omega_old**2))
     742              : 
     743              :       END DO
     744              : 
     745          244 :       CALL timestop(handle)
     746              : 
     747          244 :    END SUBROUTINE remove_scaling_factor_rpa
     748              : 
     749              : ! **************************************************************************************************
     750              : !> \brief ...
     751              : !> \param fm_mat_S ...
     752              : !> \param first_cycle ...
     753              : !> \param virtual ...
     754              : !> \param Eigenval ...
     755              : !> \param homo ...
     756              : !> \param omega ...
     757              : !> \param omega_old ...
     758              : ! **************************************************************************************************
     759        10528 :    SUBROUTINE calc_fm_mat_S_rpa(fm_mat_S, first_cycle, virtual, Eigenval, homo, &
     760              :                                 omega, omega_old)
     761              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_S
     762              :       LOGICAL, INTENT(IN)                                :: first_cycle
     763              :       INTEGER, INTENT(IN)                                :: virtual
     764              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval
     765              :       INTEGER, INTENT(IN)                                :: homo
     766              :       REAL(KIND=dp), INTENT(IN)                          :: omega, omega_old
     767              : 
     768              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'calc_fm_mat_S_rpa'
     769              : 
     770              :       INTEGER                                            :: avirt, handle, i_global, iiB, iocc, &
     771              :                                                             ncol_local
     772        10528 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices
     773              :       REAL(KIND=dp)                                      :: eigen_diff
     774              : 
     775        10528 :       CALL timeset(routineN, handle)
     776              : 
     777              :       ! get info of fm_mat_S
     778              :       CALL cp_fm_get_info(matrix=fm_mat_S, &
     779              :                           ncol_local=ncol_local, &
     780        10528 :                           col_indices=col_indices)
     781              : 
     782              :       ! update G matrix with the new value of omega
     783        10528 :       IF (first_cycle) THEN
     784              :          ! In this case just update the matrix (symmetric form) with
     785              :          ! SQRT((epsi_a-epsi_i)/((epsi_a-epsi_i)**2+omega**2))
     786              :          !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
     787          374 :          !$OMP             SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega)
     788              :          DO iiB = 1, ncol_local
     789              :             i_global = col_indices(iiB)
     790              : 
     791              :             iocc = MAX(1, i_global - 1)/virtual + 1
     792              :             avirt = i_global - (iocc - 1)*virtual
     793              :             eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
     794              : 
     795              :             fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)* &
     796              :                                           SQRT(eigen_diff/(eigen_diff**2 + omega**2))
     797              : 
     798              :          END DO
     799              :       ELSE
     800              :          ! In this case the update has to remove the old omega component thus
     801              :          ! SQRT(((epsi_a-epsi_i)**2+omega_old**2)/((epsi_a-epsi_i)**2+omega**2))
     802              :          !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(iiB,iocc,avirt,eigen_diff,i_global) &
     803        10154 :          !$OMP             SHARED(ncol_local,col_indices,Eigenval,fm_mat_S,virtual,homo,omega,omega_old)
     804              :          DO iiB = 1, ncol_local
     805              :             i_global = col_indices(iiB)
     806              : 
     807              :             iocc = MAX(1, i_global - 1)/virtual + 1
     808              :             avirt = i_global - (iocc - 1)*virtual
     809              :             eigen_diff = Eigenval(avirt + homo) - Eigenval(iocc)
     810              : 
     811              :             fm_mat_S%local_data(:, iiB) = fm_mat_S%local_data(:, iiB)* &
     812              :                                           SQRT((eigen_diff**2 + omega_old**2)/(eigen_diff**2 + omega**2))
     813              : 
     814              :          END DO
     815              :       END IF
     816              : 
     817        10528 :       CALL timestop(handle)
     818              : 
     819        10528 :    END SUBROUTINE calc_fm_mat_S_rpa
     820              : 
     821              : ! **************************************************************************************************
     822              : !> \brief ...
     823              : !> \param mm_style ...
     824              : !> \param dimen_RI ...
     825              : !> \param dimen_ia ...
     826              : !> \param alpha ...
     827              : !> \param fm_mat_S ...
     828              : !> \param fm_mat_Q_gemm ...
     829              : !> \param fm_mat_Q ...
     830              : !> \param dgemm_counter ...
     831              : ! **************************************************************************************************
     832        10566 :    SUBROUTINE contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
     833              :                               fm_mat_Q, dgemm_counter)
     834              : 
     835              :       INTEGER, INTENT(IN)                                :: mm_style, dimen_RI, dimen_ia
     836              :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     837              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_S, fm_mat_Q_gemm, fm_mat_Q
     838              :       TYPE(dgemm_counter_type), INTENT(INOUT)            :: dgemm_counter
     839              : 
     840              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'contract_S_to_Q'
     841              : 
     842              :       INTEGER                                            :: handle
     843              : 
     844        10566 :       CALL timeset(routineN, handle)
     845              : 
     846        10566 :       CALL dgemm_counter_start(dgemm_counter)
     847        21122 :       SELECT CASE (mm_style)
     848              :       CASE (wfc_mm_style_gemm)
     849              :          ! waste-fully computes the full symmetrix matrix, but maybe faster than cp_fm_syrk for optimized cp_fm_gemm !!!
     850              :          CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_ia, alpha=alpha, &
     851              :                             matrix_a=fm_mat_S, matrix_b=fm_mat_S, beta=0.0_dp, &
     852        10556 :                             matrix_c=fm_mat_Q_gemm)
     853              :       CASE (wfc_mm_style_syrk)
     854              :          ! will only compute the upper half of the matrix, which is fine, since we only use it for cholesky later
     855              :          CALL cp_fm_syrk(uplo='U', trans='N', k=dimen_ia, alpha=alpha, matrix_a=fm_mat_S, &
     856           10 :                          ia=1, ja=1, beta=0.0_dp, matrix_c=fm_mat_Q_gemm)
     857              :       CASE DEFAULT
     858        10566 :          CPABORT("")
     859              :       END SELECT
     860        10566 :       CALL dgemm_counter_stop(dgemm_counter, dimen_RI, dimen_RI, dimen_ia)
     861              : 
     862              :       ! copy/redistribute fm_mat_Q_gemm to fm_mat_Q
     863        10566 :       CALL cp_fm_set_all(matrix=fm_mat_Q, alpha=0.0_dp)
     864              :       CALL cp_fm_to_fm_submat_general(fm_mat_Q_gemm, fm_mat_Q, dimen_RI, dimen_RI, 1, 1, 1, 1, &
     865        10566 :                                       fm_mat_Q_gemm%matrix_struct%context)
     866              : 
     867        10566 :       CALL timestop(handle)
     868              : 
     869        10566 :    END SUBROUTINE contract_S_to_Q
     870              : 
     871              : ! **************************************************************************************************
     872              : !> \brief ...
     873              : !> \param dimen_RI ...
     874              : !> \param trace_Qomega ...
     875              : !> \param fm_mat_Q ...
     876              : ! **************************************************************************************************
     877        22648 :    SUBROUTINE Q_trace_and_add_unit_matrix(dimen_RI, trace_Qomega, fm_mat_Q)
     878              : 
     879              :       INTEGER, INTENT(IN)                                :: dimen_RI
     880              :       REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(OUT)    :: trace_Qomega
     881              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_Q
     882              : 
     883              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'Q_trace_and_add_unit_matrix'
     884              : 
     885              :       INTEGER                                            :: handle, i_global, iiB, j_global, jjB, &
     886              :                                                             ncol_local, nrow_local
     887        11324 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     888              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     889              : 
     890        11324 :       CALL timeset(routineN, handle)
     891              : 
     892              :       CALL cp_fm_get_info(matrix=fm_mat_Q, &
     893              :                           nrow_local=nrow_local, &
     894              :                           ncol_local=ncol_local, &
     895              :                           row_indices=row_indices, &
     896              :                           col_indices=col_indices, &
     897        11324 :                           para_env=para_env)
     898              : 
     899              :       ! calculate the trace of Q and add 1 on the diagonal
     900       938066 :       trace_Qomega = 0.0_dp
     901              : !$OMP           PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
     902        11324 : !$OMP                       SHARED(ncol_local,nrow_local,col_indices,row_indices,trace_Qomega,fm_mat_Q,dimen_RI)
     903              :       DO jjB = 1, ncol_local
     904              :          j_global = col_indices(jjB)
     905              :          DO iiB = 1, nrow_local
     906              :             i_global = row_indices(iiB)
     907              :             IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
     908              :                trace_Qomega(i_global) = fm_mat_Q%local_data(iiB, jjB)
     909              :                fm_mat_Q%local_data(iiB, jjB) = fm_mat_Q%local_data(iiB, jjB) + 1.0_dp
     910              :             END IF
     911              :          END DO
     912              :       END DO
     913        11324 :       CALL para_env%sum(trace_Qomega)
     914              : 
     915        11324 :       CALL timestop(handle)
     916              : 
     917        11324 :    END SUBROUTINE Q_trace_and_add_unit_matrix
     918              : 
     919              : ! **************************************************************************************************
     920              : !> \brief ...
     921              : !> \param dimen_RI ...
     922              : !> \param trace_Qomega ...
     923              : !> \param fm_mat_Q ...
     924              : !> \param para_env_RPA ...
     925              : !> \param Erpa ...
     926              : !> \param wjquad ...
     927              : ! **************************************************************************************************
     928        11204 :    SUBROUTINE compute_Erpa_by_freq_int(dimen_RI, trace_Qomega, fm_mat_Q, para_env_RPA, Erpa, wjquad)
     929              : 
     930              :       INTEGER, INTENT(IN)                                :: dimen_RI
     931              :       REAL(KIND=dp), DIMENSION(dimen_RI), INTENT(IN)     :: trace_Qomega
     932              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_Q
     933              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env_RPA
     934              :       REAL(KIND=dp), INTENT(INOUT)                       :: Erpa
     935              :       REAL(KIND=dp), INTENT(IN)                          :: wjquad
     936              : 
     937              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Erpa_by_freq_int'
     938              : 
     939              :       INTEGER                                            :: handle, i_global, iiB, info_chol, &
     940              :                                                             j_global, jjB, ncol_local, nrow_local
     941        11204 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     942              :       REAL(KIND=dp)                                      :: FComega
     943        11204 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Q_log
     944              : 
     945        11204 :       CALL timeset(routineN, handle)
     946              : 
     947              :       CALL cp_fm_get_info(matrix=fm_mat_Q, &
     948              :                           nrow_local=nrow_local, &
     949              :                           ncol_local=ncol_local, &
     950              :                           row_indices=row_indices, &
     951        11204 :                           col_indices=col_indices)
     952              : 
     953              :       ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition
     954        11204 :       CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q, n=dimen_RI, info_out=info_chol)
     955        11204 :       IF (info_chol .NE. 0) THEN
     956              :          CALL cp_warn(__LOCATION__, &
     957              :                       "The Cholesky decomposition before inverting the RPA matrix / dielectric "// &
     958              :                       "function failed. "// &
     959              :                       "In case of low-scaling RPA/GW, decreasing EPS_FILTER in the &LOW_SCALING "// &
     960              :                       "section might "// &
     961              :                       "increase the overall accuracy making the matrix positive definite. "// &
     962            0 :                       "Code will abort.")
     963              :       END IF
     964              : 
     965        11204 :       CPASSERT(info_chol == 0)
     966              : 
     967        33612 :       ALLOCATE (Q_log(dimen_RI))
     968       929258 :       Q_log = 0.0_dp
     969              : !$OMP             PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
     970        11204 : !$OMP                         SHARED(ncol_local,nrow_local,col_indices,row_indices,Q_log,fm_mat_Q,dimen_RI)
     971              :       DO jjB = 1, ncol_local
     972              :          j_global = col_indices(jjB)
     973              :          DO iiB = 1, nrow_local
     974              :             i_global = row_indices(iiB)
     975              :             IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
     976              :                Q_log(i_global) = 2.0_dp*LOG(fm_mat_Q%local_data(iiB, jjB))
     977              :             END IF
     978              :          END DO
     979              :       END DO
     980        11204 :       CALL para_env_RPA%sum(Q_log)
     981              : 
     982              :       ! the following frequency integration is Eq. (27) in M. Del Ben et al., JCTC 9, 2654 (2013)
     983              :       ! (https://doi.org/10.1021/ct4002202)
     984        11204 :       FComega = 0.0_dp
     985       929258 :       DO iiB = 1, dimen_RI
     986       918054 :          IF (MODULO(iiB, para_env_RPA%num_pe) /= para_env_RPA%mepos) CYCLE
     987       929258 :          FComega = FComega + (Q_log(iiB) - trace_Qomega(iiB))/2.0_dp
     988              :       END DO
     989        11204 :       Erpa = Erpa + FComega*wjquad
     990              : 
     991        11204 :       DEALLOCATE (Q_log)
     992              : 
     993        11204 :       CALL timestop(handle)
     994              : 
     995        22408 :    END SUBROUTINE compute_Erpa_by_freq_int
     996              : 
     997              : ! **************************************************************************************************
     998              : !> \brief ...
     999              : !> \param fm_struct_sub_kp ...
    1000              : !> \param para_env ...
    1001              : !> \param dimen_RI ...
    1002              : !> \param ikp_local ...
    1003              : !> \param first_ikp_local ...
    1004              : ! **************************************************************************************************
    1005           20 :    SUBROUTINE get_sub_para_kp(fm_struct_sub_kp, para_env, dimen_RI, &
    1006              :                               ikp_local, first_ikp_local)
    1007              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_sub_kp
    1008              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1009              :       INTEGER, INTENT(IN)                                :: dimen_RI
    1010              :       INTEGER, INTENT(OUT)                               :: ikp_local, first_ikp_local
    1011              : 
    1012              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'get_sub_para_kp'
    1013              : 
    1014              :       INTEGER                                            :: color_sub_kp, handle, num_proc_per_kp
    1015              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env_sub_kp
    1016              :       TYPE(mp_para_env_type), POINTER                    :: para_env_sub_kp
    1017              : 
    1018           20 :       CALL timeset(routineN, handle)
    1019              : 
    1020              :       ! we use all processors for every k-point, subgroups for cp_cfm_heevd only seems to work for
    1021              :       ! very small subgroups with 1, 2, or 3 MPI ranks. For more MPI-ranks, eigenvalues and
    1022              :       ! eigenvectors coming out of cp_cfm_heevd are totally wrong unfortunately.
    1023           20 :       num_proc_per_kp = para_env%num_pe
    1024              : 
    1025              :       ! IF(nkp > para_env%num_pe) THEN
    1026              :       !   num_proc_per_kp = para_env%num_pe
    1027              :       ! ELSE
    1028              :       !   num_proc_per_kp = para_env%num_pe/nkp
    1029              :       ! END IF
    1030              : 
    1031           20 :       color_sub_kp = para_env%mepos/num_proc_per_kp
    1032           20 :       ALLOCATE (para_env_sub_kp)
    1033           20 :       CALL para_env_sub_kp%from_split(para_env, color_sub_kp)
    1034              : 
    1035              :       ! grid_2d(1) = 1
    1036              :       ! grid_2d(2) = para_env_sub_kp%num_pe
    1037              : 
    1038           20 :       NULLIFY (blacs_env_sub_kp)
    1039              :       ! CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp, grid_2d=grid_2d)
    1040           20 :       CALL cp_blacs_env_create(blacs_env=blacs_env_sub_kp, para_env=para_env_sub_kp)
    1041              : 
    1042           20 :       NULLIFY (fm_struct_sub_kp)
    1043              :       CALL cp_fm_struct_create(fm_struct_sub_kp, context=blacs_env_sub_kp, nrow_global=dimen_RI, &
    1044           20 :                                ncol_global=dimen_RI, para_env=para_env_sub_kp)
    1045              : 
    1046           20 :       CALL cp_blacs_env_release(blacs_env_sub_kp)
    1047              : 
    1048              :       ! IF(nkp > para_env%num_pe) THEN
    1049              :       ! every processor has all ikp's
    1050           20 :       ikp_local = -1
    1051           20 :       first_ikp_local = 1
    1052              :       ! ELSE
    1053              :       !    ikp_local = 0
    1054              :       !    first_ikp_local = 1
    1055              :       !    DO ikp = 1, nkp
    1056              :       !      IF(MOD(ikp-1, para_env%num_pe/num_proc_per_kp) == color_sub_kp) THEN
    1057              :       !        ikp_local = ikp
    1058              :       !        first_ikp_local = ikp
    1059              :       !      END IF
    1060              :       !    END DO
    1061              :       ! END IF
    1062              : 
    1063           20 :       CALL mp_para_env_release(para_env_sub_kp)
    1064              : 
    1065           20 :       CALL timestop(handle)
    1066              : 
    1067           20 :    END SUBROUTINE get_sub_para_kp
    1068              : 
    1069              : ! **************************************************************************************************
    1070              : !> \brief ...
    1071              : !> \param fm_mo_coeff_occ ...
    1072              : !> \param fm_mo_coeff_virt ...
    1073              : !> \param fm_scaled_dm_occ_tau ...
    1074              : !> \param fm_scaled_dm_virt_tau ...
    1075              : !> \param index_to_cell_3c ...
    1076              : !> \param cell_to_index_3c ...
    1077              : !> \param do_ic_model ...
    1078              : !> \param do_kpoints_cubic_RPA ...
    1079              : !> \param do_kpoints_from_Gamma ...
    1080              : !> \param do_ri_Sigma_x ...
    1081              : !> \param has_mat_P_blocks ...
    1082              : !> \param wkp_W ...
    1083              : !> \param cfm_mat_Q ...
    1084              : !> \param fm_mat_Minv_L_kpoints ...
    1085              : !> \param fm_mat_L_kpoints ...
    1086              : !> \param fm_matrix_Minv ...
    1087              : !> \param fm_matrix_Minv_Vtrunc_Minv ...
    1088              : !> \param fm_mat_RI_global_work ...
    1089              : !> \param fm_mat_work ...
    1090              : !> \param fm_mo_coeff_occ_scaled ...
    1091              : !> \param fm_mo_coeff_virt_scaled ...
    1092              : !> \param mat_dm ...
    1093              : !> \param mat_L ...
    1094              : !> \param mat_MinvVMinv ...
    1095              : !> \param mat_P_omega ...
    1096              : !> \param mat_P_omega_kp ...
    1097              : !> \param t_3c_M ...
    1098              : !> \param t_3c_O ...
    1099              : !> \param t_3c_O_compressed ...
    1100              : !> \param t_3c_O_ind ...
    1101              : !> \param mat_work ...
    1102              : !> \param qs_env ...
    1103              : ! **************************************************************************************************
    1104          134 :    SUBROUTINE dealloc_im_time(fm_mo_coeff_occ, fm_mo_coeff_virt, &
    1105              :                               fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, index_to_cell_3c, &
    1106              :                               cell_to_index_3c, do_ic_model, &
    1107              :                               do_kpoints_cubic_RPA, do_kpoints_from_Gamma, do_ri_Sigma_x, &
    1108              :                               has_mat_P_blocks, &
    1109              :                               wkp_W, cfm_mat_Q, fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
    1110              :                               fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
    1111              :                               fm_mat_RI_global_work, fm_mat_work, &
    1112              :                               fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm, mat_L, &
    1113              :                               mat_MinvVMinv, mat_P_omega, mat_P_omega_kp, &
    1114              :                               t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
    1115              :                               mat_work, qs_env)
    1116              : 
    1117              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: fm_mo_coeff_occ, fm_mo_coeff_virt
    1118              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_scaled_dm_occ_tau, &
    1119              :                                                             fm_scaled_dm_virt_tau
    1120              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), &
    1121              :          INTENT(INOUT)                                   :: index_to_cell_3c
    1122              :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
    1123              :          INTENT(INOUT)                                   :: cell_to_index_3c
    1124              :       LOGICAL, INTENT(IN)                                :: do_ic_model, do_kpoints_cubic_RPA, &
    1125              :                                                             do_kpoints_from_Gamma, do_ri_Sigma_x
    1126              :       LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :), &
    1127              :          INTENT(INOUT)                                   :: has_mat_P_blocks
    1128              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
    1129              :          INTENT(INOUT)                                   :: wkp_W
    1130              :       TYPE(cp_cfm_type), INTENT(INOUT)                   :: cfm_mat_Q
    1131              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_mat_Minv_L_kpoints, fm_mat_L_kpoints, &
    1132              :                                                             fm_matrix_Minv, &
    1133              :                                                             fm_matrix_Minv_Vtrunc_Minv
    1134              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_mat_RI_global_work, fm_mat_work, &
    1135              :                                                             fm_mo_coeff_occ_scaled, &
    1136              :                                                             fm_mo_coeff_virt_scaled
    1137              :       TYPE(dbcsr_p_type), INTENT(INOUT)                  :: mat_dm, mat_L, mat_MinvVMinv
    1138              :       TYPE(dbcsr_p_type), ALLOCATABLE, &
    1139              :          DIMENSION(:, :, :), INTENT(INOUT)               :: mat_P_omega
    1140              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_P_omega_kp
    1141              :       TYPE(dbt_type)                                     :: t_3c_M
    1142              :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_O
    1143              :       TYPE(hfx_compression_type), ALLOCATABLE, &
    1144              :          DIMENSION(:, :, :), INTENT(INOUT)               :: t_3c_O_compressed
    1145              :       TYPE(block_ind_type), ALLOCATABLE, &
    1146              :          DIMENSION(:, :, :), INTENT(INOUT)               :: t_3c_O_ind
    1147              :       TYPE(dbcsr_type), POINTER                          :: mat_work
    1148              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1149              : 
    1150              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'dealloc_im_time'
    1151              : 
    1152              :       INTEGER                                            :: cut_memory, handle, i_kp, i_mem, i_size, &
    1153              :                                                             ispin, j_size, jquad, nspins, unused
    1154              :       LOGICAL                                            :: my_open_shell
    1155              : 
    1156          134 :       CALL timeset(routineN, handle)
    1157              : 
    1158          134 :       nspins = SIZE(fm_mo_coeff_occ)
    1159          134 :       my_open_shell = (nspins == 2)
    1160              : 
    1161          134 :       CALL cp_fm_release(fm_scaled_dm_occ_tau)
    1162          134 :       CALL cp_fm_release(fm_scaled_dm_virt_tau)
    1163          296 :       DO ispin = 1, SIZE(fm_mo_coeff_occ)
    1164          162 :          CALL cp_fm_release(fm_mo_coeff_occ(ispin))
    1165          296 :          CALL cp_fm_release(fm_mo_coeff_virt(ispin))
    1166              :       END DO
    1167          134 :       CALL cp_fm_release(fm_mo_coeff_occ_scaled)
    1168          134 :       CALL cp_fm_release(fm_mo_coeff_virt_scaled)
    1169              : 
    1170          134 :       CALL cp_fm_release(fm_mat_Minv_L_kpoints)
    1171          134 :       CALL cp_fm_release(fm_mat_L_kpoints)
    1172              : 
    1173          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
    1174           20 :          CALL cp_fm_release(fm_matrix_Minv_Vtrunc_Minv)
    1175           20 :          CALL cp_fm_release(fm_matrix_Minv)
    1176              :       END IF
    1177              : 
    1178          134 :       CALL cp_fm_release(fm_mat_work)
    1179              : 
    1180          134 :       IF (.NOT. (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma)) THEN
    1181          114 :          CALL dbcsr_release(mat_work)
    1182          114 :          DEALLOCATE (mat_work)
    1183              :       END IF
    1184              : 
    1185          134 :       CALL dbcsr_release(mat_L%matrix)
    1186          134 :       DEALLOCATE (mat_L%matrix)
    1187              : 
    1188          134 :       IF (do_ri_Sigma_x .OR. do_ic_model) THEN
    1189          104 :          CALL dbcsr_release(mat_MinvVMinv%matrix)
    1190          104 :          DEALLOCATE (mat_MinvVMinv%matrix)
    1191              :       END IF
    1192          134 :       IF (do_ri_Sigma_x) THEN
    1193          104 :          CALL dbcsr_release(mat_dm%matrix)
    1194          104 :          DEALLOCATE (mat_dm%matrix)
    1195              :       END IF
    1196              : 
    1197          134 :       DEALLOCATE (index_to_cell_3c, cell_to_index_3c)
    1198              : 
    1199          134 :       IF (ALLOCATED(mat_P_omega)) THEN
    1200          296 :          DO ispin = 1, SIZE(mat_P_omega, 3)
    1201          912 :             DO i_kp = 1, SIZE(mat_P_omega, 2)
    1202         4682 :                DO jquad = 1, SIZE(mat_P_omega, 1)
    1203         4520 :                   CALL dbcsr_deallocate_matrix(mat_P_omega(jquad, i_kp, ispin)%matrix)
    1204              :                END DO
    1205              :             END DO
    1206              :          END DO
    1207          134 :          DEALLOCATE (mat_P_omega)
    1208              :       END IF
    1209              : 
    1210          284 :       DO j_size = 1, SIZE(t_3c_O, 2)
    1211          514 :          DO i_size = 1, SIZE(t_3c_O, 1)
    1212          380 :             CALL dbt_destroy(t_3c_O(i_size, j_size))
    1213              :          END DO
    1214              :       END DO
    1215              : 
    1216          364 :       DEALLOCATE (t_3c_O)
    1217          134 :       CALL dbt_destroy(t_3c_M)
    1218              : 
    1219          134 :       DEALLOCATE (has_mat_P_blocks)
    1220              : 
    1221          134 :       IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
    1222           20 :          CALL cp_cfm_release(cfm_mat_Q)
    1223           20 :          CALL cp_fm_release(fm_mat_RI_global_work)
    1224           20 :          CALL dbcsr_deallocate_matrix_set(mat_P_omega_kp)
    1225           20 :          DEALLOCATE (wkp_W)
    1226              :       END IF
    1227              : 
    1228          134 :       cut_memory = SIZE(t_3c_O_compressed, 3)
    1229              : 
    1230          546 :       DEALLOCATE (t_3c_O_ind)
    1231          402 :       DO i_mem = 1, cut_memory
    1232          694 :          DO j_size = 1, SIZE(t_3c_O_compressed, 2)
    1233          972 :             DO i_size = 1, SIZE(t_3c_O_compressed, 1)
    1234          704 :                CALL dealloc_containers(t_3c_O_compressed(i_size, j_size, i_mem), unused)
    1235              :             END DO
    1236              :          END DO
    1237              :       END DO
    1238          134 :       DEALLOCATE (t_3c_O_compressed)
    1239              : 
    1240          134 :       IF (do_kpoints_from_Gamma) THEN
    1241           16 :          CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_G)
    1242           16 :          IF (qs_env%mp2_env%ri_g0w0%do_kpoints_Sigma) THEN
    1243           16 :             CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma)
    1244           16 :             CALL kpoint_release(qs_env%mp2_env%ri_rpa_im_time%kpoints_Sigma_no_xc)
    1245              :          END IF
    1246              :       END IF
    1247              : 
    1248          134 :       CALL timestop(handle)
    1249              : 
    1250          134 :    END SUBROUTINE dealloc_im_time
    1251              : 
    1252              : ! **************************************************************************************************
    1253              : !> \brief ...
    1254              : !> \param mat_P_omega ...
    1255              : !> \param mat_L ...
    1256              : !> \param mat_work ...
    1257              : !> \param eps_filter_im_time ...
    1258              : !> \param fm_mat_work ...
    1259              : !> \param dimen_RI ...
    1260              : !> \param dimen_RI_red ...
    1261              : !> \param fm_mat_L ...
    1262              : !> \param fm_mat_Q ...
    1263              : ! **************************************************************************************************
    1264         1276 :    SUBROUTINE contract_P_omega_with_mat_L(mat_P_omega, mat_L, mat_work, eps_filter_im_time, fm_mat_work, dimen_RI, &
    1265              :                                           dimen_RI_red, fm_mat_L, fm_mat_Q)
    1266              : 
    1267              :       TYPE(dbcsr_type), INTENT(IN)                       :: mat_P_omega, mat_L
    1268              :       TYPE(dbcsr_type), INTENT(INOUT)                    :: mat_work
    1269              :       REAL(KIND=dp), INTENT(IN)                          :: eps_filter_im_time
    1270              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_mat_work
    1271              :       INTEGER, INTENT(IN)                                :: dimen_RI, dimen_RI_red
    1272              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mat_L, fm_mat_Q
    1273              : 
    1274              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_P_omega_with_mat_L'
    1275              : 
    1276              :       INTEGER                                            :: handle
    1277              : 
    1278         1276 :       CALL timeset(routineN, handle)
    1279              : 
    1280              :       ! multiplication with RI metric/Coulomb operator
    1281              :       CALL dbcsr_multiply("N", "T", 1.0_dp, mat_P_omega, mat_L, &
    1282         1276 :                           0.0_dp, mat_work, filter_eps=eps_filter_im_time)
    1283              : 
    1284         1276 :       CALL copy_dbcsr_to_fm(mat_work, fm_mat_work)
    1285              : 
    1286              :       CALL parallel_gemm('N', 'N', dimen_RI_red, dimen_RI_red, dimen_RI, 1.0_dp, fm_mat_L, fm_mat_work, &
    1287         1276 :                          0.0_dp, fm_mat_Q)
    1288              : 
    1289              :       ! Reset mat_work to save memory
    1290         1276 :       CALL dbcsr_set(mat_work, 0.0_dp)
    1291         1276 :       CALL dbcsr_filter(mat_work, 1.0_dp)
    1292              : 
    1293         1276 :       CALL timestop(handle)
    1294              : 
    1295         1276 :    END SUBROUTINE contract_P_omega_with_mat_L
    1296              : 
    1297              : END MODULE rpa_util
        

Generated by: LCOV version 2.0-1