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

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines for low-scaling RPA/GW with imaginary time
      10              : !> \par History
      11              : !>      10.2015 created [Jan Wilhelm]
      12              : ! **************************************************************************************************
      13              : MODULE rpa_im_time
      14              :    USE cell_types,                      ONLY: cell_type,&
      15              :                                               get_cell
      16              :    USE cp_dbcsr_api,                    ONLY: &
      17              :         dbcsr_add, dbcsr_clear, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, &
      18              :         dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_init_p, dbcsr_p_type, &
      19              :         dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      20              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_reserve_all_blocks
      21              :    USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
      22              :                                               dbcsr_allocate_matrix_set,&
      23              :                                               dbcsr_deallocate_matrix_set
      24              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale,&
      25              :                                               cp_fm_scale
      26              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_type
      27              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      28              :                                               cp_fm_get_info,&
      29              :                                               cp_fm_release,&
      30              :                                               cp_fm_set_all,&
      31              :                                               cp_fm_to_fm,&
      32              :                                               cp_fm_type
      33              :    USE dbt_api,                         ONLY: &
      34              :         dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_contract, dbt_copy, &
      35              :         dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, dbt_filter, &
      36              :         dbt_get_info, dbt_nblks_total, dbt_nd_mp_comm, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
      37              :    USE hfx_types,                       ONLY: block_ind_type,&
      38              :                                               hfx_compression_type
      39              :    USE kinds,                           ONLY: dp,&
      40              :                                               int_8
      41              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      42              :                                               kpoint_env_type,&
      43              :                                               kpoint_type
      44              :    USE machine,                         ONLY: m_flush,&
      45              :                                               m_walltime
      46              :    USE mathconstants,                   ONLY: twopi
      47              :    USE message_passing,                 ONLY: mp_comm_type,&
      48              :                                               mp_para_env_type
      49              :    USE mp2_types,                       ONLY: mp2_type
      50              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      51              :    USE particle_types,                  ONLY: particle_type
      52              :    USE qs_environment_types,            ONLY: get_qs_env,&
      53              :                                               qs_environment_type
      54              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      55              :                                               mo_set_type
      56              :    USE qs_tensors,                      ONLY: decompress_tensor,&
      57              :                                               get_tensor_occupancy
      58              :    USE qs_tensors_types,                ONLY: create_2c_tensor
      59              :    USE rpa_gw_im_time_util,             ONLY: compute_weight_re_im,&
      60              :                                               get_atom_index_from_basis_function_index
      61              : #include "./base/base_uses.f90"
      62              : 
      63              :    IMPLICIT NONE
      64              : 
      65              :    PRIVATE
      66              : 
      67              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_im_time'
      68              : 
      69              :    PUBLIC :: compute_mat_P_omega, &
      70              :              compute_transl_dm, &
      71              :              init_cell_index_rpa, &
      72              :              zero_mat_P_omega, &
      73              :              compute_periodic_dm, &
      74              :              compute_mat_dm_global
      75              : 
      76              : CONTAINS
      77              : 
      78              : ! **************************************************************************************************
      79              : !> \brief ...
      80              : !> \param mat_P_omega ...
      81              : !> \param fm_scaled_dm_occ_tau ...
      82              : !> \param fm_scaled_dm_virt_tau ...
      83              : !> \param fm_mo_coeff_occ ...
      84              : !> \param fm_mo_coeff_virt ...
      85              : !> \param fm_mo_coeff_occ_scaled ...
      86              : !> \param fm_mo_coeff_virt_scaled ...
      87              : !> \param mat_P_global ...
      88              : !> \param matrix_s ...
      89              : !> \param ispin ...
      90              : !> \param t_3c_M ...
      91              : !> \param t_3c_O ...
      92              : !> \param t_3c_O_compressed ...
      93              : !> \param t_3c_O_ind ...
      94              : !> \param starts_array_mc ...
      95              : !> \param ends_array_mc ...
      96              : !> \param starts_array_mc_block ...
      97              : !> \param ends_array_mc_block ...
      98              : !> \param weights_cos_tf_t_to_w ...
      99              : !> \param tj ...
     100              : !> \param tau_tj ...
     101              : !> \param e_fermi ...
     102              : !> \param eps_filter ...
     103              : !> \param alpha ...
     104              : !> \param eps_filter_im_time ...
     105              : !> \param Eigenval ...
     106              : !> \param nmo ...
     107              : !> \param num_integ_points ...
     108              : !> \param cut_memory ...
     109              : !> \param unit_nr ...
     110              : !> \param mp2_env ...
     111              : !> \param para_env ...
     112              : !> \param qs_env ...
     113              : !> \param do_kpoints_from_Gamma ...
     114              : !> \param index_to_cell_3c ...
     115              : !> \param cell_to_index_3c ...
     116              : !> \param has_mat_P_blocks ...
     117              : !> \param do_ri_sos_laplace_mp2 ...
     118              : !> \param dbcsr_time ...
     119              : !> \param dbcsr_nflop ...
     120              : ! **************************************************************************************************
     121          522 :    SUBROUTINE compute_mat_P_omega(mat_P_omega, fm_scaled_dm_occ_tau, &
     122              :                                   fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
     123              :                                   fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
     124              :                                   mat_P_global, &
     125              :                                   matrix_s, &
     126              :                                   ispin, &
     127          348 :                                   t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
     128          174 :                                   starts_array_mc, ends_array_mc, &
     129          174 :                                   starts_array_mc_block, ends_array_mc_block, &
     130              :                                   weights_cos_tf_t_to_w, &
     131          174 :                                   tj, tau_tj, e_fermi, eps_filter, &
     132          174 :                                   alpha, eps_filter_im_time, Eigenval, nmo, &
     133              :                                   num_integ_points, cut_memory, unit_nr, &
     134              :                                   mp2_env, para_env, &
     135              :                                   qs_env, do_kpoints_from_Gamma, &
     136              :                                   index_to_cell_3c, cell_to_index_3c, &
     137          174 :                                   has_mat_P_blocks, do_ri_sos_laplace_mp2, &
     138              :                                   dbcsr_time, dbcsr_nflop)
     139              :       TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN)    :: mat_P_omega
     140              :       TYPE(cp_fm_type), INTENT(IN) :: fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
     141              :          fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled
     142              :       TYPE(dbcsr_p_type), INTENT(INOUT)                  :: mat_P_global
     143              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     144              :       INTEGER, INTENT(IN)                                :: ispin
     145              :       TYPE(dbt_type), INTENT(INOUT)                      :: t_3c_M
     146              :       TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT)     :: t_3c_O
     147              :       TYPE(hfx_compression_type), DIMENSION(:, :, :), &
     148              :          INTENT(INOUT)                                   :: t_3c_O_compressed
     149              :       TYPE(block_ind_type), DIMENSION(:, :, :), &
     150              :          INTENT(INOUT)                                   :: t_3c_O_ind
     151              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: starts_array_mc, ends_array_mc, &
     152              :                                                             starts_array_mc_block, &
     153              :                                                             ends_array_mc_block
     154              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
     155              :          INTENT(IN)                                      :: weights_cos_tf_t_to_w
     156              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     157              :          INTENT(IN)                                      :: tj
     158              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: tau_tj
     159              :       REAL(KIND=dp), INTENT(IN)                          :: e_fermi, eps_filter, alpha, &
     160              :                                                             eps_filter_im_time
     161              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval
     162              :       INTEGER, INTENT(IN)                                :: nmo, num_integ_points, cut_memory, &
     163              :                                                             unit_nr
     164              :       TYPE(mp2_type)                                     :: mp2_env
     165              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     166              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     167              :       LOGICAL, INTENT(IN)                                :: do_kpoints_from_Gamma
     168              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN)  :: index_to_cell_3c
     169              :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
     170              :          INTENT(IN)                                      :: cell_to_index_3c
     171              :       LOGICAL, DIMENSION(:, :, :, :, :), INTENT(INOUT)   :: has_mat_P_blocks
     172              :       LOGICAL, INTENT(IN)                                :: do_ri_sos_laplace_mp2
     173              :       REAL(dp), INTENT(INOUT)                            :: dbcsr_time
     174              :       INTEGER(int_8), INTENT(INOUT)                      :: dbcsr_nflop
     175              : 
     176              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_mat_P_omega'
     177              : 
     178              :       INTEGER :: comm_2d_handle, handle, handle2, handle3, i, i_cell, i_cell_R_1, &
     179              :          i_cell_R_1_minus_S, i_cell_R_1_minus_T, i_cell_R_2, i_cell_R_2_minus_S_minus_T, i_cell_S, &
     180              :          i_cell_T, i_mem, iquad, j, j_mem, jquad, num_3c_repl, num_cells_dm, unit_nr_dbcsr
     181              :       INTEGER(int_8)                                     :: nze, nze_dm_occ, nze_dm_virt, nze_M_occ, &
     182              :                                                             nze_M_virt, nze_O
     183              :       INTEGER(KIND=int_8)                                :: flops_1_occ, flops_1_virt, flops_2
     184          348 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: dist_1, dist_2, mc_ranges, size_dm, &
     185          174 :                                                             size_P
     186              :       INTEGER, DIMENSION(2)                              :: pdims_2d
     187              :       INTEGER, DIMENSION(2, 1)                           :: ibounds_2, jbounds_2
     188              :       INTEGER, DIMENSION(2, 2)                           :: ibounds_1, jbounds_1
     189              :       INTEGER, DIMENSION(3)                              :: bounds_3c
     190          174 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell_dm
     191              :       LOGICAL :: do_Gamma_RPA, do_kpoints_cubic_RPA, first_cycle_im_time, first_cycle_omega_loop, &
     192              :          memory_info, R_1_minus_S_needed, R_1_minus_T_needed, R_2_minus_S_minus_T_needed
     193              :       REAL(dp)                                           :: occ, occ_dm_occ, occ_dm_virt, occ_M_occ, &
     194              :                                                             occ_M_virt, occ_O, t1_flop
     195              :       REAL(KIND=dp)                                      :: omega, omega_old, t1, t2, tau, weight, &
     196              :                                                             weight_old
     197              :       TYPE(dbcsr_distribution_type)                      :: dist_P
     198          174 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_occ_global, mat_dm_virt_global
     199          522 :       TYPE(dbt_pgrid_type)                               :: pgrid_2d
     200         3306 :       TYPE(dbt_type)                                     :: t_3c_M_occ, t_3c_M_occ_tmp, t_3c_M_virt, &
     201         4872 :                                                             t_3c_M_virt_tmp, t_dm, t_dm_tmp, t_P, &
     202         1218 :                                                             t_P_tmp
     203          348 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:)          :: t_dm_occ, t_dm_virt
     204          174 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_O_occ, t_3c_O_virt
     205              :       TYPE(mp_comm_type)                                 :: comm_2d
     206              : 
     207          174 :       CALL timeset(routineN, handle)
     208              : 
     209          174 :       memory_info = mp2_env%ri_rpa_im_time%memory_info
     210          174 :       IF (memory_info) THEN
     211            0 :          unit_nr_dbcsr = unit_nr
     212              :       ELSE
     213          174 :          unit_nr_dbcsr = 0
     214              :       END IF
     215              : 
     216          174 :       do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
     217          174 :       do_Gamma_RPA = .NOT. do_kpoints_cubic_RPA
     218          744 :       num_3c_repl = MAXVAL(cell_to_index_3c)
     219              : 
     220          174 :       first_cycle_im_time = .TRUE.
     221         4052 :       ALLOCATE (t_3c_O_occ(SIZE(t_3c_O, 1), SIZE(t_3c_O, 2)), t_3c_O_virt(SIZE(t_3c_O, 1), SIZE(t_3c_O, 2)))
     222          364 :       DO i = 1, SIZE(t_3c_O, 1)
     223          634 :          DO j = 1, SIZE(t_3c_O, 2)
     224          270 :             CALL dbt_create(t_3c_O(i, j), t_3c_O_occ(i, j))
     225          460 :             CALL dbt_create(t_3c_O(i, j), t_3c_O_virt(i, j))
     226              :          END DO
     227              :       END DO
     228              : 
     229          174 :       CALL dbt_create(t_3c_M, t_3c_M_occ, name="M occ (RI | AO AO)")
     230          174 :       CALL dbt_create(t_3c_M, t_3c_M_virt, name="M virt (RI | AO AO)")
     231              : 
     232          522 :       ALLOCATE (mc_ranges(cut_memory + 1))
     233          502 :       mc_ranges(:cut_memory) = starts_array_mc_block(:)
     234          174 :       mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
     235              : 
     236         1582 :       DO jquad = 1, num_integ_points
     237              : 
     238         1408 :          CALL para_env%sync()
     239         1408 :          t1 = m_walltime()
     240              : 
     241              :          CALL compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, &
     242              :                                     fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
     243              :                                     fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, &
     244              :                                     matrix_s, ispin, &
     245              :                                     Eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
     246              :                                     jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, &
     247         1408 :                                     num_cells_dm, index_to_cell_dm, para_env)
     248              : 
     249        14128 :          ALLOCATE (t_dm_virt(num_cells_dm))
     250        12720 :          ALLOCATE (t_dm_occ(num_cells_dm))
     251         1408 :          CALL dbcsr_get_info(mat_P_global%matrix, distribution=dist_P)
     252         1408 :          CALL dbcsr_distribution_get(dist_P, group=comm_2d_handle, nprows=pdims_2d(1), npcols=pdims_2d(2))
     253         1408 :          CALL comm_2d%set_handle(comm_2d_handle)
     254              : 
     255         1408 :          pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
     256         4224 :          ALLOCATE (size_P(dbt_nblks_total(t_3c_M, 1)))
     257         1408 :          CALL dbt_get_info(t_3c_M, blk_size_1=size_P)
     258              : 
     259         4224 :          ALLOCATE (size_dm(dbt_nblks_total(t_3c_O(1, 1), 3)))
     260         1408 :          CALL dbt_get_info(t_3c_O(1, 1), blk_size_3=size_dm)
     261         1408 :          CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name="D (AO | AO)")
     262         1408 :          DEALLOCATE (size_dm)
     263         1408 :          DEALLOCATE (dist_1, dist_2)
     264              :          CALL create_2c_tensor(t_P, dist_1, dist_2, pgrid_2d, size_P, size_P, name="P (RI | RI)")
     265         1408 :          DEALLOCATE (size_P)
     266         1408 :          DEALLOCATE (dist_1, dist_2)
     267         1408 :          CALL dbt_pgrid_destroy(pgrid_2d)
     268              : 
     269         2864 :          DO i_cell = 1, num_cells_dm
     270         1456 :             CALL dbt_create(t_dm, t_dm_virt(i_cell), name="D virt (AO | AO)")
     271         1456 :             CALL dbt_create(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
     272         1456 :             CALL dbt_copy_matrix_to_tensor(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
     273         1456 :             CALL dbt_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.TRUE.)
     274         1456 :             CALL dbcsr_clear(mat_dm_virt_global(jquad, i_cell)%matrix)
     275              : 
     276         1456 :             CALL dbt_create(t_dm, t_dm_occ(i_cell), name="D occ (AO | AO)")
     277         1456 :             CALL dbt_copy_matrix_to_tensor(mat_dm_occ_global(jquad, i_cell)%matrix, t_dm_tmp)
     278         1456 :             CALL dbt_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.TRUE.)
     279         1456 :             CALL dbt_destroy(t_dm_tmp)
     280         2864 :             CALL dbcsr_clear(mat_dm_occ_global(jquad, i_cell)%matrix)
     281              :          END DO
     282              : 
     283         1408 :          CALL get_tensor_occupancy(t_dm_occ(1), nze_dm_occ, occ_dm_occ)
     284         1408 :          CALL get_tensor_occupancy(t_dm_virt(1), nze_dm_virt, occ_dm_virt)
     285              : 
     286         1408 :          CALL dbt_destroy(t_dm)
     287              : 
     288         1408 :          CALL dbt_create(t_3c_O_occ(1, 1), t_3c_M_occ_tmp, name="M (RI AO | AO)")
     289         1408 :          CALL dbt_create(t_3c_O_virt(1, 1), t_3c_M_virt_tmp, name="M (RI AO | AO)")
     290              : 
     291         1408 :          CALL timeset(routineN//"_contract", handle2)
     292              : 
     293         1408 :          CALL para_env%sync()
     294         1408 :          t1_flop = m_walltime()
     295              : 
     296         2912 :          DO i = 1, SIZE(t_3c_O_occ, 1)
     297         4896 :             DO j = 1, SIZE(t_3c_O_occ, 2)
     298         3488 :                CALL dbt_batched_contract_init(t_3c_O_occ(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
     299              :             END DO
     300              :          END DO
     301         2912 :          DO i = 1, SIZE(t_3c_O_virt, 1)
     302         4896 :             DO j = 1, SIZE(t_3c_O_virt, 2)
     303         3488 :                CALL dbt_batched_contract_init(t_3c_O_virt(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
     304              :             END DO
     305              :          END DO
     306         1408 :          CALL dbt_batched_contract_init(t_3c_M_occ_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
     307         1408 :          CALL dbt_batched_contract_init(t_3c_M_virt_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
     308         1408 :          CALL dbt_batched_contract_init(t_3c_M_occ, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
     309         1408 :          CALL dbt_batched_contract_init(t_3c_M_virt, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
     310              : 
     311         2840 :          DO i_cell_T = 1, num_cells_dm/2 + 1
     312              : 
     313         1432 :             IF (.NOT. ANY(has_mat_P_blocks(i_cell_T, :, :, :, :))) CYCLE
     314              : 
     315         1432 :             CALL dbt_batched_contract_init(t_P)
     316              : 
     317         1432 :             IF (do_Gamma_RPA) THEN
     318         1384 :                nze_O = 0
     319         1384 :                nze_M_virt = 0
     320         1384 :                nze_M_occ = 0
     321         1384 :                occ_M_virt = 0.0_dp
     322         1384 :                occ_M_occ = 0.0_dp
     323         1384 :                occ_O = 0.0_dp
     324              :             END IF
     325              : 
     326         3576 :             DO j_mem = 1, cut_memory
     327              : 
     328         2144 :                CALL dbt_get_info(t_3c_O_occ(1, 1), nfull_total=bounds_3c)
     329              : 
     330         6432 :                jbounds_1(:, 1) = [1, bounds_3c(1)]
     331         6432 :                jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
     332              : 
     333         6432 :                jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
     334              : 
     335         2144 :                IF (do_Gamma_RPA) CALL dbt_batched_contract_init(t_dm_virt(1))
     336              : 
     337         6552 :                DO i_mem = 1, cut_memory
     338              : 
     339         4568 :                   IF (.NOT. ANY(has_mat_P_blocks(i_cell_T, i_mem, j_mem, :, :))) CYCLE
     340              : 
     341        12984 :                   ibounds_1(:, 1) = [1, bounds_3c(1)]
     342        12984 :                   ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
     343              : 
     344        12984 :                   ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
     345              : 
     346         4328 :                   IF (unit_nr_dbcsr > 0) WRITE (UNIT=unit_nr_dbcsr, FMT="(T3,A,I3,1X,I3)") &
     347            0 :                      "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem
     348              : 
     349        11280 :                   DO i_cell_R_1 = 1, num_3c_repl
     350              : 
     351        16424 :                      DO i_cell_R_2 = 1, num_3c_repl
     352              : 
     353         7208 :                         IF (.NOT. has_mat_P_blocks(i_cell_T, i_mem, j_mem, i_cell_R_1, i_cell_R_2)) CYCLE
     354              : 
     355              :                         CALL get_diff_index_3c(i_cell_R_1, i_cell_T, i_cell_R_1_minus_T, &
     356              :                                                index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
     357         5258 :                                                R_1_minus_T_needed, do_kpoints_cubic_RPA)
     358              : 
     359         5258 :                         IF (do_Gamma_RPA) CALL dbt_batched_contract_init(t_dm_occ(1))
     360        12616 :                         DO i_cell_S = 1, num_cells_dm
     361              :                            CALL get_diff_index_3c(i_cell_R_1, i_cell_S, i_cell_R_1_minus_S, index_to_cell_3c, &
     362              :                                                   cell_to_index_3c, index_to_cell_dm, R_1_minus_S_needed, &
     363         7358 :                                                   do_kpoints_cubic_RPA)
     364        12616 :                            IF (R_1_minus_S_needed) THEN
     365              : 
     366         7158 :                               CALL timeset(routineN//"_calc_M_occ_t", handle3)
     367              :                               CALL decompress_tensor(t_3c_O(i_cell_R_1_minus_S, i_cell_R_2), &
     368              :                                                      t_3c_O_ind(i_cell_R_1_minus_S, i_cell_R_2, j_mem)%ind, &
     369              :                                                      t_3c_O_compressed(i_cell_R_1_minus_S, i_cell_R_2, j_mem), &
     370         7158 :                                                      qs_env%mp2_env%ri_rpa_im_time%eps_compress)
     371              : 
     372         7158 :                               IF (do_Gamma_RPA .AND. i_mem == 1) THEN
     373         2052 :                                  CALL get_tensor_occupancy(t_3c_O(1, 1), nze, occ)
     374         2052 :                                  nze_O = nze_O + nze
     375         2052 :                                  occ_O = occ_O + occ
     376              :                               END IF
     377              : 
     378              :                               CALL dbt_copy(t_3c_O(i_cell_R_1_minus_S, i_cell_R_2), &
     379         7158 :                                             t_3c_O_occ(i_cell_R_1_minus_S, i_cell_R_2), move_data=.TRUE.)
     380              : 
     381              :                               CALL dbt_contract(alpha=1.0_dp, &
     382              :                                                 tensor_1=t_3c_O_occ(i_cell_R_1_minus_S, i_cell_R_2), &
     383              :                                                 tensor_2=t_dm_occ(i_cell_S), &
     384              :                                                 beta=1.0_dp, &
     385              :                                                 tensor_3=t_3c_M_occ_tmp, &
     386              :                                                 contract_1=[3], notcontract_1=[1, 2], &
     387              :                                                 contract_2=[2], notcontract_2=[1], &
     388              :                                                 map_1=[1, 2], map_2=[3], &
     389              :                                                 bounds_2=jbounds_1, bounds_3=ibounds_2, &
     390              :                                                 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
     391         7158 :                                                 flop=flops_1_occ)
     392         7158 :                               CALL timestop(handle3)
     393              : 
     394         7158 :                               dbcsr_nflop = dbcsr_nflop + flops_1_occ
     395              : 
     396              :                            END IF
     397              :                         END DO
     398              : 
     399         5258 :                         IF (do_Gamma_RPA) CALL dbt_batched_contract_finalize(t_dm_occ(1))
     400              : 
     401              :                         ! copy matrix to optimal contraction layout - copy is done manually in order
     402              :                         ! to better control memory allocations (we can release data of previous
     403              :                         ! representation)
     404         5258 :                         CALL timeset(routineN//"_copy_M_occ_t", handle3)
     405         5258 :                         CALL dbt_copy(t_3c_M_occ_tmp, t_3c_M_occ, order=[1, 3, 2], move_data=.TRUE.)
     406         5258 :                         CALL dbt_filter(t_3c_M_occ, eps_filter)
     407         5258 :                         CALL timestop(handle3)
     408              : 
     409         5258 :                         IF (do_Gamma_RPA) THEN
     410         4208 :                            CALL get_tensor_occupancy(t_3c_M_occ, nze, occ)
     411         4208 :                            nze_M_occ = nze_M_occ + nze
     412         4208 :                            occ_M_occ = occ_M_occ + occ
     413              :                         END IF
     414              : 
     415        12616 :                         DO i_cell_S = 1, num_cells_dm
     416              :                            CALL get_diff_diff_index_3c(i_cell_R_2, i_cell_S, i_cell_T, i_cell_R_2_minus_S_minus_T, &
     417              :                                                        index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
     418         7358 :                                                        R_2_minus_S_minus_T_needed, do_kpoints_cubic_RPA)
     419              : 
     420        12616 :                            IF (R_1_minus_T_needed .AND. R_2_minus_S_minus_T_needed) THEN
     421              :                               CALL decompress_tensor(t_3c_O(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), &
     422              :                                                      t_3c_O_ind(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T, i_mem)%ind, &
     423              :                                                      t_3c_O_compressed(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T, i_mem), &
     424         6988 :                                                      qs_env%mp2_env%ri_rpa_im_time%eps_compress)
     425              : 
     426              :                               CALL dbt_copy(t_3c_O(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), &
     427         6988 :                                             t_3c_O_virt(i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), move_data=.TRUE.)
     428              : 
     429         6988 :                               CALL timeset(routineN//"_calc_M_virt_t", handle3)
     430              :                               CALL dbt_contract(alpha=alpha/2.0_dp, &
     431              :                                                 tensor_1=t_3c_O_virt( &
     432              :                                                 i_cell_R_2_minus_S_minus_T, i_cell_R_1_minus_T), &
     433              :                                                 tensor_2=t_dm_virt(i_cell_S), &
     434              :                                                 beta=1.0_dp, &
     435              :                                                 tensor_3=t_3c_M_virt_tmp, &
     436              :                                                 contract_1=[3], notcontract_1=[1, 2], &
     437              :                                                 contract_2=[2], notcontract_2=[1], &
     438              :                                                 map_1=[1, 2], map_2=[3], &
     439              :                                                 bounds_2=ibounds_1, bounds_3=jbounds_2, &
     440              :                                                 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
     441         6988 :                                                 flop=flops_1_virt)
     442         6988 :                               CALL timestop(handle3)
     443              : 
     444         6988 :                               dbcsr_nflop = dbcsr_nflop + flops_1_virt
     445              : 
     446              :                            END IF
     447              :                         END DO
     448              : 
     449         5258 :                         CALL timeset(routineN//"_copy_M_virt_t", handle3)
     450         5258 :                         CALL dbt_copy(t_3c_M_virt_tmp, t_3c_M_virt, move_data=.TRUE.)
     451         5258 :                         CALL dbt_filter(t_3c_M_virt, eps_filter)
     452         5258 :                         CALL timestop(handle3)
     453              : 
     454         5258 :                         IF (do_Gamma_RPA) THEN
     455         4208 :                            CALL get_tensor_occupancy(t_3c_M_virt, nze, occ)
     456         4208 :                            nze_M_virt = nze_M_virt + nze
     457         4208 :                            occ_M_virt = occ_M_virt + occ
     458              :                         END IF
     459              : 
     460              :                         flops_2 = 0
     461              : 
     462         5258 :                         CALL timeset(routineN//"_calc_P_t", handle3)
     463              : 
     464              :                         CALL dbt_contract(alpha=1.0_dp, tensor_1=t_3c_M_occ, &
     465              :                                           tensor_2=t_3c_M_virt, &
     466              :                                           beta=1.0_dp, &
     467              :                                           tensor_3=t_P, &
     468              :                                           contract_1=[2, 3], notcontract_1=[1], &
     469              :                                           contract_2=[2, 3], notcontract_2=[1], &
     470              :                                           map_1=[1], map_2=[2], &
     471              :                                           filter_eps=eps_filter_im_time/REAL(cut_memory**2, KIND=dp), &
     472              :                                           flop=flops_2, &
     473              :                                           move_data=.TRUE., &
     474         5258 :                                           unit_nr=unit_nr_dbcsr)
     475              : 
     476         5258 :                         CALL timestop(handle3)
     477              : 
     478         5258 :                         first_cycle_im_time = .FALSE.
     479              : 
     480         5258 :                         IF (jquad == 1 .AND. flops_2 == 0) &
     481        26246 :                            has_mat_P_blocks(i_cell_T, i_mem, j_mem, i_cell_R_1, i_cell_R_2) = .FALSE.
     482              : 
     483              :                      END DO
     484              :                   END DO
     485              :                END DO
     486         3576 :                IF (do_Gamma_RPA) CALL dbt_batched_contract_finalize(t_dm_virt(1))
     487              :             END DO
     488              : 
     489         1432 :             CALL dbt_batched_contract_finalize(t_P, unit_nr=unit_nr_dbcsr)
     490              : 
     491         1432 :             CALL dbt_create(mat_P_global%matrix, t_P_tmp)
     492         1432 :             CALL dbt_copy(t_P, t_P_tmp, move_data=.TRUE.)
     493         1432 :             CALL dbt_copy_tensor_to_matrix(t_P_tmp, mat_P_global%matrix)
     494         1432 :             CALL dbt_destroy(t_P_tmp)
     495              : 
     496         2840 :             IF (do_ri_sos_laplace_mp2) THEN
     497              :                ! For RI-SOS-Laplace-MP2 we do not perform a cosine transform,
     498              :                ! but we have to copy P_local to the output matrix
     499              : 
     500          138 :                CALL dbcsr_add(mat_P_omega(jquad, i_cell_T)%matrix, mat_P_global%matrix, 1.0_dp, 1.0_dp)
     501              :             ELSE
     502         1294 :                CALL timeset(routineN//"_Fourier_transform", handle3)
     503              : 
     504              :                ! Fourier transform of P(it) to P(iw)
     505         1294 :                first_cycle_omega_loop = .TRUE.
     506              : 
     507         1294 :                tau = tau_tj(jquad)
     508              : 
     509        24308 :                DO iquad = 1, num_integ_points
     510              : 
     511        23014 :                   omega = tj(iquad)
     512        23014 :                   weight = weights_cos_tf_t_to_w(iquad, jquad)
     513              : 
     514        23014 :                   IF (first_cycle_omega_loop) THEN
     515              :                      ! no multiplication with 2.0 as in Kresses paper (Kaltak, JCTC 10, 2498 (2014), Eq. 12)
     516              :                      ! because this factor is already absorbed in the weight w_j
     517         1294 :                      CALL dbcsr_scale(mat_P_global%matrix, COS(omega*tau)*weight)
     518              :                   ELSE
     519        21720 :                      CALL dbcsr_scale(mat_P_global%matrix, COS(omega*tau)/COS(omega_old*tau)*weight/weight_old)
     520              :                   END IF
     521              : 
     522        23014 :                   CALL dbcsr_add(mat_P_omega(iquad, i_cell_T)%matrix, mat_P_global%matrix, 1.0_dp, 1.0_dp)
     523              : 
     524        23014 :                   first_cycle_omega_loop = .FALSE.
     525              : 
     526        23014 :                   omega_old = omega
     527        24308 :                   weight_old = weight
     528              : 
     529              :                END DO
     530              : 
     531         1294 :                CALL timestop(handle3)
     532              :             END IF
     533              : 
     534              :          END DO
     535              : 
     536         1408 :          CALL timestop(handle2)
     537              : 
     538         1408 :          CALL dbt_batched_contract_finalize(t_3c_M_occ_tmp)
     539         1408 :          CALL dbt_batched_contract_finalize(t_3c_M_virt_tmp)
     540         1408 :          CALL dbt_batched_contract_finalize(t_3c_M_occ)
     541         1408 :          CALL dbt_batched_contract_finalize(t_3c_M_virt)
     542              : 
     543         2912 :          DO i = 1, SIZE(t_3c_O_occ, 1)
     544         4896 :             DO j = 1, SIZE(t_3c_O_occ, 2)
     545         3488 :                CALL dbt_batched_contract_finalize(t_3c_O_occ(i, j))
     546              :             END DO
     547              :          END DO
     548              : 
     549         2912 :          DO i = 1, SIZE(t_3c_O_virt, 1)
     550         4896 :             DO j = 1, SIZE(t_3c_O_virt, 2)
     551         3488 :                CALL dbt_batched_contract_finalize(t_3c_O_virt(i, j))
     552              :             END DO
     553              :          END DO
     554              : 
     555         1408 :          CALL dbt_destroy(t_P)
     556         2864 :          DO i_cell = 1, num_cells_dm
     557         1456 :             CALL dbt_destroy(t_dm_virt(i_cell))
     558         2864 :             CALL dbt_destroy(t_dm_occ(i_cell))
     559              :          END DO
     560              : 
     561         1408 :          CALL dbt_destroy(t_3c_M_occ_tmp)
     562         1408 :          CALL dbt_destroy(t_3c_M_virt_tmp)
     563         2864 :          DEALLOCATE (t_dm_virt)
     564         2864 :          DEALLOCATE (t_dm_occ)
     565              : 
     566         1408 :          CALL para_env%sync()
     567         1408 :          t2 = m_walltime()
     568              : 
     569         1408 :          dbcsr_time = dbcsr_time + t2 - t1_flop
     570              : 
     571         8622 :          IF (unit_nr > 0) THEN
     572              :             WRITE (unit_nr, '(/T3,A,1X,I3)') &
     573          704 :                'RPA_LOW_SCALING_INFO| Info for time point', jquad
     574              :             WRITE (unit_nr, '(T6,A,T56,F25.1)') &
     575          704 :                'Execution time (s):', t2 - t1
     576              :             WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
     577          704 :                'Occupancy of D occ:', REAL(nze_dm_occ, dp), '/', occ_dm_occ*100, '%'
     578              :             WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
     579          704 :                'Occupancy of D virt:', REAL(nze_dm_virt, dp), '/', occ_dm_virt*100, '%'
     580          704 :             IF (do_Gamma_RPA) THEN
     581              :                WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
     582          692 :                   'Occupancy of 3c ints:', REAL(nze_O, dp), '/', occ_O*100, '%'
     583              :                WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
     584          692 :                   'Occupancy of M occ:', REAL(nze_M_occ, dp), '/', occ_M_occ*100, '%'
     585              :                WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
     586          692 :                   'Occupancy of M virt:', REAL(nze_M_virt, dp), '/', occ_M_virt*100, '%'
     587              :             END IF
     588          704 :             WRITE (unit_nr, *)
     589          704 :             CALL m_flush(unit_nr)
     590              :          END IF
     591              : 
     592              :       END DO ! time points
     593              : 
     594          174 :       CALL dbt_destroy(t_3c_M_occ)
     595          174 :       CALL dbt_destroy(t_3c_M_virt)
     596              : 
     597          364 :       DO i = 1, SIZE(t_3c_O, 1)
     598          634 :          DO j = 1, SIZE(t_3c_O, 2)
     599          270 :             CALL dbt_destroy(t_3c_O_occ(i, j))
     600          460 :             CALL dbt_destroy(t_3c_O_virt(i, j))
     601              :          END DO
     602              :       END DO
     603              : 
     604          174 :       CALL clean_up(mat_dm_occ_global, mat_dm_virt_global)
     605              : 
     606          174 :       CALL timestop(handle)
     607              : 
     608         1062 :    END SUBROUTINE compute_mat_P_omega
     609              : 
     610              : ! **************************************************************************************************
     611              : !> \brief ...
     612              : !> \param mat_P_omega ...
     613              : ! **************************************************************************************************
     614          174 :    SUBROUTINE zero_mat_P_omega(mat_P_omega)
     615              :       TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN)    :: mat_P_omega
     616              : 
     617              :       INTEGER                                            :: i_kp, jquad
     618              : 
     619         1582 :       DO jquad = 1, SIZE(mat_P_omega, 1)
     620         5714 :          DO i_kp = 1, SIZE(mat_P_omega, 2)
     621              : 
     622         5540 :             CALL dbcsr_set(mat_P_omega(jquad, i_kp)%matrix, 0.0_dp)
     623              : 
     624              :          END DO
     625              :       END DO
     626              : 
     627          174 :    END SUBROUTINE zero_mat_P_omega
     628              : 
     629              : ! **************************************************************************************************
     630              : !> \brief ...
     631              : !> \param fm_scaled_dm_occ_tau ...
     632              : !> \param fm_scaled_dm_virt_tau ...
     633              : !> \param tau_tj ...
     634              : !> \param num_integ_points ...
     635              : !> \param nmo ...
     636              : !> \param fm_mo_coeff_occ ...
     637              : !> \param fm_mo_coeff_virt ...
     638              : !> \param fm_mo_coeff_occ_scaled ...
     639              : !> \param fm_mo_coeff_virt_scaled ...
     640              : !> \param mat_dm_occ_global ...
     641              : !> \param mat_dm_virt_global ...
     642              : !> \param matrix_s ...
     643              : !> \param ispin ...
     644              : !> \param Eigenval ...
     645              : !> \param e_fermi ...
     646              : !> \param eps_filter ...
     647              : !> \param memory_info ...
     648              : !> \param unit_nr ...
     649              : !> \param jquad ...
     650              : !> \param do_kpoints_cubic_RPA ...
     651              : !> \param do_kpoints_from_Gamma ...
     652              : !> \param qs_env ...
     653              : !> \param num_cells_dm ...
     654              : !> \param index_to_cell_dm ...
     655              : !> \param para_env ...
     656              : ! **************************************************************************************************
     657         1578 :    SUBROUTINE compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, &
     658              :                                     fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
     659              :                                     fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, &
     660         1578 :                                     matrix_s, ispin, &
     661         1578 :                                     Eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
     662              :                                     jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, &
     663              :                                     num_cells_dm, index_to_cell_dm, para_env)
     664              : 
     665              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_scaled_dm_occ_tau, &
     666              :                                                             fm_scaled_dm_virt_tau
     667              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: tau_tj
     668              :       INTEGER, INTENT(IN)                                :: num_integ_points, nmo
     669              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_mo_coeff_occ, fm_mo_coeff_virt, &
     670              :                                                             fm_mo_coeff_occ_scaled, &
     671              :                                                             fm_mo_coeff_virt_scaled
     672              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_occ_global, mat_dm_virt_global
     673              :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN)       :: matrix_s
     674              :       INTEGER, INTENT(IN)                                :: ispin
     675              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval
     676              :       REAL(KIND=dp), INTENT(IN)                          :: e_fermi, eps_filter
     677              :       LOGICAL, INTENT(IN)                                :: memory_info
     678              :       INTEGER, INTENT(IN)                                :: unit_nr, jquad
     679              :       LOGICAL, INTENT(IN)                                :: do_kpoints_cubic_RPA, &
     680              :                                                             do_kpoints_from_Gamma
     681              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     682              :       INTEGER, INTENT(OUT)                               :: num_cells_dm
     683              :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell_dm
     684              :       TYPE(mp_para_env_type), INTENT(IN)                 :: para_env
     685              : 
     686              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_mat_dm_global'
     687              :       REAL(KIND=dp), PARAMETER                           :: stabilize_exp = 70.0_dp
     688              : 
     689              :       INTEGER                                            :: handle, i_global, iiB, iquad, jjB, &
     690              :                                                             ncol_local, nrow_local, size_dm_occ, &
     691              :                                                             size_dm_virt
     692         1578 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     693              :       REAL(KIND=dp)                                      :: tau
     694              : 
     695         1578 :       CALL timeset(routineN, handle)
     696              : 
     697         1578 :       IF (memory_info .AND. unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
     698            0 :          "RPA_LOW_SCALING_INFO| Started with time point: ", jquad
     699              : 
     700         1578 :       tau = tau_tj(jquad)
     701              : 
     702         1578 :       IF (do_kpoints_cubic_RPA) THEN
     703              : 
     704              :          CALL compute_transl_dm(mat_dm_occ_global, qs_env, &
     705              :                                 ispin, num_integ_points, jquad, e_fermi, tau, &
     706              :                                 eps_filter, num_cells_dm, index_to_cell_dm, &
     707           24 :                                 remove_occ=.FALSE., remove_virt=.TRUE., first_jquad=1)
     708              : 
     709              :          CALL compute_transl_dm(mat_dm_virt_global, qs_env, &
     710              :                                 ispin, num_integ_points, jquad, e_fermi, tau, &
     711              :                                 eps_filter, num_cells_dm, index_to_cell_dm, &
     712           24 :                                 remove_occ=.TRUE., remove_virt=.FALSE., first_jquad=1)
     713              : 
     714         1554 :       ELSE IF (do_kpoints_from_Gamma) THEN
     715              : 
     716              :          CALL compute_periodic_dm(mat_dm_occ_global, qs_env, &
     717              :                                   ispin, num_integ_points, jquad, e_fermi, tau, &
     718              :                                   remove_occ=.FALSE., remove_virt=.TRUE., &
     719          108 :                                   alloc_dm=(jquad == 1))
     720              : 
     721              :          CALL compute_periodic_dm(mat_dm_virt_global, qs_env, &
     722              :                                   ispin, num_integ_points, jquad, e_fermi, tau, &
     723              :                                   remove_occ=.TRUE., remove_virt=.FALSE., &
     724          108 :                                   alloc_dm=(jquad == 1))
     725              : 
     726          108 :          num_cells_dm = 1
     727              : 
     728              :       ELSE
     729              : 
     730         1446 :          num_cells_dm = 1
     731              : 
     732         1446 :          CALL para_env%sync()
     733              : 
     734              :          ! get info of fm_mo_coeff_occ
     735              :          CALL cp_fm_get_info(matrix=fm_mo_coeff_occ, &
     736              :                              nrow_local=nrow_local, &
     737              :                              ncol_local=ncol_local, &
     738              :                              row_indices=row_indices, &
     739         1446 :                              col_indices=col_indices)
     740              : 
     741              :          ! Multiply the occupied and the virtual MO coefficients with the factor exp((-e_i-e_F)*tau/2).
     742              :          ! Then, we simply get the sum over all occ states and virt. states by a simple matrix-matrix
     743              :          ! multiplication.
     744              : 
     745              :          ! first, the occ
     746        19689 :          DO jjB = 1, nrow_local
     747       553018 :             DO iiB = 1, ncol_local
     748       533329 :                i_global = col_indices(iiB)
     749              : 
     750              :                ! hard coded: exponential function gets NaN if argument is negative with large absolute value
     751              :                ! use 69, since e^(-69) = 10^(-30) which should be sufficiently small that it does not matter
     752       551572 :                IF (ABS(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) < stabilize_exp) THEN
     753              :                   fm_mo_coeff_occ_scaled%local_data(jjB, iiB) = &
     754       434957 :                      fm_mo_coeff_occ%local_data(jjB, iiB)*EXP(tau*0.5_dp*(Eigenval(i_global) - e_fermi))
     755              :                ELSE
     756        98372 :                   fm_mo_coeff_occ_scaled%local_data(jjB, iiB) = 0.0_dp
     757              :                END IF
     758              : 
     759              :             END DO
     760              :          END DO
     761              : 
     762              :          ! get info of fm_mo_coeff_virt
     763              :          CALL cp_fm_get_info(matrix=fm_mo_coeff_virt, &
     764              :                              nrow_local=nrow_local, &
     765              :                              ncol_local=ncol_local, &
     766              :                              row_indices=row_indices, &
     767         1446 :                              col_indices=col_indices)
     768              : 
     769              :          ! the same for virt
     770        19689 :          DO jjB = 1, nrow_local
     771       553018 :             DO iiB = 1, ncol_local
     772       533329 :                i_global = col_indices(iiB)
     773              : 
     774       551572 :                IF (ABS(tau*0.5_dp*(Eigenval(i_global) - e_fermi)) < stabilize_exp) THEN
     775              :                   fm_mo_coeff_virt_scaled%local_data(jjB, iiB) = &
     776       434957 :                      fm_mo_coeff_virt%local_data(jjB, iiB)*EXP(-tau*0.5_dp*(Eigenval(i_global) - e_fermi))
     777              :                ELSE
     778        98372 :                   fm_mo_coeff_virt_scaled%local_data(jjB, iiB) = 0.0_dp
     779              :                END IF
     780              : 
     781              :             END DO
     782              :          END DO
     783              : 
     784         1446 :          CALL para_env%sync()
     785              : 
     786         1446 :          size_dm_occ = nmo
     787         1446 :          size_dm_virt = nmo
     788              : 
     789              :          CALL parallel_gemm(transa="N", transb="T", m=size_dm_occ, n=size_dm_occ, k=nmo, alpha=1.0_dp, &
     790              :                             matrix_a=fm_mo_coeff_occ_scaled, matrix_b=fm_mo_coeff_occ_scaled, beta=0.0_dp, &
     791         1446 :                             matrix_c=fm_scaled_dm_occ_tau)
     792              : 
     793              :          CALL parallel_gemm(transa="N", transb="T", m=size_dm_virt, n=size_dm_virt, k=nmo, alpha=1.0_dp, &
     794              :                             matrix_a=fm_mo_coeff_virt_scaled, matrix_b=fm_mo_coeff_virt_scaled, beta=0.0_dp, &
     795         1446 :                             matrix_c=fm_scaled_dm_virt_tau)
     796              : 
     797         1446 :          IF (jquad == 1) THEN
     798              : 
     799              :             ! transfer fm density matrices to dbcsr matrix
     800          214 :             NULLIFY (mat_dm_occ_global)
     801          214 :             CALL dbcsr_allocate_matrix_set(mat_dm_occ_global, num_integ_points, 1)
     802              : 
     803         1660 :             DO iquad = 1, num_integ_points
     804              : 
     805         1446 :                ALLOCATE (mat_dm_occ_global(iquad, 1)%matrix)
     806              :                CALL dbcsr_create(matrix=mat_dm_occ_global(iquad, 1)%matrix, &
     807              :                                  template=matrix_s(1)%matrix, &
     808         1660 :                                  matrix_type=dbcsr_type_no_symmetry)
     809              : 
     810              :             END DO
     811              : 
     812              :          END IF
     813              : 
     814              :          CALL copy_fm_to_dbcsr(fm_scaled_dm_occ_tau, &
     815              :                                mat_dm_occ_global(jquad, 1)%matrix, &
     816         1446 :                                keep_sparsity=.FALSE.)
     817              : 
     818         1446 :          CALL dbcsr_filter(mat_dm_occ_global(jquad, 1)%matrix, eps_filter)
     819              : 
     820         1446 :          IF (jquad == 1) THEN
     821              : 
     822          214 :             NULLIFY (mat_dm_virt_global)
     823          214 :             CALL dbcsr_allocate_matrix_set(mat_dm_virt_global, num_integ_points, 1)
     824              : 
     825              :          END IF
     826              : 
     827         1446 :          ALLOCATE (mat_dm_virt_global(jquad, 1)%matrix)
     828              :          CALL dbcsr_create(matrix=mat_dm_virt_global(jquad, 1)%matrix, &
     829              :                            template=matrix_s(1)%matrix, &
     830         1446 :                            matrix_type=dbcsr_type_no_symmetry)
     831              :          CALL copy_fm_to_dbcsr(fm_scaled_dm_virt_tau, &
     832              :                                mat_dm_virt_global(jquad, 1)%matrix, &
     833         1446 :                                keep_sparsity=.FALSE.)
     834              : 
     835         1446 :          CALL dbcsr_filter(mat_dm_virt_global(jquad, 1)%matrix, eps_filter)
     836              : 
     837              :          ! release memory
     838         1446 :          IF (jquad > 1) THEN
     839         1232 :             CALL dbcsr_set(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
     840         1232 :             CALL dbcsr_set(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp)
     841         1232 :             CALL dbcsr_filter(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
     842         1232 :             CALL dbcsr_filter(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp)
     843              :          END IF
     844              : 
     845              :       END IF ! do kpoints
     846              : 
     847         1578 :       CALL timestop(handle)
     848              : 
     849         1578 :    END SUBROUTINE compute_mat_dm_global
     850              : 
     851              : ! **************************************************************************************************
     852              : !> \brief ...
     853              : !> \param mat_dm_occ_global ...
     854              : !> \param mat_dm_virt_global ...
     855              : ! **************************************************************************************************
     856          174 :    SUBROUTINE clean_up(mat_dm_occ_global, mat_dm_virt_global)
     857              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_occ_global, mat_dm_virt_global
     858              : 
     859          174 :       CALL dbcsr_deallocate_matrix_set(mat_dm_occ_global)
     860          174 :       CALL dbcsr_deallocate_matrix_set(mat_dm_virt_global)
     861              : 
     862          174 :    END SUBROUTINE clean_up
     863              : 
     864              : ! **************************************************************************************************
     865              : !> \brief Calculate kpoint density matrices (rho(k), owned by kpoint groups)
     866              : !> \param kpoint    kpoint environment
     867              : !> \param tau ...
     868              : !> \param e_fermi ...
     869              : !> \param remove_occ ...
     870              : !> \param remove_virt ...
     871              : ! **************************************************************************************************
     872          508 :    SUBROUTINE kpoint_density_matrices_rpa(kpoint, tau, e_fermi, remove_occ, remove_virt)
     873              : 
     874              :       TYPE(kpoint_type), POINTER                         :: kpoint
     875              :       REAL(KIND=dp), INTENT(IN)                          :: tau, e_fermi
     876              :       LOGICAL, INTENT(IN)                                :: remove_occ, remove_virt
     877              : 
     878              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_density_matrices_rpa'
     879              :       REAL(KIND=dp), PARAMETER                           :: stabilize_exp = 70.0_dp
     880              : 
     881              :       INTEGER                                            :: handle, i_mo, ikpgr, ispin, kplocal, &
     882              :                                                             nao, nmo, nspin
     883              :       INTEGER, DIMENSION(2)                              :: kp_range
     884          508 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eigenvalues, exp_scaling, occupation
     885              :       TYPE(cp_fm_struct_type), POINTER                   :: matrix_struct
     886              :       TYPE(cp_fm_type)                                   :: fwork
     887              :       TYPE(cp_fm_type), POINTER                          :: cpmat, rpmat
     888              :       TYPE(kpoint_env_type), POINTER                     :: kp
     889              :       TYPE(mo_set_type), POINTER                         :: mo_set
     890              : 
     891          508 :       CALL timeset(routineN, handle)
     892              : 
     893              :       ! only imaginary wavefunctions supported in kpoint cubic scaling RPA
     894          508 :       CPASSERT(kpoint%use_real_wfn .EQV. .FALSE.)
     895              : 
     896              :       ! work matrix
     897          508 :       mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
     898          508 :       CALL get_mo_set(mo_set, nao=nao, nmo=nmo)
     899              : 
     900              :       ! if this CPASSERT is triggered, please add all virtual MOs to SCF section,
     901              :       ! e.g. ADDED_MOS 1000000
     902          508 :       CPASSERT(nao == nmo)
     903              : 
     904         1524 :       ALLOCATE (exp_scaling(nmo))
     905              : 
     906          508 :       CALL cp_fm_get_info(mo_set%mo_coeff, matrix_struct=matrix_struct)
     907          508 :       CALL cp_fm_create(fwork, matrix_struct)
     908              : 
     909          508 :       CALL get_kpoint_info(kpoint, kp_range=kp_range)
     910          508 :       kplocal = kp_range(2) - kp_range(1) + 1
     911              : 
     912         1064 :       DO ikpgr = 1, kplocal
     913          556 :          kp => kpoint%kp_env(ikpgr)%kpoint_env
     914          556 :          nspin = SIZE(kp%mos, 2)
     915         1724 :          DO ispin = 1, nspin
     916          660 :             mo_set => kp%mos(1, ispin)
     917          660 :             CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
     918          660 :             rpmat => kp%wmat(1, ispin)
     919          660 :             cpmat => kp%wmat(2, ispin)
     920          660 :             CALL get_mo_set(mo_set, occupation_numbers=occupation)
     921          660 :             CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
     922              : 
     923          660 :             IF (remove_virt) THEN
     924          330 :                CALL cp_fm_column_scale(fwork, occupation)
     925              :             END IF
     926          660 :             IF (remove_occ) THEN
     927         6724 :                CALL cp_fm_column_scale(fwork, 2.0_dp/REAL(nspin, KIND=dp) - occupation)
     928              :             END IF
     929              : 
     930              :             ! proper spin
     931          660 :             IF (nspin == 1) THEN
     932          452 :                CALL cp_fm_scale(0.5_dp, fwork)
     933              :             END IF
     934              : 
     935        13448 :             DO i_mo = 1, nmo
     936              : 
     937        13448 :                IF (ABS(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
     938        12788 :                   exp_scaling(i_mo) = EXP(-ABS(tau*(eigenvalues(i_mo) - e_fermi)))
     939              :                ELSE
     940            0 :                   exp_scaling(i_mo) = 0.0_dp
     941              :                END IF
     942              :             END DO
     943              : 
     944          660 :             CALL cp_fm_column_scale(fwork, exp_scaling)
     945              : 
     946              :             ! Re(c)*Re(c)
     947          660 :             CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
     948          660 :             mo_set => kp%mos(2, ispin)
     949              :             ! Im(c)*Re(c)
     950          660 :             CALL parallel_gemm("N", "T", nao, nao, nmo, -1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
     951              :             ! Re(c)*Im(c)
     952          660 :             CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
     953              : 
     954          660 :             CALL cp_fm_to_fm(mo_set%mo_coeff, fwork)
     955              : 
     956          660 :             IF (remove_virt) THEN
     957          330 :                CALL cp_fm_column_scale(fwork, occupation)
     958              :             END IF
     959          660 :             IF (remove_occ) THEN
     960         6724 :                CALL cp_fm_column_scale(fwork, 2.0_dp/REAL(nspin, KIND=dp) - occupation)
     961              :             END IF
     962              : 
     963              :             ! proper spin
     964          660 :             IF (nspin == 1) THEN
     965          452 :                CALL cp_fm_scale(0.5_dp, fwork)
     966              :             END IF
     967              : 
     968        13448 :             DO i_mo = 1, nmo
     969        13448 :                IF (ABS(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp) THEN
     970        12788 :                   exp_scaling(i_mo) = EXP(-ABS(tau*(eigenvalues(i_mo) - e_fermi)))
     971              :                ELSE
     972            0 :                   exp_scaling(i_mo) = 0.0_dp
     973              :                END IF
     974              :             END DO
     975              : 
     976          660 :             CALL cp_fm_column_scale(fwork, exp_scaling)
     977              :             ! Im(c)*Im(c)
     978         1216 :             CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
     979              : 
     980              :          END DO
     981              : 
     982              :       END DO
     983              : 
     984          508 :       CALL cp_fm_release(fwork)
     985          508 :       DEALLOCATE (exp_scaling)
     986              : 
     987          508 :       CALL timestop(handle)
     988              : 
     989         1524 :    END SUBROUTINE kpoint_density_matrices_rpa
     990              : 
     991              : ! **************************************************************************************************
     992              : !> \brief ...
     993              : !> \param mat_dm_global ...
     994              : !> \param qs_env ...
     995              : !> \param ispin ...
     996              : !> \param num_integ_points ...
     997              : !> \param jquad ...
     998              : !> \param e_fermi ...
     999              : !> \param tau ...
    1000              : !> \param eps_filter ...
    1001              : !> \param num_cells_dm ...
    1002              : !> \param index_to_cell_dm ...
    1003              : !> \param remove_occ ...
    1004              : !> \param remove_virt ...
    1005              : !> \param first_jquad ...
    1006              : ! **************************************************************************************************
    1007           48 :    SUBROUTINE compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, &
    1008              :                                 eps_filter, num_cells_dm, index_to_cell_dm, remove_occ, remove_virt, &
    1009              :                                 first_jquad)
    1010              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_global
    1011              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1012              :       INTEGER, INTENT(IN)                                :: ispin, num_integ_points, jquad
    1013              :       REAL(KIND=dp), INTENT(IN)                          :: e_fermi, tau, eps_filter
    1014              :       INTEGER, INTENT(OUT)                               :: num_cells_dm
    1015              :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell_dm
    1016              :       LOGICAL, INTENT(IN)                                :: remove_occ, remove_virt
    1017              :       INTEGER, INTENT(IN)                                :: first_jquad
    1018              : 
    1019              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_transl_dm'
    1020              : 
    1021              :       INTEGER                                            :: handle, i_dim, i_img, iquad, jspin, nspin
    1022              :       INTEGER, DIMENSION(3)                              :: cell_grid_dm
    1023              :       TYPE(cell_type), POINTER                           :: cell
    1024           48 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_global_work, matrix_s_kp
    1025              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1026           48 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1027              : 
    1028           48 :       CALL timeset(routineN, handle)
    1029              : 
    1030              :       CALL get_qs_env(qs_env, &
    1031              :                       matrix_s_kp=matrix_s_kp, &
    1032              :                       mos=mos, &
    1033              :                       kpoints=kpoints, &
    1034           48 :                       cell=cell)
    1035              : 
    1036           48 :       nspin = SIZE(mos)
    1037              : 
    1038              :       ! we always use an odd number of image cells
    1039              :       ! CAUTION: also at another point, cell_grid_dm is defined, these definitions have to be identical
    1040          192 :       DO i_dim = 1, 3
    1041          192 :          cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
    1042              :       END DO
    1043              : 
    1044           48 :       num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
    1045              : 
    1046           48 :       NULLIFY (mat_dm_global_work)
    1047           48 :       CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm)
    1048              : 
    1049           96 :       DO jspin = 1, nspin
    1050              : 
    1051          240 :          DO i_img = 1, num_cells_dm
    1052              : 
    1053          144 :             ALLOCATE (mat_dm_global_work(jspin, i_img)%matrix)
    1054              :             CALL dbcsr_create(matrix=mat_dm_global_work(jspin, i_img)%matrix, &
    1055              :                               template=matrix_s_kp(1, 1)%matrix, &
    1056              :                               !                              matrix_type=dbcsr_type_symmetric)
    1057          144 :                               matrix_type=dbcsr_type_no_symmetry)
    1058              : 
    1059          144 :             CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, i_img)%matrix)
    1060              : 
    1061          192 :             CALL dbcsr_set(mat_dm_global_work(ispin, i_img)%matrix, 0.0_dp)
    1062              : 
    1063              :          END DO
    1064              : 
    1065              :       END DO
    1066              : 
    1067              :       ! density matrices in k-space weighted with EXP(-|e_i-e_F|*t) for occupied orbitals
    1068              :       CALL kpoint_density_matrices_rpa(kpoints, tau, e_fermi, &
    1069           48 :                                        remove_occ=remove_occ, remove_virt=remove_virt)
    1070              : 
    1071              :       ! overwrite the cell indices in kpoints
    1072           48 :       CALL init_cell_index_rpa(cell_grid_dm, kpoints%cell_to_index, kpoints%index_to_cell, cell)
    1073              : 
    1074              :       ! density matrices in real space, the cell vectors T for transforming are taken from kpoints%index_to_cell
    1075              :       ! (custom made for RPA) and not from sab_nl (which is symmetric and from SCF)
    1076           48 :       CALL density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, kpoints%index_to_cell)
    1077              : 
    1078              :       ! we need the index to cell for the density matrices later
    1079           48 :       index_to_cell_dm => kpoints%index_to_cell
    1080              : 
    1081              :       ! normally, jquad = 1 to allocate the matrix set, but for GW jquad = 0 is the exchange self-energy
    1082           48 :       IF (jquad == first_jquad) THEN
    1083              : 
    1084            8 :          NULLIFY (mat_dm_global)
    1085          200 :          ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm))
    1086              : 
    1087           56 :          DO iquad = first_jquad, num_integ_points
    1088          200 :             DO i_img = 1, num_cells_dm
    1089          144 :                NULLIFY (mat_dm_global(iquad, i_img)%matrix)
    1090          144 :                ALLOCATE (mat_dm_global(iquad, i_img)%matrix)
    1091              :                CALL dbcsr_create(matrix=mat_dm_global(iquad, i_img)%matrix, &
    1092              :                                  template=matrix_s_kp(1, 1)%matrix, &
    1093          192 :                                  matrix_type=dbcsr_type_no_symmetry)
    1094              : 
    1095              :             END DO
    1096              :          END DO
    1097              : 
    1098              :       END IF
    1099              : 
    1100          192 :       DO i_img = 1, num_cells_dm
    1101              : 
    1102              :          ! filter to get rid of the blocks full with zeros on the lower half, otherwise blocks doubled
    1103          144 :          CALL dbcsr_filter(mat_dm_global_work(ispin, i_img)%matrix, eps_filter)
    1104              : 
    1105              :          CALL dbcsr_copy(mat_dm_global(jquad, i_img)%matrix, &
    1106          192 :                          mat_dm_global_work(ispin, i_img)%matrix)
    1107              : 
    1108              :       END DO
    1109              : 
    1110           48 :       CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
    1111              : 
    1112           48 :       CALL timestop(handle)
    1113              : 
    1114           48 :    END SUBROUTINE compute_transl_dm
    1115              : 
    1116              : ! **************************************************************************************************
    1117              : !> \brief ...
    1118              : !> \param mat_dm_global ...
    1119              : !> \param qs_env ...
    1120              : !> \param ispin ...
    1121              : !> \param num_integ_points ...
    1122              : !> \param jquad ...
    1123              : !> \param e_fermi ...
    1124              : !> \param tau ...
    1125              : !> \param remove_occ ...
    1126              : !> \param remove_virt ...
    1127              : !> \param alloc_dm ...
    1128              : ! **************************************************************************************************
    1129          460 :    SUBROUTINE compute_periodic_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, &
    1130              :                                   remove_occ, remove_virt, alloc_dm)
    1131              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_global
    1132              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1133              :       INTEGER, INTENT(IN)                                :: ispin, num_integ_points, jquad
    1134              :       REAL(KIND=dp), INTENT(IN)                          :: e_fermi, tau
    1135              :       LOGICAL, INTENT(IN)                                :: remove_occ, remove_virt, alloc_dm
    1136              : 
    1137              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_periodic_dm'
    1138              : 
    1139              :       INTEGER                                            :: handle, iquad, jspin, nspin, num_cells_dm
    1140          460 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_global_work, matrix_s_kp
    1141              :       TYPE(kpoint_type), POINTER                         :: kpoints_G
    1142          460 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1143              : 
    1144          460 :       CALL timeset(routineN, handle)
    1145              : 
    1146          460 :       NULLIFY (matrix_s_kp, mos)
    1147              : 
    1148              :       CALL get_qs_env(qs_env, &
    1149              :                       matrix_s_kp=matrix_s_kp, &
    1150          460 :                       mos=mos)
    1151              : 
    1152          460 :       kpoints_G => qs_env%mp2_env%ri_rpa_im_time%kpoints_G
    1153              : 
    1154          460 :       nspin = SIZE(mos)
    1155              : 
    1156          460 :       num_cells_dm = 1
    1157              : 
    1158          460 :       NULLIFY (mat_dm_global_work)
    1159          460 :       CALL dbcsr_allocate_matrix_set(mat_dm_global_work, nspin, num_cells_dm)
    1160              : 
    1161              :       ! if necessaray, allocate mat_dm_global
    1162          460 :       IF (alloc_dm) THEN
    1163              : 
    1164           68 :          NULLIFY (mat_dm_global)
    1165          704 :          ALLOCATE (mat_dm_global(1:num_integ_points, num_cells_dm))
    1166              : 
    1167          500 :          DO iquad = 1, num_integ_points
    1168          432 :             NULLIFY (mat_dm_global(iquad, 1)%matrix)
    1169          432 :             ALLOCATE (mat_dm_global(iquad, 1)%matrix)
    1170              :             CALL dbcsr_create(matrix=mat_dm_global(iquad, 1)%matrix, &
    1171              :                               template=matrix_s_kp(1, 1)%matrix, &
    1172          500 :                               matrix_type=dbcsr_type_no_symmetry)
    1173              : 
    1174              :          END DO
    1175              : 
    1176              :       END IF
    1177              : 
    1178         1024 :       DO jspin = 1, nspin
    1179              : 
    1180          564 :          ALLOCATE (mat_dm_global_work(jspin, 1)%matrix)
    1181              :          CALL dbcsr_create(matrix=mat_dm_global_work(jspin, 1)%matrix, &
    1182              :                            template=matrix_s_kp(1, 1)%matrix, &
    1183          564 :                            matrix_type=dbcsr_type_no_symmetry)
    1184              : 
    1185          564 :          CALL dbcsr_reserve_all_blocks(mat_dm_global_work(jspin, 1)%matrix)
    1186              : 
    1187         1024 :          CALL dbcsr_set(mat_dm_global_work(jspin, 1)%matrix, 0.0_dp)
    1188              : 
    1189              :       END DO
    1190              : 
    1191              :       ! density matrices in k-space weighted with EXP(-|e_i-e_F|*t) for occupied orbitals
    1192              :       CALL kpoint_density_matrices_rpa(kpoints_G, tau, e_fermi, &
    1193          460 :                                        remove_occ=remove_occ, remove_virt=remove_virt)
    1194              : 
    1195          460 :       CALL density_matrix_from_kp_to_mic(kpoints_G, mat_dm_global_work, qs_env)
    1196              : 
    1197              :       CALL dbcsr_copy(mat_dm_global(jquad, 1)%matrix, &
    1198          460 :                       mat_dm_global_work(ispin, 1)%matrix)
    1199              : 
    1200          460 :       CALL dbcsr_deallocate_matrix_set(mat_dm_global_work)
    1201              : 
    1202          460 :       CALL timestop(handle)
    1203              : 
    1204          460 :    END SUBROUTINE compute_periodic_dm
    1205              : 
    1206              :    ! **************************************************************************************************
    1207              : !> \brief ...
    1208              : !> \param kpoints_G ...
    1209              : !> \param mat_dm_global_work ...
    1210              : !> \param qs_env ...
    1211              : ! **************************************************************************************************
    1212          460 :    SUBROUTINE density_matrix_from_kp_to_mic(kpoints_G, mat_dm_global_work, qs_env)
    1213              : 
    1214              :       TYPE(kpoint_type), POINTER                         :: kpoints_G
    1215              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_dm_global_work
    1216              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1217              : 
    1218              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_from_kp_to_mic'
    1219              : 
    1220              :       INTEGER                                            :: handle, iatom, iatom_old, ik, irow, &
    1221              :                                                             ispin, jatom, jatom_old, jcol, nao, &
    1222              :                                                             ncol_local, nkp, nrow_local, nspin, &
    1223              :                                                             num_cells
    1224              :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_from_ao_index
    1225          460 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1226          460 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
    1227              :       REAL(KIND=dp)                                      :: contribution, weight_im, weight_re
    1228              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
    1229          460 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: wkp
    1230          460 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
    1231              :       TYPE(cell_type), POINTER                           :: cell
    1232              :       TYPE(cp_fm_type)                                   :: fm_mat_work
    1233              :       TYPE(cp_fm_type), POINTER                          :: cpmat, rpmat
    1234              :       TYPE(kpoint_env_type), POINTER                     :: kp
    1235              :       TYPE(mo_set_type), POINTER                         :: mo_set
    1236          460 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1237              : 
    1238          460 :       CALL timeset(routineN, handle)
    1239              : 
    1240          460 :       NULLIFY (xkp, wkp)
    1241              : 
    1242          460 :       CALL cp_fm_create(fm_mat_work, kpoints_G%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct)
    1243          460 :       CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
    1244              : 
    1245          460 :       CALL get_kpoint_info(kpoints_G, nkp=nkp, xkp=xkp, wkp=wkp)
    1246          460 :       index_to_cell => kpoints_G%index_to_cell
    1247          460 :       num_cells = SIZE(index_to_cell, 2)
    1248              : 
    1249          460 :       nspin = SIZE(mat_dm_global_work, 1)
    1250              : 
    1251          460 :       mo_set => kpoints_G%kp_env(1)%kpoint_env%mos(1, 1)
    1252          460 :       CALL get_mo_set(mo_set, nao=nao)
    1253              : 
    1254         1380 :       ALLOCATE (atom_from_ao_index(nao))
    1255              : 
    1256          460 :       CALL get_atom_index_from_basis_function_index(qs_env, atom_from_ao_index, nao, "ORB")
    1257              : 
    1258              :       CALL cp_fm_get_info(matrix=kpoints_G%kp_env(1)%kpoint_env%wmat(1, 1), &
    1259              :                           nrow_local=nrow_local, &
    1260              :                           ncol_local=ncol_local, &
    1261              :                           row_indices=row_indices, &
    1262          460 :                           col_indices=col_indices)
    1263              : 
    1264          460 :       NULLIFY (cell, particle_set)
    1265          460 :       CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
    1266          460 :       CALL get_cell(cell=cell, h=hmat)
    1267              : 
    1268          460 :       iatom_old = 0
    1269          460 :       jatom_old = 0
    1270              : 
    1271         1024 :       DO ispin = 1, nspin
    1272              : 
    1273          564 :          CALL dbcsr_set(mat_dm_global_work(ispin, 1)%matrix, 0.0_dp)
    1274              : 
    1275         1128 :          DO ik = 1, nkp
    1276              : 
    1277          564 :             kp => kpoints_G%kp_env(ik)%kpoint_env
    1278          564 :             rpmat => kp%wmat(1, ispin)
    1279          564 :             cpmat => kp%wmat(2, ispin)
    1280              : 
    1281         6418 :             DO irow = 1, nrow_local
    1282       111768 :                DO jcol = 1, ncol_local
    1283              : 
    1284       105914 :                   iatom = atom_from_ao_index(row_indices(irow))
    1285       105914 :                   jatom = atom_from_ao_index(col_indices(jcol))
    1286              : 
    1287       105914 :                   IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN
    1288              : 
    1289              :                      CALL compute_weight_re_im(weight_re, weight_im, &
    1290              :                                                num_cells, iatom, jatom, xkp(1:3, ik), wkp(ik), &
    1291        15636 :                                                cell, index_to_cell, hmat, particle_set)
    1292              : 
    1293        15636 :                      iatom_old = iatom
    1294        15636 :                      jatom_old = jatom
    1295              : 
    1296              :                   END IF
    1297              : 
    1298              :                   ! minus sign because of i^2 = -1
    1299              :                   contribution = weight_re*rpmat%local_data(irow, jcol) - &
    1300       105914 :                                  weight_im*cpmat%local_data(irow, jcol)
    1301              : 
    1302       111204 :                   fm_mat_work%local_data(irow, jcol) = fm_mat_work%local_data(irow, jcol) + contribution
    1303              : 
    1304              :                END DO
    1305              :             END DO
    1306              : 
    1307              :          END DO ! ik
    1308              : 
    1309          564 :          CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(ispin, 1)%matrix, keep_sparsity=.FALSE.)
    1310         1024 :          CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
    1311              : 
    1312              :       END DO
    1313              : 
    1314          460 :       CALL cp_fm_release(fm_mat_work)
    1315          460 :       DEALLOCATE (atom_from_ao_index)
    1316              : 
    1317          460 :       CALL timestop(handle)
    1318              : 
    1319         1380 :    END SUBROUTINE density_matrix_from_kp_to_mic
    1320              : 
    1321              : ! **************************************************************************************************
    1322              : !> \brief ...
    1323              : !> \param kpoints ...
    1324              : !> \param mat_dm_global_work ...
    1325              : !> \param index_to_cell ...
    1326              : ! **************************************************************************************************
    1327           48 :    SUBROUTINE density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, index_to_cell)
    1328              : 
    1329              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1330              :       TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN)    :: mat_dm_global_work
    1331              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: index_to_cell
    1332              : 
    1333              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_from_kp_to_transl'
    1334              : 
    1335              :       INTEGER                                            :: handle, icell, ik, ispin, nkp, nspin, &
    1336              :                                                             xcell, ycell, zcell
    1337              :       REAL(KIND=dp)                                      :: arg, coskl, sinkl
    1338           48 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: wkp
    1339           48 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
    1340              :       TYPE(cp_fm_type), POINTER                          :: cpmat, rpmat
    1341              :       TYPE(dbcsr_type), POINTER                          :: mat_work_im, mat_work_re
    1342              :       TYPE(kpoint_env_type), POINTER                     :: kp
    1343              : 
    1344           48 :       CALL timeset(routineN, handle)
    1345              : 
    1346           48 :       NULLIFY (xkp, wkp)
    1347              : 
    1348           48 :       NULLIFY (mat_work_re)
    1349           48 :       CALL dbcsr_init_p(mat_work_re)
    1350              :       CALL dbcsr_create(matrix=mat_work_re, &
    1351              :                         template=mat_dm_global_work(1, 1)%matrix, &
    1352           48 :                         matrix_type=dbcsr_type_no_symmetry)
    1353              : 
    1354           48 :       NULLIFY (mat_work_im)
    1355           48 :       CALL dbcsr_init_p(mat_work_im)
    1356              :       CALL dbcsr_create(matrix=mat_work_im, &
    1357              :                         template=mat_dm_global_work(1, 1)%matrix, &
    1358           48 :                         matrix_type=dbcsr_type_no_symmetry)
    1359              : 
    1360           48 :       CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, wkp=wkp)
    1361              : 
    1362           48 :       nspin = SIZE(mat_dm_global_work, 1)
    1363              : 
    1364           48 :       CPASSERT(SIZE(mat_dm_global_work, 2) == SIZE(index_to_cell, 2))
    1365              : 
    1366           96 :       DO ispin = 1, nspin
    1367              : 
    1368          240 :          DO icell = 1, SIZE(mat_dm_global_work, 2)
    1369              : 
    1370          192 :             CALL dbcsr_set(mat_dm_global_work(ispin, icell)%matrix, 0.0_dp)
    1371              : 
    1372              :          END DO
    1373              : 
    1374              :       END DO
    1375              : 
    1376           96 :       DO ispin = 1, nspin
    1377              : 
    1378          192 :          DO ik = 1, nkp
    1379              : 
    1380           96 :             kp => kpoints%kp_env(ik)%kpoint_env
    1381           96 :             rpmat => kp%wmat(1, ispin)
    1382           96 :             cpmat => kp%wmat(2, ispin)
    1383              : 
    1384           96 :             CALL copy_fm_to_dbcsr(rpmat, mat_work_re, keep_sparsity=.FALSE.)
    1385           96 :             CALL copy_fm_to_dbcsr(cpmat, mat_work_im, keep_sparsity=.FALSE.)
    1386              : 
    1387          432 :             DO icell = 1, SIZE(mat_dm_global_work, 2)
    1388              : 
    1389          288 :                xcell = index_to_cell(1, icell)
    1390          288 :                ycell = index_to_cell(2, icell)
    1391          288 :                zcell = index_to_cell(3, icell)
    1392              : 
    1393          288 :                arg = REAL(xcell, dp)*xkp(1, ik) + REAL(ycell, dp)*xkp(2, ik) + REAL(zcell, dp)*xkp(3, ik)
    1394          288 :                coskl = wkp(ik)*COS(twopi*arg)
    1395          288 :                sinkl = wkp(ik)*SIN(twopi*arg)
    1396              : 
    1397          288 :                CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_re, 1.0_dp, coskl)
    1398          384 :                CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_im, 1.0_dp, sinkl)
    1399              : 
    1400              :             END DO
    1401              : 
    1402              :          END DO
    1403              :       END DO
    1404              : 
    1405           48 :       CALL dbcsr_release_p(mat_work_re)
    1406           48 :       CALL dbcsr_release_p(mat_work_im)
    1407              : 
    1408           48 :       CALL timestop(handle)
    1409              : 
    1410           48 :    END SUBROUTINE density_matrix_from_kp_to_transl
    1411              : 
    1412              : ! **************************************************************************************************
    1413              : !> \brief ...
    1414              : !> \param cell_grid ...
    1415              : !> \param cell_to_index ...
    1416              : !> \param index_to_cell ...
    1417              : !> \param cell ...
    1418              : ! **************************************************************************************************
    1419          536 :    SUBROUTINE init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
    1420              :       INTEGER, DIMENSION(3), INTENT(IN)                  :: cell_grid
    1421              :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    1422              :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
    1423              :       TYPE(cell_type), INTENT(IN), POINTER               :: cell
    1424              : 
    1425              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'init_cell_index_rpa'
    1426              : 
    1427              :       INTEGER                                            :: cell_counter, handle, i_cell, &
    1428              :                                                             index_min_dist, num_cells, xcell, &
    1429              :                                                             ycell, zcell
    1430              :       INTEGER, DIMENSION(3)                              :: itm
    1431          536 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell_unsorted
    1432          536 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index_unsorted
    1433          536 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: abs_cell_vectors
    1434              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_vector
    1435              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
    1436              : 
    1437          536 :       CALL timeset(routineN, handle)
    1438              : 
    1439          536 :       CALL get_cell(cell=cell, h=hmat)
    1440              : 
    1441          536 :       num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3)
    1442         2144 :       itm(:) = cell_grid(:)/2
    1443              : 
    1444              :       ! check that real space super lattice is a (2n+1)x(2m+1)x(2k+1) super lattice with the unit cell
    1445              :       ! in the middle
    1446          536 :       CPASSERT(cell_grid(1) .NE. itm(1)*2)
    1447          536 :       CPASSERT(cell_grid(2) .NE. itm(2)*2)
    1448          536 :       CPASSERT(cell_grid(3) .NE. itm(3)*2)
    1449              : 
    1450          536 :       IF (ASSOCIATED(cell_to_index)) DEALLOCATE (cell_to_index)
    1451          536 :       IF (ASSOCIATED(index_to_cell)) DEALLOCATE (index_to_cell)
    1452              : 
    1453         2680 :       ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
    1454        10568 :       cell_to_index_unsorted(:, :, :) = 0
    1455              : 
    1456         1608 :       ALLOCATE (index_to_cell_unsorted(3, num_cells))
    1457        18680 :       index_to_cell_unsorted(:, :) = 0
    1458              : 
    1459         2144 :       ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
    1460        10568 :       cell_to_index(:, :, :) = 0
    1461              : 
    1462         1072 :       ALLOCATE (index_to_cell(3, num_cells))
    1463        18680 :       index_to_cell(:, :) = 0
    1464              : 
    1465         1608 :       ALLOCATE (abs_cell_vectors(1:num_cells))
    1466              : 
    1467         1288 :       cell_counter = 0
    1468              : 
    1469         1288 :       DO xcell = -itm(1), itm(1)
    1470         2800 :          DO ycell = -itm(2), itm(2)
    1471         6800 :             DO zcell = -itm(3), itm(3)
    1472              : 
    1473         4536 :                cell_counter = cell_counter + 1
    1474         4536 :                cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter
    1475              : 
    1476         4536 :                index_to_cell_unsorted(1, cell_counter) = xcell
    1477         4536 :                index_to_cell_unsorted(2, cell_counter) = ycell
    1478         4536 :                index_to_cell_unsorted(3, cell_counter) = zcell
    1479              : 
    1480        72576 :                cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell_unsorted(1:3, cell_counter), dp))
    1481              : 
    1482         6048 :                abs_cell_vectors(cell_counter) = SQRT(cell_vector(1)**2 + cell_vector(2)**2 + cell_vector(3)**2)
    1483              : 
    1484              :             END DO
    1485              :          END DO
    1486              :       END DO
    1487              : 
    1488              :       ! first only do all symmetry non-equivalent cells, we need that because chi^T is computed for
    1489              :       ! cell indices T from index_to_cell(:,1:num_cells/2+1)
    1490         3072 :       DO i_cell = 1, num_cells/2 + 1
    1491              : 
    1492        14928 :          index_min_dist = MINLOC(abs_cell_vectors(1:num_cells/2 + 1), DIM=1)
    1493              : 
    1494         2536 :          xcell = index_to_cell_unsorted(1, index_min_dist)
    1495         2536 :          ycell = index_to_cell_unsorted(2, index_min_dist)
    1496         2536 :          zcell = index_to_cell_unsorted(3, index_min_dist)
    1497              : 
    1498         2536 :          index_to_cell(1, i_cell) = xcell
    1499         2536 :          index_to_cell(2, i_cell) = ycell
    1500         2536 :          index_to_cell(3, i_cell) = zcell
    1501              : 
    1502         2536 :          cell_to_index(xcell, ycell, zcell) = i_cell
    1503              : 
    1504         3072 :          abs_cell_vectors(index_min_dist) = 10000000000.0_dp
    1505              : 
    1506              :       END DO
    1507              : 
    1508              :       ! now all the remaining cells
    1509         2536 :       DO i_cell = num_cells/2 + 2, num_cells
    1510              : 
    1511        19712 :          index_min_dist = MINLOC(abs_cell_vectors(1:num_cells), DIM=1)
    1512              : 
    1513         2000 :          xcell = index_to_cell_unsorted(1, index_min_dist)
    1514         2000 :          ycell = index_to_cell_unsorted(2, index_min_dist)
    1515         2000 :          zcell = index_to_cell_unsorted(3, index_min_dist)
    1516              : 
    1517         2000 :          index_to_cell(1, i_cell) = xcell
    1518         2000 :          index_to_cell(2, i_cell) = ycell
    1519         2000 :          index_to_cell(3, i_cell) = zcell
    1520              : 
    1521         2000 :          cell_to_index(xcell, ycell, zcell) = i_cell
    1522              : 
    1523         2536 :          abs_cell_vectors(index_min_dist) = 10000000000.0_dp
    1524              : 
    1525              :       END DO
    1526              : 
    1527          536 :       DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors)
    1528              : 
    1529          536 :       CALL timestop(handle)
    1530              : 
    1531          536 :    END SUBROUTINE init_cell_index_rpa
    1532              : 
    1533              : ! **************************************************************************************************
    1534              : !> \brief ...
    1535              : !> \param i_cell_R ...
    1536              : !> \param i_cell_S ...
    1537              : !> \param i_cell_R_minus_S ...
    1538              : !> \param index_to_cell_3c ...
    1539              : !> \param cell_to_index_3c ...
    1540              : !> \param index_to_cell_dm ...
    1541              : !> \param R_minus_S_needed ...
    1542              : !> \param do_kpoints_cubic_RPA ...
    1543              : ! **************************************************************************************************
    1544        12616 :    SUBROUTINE get_diff_index_3c(i_cell_R, i_cell_S, i_cell_R_minus_S, index_to_cell_3c, &
    1545              :                                 cell_to_index_3c, index_to_cell_dm, R_minus_S_needed, &
    1546              :                                 do_kpoints_cubic_RPA)
    1547              : 
    1548              :       INTEGER, INTENT(IN)                                :: i_cell_R, i_cell_S
    1549              :       INTEGER, INTENT(OUT)                               :: i_cell_R_minus_S
    1550              :       INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN)  :: index_to_cell_3c
    1551              :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
    1552              :          INTENT(IN)                                      :: cell_to_index_3c
    1553              :       INTEGER, DIMENSION(:, :), INTENT(IN), POINTER      :: index_to_cell_dm
    1554              :       LOGICAL, INTENT(OUT)                               :: R_minus_S_needed
    1555              :       LOGICAL, INTENT(IN)                                :: do_kpoints_cubic_RPA
    1556              : 
    1557              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_diff_index_3c'
    1558              : 
    1559              :       INTEGER :: handle, x_cell_R, x_cell_R_minus_S, x_cell_S, y_cell_R, y_cell_R_minus_S, &
    1560              :          y_cell_S, z_cell_R, z_cell_R_minus_S, z_cell_S
    1561              : 
    1562        12616 :       CALL timeset(routineN, handle)
    1563              : 
    1564        12616 :       IF (do_kpoints_cubic_RPA) THEN
    1565              : 
    1566         4200 :          x_cell_R = index_to_cell_3c(1, i_cell_R)
    1567         4200 :          y_cell_R = index_to_cell_3c(2, i_cell_R)
    1568         4200 :          z_cell_R = index_to_cell_3c(3, i_cell_R)
    1569              : 
    1570         4200 :          x_cell_S = index_to_cell_dm(1, i_cell_S)
    1571         4200 :          y_cell_S = index_to_cell_dm(2, i_cell_S)
    1572         4200 :          z_cell_S = index_to_cell_dm(3, i_cell_S)
    1573              : 
    1574         4200 :          x_cell_R_minus_S = x_cell_R - x_cell_S
    1575         4200 :          y_cell_R_minus_S = y_cell_R - y_cell_S
    1576         4200 :          z_cell_R_minus_S = z_cell_R - z_cell_S
    1577              : 
    1578              :          IF (x_cell_R_minus_S .GE. LBOUND(cell_to_index_3c, 1) .AND. &
    1579              :              x_cell_R_minus_S .LE. UBOUND(cell_to_index_3c, 1) .AND. &
    1580              :              y_cell_R_minus_S .GE. LBOUND(cell_to_index_3c, 2) .AND. &
    1581              :              y_cell_R_minus_S .LE. UBOUND(cell_to_index_3c, 2) .AND. &
    1582        29300 :              z_cell_R_minus_S .GE. LBOUND(cell_to_index_3c, 3) .AND. &
    1583              :              z_cell_R_minus_S .LE. UBOUND(cell_to_index_3c, 3)) THEN
    1584              : 
    1585         3950 :             i_cell_R_minus_S = cell_to_index_3c(x_cell_R_minus_S, y_cell_R_minus_S, z_cell_R_minus_S)
    1586              : 
    1587              :             ! 0 means that there is no 3c index with this R-S vector because R-S is too big and the 3c integral is 0
    1588         3950 :             IF (i_cell_R_minus_S == 0) THEN
    1589              : 
    1590            0 :                R_minus_S_needed = .FALSE.
    1591            0 :                i_cell_R_minus_S = 0
    1592              : 
    1593              :             ELSE
    1594              : 
    1595         3950 :                R_minus_S_needed = .TRUE.
    1596              : 
    1597              :             END IF
    1598              : 
    1599              :          ELSE
    1600              : 
    1601          250 :             i_cell_R_minus_S = 0
    1602          250 :             R_minus_S_needed = .FALSE.
    1603              : 
    1604              :          END IF
    1605              : 
    1606              :       ELSE ! no k-points
    1607              : 
    1608         8416 :          R_minus_S_needed = .TRUE.
    1609         8416 :          i_cell_R_minus_S = 1
    1610              : 
    1611              :       END IF
    1612              : 
    1613        12616 :       CALL timestop(handle)
    1614              : 
    1615        12616 :    END SUBROUTINE get_diff_index_3c
    1616              : 
    1617              : ! **************************************************************************************************
    1618              : !> \brief ...
    1619              : !> \param i_cell_R ...
    1620              : !> \param i_cell_S ...
    1621              : !> \param i_cell_T ...
    1622              : !> \param i_cell_R_minus_S_minus_T ...
    1623              : !> \param index_to_cell_3c ...
    1624              : !> \param cell_to_index_3c ...
    1625              : !> \param index_to_cell_dm ...
    1626              : !> \param R_minus_S_minus_T_needed ...
    1627              : !> \param do_kpoints_cubic_RPA ...
    1628              : ! **************************************************************************************************
    1629        14716 :    SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, &
    1630         7358 :                                      index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
    1631              :                                      R_minus_S_minus_T_needed, &
    1632              :                                      do_kpoints_cubic_RPA)
    1633              : 
    1634              :       INTEGER, INTENT(IN)                                :: i_cell_R, i_cell_S, i_cell_T
    1635              :       INTEGER, INTENT(OUT)                               :: i_cell_R_minus_S_minus_T
    1636              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: index_to_cell_3c
    1637              :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :), &
    1638              :          INTENT(IN)                                      :: cell_to_index_3c
    1639              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: index_to_cell_dm
    1640              :       LOGICAL, INTENT(OUT)                               :: R_minus_S_minus_T_needed
    1641              :       LOGICAL, INTENT(IN)                                :: do_kpoints_cubic_RPA
    1642              : 
    1643              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_diff_diff_index_3c'
    1644              : 
    1645              :       INTEGER :: handle, x_cell_R, x_cell_R_minus_S_minus_T, x_cell_S, x_cell_T, y_cell_R, &
    1646              :          y_cell_R_minus_S_minus_T, y_cell_S, y_cell_T, z_cell_R, z_cell_R_minus_S_minus_T, &
    1647              :          z_cell_S, z_cell_T
    1648              : 
    1649         7358 :       CALL timeset(routineN, handle)
    1650              : 
    1651         7358 :       IF (do_kpoints_cubic_RPA) THEN
    1652              : 
    1653         3150 :          x_cell_R = index_to_cell_3c(1, i_cell_R)
    1654         3150 :          y_cell_R = index_to_cell_3c(2, i_cell_R)
    1655         3150 :          z_cell_R = index_to_cell_3c(3, i_cell_R)
    1656              : 
    1657         3150 :          x_cell_S = index_to_cell_dm(1, i_cell_S)
    1658         3150 :          y_cell_S = index_to_cell_dm(2, i_cell_S)
    1659         3150 :          z_cell_S = index_to_cell_dm(3, i_cell_S)
    1660              : 
    1661         3150 :          x_cell_T = index_to_cell_dm(1, i_cell_T)
    1662         3150 :          y_cell_T = index_to_cell_dm(2, i_cell_T)
    1663         3150 :          z_cell_T = index_to_cell_dm(3, i_cell_T)
    1664              : 
    1665         3150 :          x_cell_R_minus_S_minus_T = x_cell_R - x_cell_S - x_cell_T
    1666         3150 :          y_cell_R_minus_S_minus_T = y_cell_R - y_cell_S - y_cell_T
    1667         3150 :          z_cell_R_minus_S_minus_T = z_cell_R - z_cell_S - z_cell_T
    1668              : 
    1669              :          IF (x_cell_R_minus_S_minus_T .GE. LBOUND(cell_to_index_3c, 1) .AND. &
    1670              :              x_cell_R_minus_S_minus_T .LE. UBOUND(cell_to_index_3c, 1) .AND. &
    1671              :              y_cell_R_minus_S_minus_T .GE. LBOUND(cell_to_index_3c, 2) .AND. &
    1672              :              y_cell_R_minus_S_minus_T .LE. UBOUND(cell_to_index_3c, 2) .AND. &
    1673        22000 :              z_cell_R_minus_S_minus_T .GE. LBOUND(cell_to_index_3c, 3) .AND. &
    1674              :              z_cell_R_minus_S_minus_T .LE. UBOUND(cell_to_index_3c, 3)) THEN
    1675              : 
    1676              :             i_cell_R_minus_S_minus_T = cell_to_index_3c(x_cell_R_minus_S_minus_T, &
    1677              :                                                         y_cell_R_minus_S_minus_T, &
    1678         2900 :                                                         z_cell_R_minus_S_minus_T)
    1679              : 
    1680              :             ! index 0 means that there are only no 3c matrix elements because R-S-T is too big
    1681         2900 :             IF (i_cell_R_minus_S_minus_T == 0) THEN
    1682              : 
    1683            0 :                R_minus_S_minus_T_needed = .FALSE.
    1684              : 
    1685              :             ELSE
    1686              : 
    1687         2900 :                R_minus_S_minus_T_needed = .TRUE.
    1688              : 
    1689              :             END IF
    1690              : 
    1691              :          ELSE
    1692              : 
    1693          250 :             i_cell_R_minus_S_minus_T = 0
    1694          250 :             R_minus_S_minus_T_needed = .FALSE.
    1695              : 
    1696              :          END IF
    1697              : 
    1698              :          !  no k-kpoints
    1699              :       ELSE
    1700              : 
    1701         4208 :          R_minus_S_minus_T_needed = .TRUE.
    1702         4208 :          i_cell_R_minus_S_minus_T = 1
    1703              : 
    1704              :       END IF
    1705              : 
    1706         7358 :       CALL timestop(handle)
    1707              : 
    1708         7358 :    END SUBROUTINE get_diff_diff_index_3c
    1709              : 
    1710         1408 : END MODULE rpa_im_time
        

Generated by: LCOV version 2.0-1