LCOV - code coverage report
Current view: top level - src - gw_large_cell_gamma.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.7 % 740 686
Test Date: 2025-07-25 12:55:17 Functions: 97.5 % 40 39

            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 from paper [Graml2024]
      10              : !> \author Jan Wilhelm
      11              : !> \date 07.2023
      12              : ! **************************************************************************************************
      13              : MODULE gw_large_cell_gamma
      14              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      15              :    USE cell_types,                      ONLY: cell_type,&
      16              :                                               get_cell,&
      17              :                                               pbc
      18              :    USE constants_operator,              ONLY: operator_coulomb
      19              :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_uplo_to_full
      20              :    USE cp_cfm_cholesky,                 ONLY: cp_cfm_cholesky_decompose,&
      21              :                                               cp_cfm_cholesky_invert
      22              :    USE cp_cfm_diag,                     ONLY: cp_cfm_geeig
      23              :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      24              :                                               cp_cfm_get_info,&
      25              :                                               cp_cfm_release,&
      26              :                                               cp_cfm_to_cfm,&
      27              :                                               cp_cfm_to_fm,&
      28              :                                               cp_cfm_type,&
      29              :                                               cp_fm_to_cfm
      30              :    USE cp_dbcsr_api,                    ONLY: &
      31              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_block_p, &
      32              :         dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
      33              :         dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_release, dbcsr_set, &
      34              :         dbcsr_type
      35              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_reserve_all_blocks
      36              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      37              :                                               copy_fm_to_dbcsr,&
      38              :                                               dbcsr_deallocate_matrix_set
      39              :    USE cp_files,                        ONLY: close_file,&
      40              :                                               open_file
      41              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
      42              :    USE cp_fm_diag,                      ONLY: cp_fm_power
      43              :    USE cp_fm_types,                     ONLY: &
      44              :         cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_read_unformatted, cp_fm_release, &
      45              :         cp_fm_set_all, cp_fm_to_fm, cp_fm_type, cp_fm_write_unformatted
      46              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      47              :                                               cp_logger_type
      48              :    USE cp_output_handling,              ONLY: cp_p_file,&
      49              :                                               cp_print_key_should_output,&
      50              :                                               cp_print_key_unit_nr
      51              :    USE dbt_api,                         ONLY: dbt_clear,&
      52              :                                               dbt_contract,&
      53              :                                               dbt_copy,&
      54              :                                               dbt_create,&
      55              :                                               dbt_destroy,&
      56              :                                               dbt_type
      57              :    USE gw_communication,                ONLY: fm_to_local_tensor,&
      58              :                                               local_dbt_to_global_mat
      59              :    USE gw_utils,                        ONLY: analyt_conti_and_print,&
      60              :                                               de_init_bs_env,&
      61              :                                               time_to_freq
      62              :    USE input_constants,                 ONLY: rtp_method_bse
      63              :    USE input_section_types,             ONLY: section_vals_type
      64              :    USE kinds,                           ONLY: default_string_length,&
      65              :                                               dp,&
      66              :                                               int_8
      67              :    USE kpoint_coulomb_2c,               ONLY: build_2c_coulomb_matrix_kp
      68              :    USE kpoint_types,                    ONLY: kpoint_type
      69              :    USE machine,                         ONLY: m_walltime
      70              :    USE mathconstants,                   ONLY: twopi,&
      71              :                                               z_one,&
      72              :                                               z_zero
      73              :    USE message_passing,                 ONLY: mp_file_delete
      74              :    USE mp2_ri_2c,                       ONLY: RI_2c_integral_mat
      75              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      76              :    USE particle_types,                  ONLY: particle_type
      77              :    USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
      78              :    USE post_scf_bandstructure_utils,    ONLY: MIC_contribution_from_ikp,&
      79              :                                               cfm_ikp_from_fm_Gamma,&
      80              :                                               get_all_VBM_CBM_bandgaps
      81              :    USE qs_environment_types,            ONLY: get_qs_env,&
      82              :                                               qs_environment_type
      83              :    USE qs_kind_types,                   ONLY: qs_kind_type
      84              :    USE qs_tensors,                      ONLY: build_3c_integrals
      85              :    USE rpa_gw_kpoints_util,             ONLY: cp_cfm_power
      86              : #include "./base/base_uses.f90"
      87              : 
      88              :    IMPLICIT NONE
      89              : 
      90              :    PRIVATE
      91              : 
      92              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_large_cell_gamma'
      93              : 
      94              :    PUBLIC :: gw_calc_large_cell_Gamma, &
      95              :              compute_3c_integrals
      96              : 
      97              : CONTAINS
      98              : 
      99              : ! **************************************************************************************************
     100              : !> \brief Perform GW band structure calculation
     101              : !> \param qs_env ...
     102              : !> \param bs_env ...
     103              : !> \par History
     104              : !>    * 07.2023 created [Jan Wilhelm]
     105              : ! **************************************************************************************************
     106           22 :    SUBROUTINE gw_calc_large_cell_Gamma(qs_env, bs_env)
     107              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     108              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     109              : 
     110              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_large_cell_Gamma'
     111              : 
     112              :       INTEGER                                            :: handle
     113           22 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_Sigma_x_Gamma, fm_W_MIC_time
     114           22 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :)  :: fm_Sigma_c_Gamma_time
     115              : 
     116           22 :       CALL timeset(routineN, handle)
     117              : 
     118              :       ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
     119              :       ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
     120              :       ! χ_PQ(iτ,k=0) = sum_λν [sum_µ (µν|P) G^occ_µλ(i|τ|)] [sum_σ (σλ|Q) G^vir_σν(i|τ|)]
     121           22 :       CALL get_mat_chi_Gamma_tau(bs_env, qs_env, bs_env%mat_chi_Gamma_tau)
     122              : 
     123              :       ! χ_PQ(iτ,k=0) -> χ_PQ(iω,k) -> ε_PQ(iω,k) -> W_PQ(iω,k) -> W^MIC_PQ(iτ) -> M^-1*W^MIC*M^-1
     124           22 :       CALL get_W_MIC(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_MIC_time)
     125              : 
     126              :       ! D_µν = sum_n^occ C_µn(k=0) C_νn(k=0), V^trunc_PQ = sum_cell_R <phi_P,0|V^trunc|phi_Q,R>
     127              :       ! Σ^x_λσ(k=0) = sum_νQ [sum_P (νσ|P) V^trunc_PQ] [sum_µ (λµ|Q) D_µν)]
     128           22 :       CALL get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
     129              : 
     130              :       ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^occ_µν(i|τ|)], τ < 0
     131              :       ! Σ^c_λσ(iτ,k=0) = sum_νQ [sum_P (νσ|P) W^MIC_PQ(iτ)] [sum_µ (λµ|Q) G^vir_µν(i|τ|)], τ > 0
     132           22 :       CALL get_Sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
     133              : 
     134              :       ! Σ^c_λσ(iτ,k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
     135           22 :       CALL compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
     136              : 
     137           22 :       CALL de_init_bs_env(bs_env)
     138              : 
     139           22 :       CALL timestop(handle)
     140              : 
     141           22 :    END SUBROUTINE gw_calc_large_cell_Gamma
     142              : 
     143              : ! **************************************************************************************************
     144              : !> \brief ...
     145              : !> \param bs_env ...
     146              : !> \param qs_env ...
     147              : !> \param mat_chi_Gamma_tau ...
     148              : ! **************************************************************************************************
     149           22 :    SUBROUTINE get_mat_chi_Gamma_tau(bs_env, qs_env, mat_chi_Gamma_tau)
     150              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     151              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     152              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_Gamma_tau
     153              : 
     154              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
     155              : 
     156              :       INTEGER :: handle, i_intval_idx, i_t, inner_loop_atoms_interval_index, ispin, j_intval_idx
     157              :       INTEGER, DIMENSION(2)                              :: i_atoms, IL_atoms, j_atoms
     158              :       LOGICAL                                            :: dist_too_long_i, dist_too_long_j
     159              :       REAL(KIND=dp)                                      :: t1, tau
     160          550 :       TYPE(dbt_type)                                     :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
     161          374 :                                                             t_3c_for_Gvir, t_3c_x_Gocc, &
     162          374 :                                                             t_3c_x_Gocc_2, t_3c_x_Gvir, &
     163          198 :                                                             t_3c_x_Gvir_2
     164              : 
     165           22 :       CALL timeset(routineN, handle)
     166              : 
     167          346 :       DO i_t = 1, bs_env%num_time_freq_points
     168              : 
     169          324 :          t1 = m_walltime()
     170              : 
     171          324 :          IF (bs_env%read_chi(i_t)) THEN
     172              : 
     173            0 :             CALL fm_read(bs_env%fm_RI_RI, bs_env, bs_env%chi_name, i_t)
     174              : 
     175              :             CALL copy_fm_to_dbcsr(bs_env%fm_RI_RI, mat_chi_Gamma_tau(i_t)%matrix, &
     176            0 :                                   keep_sparsity=.FALSE.)
     177              : 
     178            0 :             IF (bs_env%unit_nr > 0) THEN
     179              :                WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') &
     180            0 :                   'Read χ(iτ,k=0) from file for time point  ', i_t, ' /', &
     181            0 :                   bs_env%num_time_freq_points, &
     182            0 :                   ',    Execution time', m_walltime() - t1, ' s'
     183              :             END IF
     184              : 
     185              :             CYCLE
     186              : 
     187              :          END IF
     188              : 
     189          324 :          IF (.NOT. bs_env%calc_chi(i_t)) CYCLE
     190              : 
     191              :          CALL create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
     192          224 :                                  t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
     193              : 
     194              :          ! 1. compute G^occ and G^vir
     195              :          !    Background: G^σ(iτ) = G^occ,σ(iτ) * Θ(-τ) + G^vir,σ(iτ) * Θ(τ), σ ∈ {↑,↓}
     196              :          !    G^occ,σ_µλ(i|τ|,k=0) = sum_n^occ C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
     197              :          !    G^vir,σ_µλ(i|τ|,k=0) = sum_n^vir C^σ_µn(k=0) e^(-|(ϵ^σ_nk=0-ϵ_F)τ|) C^σ_λn(k=0)
     198          224 :          tau = bs_env%imag_time_points(i_t)
     199              : 
     200          468 :          DO ispin = 1, bs_env%n_spin
     201          244 :             CALL G_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.TRUE., vir=.FALSE.)
     202          244 :             CALL G_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.FALSE., vir=.TRUE.)
     203              : 
     204              :             CALL fm_to_local_tensor(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, &
     205              :                                     bs_env%mat_ao_ao_tensor%matrix, t_2c_Gocc, bs_env, &
     206          244 :                                     bs_env%atoms_j_t_group)
     207              :             CALL fm_to_local_tensor(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, &
     208              :                                     bs_env%mat_ao_ao_tensor%matrix, t_2c_Gvir, bs_env, &
     209          244 :                                     bs_env%atoms_i_t_group)
     210              : 
     211              :             ! every group has its own range of i_atoms and j_atoms; only deal with a
     212              :             ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
     213          712 :             DO i_intval_idx = 1, bs_env%n_intervals_i
     214          732 :                DO j_intval_idx = 1, bs_env%n_intervals_j
     215          732 :                   i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
     216          732 :                   j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
     217              : 
     218          488 :                   DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
     219              : 
     220          732 :                      IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
     221              : 
     222          244 :                      CALL check_dist(i_atoms, IL_atoms, qs_env, bs_env, dist_too_long_i)
     223          244 :                      CALL check_dist(j_atoms, IL_atoms, qs_env, bs_env, dist_too_long_j)
     224          244 :                      IF (dist_too_long_i .OR. dist_too_long_j) CYCLE
     225              : 
     226              :                      ! 2. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
     227          244 :                      CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gocc, i_atoms, IL_atoms)
     228              : 
     229              :                      ! 3. tensor operation M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
     230              :                      CALL G_times_3c(t_3c_for_Gocc, t_2c_Gocc, t_3c_x_Gocc, bs_env, &
     231          244 :                                      j_atoms, i_atoms, IL_atoms)
     232              : 
     233              :                      ! 4. compute 3-center integrals (σλ|Q) ("|": truncated Coulomb operator)
     234          244 :                      CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gvir, j_atoms, IL_atoms)
     235              : 
     236              :                      ! 5. tensor operation N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
     237              :                      CALL G_times_3c(t_3c_for_Gvir, t_2c_Gvir, t_3c_x_Gvir, bs_env, &
     238          488 :                                      i_atoms, j_atoms, IL_atoms)
     239              : 
     240              :                   END DO ! IL_atoms
     241              : 
     242              :                   ! 6. reorder tensors
     243          244 :                   CALL dbt_copy(t_3c_x_Gocc, t_3c_x_Gocc_2, move_data=.TRUE., order=[1, 3, 2])
     244          244 :                   CALL dbt_copy(t_3c_x_Gvir, t_3c_x_Gvir_2, move_data=.TRUE.)
     245              : 
     246              :                   ! 7. tensor operation χ_PQ(iτ,k=0) = sum_λν M_λνP(iτ) N_νλQ(iτ),
     247              :                   CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
     248              :                                     tensor_1=t_3c_x_Gocc_2, tensor_2=t_3c_x_Gvir_2, &
     249              :                                     beta=1.0_dp, tensor_3=bs_env%t_chi, &
     250              :                                     contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
     251              :                                     contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
     252          488 :                                     filter_eps=bs_env%eps_filter, move_data=.TRUE.)
     253              : 
     254              :                END DO ! j_atoms
     255              :             END DO ! i_atoms
     256              :          END DO ! ispin
     257              : 
     258              :          ! 8. communicate data of χ_PQ(iτ,k=0) in tensor bs_env%t_chi (which local in the
     259              :          !    subgroup) to the global dbcsr matrix mat_chi_Gamma_tau (which stores
     260              :          !    χ_PQ(iτ,k=0) for all time points)
     261              :          CALL local_dbt_to_global_mat(bs_env%t_chi, bs_env%mat_RI_RI_tensor%matrix, &
     262          224 :                                       mat_chi_Gamma_tau(i_t)%matrix, bs_env%para_env)
     263              : 
     264              :          CALL write_matrix(mat_chi_Gamma_tau(i_t)%matrix, i_t, bs_env%chi_name, &
     265          224 :                            bs_env%fm_RI_RI, qs_env)
     266              : 
     267              :          CALL destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
     268          224 :                                   t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
     269              : 
     270          246 :          IF (bs_env%unit_nr > 0) THEN
     271              :             WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
     272          112 :                'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
     273          224 :                ',    Execution time', m_walltime() - t1, ' s'
     274              :          END IF
     275              : 
     276              :       END DO ! i_t
     277              : 
     278           22 :       IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
     279              : 
     280           22 :       CALL timestop(handle)
     281              : 
     282           22 :    END SUBROUTINE get_mat_chi_Gamma_tau
     283              : 
     284              : ! **************************************************************************************************
     285              : !> \brief ...
     286              : !> \param fm ...
     287              : !> \param bs_env ...
     288              : !> \param mat_name ...
     289              : !> \param idx ...
     290              : ! **************************************************************************************************
     291          352 :    SUBROUTINE fm_read(fm, bs_env, mat_name, idx)
     292              :       TYPE(cp_fm_type)                                   :: fm
     293              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     294              :       CHARACTER(LEN=*)                                   :: mat_name
     295              :       INTEGER                                            :: idx
     296              : 
     297              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'fm_read'
     298              : 
     299              :       CHARACTER(LEN=default_string_length)               :: f_chi
     300              :       INTEGER                                            :: handle, unit_nr
     301              : 
     302          352 :       CALL timeset(routineN, handle)
     303              : 
     304          352 :       unit_nr = -1
     305          352 :       IF (bs_env%para_env%is_source()) THEN
     306              : 
     307          176 :          IF (idx < 10) THEN
     308           87 :             WRITE (f_chi, '(3A,I1,A)') TRIM(bs_env%prefix), TRIM(mat_name), "_0", idx, ".matrix"
     309           89 :          ELSE IF (idx < 100) THEN
     310           89 :             WRITE (f_chi, '(3A,I2,A)') TRIM(bs_env%prefix), TRIM(mat_name), "_", idx, ".matrix"
     311              :          ELSE
     312            0 :             CPABORT('Please implement more than 99 time/frequency points.')
     313              :          END IF
     314              : 
     315              :          CALL open_file(file_name=TRIM(f_chi), file_action="READ", file_form="UNFORMATTED", &
     316          176 :                         file_position="REWIND", file_status="OLD", unit_number=unit_nr)
     317              : 
     318              :       END IF
     319              : 
     320          352 :       CALL cp_fm_read_unformatted(fm, unit_nr)
     321              : 
     322          352 :       IF (bs_env%para_env%is_source()) CALL close_file(unit_number=unit_nr)
     323              : 
     324          352 :       CALL timestop(handle)
     325              : 
     326          352 :    END SUBROUTINE fm_read
     327              : 
     328              : ! **************************************************************************************************
     329              : !> \brief ...
     330              : !> \param t_2c_Gocc ...
     331              : !> \param t_2c_Gvir ...
     332              : !> \param t_3c_for_Gocc ...
     333              : !> \param t_3c_for_Gvir ...
     334              : !> \param t_3c_x_Gocc ...
     335              : !> \param t_3c_x_Gvir ...
     336              : !> \param t_3c_x_Gocc_2 ...
     337              : !> \param t_3c_x_Gvir_2 ...
     338              : !> \param bs_env ...
     339              : ! **************************************************************************************************
     340          224 :    SUBROUTINE create_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
     341              :                                  t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2, bs_env)
     342              : 
     343              :       TYPE(dbt_type)                                     :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
     344              :                                                             t_3c_for_Gvir, t_3c_x_Gocc, &
     345              :                                                             t_3c_x_Gvir, t_3c_x_Gocc_2, &
     346              :                                                             t_3c_x_Gvir_2
     347              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     348              : 
     349              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_tensors_chi'
     350              : 
     351              :       INTEGER                                            :: handle
     352              : 
     353          224 :       CALL timeset(routineN, handle)
     354              : 
     355          224 :       CALL dbt_create(bs_env%t_G, t_2c_Gocc)
     356          224 :       CALL dbt_create(bs_env%t_G, t_2c_Gvir)
     357          224 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gocc)
     358          224 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gvir)
     359          224 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gocc)
     360          224 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gvir)
     361          224 :       CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gocc_2)
     362          224 :       CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gvir_2)
     363              : 
     364          224 :       CALL timestop(handle)
     365              : 
     366          224 :    END SUBROUTINE create_tensors_chi
     367              : 
     368              : ! **************************************************************************************************
     369              : !> \brief ...
     370              : !> \param t_2c_Gocc ...
     371              : !> \param t_2c_Gvir ...
     372              : !> \param t_3c_for_Gocc ...
     373              : !> \param t_3c_for_Gvir ...
     374              : !> \param t_3c_x_Gocc ...
     375              : !> \param t_3c_x_Gvir ...
     376              : !> \param t_3c_x_Gocc_2 ...
     377              : !> \param t_3c_x_Gvir_2 ...
     378              : ! **************************************************************************************************
     379          224 :    SUBROUTINE destroy_tensors_chi(t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, t_3c_for_Gvir, &
     380              :                                   t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
     381              :       TYPE(dbt_type)                                     :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
     382              :                                                             t_3c_for_Gvir, t_3c_x_Gocc, &
     383              :                                                             t_3c_x_Gvir, t_3c_x_Gocc_2, &
     384              :                                                             t_3c_x_Gvir_2
     385              : 
     386              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_tensors_chi'
     387              : 
     388              :       INTEGER                                            :: handle
     389              : 
     390          224 :       CALL timeset(routineN, handle)
     391              : 
     392          224 :       CALL dbt_destroy(t_2c_Gocc)
     393          224 :       CALL dbt_destroy(t_2c_Gvir)
     394          224 :       CALL dbt_destroy(t_3c_for_Gocc)
     395          224 :       CALL dbt_destroy(t_3c_for_Gvir)
     396          224 :       CALL dbt_destroy(t_3c_x_Gocc)
     397          224 :       CALL dbt_destroy(t_3c_x_Gvir)
     398          224 :       CALL dbt_destroy(t_3c_x_Gocc_2)
     399          224 :       CALL dbt_destroy(t_3c_x_Gvir_2)
     400              : 
     401          224 :       CALL timestop(handle)
     402              : 
     403          224 :    END SUBROUTINE destroy_tensors_chi
     404              : 
     405              : ! **************************************************************************************************
     406              : !> \brief ...
     407              : !> \param matrix ...
     408              : !> \param matrix_index ...
     409              : !> \param matrix_name ...
     410              : !> \param fm ...
     411              : !> \param qs_env ...
     412              : ! **************************************************************************************************
     413          730 :    SUBROUTINE write_matrix(matrix, matrix_index, matrix_name, fm, qs_env)
     414              :       TYPE(dbcsr_type)                                   :: matrix
     415              :       INTEGER                                            :: matrix_index
     416              :       CHARACTER(LEN=*)                                   :: matrix_name
     417              :       TYPE(cp_fm_type), INTENT(IN), POINTER              :: fm
     418              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     419              : 
     420              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'write_matrix'
     421              : 
     422              :       INTEGER                                            :: handle
     423              : 
     424          730 :       CALL timeset(routineN, handle)
     425              : 
     426          730 :       CALL cp_fm_set_all(fm, 0.0_dp)
     427              : 
     428          730 :       CALL copy_dbcsr_to_fm(matrix, fm)
     429              : 
     430          730 :       CALL fm_write(fm, matrix_index, matrix_name, qs_env)
     431              : 
     432          730 :       CALL timestop(handle)
     433              : 
     434          730 :    END SUBROUTINE write_matrix
     435              : 
     436              : ! **************************************************************************************************
     437              : !> \brief ...
     438              : !> \param fm ...
     439              : !> \param matrix_index ...
     440              : !> \param matrix_name ...
     441              : !> \param qs_env ...
     442              : ! **************************************************************************************************
     443          962 :    SUBROUTINE fm_write(fm, matrix_index, matrix_name, qs_env)
     444              :       TYPE(cp_fm_type)                                   :: fm
     445              :       INTEGER                                            :: matrix_index
     446              :       CHARACTER(LEN=*)                                   :: matrix_name
     447              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     448              : 
     449              :       CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
     450              :          routineN = 'fm_write'
     451              : 
     452              :       CHARACTER(LEN=default_string_length)               :: filename
     453              :       INTEGER                                            :: handle, unit_nr
     454              :       TYPE(cp_logger_type), POINTER                      :: logger
     455              :       TYPE(section_vals_type), POINTER                   :: input
     456              : 
     457          962 :       CALL timeset(routineN, handle)
     458              : 
     459          962 :       CALL get_qs_env(qs_env, input=input)
     460              : 
     461          962 :       logger => cp_get_default_logger()
     462              : 
     463          962 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
     464              : 
     465          780 :          IF (matrix_index < 10) THEN
     466          380 :             WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
     467          400 :          ELSE IF (matrix_index < 100) THEN
     468          400 :             WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
     469              :          ELSE
     470            0 :             CPABORT('Please implement more than 99 time/frequency points.')
     471              :          END IF
     472              : 
     473              :          unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
     474              :                                         file_form="UNFORMATTED", middle_name=TRIM(filename), &
     475          780 :                                         file_position="REWIND", file_action="WRITE")
     476              : 
     477          780 :          CALL cp_fm_write_unformatted(fm, unit_nr)
     478          780 :          IF (unit_nr > 0) THEN
     479          390 :             CALL close_file(unit_nr)
     480              :          END IF
     481              :       END IF
     482              : 
     483          962 :       CALL timestop(handle)
     484              : 
     485          962 :    END SUBROUTINE fm_write
     486              : 
     487              : ! **************************************************************************************************
     488              : !> \brief ...
     489              : !> \param bs_env ...
     490              : !> \param tau ...
     491              : !> \param fm_G_Gamma ...
     492              : !> \param ispin ...
     493              : !> \param occ ...
     494              : !> \param vir ...
     495              : ! **************************************************************************************************
     496         1988 :    SUBROUTINE G_occ_vir(bs_env, tau, fm_G_Gamma, ispin, occ, vir)
     497              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     498              :       REAL(KIND=dp)                                      :: tau
     499              :       TYPE(cp_fm_type)                                   :: fm_G_Gamma
     500              :       INTEGER                                            :: ispin
     501              :       LOGICAL                                            :: occ, vir
     502              : 
     503              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'G_occ_vir'
     504              : 
     505              :       INTEGER                                            :: handle, homo, i_row_local, j_col, &
     506              :                                                             j_col_local, n_mo, ncol_local, &
     507              :                                                             nrow_local
     508          994 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices
     509              :       REAL(KIND=dp)                                      :: tau_E
     510              : 
     511          994 :       CALL timeset(routineN, handle)
     512              : 
     513          994 :       CPASSERT(occ .NEQV. vir)
     514              : 
     515              :       CALL cp_fm_get_info(matrix=bs_env%fm_work_mo(1), &
     516              :                           nrow_local=nrow_local, &
     517              :                           ncol_local=ncol_local, &
     518          994 :                           col_indices=col_indices)
     519              : 
     520          994 :       n_mo = bs_env%n_ao
     521          994 :       homo = bs_env%n_occ(ispin)
     522              : 
     523          994 :       CALL cp_fm_to_fm(bs_env%fm_mo_coeff_Gamma(ispin), bs_env%fm_work_mo(1))
     524              : 
     525         3899 :       DO i_row_local = 1, nrow_local
     526        41608 :          DO j_col_local = 1, ncol_local
     527              : 
     528        37709 :             j_col = col_indices(j_col_local)
     529              : 
     530        37709 :             tau_E = ABS(tau*0.5_dp*(bs_env%eigenval_scf_Gamma(j_col, ispin) - bs_env%e_fermi(ispin)))
     531              : 
     532        37709 :             IF (tau_E < bs_env%stabilize_exp) THEN
     533              :                bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = &
     534        36917 :                   bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local)*EXP(-tau_E)
     535              :             ELSE
     536          792 :                bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
     537              :             END IF
     538              : 
     539        40614 :             IF ((occ .AND. j_col > homo) .OR. (vir .AND. j_col <= homo)) THEN
     540        19222 :                bs_env%fm_work_mo(1)%local_data(i_row_local, j_col_local) = 0.0_dp
     541              :             END IF
     542              : 
     543              :          END DO
     544              :       END DO
     545              : 
     546              :       CALL parallel_gemm(transa="N", transb="T", m=n_mo, n=n_mo, k=n_mo, alpha=1.0_dp, &
     547              :                          matrix_a=bs_env%fm_work_mo(1), matrix_b=bs_env%fm_work_mo(1), &
     548          994 :                          beta=0.0_dp, matrix_c=fm_G_Gamma)
     549              : 
     550          994 :       CALL timestop(handle)
     551              : 
     552          994 :    END SUBROUTINE G_occ_vir
     553              : 
     554              : ! **************************************************************************************************
     555              : !> \brief ...
     556              : !> \param qs_env ...
     557              : !> \param bs_env ...
     558              : !> \param t_3c ...
     559              : !> \param atoms_AO_1 ...
     560              : !> \param atoms_AO_2 ...
     561              : !> \param atoms_RI ...
     562              : ! **************************************************************************************************
     563         1214 :    SUBROUTINE compute_3c_integrals(qs_env, bs_env, t_3c, atoms_AO_1, atoms_AO_2, atoms_RI)
     564              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     565              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     566              :       TYPE(dbt_type)                                     :: t_3c
     567              :       INTEGER, DIMENSION(2), OPTIONAL                    :: atoms_AO_1, atoms_AO_2, atoms_RI
     568              : 
     569              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_integrals'
     570              : 
     571              :       INTEGER                                            :: handle
     572              :       INTEGER, DIMENSION(2)                              :: my_atoms_AO_1, my_atoms_AO_2, my_atoms_RI
     573         1214 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_array
     574              : 
     575         1214 :       CALL timeset(routineN, handle)
     576              : 
     577              :       ! free memory (not clear whether memory has been freed previously)
     578         1214 :       CALL dbt_clear(t_3c)
     579              : 
     580        13354 :       ALLOCATE (t_3c_array(1, 1))
     581         1214 :       CALL dbt_create(t_3c, t_3c_array(1, 1))
     582              : 
     583         1214 :       IF (PRESENT(atoms_AO_1)) THEN
     584              :          my_atoms_AO_1 = atoms_AO_1
     585              :       ELSE
     586         1446 :          my_atoms_AO_1 = [1, bs_env%n_atom]
     587              :       END IF
     588         1214 :       IF (PRESENT(atoms_AO_2)) THEN
     589              :          my_atoms_AO_2 = atoms_AO_2
     590              :       ELSE
     591          768 :          my_atoms_AO_2 = [1, bs_env%n_atom]
     592              :       END IF
     593         1214 :       IF (PRESENT(atoms_RI)) THEN
     594              :          my_atoms_RI = atoms_RI
     595              :       ELSE
     596         1500 :          my_atoms_RI = [1, bs_env%n_atom]
     597              :       END IF
     598              : 
     599              :       CALL build_3c_integrals(t_3c_array, &
     600              :                               bs_env%eps_filter, &
     601              :                               qs_env, &
     602              :                               bs_env%nl_3c, &
     603              :                               int_eps=bs_env%eps_filter, &
     604              :                               basis_i=bs_env%basis_set_RI, &
     605              :                               basis_j=bs_env%basis_set_AO, &
     606              :                               basis_k=bs_env%basis_set_AO, &
     607              :                               potential_parameter=bs_env%ri_metric, &
     608              :                               bounds_i=atoms_RI, &
     609              :                               bounds_j=atoms_AO_1, &
     610              :                               bounds_k=atoms_AO_2, &
     611         1214 :                               desymmetrize=.FALSE.)
     612              : 
     613         1214 :       CALL dbt_copy(t_3c_array(1, 1), t_3c, move_data=.TRUE.)
     614              : 
     615         1214 :       CALL dbt_destroy(t_3c_array(1, 1))
     616         2428 :       DEALLOCATE (t_3c_array)
     617              : 
     618         1214 :       CALL timestop(handle)
     619              : 
     620         2428 :    END SUBROUTINE compute_3c_integrals
     621              : 
     622              : ! **************************************************************************************************
     623              : !> \brief ...
     624              : !> \param t_3c_for_G ...
     625              : !> \param t_G ...
     626              : !> \param t_M ...
     627              : !> \param bs_env ...
     628              : !> \param atoms_AO_1 ...
     629              : !> \param atoms_AO_2 ...
     630              : !> \param atoms_IL ...
     631              : ! **************************************************************************************************
     632          488 :    SUBROUTINE G_times_3c(t_3c_for_G, t_G, t_M, bs_env, atoms_AO_1, atoms_AO_2, atoms_IL)
     633              :       TYPE(dbt_type)                                     :: t_3c_for_G, t_G, t_M
     634              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     635              :       INTEGER, DIMENSION(2)                              :: atoms_AO_1, atoms_AO_2, atoms_IL
     636              : 
     637              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'G_times_3c'
     638              : 
     639              :       INTEGER                                            :: handle
     640              :       INTEGER, DIMENSION(2)                              :: bounds_IL, bounds_l
     641              :       INTEGER, DIMENSION(2, 2)                           :: bounds_k
     642              : 
     643          488 :       CALL timeset(routineN, handle)
     644              : 
     645              :       ! JW bounds_IL and bounds_k do not safe any operations, but maybe communication
     646              :       !    maybe remove "bounds_1=bounds_IL, &" and "bounds_2=bounds_k, &" later and
     647              :       !    check whether performance improves
     648              : 
     649              :       bounds_IL(1:2) = [bs_env%i_ao_start_from_atom(atoms_IL(1)), &
     650         1464 :                         bs_env%i_ao_end_from_atom(atoms_IL(2))]
     651         1464 :       bounds_k(1:2, 1) = [1, bs_env%n_RI]
     652              :       bounds_k(1:2, 2) = [bs_env%i_ao_start_from_atom(atoms_AO_2(1)), &
     653         1464 :                           bs_env%i_ao_end_from_atom(atoms_AO_2(2))]
     654              :       bounds_l(1:2) = [bs_env%i_ao_start_from_atom(atoms_AO_1(1)), &
     655         1464 :                        bs_env%i_ao_end_from_atom(atoms_AO_1(2))]
     656              : 
     657              :       CALL dbt_contract(alpha=1.0_dp, &
     658              :                         tensor_1=t_3c_for_G, &
     659              :                         tensor_2=t_G, &
     660              :                         beta=1.0_dp, &
     661              :                         tensor_3=t_M, &
     662              :                         contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
     663              :                         contract_2=[2], notcontract_2=[1], map_2=[3], &
     664              :                         bounds_1=bounds_IL, &
     665              :                         bounds_2=bounds_k, &
     666              :                         bounds_3=bounds_l, &
     667          488 :                         filter_eps=bs_env%eps_filter)
     668              : 
     669          488 :       CALL dbt_clear(t_3c_for_G)
     670              : 
     671          488 :       CALL timestop(handle)
     672              : 
     673          488 :    END SUBROUTINE G_times_3c
     674              : 
     675              : ! **************************************************************************************************
     676              : !> \brief ...
     677              : !> \param atoms_1 ...
     678              : !> \param atoms_2 ...
     679              : !> \param qs_env ...
     680              : !> \param bs_env ...
     681              : !> \param dist_too_long ...
     682              : ! **************************************************************************************************
     683          488 :    SUBROUTINE check_dist(atoms_1, atoms_2, qs_env, bs_env, dist_too_long)
     684              :       INTEGER, DIMENSION(2)                              :: atoms_1, atoms_2
     685              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     686              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     687              :       LOGICAL                                            :: dist_too_long
     688              : 
     689              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'check_dist'
     690              : 
     691              :       INTEGER                                            :: atom_1, atom_2, handle
     692              :       REAL(dp)                                           :: abs_rab, min_dist_AO_atoms
     693              :       REAL(KIND=dp), DIMENSION(3)                        :: rab
     694              :       TYPE(cell_type), POINTER                           :: cell
     695          488 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     696              : 
     697          488 :       CALL timeset(routineN, handle)
     698              : 
     699          488 :       CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
     700              : 
     701          488 :       min_dist_AO_atoms = 1.0E5_dp
     702         1512 :       DO atom_1 = atoms_1(1), atoms_1(2)
     703         3704 :          DO atom_2 = atoms_2(1), atoms_2(2)
     704              : 
     705         2192 :             rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
     706              : 
     707         2192 :             abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
     708              : 
     709         3216 :             min_dist_AO_atoms = MIN(min_dist_AO_atoms, abs_rab)
     710              : 
     711              :          END DO
     712              :       END DO
     713              : 
     714          488 :       dist_too_long = (min_dist_AO_atoms > bs_env%max_dist_AO_atoms)
     715              : 
     716          488 :       CALL timestop(handle)
     717              : 
     718          488 :    END SUBROUTINE check_dist
     719              : 
     720              : ! **************************************************************************************************
     721              : !> \brief ...
     722              : !> \param bs_env ...
     723              : !> \param qs_env ...
     724              : !> \param mat_chi_Gamma_tau ...
     725              : !> \param fm_W_MIC_time ...
     726              : ! **************************************************************************************************
     727           22 :    SUBROUTINE get_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
     728              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     729              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     730              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_Gamma_tau
     731              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
     732              : 
     733              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_W_MIC'
     734              : 
     735              :       INTEGER                                            :: handle
     736              : 
     737           22 :       CALL timeset(routineN, handle)
     738              : 
     739           22 :       IF (bs_env%all_W_exist) THEN
     740            6 :          CALL read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
     741              :       ELSE
     742           16 :          CALL compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
     743              :       END IF
     744              : 
     745           22 :       CALL timestop(handle)
     746              : 
     747           22 :    END SUBROUTINE get_W_MIC
     748              : 
     749              : ! **************************************************************************************************
     750              : !> \brief ...
     751              : !> \param bs_env ...
     752              : !> \param qs_env ...
     753              : !> \param fm_V_kp ...
     754              : !> \param ikp_batch ...
     755              : ! **************************************************************************************************
     756           64 :    SUBROUTINE compute_V_k_by_lattice_sum(bs_env, qs_env, fm_V_kp, ikp_batch)
     757              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     758              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     759              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_V_kp
     760              :       INTEGER                                            :: ikp_batch
     761              : 
     762              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_k_by_lattice_sum'
     763              : 
     764              :       INTEGER                                            :: handle, ikp, ikp_end, ikp_start, &
     765              :                                                             nkp_chi_eps_W_batch, re_im
     766           64 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     767              :       TYPE(cell_type), POINTER                           :: cell
     768           64 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_V_kp
     769           64 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     770           64 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     771              : 
     772           64 :       CALL timeset(routineN, handle)
     773              : 
     774           64 :       nkp_chi_eps_W_batch = bs_env%nkp_chi_eps_W_batch
     775              : 
     776           64 :       ikp_start = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + 1
     777           64 :       ikp_end = MIN(ikp_batch*bs_env%nkp_chi_eps_W_batch, bs_env%kpoints_chi_eps_W%nkp)
     778              : 
     779           64 :       NULLIFY (mat_V_kp)
     780          816 :       ALLOCATE (mat_V_kp(ikp_start:ikp_end, 2))
     781              : 
     782          280 :       DO ikp = ikp_start, ikp_end
     783          712 :          DO re_im = 1, 2
     784          432 :             NULLIFY (mat_V_kp(ikp, re_im)%matrix)
     785          432 :             ALLOCATE (mat_V_kp(ikp, re_im)%matrix)
     786          432 :             CALL dbcsr_create(mat_V_kp(ikp, re_im)%matrix, template=bs_env%mat_RI_RI%matrix)
     787          432 :             CALL dbcsr_reserve_all_blocks(mat_V_kp(ikp, re_im)%matrix)
     788          648 :             CALL dbcsr_set(mat_V_kp(ikp, re_im)%matrix, 0.0_dp)
     789              : 
     790              :          END DO ! re_im
     791              :       END DO ! ikp
     792              : 
     793              :       CALL get_qs_env(qs_env=qs_env, &
     794              :                       particle_set=particle_set, &
     795              :                       cell=cell, &
     796              :                       qs_kind_set=qs_kind_set, &
     797           64 :                       atomic_kind_set=atomic_kind_set)
     798              : 
     799           64 :       IF (ikp_end .LE. bs_env%nkp_chi_eps_W_orig) THEN
     800              : 
     801              :          ! 1. 2c Coulomb integrals for the first "original" k-point grid
     802           96 :          bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
     803              : 
     804           40 :       ELSE IF (ikp_start > bs_env%nkp_chi_eps_W_orig .AND. &
     805              :                ikp_end .LE. bs_env%nkp_chi_eps_W_orig_plus_extra) THEN
     806              : 
     807              :          ! 2. 2c Coulomb integrals for the second "extrapolation" k-point grid
     808          160 :          bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_extra
     809              : 
     810              :       ELSE
     811              : 
     812            0 :          CPABORT("Error with k-point parallelization.")
     813              : 
     814              :       END IF
     815              : 
     816              :       CALL build_2c_coulomb_matrix_kp(mat_V_kp, &
     817              :                                       bs_env%kpoints_chi_eps_W, &
     818              :                                       basis_type="RI_AUX", &
     819              :                                       cell=cell, &
     820              :                                       particle_set=particle_set, &
     821              :                                       qs_kind_set=qs_kind_set, &
     822              :                                       atomic_kind_set=atomic_kind_set, &
     823              :                                       size_lattice_sum=bs_env%size_lattice_sum_V, &
     824              :                                       operator_type=operator_coulomb, &
     825              :                                       ikp_start=ikp_start, &
     826           64 :                                       ikp_end=ikp_end)
     827              : 
     828          256 :       bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
     829              : 
     830          816 :       ALLOCATE (fm_V_kp(ikp_start:ikp_end, 2))
     831          280 :       DO ikp = ikp_start, ikp_end
     832          712 :          DO re_im = 1, 2
     833          432 :             CALL cp_fm_create(fm_V_kp(ikp, re_im), bs_env%fm_RI_RI%matrix_struct)
     834          432 :             CALL copy_dbcsr_to_fm(mat_V_kp(ikp, re_im)%matrix, fm_V_kp(ikp, re_im))
     835          648 :             CALL dbcsr_deallocate_matrix(mat_V_kp(ikp, re_im)%matrix)
     836              :          END DO
     837              :       END DO
     838           64 :       DEALLOCATE (mat_V_kp)
     839              : 
     840           64 :       CALL timestop(handle)
     841              : 
     842           64 :    END SUBROUTINE compute_V_k_by_lattice_sum
     843              : 
     844              : ! **************************************************************************************************
     845              : !> \brief ...
     846              : !> \param bs_env ...
     847              : !> \param qs_env ...
     848              : !> \param fm_V_kp ...
     849              : !> \param cfm_V_sqrt_ikp ...
     850              : !> \param cfm_M_inv_V_sqrt_ikp ...
     851              : !> \param ikp ...
     852              : ! **************************************************************************************************
     853          216 :    SUBROUTINE compute_MinvVsqrt_Vsqrt(bs_env, qs_env, fm_V_kp, cfm_V_sqrt_ikp, &
     854              :                                       cfm_M_inv_V_sqrt_ikp, ikp)
     855              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     856              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     857              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_V_kp
     858              :       TYPE(cp_cfm_type)                                  :: cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp
     859              :       INTEGER                                            :: ikp
     860              : 
     861              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_MinvVsqrt_Vsqrt'
     862              : 
     863              :       INTEGER                                            :: handle, info, n_RI
     864              :       TYPE(cp_cfm_type)                                  :: cfm_M_inv_ikp, cfm_work
     865          216 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_M_ikp
     866              : 
     867          216 :       CALL timeset(routineN, handle)
     868              : 
     869          216 :       n_RI = bs_env%n_RI
     870              : 
     871              :       ! get here M(k) and write it to fm_M_ikp
     872              :       CALL RI_2c_integral_mat(qs_env, fm_M_ikp, fm_V_kp(ikp, 1), &
     873              :                               n_RI, bs_env%ri_metric, do_kpoints=.TRUE., &
     874              :                               kpoints=bs_env%kpoints_chi_eps_W, &
     875              :                               regularization_RI=bs_env%regularization_RI, ikp_ext=ikp, &
     876          216 :                               do_build_cell_index=(ikp == 1))
     877              : 
     878          216 :       IF (ikp == 1) THEN
     879           16 :          CALL cp_cfm_create(cfm_V_sqrt_ikp, fm_V_kp(ikp, 1)%matrix_struct)
     880           16 :          CALL cp_cfm_create(cfm_M_inv_V_sqrt_ikp, fm_V_kp(ikp, 1)%matrix_struct)
     881              :       END IF
     882          216 :       CALL cp_cfm_create(cfm_M_inv_ikp, fm_V_kp(ikp, 1)%matrix_struct)
     883              : 
     884          216 :       CALL cp_fm_to_cfm(fm_M_ikp(1, 1), fm_M_ikp(1, 2), cfm_M_inv_ikp)
     885          216 :       CALL cp_fm_to_cfm(fm_V_kp(ikp, 1), fm_V_kp(ikp, 2), cfm_V_sqrt_ikp)
     886              : 
     887          216 :       CALL cp_fm_release(fm_M_ikp)
     888              : 
     889          216 :       CALL cp_cfm_create(cfm_work, fm_V_kp(ikp, 1)%matrix_struct)
     890              : 
     891              :       ! M(k) -> M^-1(k)
     892          216 :       CALL cp_cfm_to_cfm(cfm_M_inv_ikp, cfm_work)
     893          216 :       CALL cp_cfm_cholesky_decompose(matrix=cfm_M_inv_ikp, n=n_RI, info_out=info)
     894          216 :       IF (info == 0) THEN
     895              :          ! successful Cholesky decomposition
     896          216 :          CALL cp_cfm_cholesky_invert(cfm_M_inv_ikp)
     897              :          ! symmetrize the result
     898          216 :          CALL cp_cfm_uplo_to_full(cfm_M_inv_ikp)
     899              :       ELSE
     900              :          ! Cholesky decomposition not successful: use expensive diagonalization
     901            0 :          CALL cp_cfm_power(cfm_work, threshold=bs_env%eps_eigval_mat_RI, exponent=-1.0_dp)
     902            0 :          CALL cp_cfm_to_cfm(cfm_work, cfm_M_inv_ikp)
     903              :       END IF
     904              : 
     905              :       ! V(k) -> L(k) with L^H(k)*L(k) = V(k) [L(k) can be just considered to be V^0.5(k)]
     906          216 :       CALL cp_cfm_to_cfm(cfm_V_sqrt_ikp, cfm_work)
     907          216 :       CALL cp_cfm_cholesky_decompose(matrix=cfm_V_sqrt_ikp, n=n_RI, info_out=info)
     908          216 :       IF (info == 0) THEN
     909              :          ! successful Cholesky decomposition
     910          216 :          CALL clean_lower_part(cfm_V_sqrt_ikp)
     911              :       ELSE
     912              :          ! Cholesky decomposition not successful: use expensive diagonalization
     913            0 :          CALL cp_cfm_power(cfm_work, threshold=0.0_dp, exponent=0.5_dp)
     914            0 :          CALL cp_cfm_to_cfm(cfm_work, cfm_V_sqrt_ikp)
     915              :       END IF
     916          216 :       CALL cp_cfm_release(cfm_work)
     917              : 
     918              :       ! get M^-1(k)*V^0.5(k)
     919              :       CALL parallel_gemm("N", "C", n_RI, n_RI, n_RI, z_one, cfm_M_inv_ikp, cfm_V_sqrt_ikp, &
     920          216 :                          z_zero, cfm_M_inv_V_sqrt_ikp)
     921              : 
     922          216 :       CALL cp_cfm_release(cfm_M_inv_ikp)
     923              : 
     924          216 :       CALL timestop(handle)
     925              : 
     926          432 :    END SUBROUTINE compute_MinvVsqrt_Vsqrt
     927              : 
     928              : ! **************************************************************************************************
     929              : !> \brief ...
     930              : !> \param bs_env ...
     931              : !> \param mat_chi_Gamma_tau ...
     932              : !> \param fm_W_MIC_time ...
     933              : ! **************************************************************************************************
     934            6 :    SUBROUTINE read_W_MIC_time(bs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
     935              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     936              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_Gamma_tau
     937              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
     938              : 
     939              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'read_W_MIC_time'
     940              : 
     941              :       INTEGER                                            :: handle, i_t
     942              :       REAL(KIND=dp)                                      :: t1
     943              : 
     944            6 :       CALL timeset(routineN, handle)
     945              : 
     946            6 :       CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
     947            6 :       CALL create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
     948              : 
     949          106 :       DO i_t = 1, bs_env%num_time_freq_points
     950              : 
     951          100 :          t1 = m_walltime()
     952              : 
     953          100 :          CALL fm_read(fm_W_MIC_time(i_t), bs_env, bs_env%W_time_name, i_t)
     954              : 
     955          106 :          IF (bs_env%unit_nr > 0) THEN
     956              :             WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') &
     957           50 :                'Read W^MIC(iτ) from file for time point  ', i_t, ' /', bs_env%num_time_freq_points, &
     958          100 :                ',    Execution time', m_walltime() - t1, ' s'
     959              :          END IF
     960              : 
     961              :       END DO
     962              : 
     963            6 :       IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
     964              : 
     965              :       ! Marek : Reading of the W(w=0) potential for RTP
     966              :       ! TODO : is the condition bs_env%all_W_exist sufficient for reading?
     967            6 :       IF (bs_env%rtp_method == rtp_method_bse) THEN
     968            4 :          CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
     969            4 :          t1 = m_walltime()
     970            4 :          CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0)
     971            4 :          IF (bs_env%unit_nr > 0) THEN
     972              :             WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F7.1,A)') &
     973            2 :                'Read W^MIC(f=0) from file for freq. point  ', 1, ' /', 1, &
     974            4 :                ',    Execution time', m_walltime() - t1, ' s'
     975              :          END IF
     976              :       END IF
     977              : 
     978            6 :       CALL timestop(handle)
     979              : 
     980            6 :    END SUBROUTINE read_W_MIC_time
     981              : 
     982              : ! **************************************************************************************************
     983              : !> \brief ...
     984              : !> \param bs_env ...
     985              : !> \param qs_env ...
     986              : !> \param mat_chi_Gamma_tau ...
     987              : !> \param fm_W_MIC_time ...
     988              : ! **************************************************************************************************
     989           16 :    SUBROUTINE compute_W_MIC(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_time)
     990              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     991              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     992              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_Gamma_tau
     993              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
     994              : 
     995              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_W_MIC'
     996              : 
     997              :       INTEGER                                            :: handle, i_t, ikp, ikp_batch, &
     998              :                                                             ikp_in_batch, j_w
     999              :       REAL(KIND=dp)                                      :: t1
    1000              :       TYPE(cp_cfm_type)                                  :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
    1001           16 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_V_kp
    1002              : 
    1003           16 :       CALL timeset(routineN, handle)
    1004              : 
    1005           16 :       CALL create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
    1006              : 
    1007           80 :       DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
    1008              : 
    1009           64 :          t1 = m_walltime()
    1010              : 
    1011              :          ! Compute V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
    1012           64 :          CALL compute_V_k_by_lattice_sum(bs_env, qs_env, fm_V_kp, ikp_batch)
    1013              : 
    1014          320 :          DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
    1015              : 
    1016          256 :             ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
    1017              : 
    1018          256 :             IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) CYCLE
    1019              : 
    1020              :             CALL compute_MinvVsqrt_Vsqrt(bs_env, qs_env, fm_V_kp, &
    1021          216 :                                          cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp, ikp)
    1022              : 
    1023          216 :             CALL bs_env%para_env%sync()
    1024              : 
    1025         2040 :             DO j_w = 1, bs_env%num_time_freq_points
    1026              : 
    1027              :                ! check if we need this (ikp, ω_j) combination for approximate k-point extrapolation
    1028         1824 :                IF (bs_env%approx_kp_extrapol .AND. j_w > 1 .AND. &
    1029              :                    ikp > bs_env%nkp_chi_eps_W_orig) CYCLE
    1030              : 
    1031              :                CALL compute_fm_W_MIC_freq_j(bs_env, qs_env, bs_env%fm_W_MIC_freq, j_w, ikp, &
    1032              :                                             mat_chi_Gamma_tau, cfm_M_inv_V_sqrt_ikp, &
    1033         1500 :                                             cfm_V_sqrt_ikp)
    1034              : 
    1035              :                ! Fourier trafo from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
    1036         2040 :                CALL Fourier_transform_w_to_t(bs_env, fm_W_MIC_time, bs_env%fm_W_MIC_freq, j_w)
    1037              : 
    1038              :             END DO ! ω_j
    1039              : 
    1040          216 :             CALL cp_fm_release(fm_V_kp(ikp, 1))
    1041          320 :             CALL cp_fm_release(fm_V_kp(ikp, 2))
    1042              : 
    1043              :          END DO ! ikp_in_batch
    1044              : 
    1045           64 :          DEALLOCATE (fm_V_kp)
    1046              : 
    1047           80 :          IF (bs_env%unit_nr > 0) THEN
    1048              :             WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F7.1,A)') &
    1049           32 :                'Computed W(iτ,k) for k-point batch', &
    1050           32 :                ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
    1051           64 :                ',    Execution time', m_walltime() - t1, ' s'
    1052              :          END IF
    1053              : 
    1054              :       END DO ! ikp_batch
    1055              : 
    1056           16 :       IF (bs_env%approx_kp_extrapol) THEN
    1057            2 :          CALL apply_extrapol_factor(bs_env, fm_W_MIC_time)
    1058              :       END IF
    1059              : 
    1060              :       ! M^-1(k=0)*W^MIC(iτ)*M^-1(k=0)
    1061           16 :       CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_W_MIC_time)
    1062              : 
    1063          240 :       DO i_t = 1, bs_env%num_time_freq_points
    1064          240 :          CALL fm_write(fm_W_MIC_time(i_t), i_t, bs_env%W_time_name, qs_env)
    1065              :       END DO
    1066              : 
    1067           16 :       CALL cp_cfm_release(cfm_M_inv_V_sqrt_ikp)
    1068           16 :       CALL cp_cfm_release(cfm_V_sqrt_ikp)
    1069           16 :       CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
    1070              : 
    1071              :       ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
    1072           16 :       IF (bs_env%rtp_method == rtp_method_bse) THEN
    1073            8 :          t1 = m_walltime()
    1074            8 :          CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
    1075              :          ! Set to zero
    1076            8 :          CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
    1077              :          ! Sum over all times
    1078          168 :          DO i_t = 1, bs_env%num_time_freq_points
    1079              :             ! Add the relevant structure with correct weight
    1080              :             CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
    1081          168 :                                      bs_env%imag_time_weights_freq_zero(i_t), fm_W_MIC_time(i_t))
    1082              :          END DO
    1083              :          ! Done, save to file
    1084            8 :          CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
    1085              :          ! Report calculation
    1086            8 :          IF (bs_env%unit_nr > 0) THEN
    1087              :             WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F7.1,A)') &
    1088            4 :                'Computed W(f=0,k) for k-point batch', &
    1089            4 :                1, ' /', 1, &
    1090            8 :                ',    Execution time', m_walltime() - t1, ' s'
    1091              :          END IF
    1092              :       END IF
    1093              : 
    1094           16 :       IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
    1095              : 
    1096           16 :       CALL timestop(handle)
    1097              : 
    1098           32 :    END SUBROUTINE compute_W_MIC
    1099              : 
    1100              : ! **************************************************************************************************
    1101              : !> \brief ...
    1102              : !> \param bs_env ...
    1103              : !> \param qs_env ...
    1104              : !> \param fm_W_MIC_freq_j ...
    1105              : !> \param j_w ...
    1106              : !> \param ikp ...
    1107              : !> \param mat_chi_Gamma_tau ...
    1108              : !> \param cfm_M_inv_V_sqrt_ikp ...
    1109              : !> \param cfm_V_sqrt_ikp ...
    1110              : ! **************************************************************************************************
    1111         1500 :    SUBROUTINE compute_fm_W_MIC_freq_j(bs_env, qs_env, fm_W_MIC_freq_j, j_w, ikp, mat_chi_Gamma_tau, &
    1112              :                                       cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
    1113              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1114              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1115              :       TYPE(cp_fm_type)                                   :: fm_W_MIC_freq_j
    1116              :       INTEGER                                            :: j_w, ikp
    1117              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_Gamma_tau
    1118              :       TYPE(cp_cfm_type)                                  :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
    1119              : 
    1120              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_MIC_freq_j'
    1121              : 
    1122              :       INTEGER                                            :: handle
    1123              :       TYPE(cp_cfm_type)                                  :: cfm_chi_ikp_freq_j, cfm_W_ikp_freq_j
    1124              : 
    1125         1500 :       CALL timeset(routineN, handle)
    1126              : 
    1127              :       ! 1. Fourier transformation of χ_PQ(iτ,k=0) to χ_PQ(iω_j,k=0)
    1128         1500 :       CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
    1129              : 
    1130         1500 :       CALL cp_fm_set_all(fm_W_MIC_freq_j, 0.0_dp)
    1131              : 
    1132              :       ! 2. Get χ_PQ(iω_j,k_i) from χ_PQ(iω_j,k=0) using the minimum image convention
    1133              :       CALL cfm_ikp_from_fm_Gamma(cfm_chi_ikp_freq_j, bs_env%fm_chi_Gamma_freq, &
    1134         1500 :                                  ikp, qs_env, bs_env%kpoints_chi_eps_W, "RI_AUX")
    1135              : 
    1136              :       ! 3. Remove all negative eigenvalues from χ_PQ(iω_j,k_i)
    1137         1500 :       CALL cp_cfm_power(cfm_chi_ikp_freq_j, threshold=0.0_dp, exponent=1.0_dp)
    1138              : 
    1139              :       ! 4. ε(iω_j,k_i) = Id - V^0.5(k_i)*M^-1(k_i)*χ(iω_j,k_i)*M^-1(k_i)*V^0.5(k_i)
    1140              :       !    W(iω_j,k_i) = V^0.5(k_i)*(ε^-1(iω_j,k_i)-Id)*V^0.5(k_i)
    1141              :       CALL compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
    1142         1500 :                                     cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j)
    1143              : 
    1144              :       ! 5. k-point integration W_PQ(iω_j, k_i) to W_PQ^MIC(iω_j)
    1145         1500 :       SELECT CASE (bs_env%approx_kp_extrapol)
    1146              :       CASE (.FALSE.)
    1147              :          ! default: standard k-point extrapolation
    1148              :          CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, cfm_W_ikp_freq_j, ikp, &
    1149         1500 :                                         bs_env%kpoints_chi_eps_W, "RI_AUX")
    1150              :       CASE (.TRUE.)
    1151              :          ! for approximate kpoint extrapolation: get W_PQ^MIC(iω_1) with and without k-point
    1152              :          ! extrapolation to compute the extrapolation factor f_PQ for every PQ-matrix element,
    1153              :          ! f_PQ = (W_PQ^MIC(iω_1) with extrapolation) / (W_PQ^MIC(iω_1) without extrapolation)
    1154              : 
    1155              :          ! for ω_1, we compute the k-point extrapolated result using all k-points
    1156          196 :          IF (j_w == 1) THEN
    1157              : 
    1158              :             ! k-point extrapolated
    1159              :             CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_extra, &
    1160              :                                            cfm_W_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, &
    1161           52 :                                            "RI_AUX")
    1162              :             ! non-kpoint extrapolated
    1163           52 :             IF (ikp .LE. bs_env%nkp_chi_eps_W_orig) THEN
    1164              :                CALL MIC_contribution_from_ikp(bs_env, qs_env, bs_env%fm_W_MIC_freq_1_no_extra, &
    1165              :                                               cfm_W_ikp_freq_j, ikp, bs_env%kpoints_chi_eps_W, &
    1166           16 :                                               "RI_AUX", wkp_ext=bs_env%wkp_orig)
    1167              :             END IF
    1168              : 
    1169              :          END IF
    1170              : 
    1171              :          ! for all ω_j, we need to compute W^MIC without k-point extrpolation
    1172          196 :          IF (ikp .LE. bs_env%nkp_chi_eps_W_orig) THEN
    1173              :             CALL MIC_contribution_from_ikp(bs_env, qs_env, fm_W_MIC_freq_j, cfm_W_ikp_freq_j, &
    1174              :                                            ikp, bs_env%kpoints_chi_eps_W, "RI_AUX", &
    1175          160 :                                            wkp_ext=bs_env%wkp_orig)
    1176              :          END IF
    1177              :       END SELECT
    1178              : 
    1179         1500 :       CALL cp_cfm_release(cfm_W_ikp_freq_j)
    1180              : 
    1181         1500 :       CALL timestop(handle)
    1182              : 
    1183         1500 :    END SUBROUTINE compute_fm_W_MIC_freq_j
    1184              : 
    1185              : ! **************************************************************************************************
    1186              : !> \brief ...
    1187              : !> \param cfm_mat ...
    1188              : ! **************************************************************************************************
    1189          432 :    SUBROUTINE clean_lower_part(cfm_mat)
    1190              :       TYPE(cp_cfm_type)                                  :: cfm_mat
    1191              : 
    1192              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'clean_lower_part'
    1193              : 
    1194              :       INTEGER                                            :: handle, i_global, i_row, j_col, &
    1195              :                                                             j_global, ncol_local, nrow_local
    1196          216 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1197              : 
    1198          216 :       CALL timeset(routineN, handle)
    1199              : 
    1200              :       CALL cp_cfm_get_info(matrix=cfm_mat, &
    1201              :                            nrow_local=nrow_local, ncol_local=ncol_local, &
    1202          216 :                            row_indices=row_indices, col_indices=col_indices)
    1203              : 
    1204          980 :       DO i_row = 1, nrow_local
    1205         7544 :          DO j_col = 1, ncol_local
    1206         6564 :             i_global = row_indices(i_row)
    1207         6564 :             j_global = col_indices(j_col)
    1208         7328 :             IF (j_global < i_global) cfm_mat%local_data(i_row, j_col) = z_zero
    1209              :          END DO
    1210              :       END DO
    1211              : 
    1212          216 :       CALL timestop(handle)
    1213              : 
    1214          216 :    END SUBROUTINE clean_lower_part
    1215              : 
    1216              : ! **************************************************************************************************
    1217              : !> \brief ...
    1218              : !> \param bs_env ...
    1219              : !> \param fm_W_MIC_time ...
    1220              : ! **************************************************************************************************
    1221            4 :    SUBROUTINE apply_extrapol_factor(bs_env, fm_W_MIC_time)
    1222              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1223              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
    1224              : 
    1225              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_extrapol_factor'
    1226              : 
    1227              :       INTEGER                                            :: handle, i, i_t, j, ncol_local, nrow_local
    1228              :       REAL(KIND=dp)                                      :: extrapol_factor, W_extra_1, W_no_extra_1
    1229              : 
    1230            2 :       CALL timeset(routineN, handle)
    1231              : 
    1232            2 :       CALL cp_fm_get_info(matrix=fm_W_MIC_time(1), nrow_local=nrow_local, ncol_local=ncol_local)
    1233              : 
    1234           22 :       DO i_t = 1, bs_env%num_time_freq_points
    1235           72 :          DO i = 1, nrow_local
    1236          320 :             DO j = 1, ncol_local
    1237              : 
    1238          250 :                W_extra_1 = bs_env%fm_W_MIC_freq_1_extra%local_data(i, j)
    1239          250 :                W_no_extra_1 = bs_env%fm_W_MIC_freq_1_no_extra%local_data(i, j)
    1240              : 
    1241          250 :                IF (ABS(W_no_extra_1) > 1.0E-13) THEN
    1242          190 :                   extrapol_factor = W_extra_1/W_no_extra_1
    1243              :                ELSE
    1244              :                   extrapol_factor = 1.0_dp
    1245              :                END IF
    1246              : 
    1247              :                ! reset extrapolation factor if it is very large
    1248          250 :                IF (ABS(extrapol_factor) > 10.0_dp) extrapol_factor = 1.0_dp
    1249              : 
    1250              :                fm_W_MIC_time(i_t)%local_data(i, j) = fm_W_MIC_time(i_t)%local_data(i, j) &
    1251          300 :                                                      *extrapol_factor
    1252              :             END DO
    1253              :          END DO
    1254              :       END DO
    1255              : 
    1256            2 :       CALL timestop(handle)
    1257              : 
    1258            2 :    END SUBROUTINE apply_extrapol_factor
    1259              : 
    1260              : ! **************************************************************************************************
    1261              : !> \brief ...
    1262              : !> \param bs_env ...
    1263              : !> \param fm_chi_Gamma_freq ...
    1264              : !> \param j_w ...
    1265              : !> \param mat_chi_Gamma_tau ...
    1266              : ! **************************************************************************************************
    1267         1500 :    SUBROUTINE compute_fm_chi_Gamma_freq(bs_env, fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
    1268              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1269              :       TYPE(cp_fm_type)                                   :: fm_chi_Gamma_freq
    1270              :       INTEGER                                            :: j_w
    1271              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_Gamma_tau
    1272              : 
    1273              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_chi_Gamma_freq'
    1274              : 
    1275              :       INTEGER                                            :: handle, i_t
    1276              :       REAL(KIND=dp)                                      :: freq_j, time_i, weight_ij
    1277              : 
    1278         1500 :       CALL timeset(routineN, handle)
    1279              : 
    1280         1500 :       CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
    1281              : 
    1282         1500 :       freq_j = bs_env%imag_freq_points(j_w)
    1283              : 
    1284        15604 :       DO i_t = 1, bs_env%num_time_freq_points
    1285              : 
    1286        14104 :          time_i = bs_env%imag_time_points(i_t)
    1287        14104 :          weight_ij = bs_env%weights_cos_t_to_w(j_w, i_t)
    1288              : 
    1289              :          ! actual Fourier transform
    1290              :          CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_Gamma_tau(i_t)%matrix, &
    1291        15604 :                         1.0_dp, COS(time_i*freq_j)*weight_ij)
    1292              : 
    1293              :       END DO
    1294              : 
    1295         1500 :       CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, fm_chi_Gamma_freq)
    1296              : 
    1297         1500 :       CALL timestop(handle)
    1298              : 
    1299         1500 :    END SUBROUTINE compute_fm_chi_Gamma_freq
    1300              : 
    1301              : ! **************************************************************************************************
    1302              : !> \brief ...
    1303              : !> \param mat_ikp_re ...
    1304              : !> \param mat_ikp_im ...
    1305              : !> \param mat_Gamma ...
    1306              : !> \param kpoints ...
    1307              : !> \param ikp ...
    1308              : !> \param qs_env ...
    1309              : ! **************************************************************************************************
    1310            0 :    SUBROUTINE mat_ikp_from_mat_Gamma(mat_ikp_re, mat_ikp_im, mat_Gamma, kpoints, ikp, qs_env)
    1311              :       TYPE(dbcsr_type)                                   :: mat_ikp_re, mat_ikp_im, mat_Gamma
    1312              :       TYPE(kpoint_type), POINTER                         :: kpoints
    1313              :       INTEGER                                            :: ikp
    1314              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1315              : 
    1316              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'mat_ikp_from_mat_Gamma'
    1317              : 
    1318              :       INTEGER                                            :: col, handle, i_cell, j_cell, num_cells, &
    1319              :                                                             row
    1320            0 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
    1321              :       LOGICAL :: f, i_cell_is_the_minimum_image_cell
    1322              :       REAL(KIND=dp)                                      :: abs_rab_cell_i, abs_rab_cell_j, arg
    1323              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_vector, cell_vector_j, rab_cell_i, &
    1324              :                                                             rab_cell_j
    1325              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
    1326            0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: block_im, block_re, data_block
    1327              :       TYPE(cell_type), POINTER                           :: cell
    1328              :       TYPE(dbcsr_iterator_type)                          :: iter
    1329            0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1330              : 
    1331            0 :       CALL timeset(routineN, handle)
    1332              : 
    1333              :       ! get the same blocks in mat_ikp_re and mat_ikp_im as in mat_Gamma
    1334            0 :       CALL dbcsr_copy(mat_ikp_re, mat_Gamma)
    1335            0 :       CALL dbcsr_copy(mat_ikp_im, mat_Gamma)
    1336            0 :       CALL dbcsr_set(mat_ikp_re, 0.0_dp)
    1337            0 :       CALL dbcsr_set(mat_ikp_im, 0.0_dp)
    1338              : 
    1339            0 :       NULLIFY (cell, particle_set)
    1340            0 :       CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
    1341            0 :       CALL get_cell(cell=cell, h=hmat)
    1342              : 
    1343            0 :       index_to_cell => kpoints%index_to_cell
    1344              : 
    1345            0 :       num_cells = SIZE(index_to_cell, 2)
    1346              : 
    1347            0 :       DO i_cell = 1, num_cells
    1348              : 
    1349            0 :          CALL dbcsr_iterator_start(iter, mat_Gamma)
    1350            0 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
    1351            0 :             CALL dbcsr_iterator_next_block(iter, row, col, data_block)
    1352              : 
    1353            0 :             cell_vector(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, i_cell), dp))
    1354              : 
    1355              :             rab_cell_i(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
    1356            0 :                               (pbc(particle_set(col)%r(1:3), cell) + cell_vector(1:3))
    1357            0 :             abs_rab_cell_i = SQRT(rab_cell_i(1)**2 + rab_cell_i(2)**2 + rab_cell_i(3)**2)
    1358              : 
    1359              :             ! minimum image convention
    1360            0 :             i_cell_is_the_minimum_image_cell = .TRUE.
    1361            0 :             DO j_cell = 1, num_cells
    1362            0 :                cell_vector_j(1:3) = MATMUL(hmat, REAL(index_to_cell(1:3, j_cell), dp))
    1363              :                rab_cell_j(1:3) = pbc(particle_set(row)%r(1:3), cell) - &
    1364            0 :                                  (pbc(particle_set(col)%r(1:3), cell) + cell_vector_j(1:3))
    1365            0 :                abs_rab_cell_j = SQRT(rab_cell_j(1)**2 + rab_cell_j(2)**2 + rab_cell_j(3)**2)
    1366              : 
    1367            0 :                IF (abs_rab_cell_i > abs_rab_cell_j + 1.0E-6_dp) THEN
    1368            0 :                   i_cell_is_the_minimum_image_cell = .FALSE.
    1369              :                END IF
    1370              :             END DO
    1371              : 
    1372            0 :             IF (i_cell_is_the_minimum_image_cell) THEN
    1373            0 :                NULLIFY (block_re, block_im)
    1374            0 :                CALL dbcsr_get_block_p(matrix=mat_ikp_re, row=row, col=col, block=block_re, found=f)
    1375            0 :                CALL dbcsr_get_block_p(matrix=mat_ikp_im, row=row, col=col, block=block_im, found=f)
    1376            0 :                CPASSERT(ALL(ABS(block_re) < 1.0E-10_dp))
    1377            0 :                CPASSERT(ALL(ABS(block_im) < 1.0E-10_dp))
    1378              : 
    1379              :                arg = REAL(index_to_cell(1, i_cell), dp)*kpoints%xkp(1, ikp) + &
    1380              :                      REAL(index_to_cell(2, i_cell), dp)*kpoints%xkp(2, ikp) + &
    1381            0 :                      REAL(index_to_cell(3, i_cell), dp)*kpoints%xkp(3, ikp)
    1382              : 
    1383            0 :                block_re(:, :) = COS(twopi*arg)*data_block(:, :)
    1384            0 :                block_im(:, :) = SIN(twopi*arg)*data_block(:, :)
    1385              :             END IF
    1386              : 
    1387              :          END DO
    1388            0 :          CALL dbcsr_iterator_stop(iter)
    1389              : 
    1390              :       END DO
    1391              : 
    1392            0 :       CALL timestop(handle)
    1393              : 
    1394            0 :    END SUBROUTINE mat_ikp_from_mat_Gamma
    1395              : 
    1396              : ! **************************************************************************************************
    1397              : !> \brief ...
    1398              : !> \param bs_env ...
    1399              : !> \param cfm_chi_ikp_freq_j ...
    1400              : !> \param cfm_V_sqrt_ikp ...
    1401              : !> \param cfm_M_inv_V_sqrt_ikp ...
    1402              : !> \param cfm_W_ikp_freq_j ...
    1403              : ! **************************************************************************************************
    1404         7500 :    SUBROUTINE compute_cfm_W_ikp_freq_j(bs_env, cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
    1405              :                                        cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j)
    1406              : 
    1407              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1408              :       TYPE(cp_cfm_type)                                  :: cfm_chi_ikp_freq_j, cfm_V_sqrt_ikp, &
    1409              :                                                             cfm_M_inv_V_sqrt_ikp, cfm_W_ikp_freq_j
    1410              : 
    1411              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_cfm_W_ikp_freq_j'
    1412              : 
    1413              :       INTEGER                                            :: handle, info, n_RI
    1414              :       TYPE(cp_cfm_type)                                  :: cfm_eps_ikp_freq_j, cfm_work
    1415              : 
    1416         1500 :       CALL timeset(routineN, handle)
    1417              : 
    1418         1500 :       CALL cp_cfm_create(cfm_work, cfm_chi_ikp_freq_j%matrix_struct)
    1419         1500 :       n_RI = bs_env%n_RI
    1420              : 
    1421              :       ! 1. ε(iω_j,k) = Id - V^0.5(k)*M^-1(k)*χ(iω_j,k)*M^-1(k)*V^0.5(k)
    1422              : 
    1423              :       ! 1. a) work = χ(iω_j,k)*M^-1(k)*V^0.5(k)
    1424              :       CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, z_one, &
    1425         1500 :                          cfm_chi_ikp_freq_j, cfm_M_inv_V_sqrt_ikp, z_zero, cfm_work)
    1426         1500 :       CALL cp_cfm_release(cfm_chi_ikp_freq_j)
    1427              : 
    1428              :       ! 1. b) eps_work = V^0.5(k)*M^-1(k)*work
    1429         1500 :       CALL cp_cfm_create(cfm_eps_ikp_freq_j, cfm_work%matrix_struct)
    1430              :       CALL parallel_gemm('C', 'N', n_RI, n_RI, n_RI, z_one, &
    1431         1500 :                          cfm_M_inv_V_sqrt_ikp, cfm_work, z_zero, cfm_eps_ikp_freq_j)
    1432              : 
    1433              :       ! 1. c) ε(iω_j,k) = eps_work - Id
    1434         1500 :       CALL cfm_add_on_diag(cfm_eps_ikp_freq_j, z_one)
    1435              : 
    1436              :       ! 2. W(iω_j,k) = V^0.5(k)*(ε^-1(iω_j,k)-Id)*V^0.5(k)
    1437              : 
    1438              :       ! 2. a) Cholesky decomposition of ε(iω_j,k) as preparation for inversion
    1439         1500 :       CALL cp_cfm_cholesky_decompose(matrix=cfm_eps_ikp_freq_j, n=n_RI, info_out=info)
    1440         1500 :       CPASSERT(info == 0)
    1441              : 
    1442              :       ! 2. b) Inversion of ε(iω_j,k) using its Cholesky decomposition
    1443         1500 :       CALL cp_cfm_cholesky_invert(cfm_eps_ikp_freq_j)
    1444         1500 :       CALL cp_cfm_uplo_to_full(cfm_eps_ikp_freq_j)
    1445              : 
    1446              :       ! 2. c) ε^-1(iω_j,k)-Id
    1447         1500 :       CALL cfm_add_on_diag(cfm_eps_ikp_freq_j, -z_one)
    1448              : 
    1449              :       ! 2. d) work = (ε^-1(iω_j,k)-Id)*V^0.5(k)
    1450              :       CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, z_one, cfm_eps_ikp_freq_j, cfm_V_sqrt_ikp, &
    1451         1500 :                          z_zero, cfm_work)
    1452              : 
    1453              :       ! 2. e) W(iw,k) = V^0.5(k)*work
    1454         1500 :       CALL cp_cfm_create(cfm_W_ikp_freq_j, cfm_work%matrix_struct)
    1455              :       CALL parallel_gemm('C', 'N', n_RI, n_RI, n_RI, z_one, cfm_V_sqrt_ikp, cfm_work, &
    1456         1500 :                          z_zero, cfm_W_ikp_freq_j)
    1457              : 
    1458         1500 :       CALL cp_cfm_release(cfm_work)
    1459         1500 :       CALL cp_cfm_release(cfm_eps_ikp_freq_j)
    1460              : 
    1461         1500 :       CALL timestop(handle)
    1462              : 
    1463         1500 :    END SUBROUTINE compute_cfm_W_ikp_freq_j
    1464              : 
    1465              : ! **************************************************************************************************
    1466              : !> \brief ...
    1467              : !> \param cfm ...
    1468              : !> \param alpha ...
    1469              : ! **************************************************************************************************
    1470         6000 :    SUBROUTINE cfm_add_on_diag(cfm, alpha)
    1471              : 
    1472              :       TYPE(cp_cfm_type)                                  :: cfm
    1473              :       COMPLEX(KIND=dp)                                   :: alpha
    1474              : 
    1475              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cfm_add_on_diag'
    1476              : 
    1477              :       INTEGER                                            :: handle, i_global, i_row, j_col, &
    1478              :                                                             j_global, ncol_local, nrow_local
    1479         3000 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1480              : 
    1481         3000 :       CALL timeset(routineN, handle)
    1482              : 
    1483              :       CALL cp_cfm_get_info(matrix=cfm, &
    1484              :                            nrow_local=nrow_local, &
    1485              :                            ncol_local=ncol_local, &
    1486              :                            row_indices=row_indices, &
    1487         3000 :                            col_indices=col_indices)
    1488              : 
    1489              :       ! add 1 on the diagonal
    1490        27184 :       DO j_col = 1, ncol_local
    1491        24184 :          j_global = col_indices(j_col)
    1492       162460 :          DO i_row = 1, nrow_local
    1493       135276 :             i_global = row_indices(i_row)
    1494       159460 :             IF (j_global == i_global) THEN
    1495        12092 :                cfm%local_data(i_row, j_col) = cfm%local_data(i_row, j_col) + alpha
    1496              :             END IF
    1497              :          END DO
    1498              :       END DO
    1499              : 
    1500         3000 :       CALL timestop(handle)
    1501              : 
    1502         3000 :    END SUBROUTINE cfm_add_on_diag
    1503              : 
    1504              : ! **************************************************************************************************
    1505              : !> \brief ...
    1506              : !> \param bs_env ...
    1507              : !> \param fm_W_MIC_time ...
    1508              : ! **************************************************************************************************
    1509           22 :    SUBROUTINE create_fm_W_MIC_time(bs_env, fm_W_MIC_time)
    1510              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1511              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
    1512              : 
    1513              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_fm_W_MIC_time'
    1514              : 
    1515              :       INTEGER                                            :: handle, i_t
    1516              : 
    1517           22 :       CALL timeset(routineN, handle)
    1518              : 
    1519          390 :       ALLOCATE (fm_W_MIC_time(bs_env%num_time_freq_points))
    1520          346 :       DO i_t = 1, bs_env%num_time_freq_points
    1521          346 :          CALL cp_fm_create(fm_W_MIC_time(i_t), bs_env%fm_RI_RI%matrix_struct)
    1522              :       END DO
    1523              : 
    1524           22 :       CALL timestop(handle)
    1525              : 
    1526           22 :    END SUBROUTINE create_fm_W_MIC_time
    1527              : 
    1528              : ! **************************************************************************************************
    1529              : !> \brief ...
    1530              : !> \param bs_env ...
    1531              : !> \param fm_W_MIC_time ...
    1532              : !> \param fm_W_MIC_freq_j ...
    1533              : !> \param j_w ...
    1534              : ! **************************************************************************************************
    1535         1500 :    SUBROUTINE Fourier_transform_w_to_t(bs_env, fm_W_MIC_time, fm_W_MIC_freq_j, j_w)
    1536              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1537              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
    1538              :       TYPE(cp_fm_type)                                   :: fm_W_MIC_freq_j
    1539              :       INTEGER                                            :: j_w
    1540              : 
    1541              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'Fourier_transform_w_to_t'
    1542              : 
    1543              :       INTEGER                                            :: handle, i_t
    1544              :       REAL(KIND=dp)                                      :: freq_j, time_i, weight_ij
    1545              : 
    1546         1500 :       CALL timeset(routineN, handle)
    1547              : 
    1548         1500 :       freq_j = bs_env%imag_freq_points(j_w)
    1549              : 
    1550        15604 :       DO i_t = 1, bs_env%num_time_freq_points
    1551              : 
    1552        14104 :          time_i = bs_env%imag_time_points(i_t)
    1553        14104 :          weight_ij = bs_env%weights_cos_w_to_t(i_t, j_w)
    1554              : 
    1555              :          ! actual Fourier transform
    1556              :          CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=fm_W_MIC_time(i_t), &
    1557        15604 :                                   beta=weight_ij*COS(time_i*freq_j), matrix_b=fm_W_MIC_freq_j)
    1558              : 
    1559              :       END DO
    1560              : 
    1561         1500 :       CALL timestop(handle)
    1562              : 
    1563         1500 :    END SUBROUTINE Fourier_transform_w_to_t
    1564              : 
    1565              : ! **************************************************************************************************
    1566              : !> \brief ...
    1567              : !> \param bs_env ...
    1568              : !> \param qs_env ...
    1569              : !> \param fm_W_MIC_time ...
    1570              : ! **************************************************************************************************
    1571           32 :    SUBROUTINE multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_W_MIC_time)
    1572              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1573              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1574              :       TYPE(cp_fm_type), DIMENSION(:)                     :: fm_W_MIC_time
    1575              : 
    1576              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'multiply_fm_W_MIC_time_with_Minv_Gamma'
    1577              : 
    1578              :       INTEGER                                            :: handle, i_t, n_RI, ndep
    1579              :       TYPE(cp_fm_type)                                   :: fm_work
    1580           32 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_Minv_Gamma
    1581              : 
    1582           32 :       CALL timeset(routineN, handle)
    1583              : 
    1584           32 :       n_RI = bs_env%n_RI
    1585              : 
    1586           32 :       CALL cp_fm_create(fm_work, fm_W_MIC_time(1)%matrix_struct)
    1587              : 
    1588              :       ! compute Gamma-only RI-metric matrix M(k=0); no regularization
    1589              :       CALL RI_2c_integral_mat(qs_env, fm_Minv_Gamma, fm_W_MIC_time(1), n_RI, &
    1590           32 :                               bs_env%ri_metric, do_kpoints=.FALSE.)
    1591              : 
    1592           32 :       CALL cp_fm_power(fm_Minv_Gamma(1, 1), fm_work, -1.0_dp, 0.0_dp, ndep)
    1593              : 
    1594              :       ! M^-1(k=0)*W^MIC(iτ)*M^-1(k=0)
    1595          272 :       DO i_t = 1, SIZE(fm_W_MIC_time)
    1596              : 
    1597              :          CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_Minv_Gamma(1, 1), &
    1598          240 :                             fm_W_MIC_time(i_t), 0.0_dp, fm_work)
    1599              : 
    1600              :          CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_work, &
    1601          272 :                             fm_Minv_Gamma(1, 1), 0.0_dp, fm_W_MIC_time(i_t))
    1602              : 
    1603              :       END DO
    1604              : 
    1605           32 :       CALL cp_fm_release(fm_work)
    1606           32 :       CALL cp_fm_release(fm_Minv_Gamma)
    1607              : 
    1608           32 :       CALL timestop(handle)
    1609              : 
    1610           64 :    END SUBROUTINE multiply_fm_W_MIC_time_with_Minv_Gamma
    1611              : 
    1612              : ! **************************************************************************************************
    1613              : !> \brief ...
    1614              : !> \param bs_env ...
    1615              : !> \param qs_env ...
    1616              : !> \param fm_Sigma_x_Gamma ...
    1617              : ! **************************************************************************************************
    1618           22 :    SUBROUTINE get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
    1619              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1620              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1621              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_Sigma_x_Gamma
    1622              : 
    1623              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_Sigma_x'
    1624              : 
    1625              :       INTEGER                                            :: handle, ispin
    1626              : 
    1627           22 :       CALL timeset(routineN, handle)
    1628              : 
    1629           92 :       ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
    1630           48 :       DO ispin = 1, bs_env%n_spin
    1631           48 :          CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
    1632              :       END DO
    1633              : 
    1634           22 :       IF (bs_env%Sigma_x_exists) THEN
    1635           14 :          DO ispin = 1, bs_env%n_spin
    1636           14 :             CALL fm_read(fm_Sigma_x_Gamma(ispin), bs_env, bs_env%Sigma_x_name, ispin)
    1637              :          END DO
    1638              :       ELSE
    1639           16 :          CALL compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
    1640              :       END IF
    1641              : 
    1642           22 :       CALL timestop(handle)
    1643              : 
    1644           22 :    END SUBROUTINE get_Sigma_x
    1645              : 
    1646              : ! **************************************************************************************************
    1647              : !> \brief ...
    1648              : !> \param bs_env ...
    1649              : !> \param qs_env ...
    1650              : !> \param fm_Sigma_x_Gamma ...
    1651              : ! **************************************************************************************************
    1652           16 :    SUBROUTINE compute_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
    1653              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1654              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1655              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_Sigma_x_Gamma
    1656              : 
    1657              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_Sigma_x'
    1658              : 
    1659              :       INTEGER                                            :: handle, i_intval_idx, ispin, j_intval_idx
    1660              :       INTEGER, DIMENSION(2)                              :: i_atoms, j_atoms
    1661              :       REAL(KIND=dp)                                      :: t1
    1662           16 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: fm_Vtr_Gamma
    1663              :       TYPE(dbcsr_type)                                   :: mat_Sigma_x_Gamma
    1664          528 :       TYPE(dbt_type)                                     :: t_2c_D, t_2c_Sigma_x, t_2c_V, t_3c_x_V
    1665              : 
    1666           16 :       CALL timeset(routineN, handle)
    1667              : 
    1668           16 :       t1 = m_walltime()
    1669              : 
    1670           16 :       CALL dbt_create(bs_env%t_G, t_2c_D)
    1671           16 :       CALL dbt_create(bs_env%t_W, t_2c_V)
    1672           16 :       CALL dbt_create(bs_env%t_G, t_2c_Sigma_x)
    1673           16 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_V)
    1674           16 :       CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
    1675              : 
    1676              :       ! 1. Compute truncated Coulomb operator matrix V^tr(k=0) (cutoff rad: cellsize/2)
    1677              :       CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
    1678           16 :                               bs_env%trunc_coulomb, do_kpoints=.FALSE.)
    1679              : 
    1680              :       ! 2. Compute M^-1(k=0) and get M^-1(k=0)*V^tr(k=0)*M^-1(k=0)
    1681           16 :       CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_Vtr_Gamma(:, 1))
    1682              : 
    1683           34 :       DO ispin = 1, bs_env%n_spin
    1684              : 
    1685              :          ! 3. Compute density matrix D_µν
    1686           18 :          CALL G_occ_vir(bs_env, 0.0_dp, bs_env%fm_work_mo(2), ispin, occ=.TRUE., vir=.FALSE.)
    1687              : 
    1688              :          CALL fm_to_local_tensor(bs_env%fm_work_mo(2), bs_env%mat_ao_ao%matrix, &
    1689              :                                  bs_env%mat_ao_ao_tensor%matrix, t_2c_D, bs_env, &
    1690           18 :                                  bs_env%atoms_i_t_group)
    1691              : 
    1692              :          CALL fm_to_local_tensor(fm_Vtr_Gamma(1, 1), bs_env%mat_RI_RI%matrix, &
    1693              :                                  bs_env%mat_RI_RI_tensor%matrix, t_2c_V, bs_env, &
    1694           18 :                                  bs_env%atoms_j_t_group)
    1695              : 
    1696              :          ! every group has its own range of i_atoms and j_atoms; only deal with a
    1697              :          ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
    1698           36 :          DO i_intval_idx = 1, bs_env%n_intervals_i
    1699           54 :             DO j_intval_idx = 1, bs_env%n_intervals_j
    1700           54 :                i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
    1701           54 :                j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
    1702              : 
    1703              :                ! 4. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
    1704              :                ! 5. M_νσQ(iτ) = sum_P (νσ|P) (M^-1(k=0)*V^tr(k=0)*M^-1(k=0))_PQ(iτ)
    1705           18 :                CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_V, t_2c_V)
    1706              : 
    1707              :                ! 6. tensor operations with D and computation of Σ^x
    1708              :                !    Σ^x_λσ(k=0) = sum_νQ M_νσQ(iτ) sum_µ (λµ|Q) D_µν
    1709              :                CALL contract_to_Sigma(t_2c_D, t_3c_x_V, t_2c_Sigma_x, i_atoms, j_atoms, &
    1710           36 :                                       qs_env, bs_env, occ=.TRUE., vir=.FALSE., clear_W=.TRUE.)
    1711              : 
    1712              :             END DO ! j_atoms
    1713              :          END DO ! i_atoms
    1714              : 
    1715              :          CALL local_dbt_to_global_mat(t_2c_Sigma_x, bs_env%mat_ao_ao_tensor%matrix, &
    1716           18 :                                       mat_Sigma_x_Gamma, bs_env%para_env)
    1717              : 
    1718              :          CALL write_matrix(mat_Sigma_x_Gamma, ispin, bs_env%Sigma_x_name, &
    1719           18 :                            bs_env%fm_work_mo(1), qs_env)
    1720              : 
    1721           34 :          CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
    1722              : 
    1723              :       END DO ! ispin
    1724              : 
    1725           16 :       IF (bs_env%unit_nr > 0) THEN
    1726              :          WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
    1727            8 :             'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
    1728            8 :          WRITE (bs_env%unit_nr, '(A)') ' '
    1729              :       END IF
    1730              : 
    1731           16 :       CALL dbcsr_release(mat_Sigma_x_Gamma)
    1732           16 :       CALL dbt_destroy(t_2c_D)
    1733           16 :       CALL dbt_destroy(t_2c_V)
    1734           16 :       CALL dbt_destroy(t_2c_Sigma_x)
    1735           16 :       CALL dbt_destroy(t_3c_x_V)
    1736           16 :       CALL cp_fm_release(fm_Vtr_Gamma)
    1737              : 
    1738           16 :       CALL timestop(handle)
    1739              : 
    1740           32 :    END SUBROUTINE compute_Sigma_x
    1741              : 
    1742              : ! **************************************************************************************************
    1743              : !> \brief ...
    1744              : !> \param bs_env ...
    1745              : !> \param qs_env ...
    1746              : !> \param fm_W_MIC_time ...
    1747              : !> \param fm_Sigma_c_Gamma_time ...
    1748              : ! **************************************************************************************************
    1749           22 :    SUBROUTINE get_Sigma_c(bs_env, qs_env, fm_W_MIC_time, fm_Sigma_c_Gamma_time)
    1750              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1751              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1752              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
    1753              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :)  :: fm_Sigma_c_Gamma_time
    1754              : 
    1755              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_Sigma_c'
    1756              : 
    1757              :       INTEGER                                            :: handle, i_intval_idx, i_t, ispin, &
    1758              :                                                             j_intval_idx, read_write_index
    1759              :       INTEGER, DIMENSION(2)                              :: i_atoms, j_atoms
    1760              :       REAL(KIND=dp)                                      :: t1, tau
    1761           22 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
    1762          374 :       TYPE(dbt_type)                                     :: t_2c_Gocc, t_2c_Gvir, &
    1763          198 :                                                             t_2c_Sigma_neg_tau, &
    1764          550 :                                                             t_2c_Sigma_pos_tau, t_2c_W, t_3c_x_W
    1765              : 
    1766           22 :       CALL timeset(routineN, handle)
    1767              : 
    1768              :       CALL create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
    1769              :                                   t_2c_Sigma_pos_tau, t_3c_x_W, &
    1770           22 :                                   mat_Sigma_neg_tau, mat_Sigma_pos_tau)
    1771              : 
    1772          346 :       DO i_t = 1, bs_env%num_time_freq_points
    1773              : 
    1774          710 :          DO ispin = 1, bs_env%n_spin
    1775              : 
    1776          364 :             t1 = m_walltime()
    1777              : 
    1778          364 :             read_write_index = i_t + (ispin - 1)*bs_env%num_time_freq_points
    1779              : 
    1780              :             ! read self-energy from restart
    1781          364 :             IF (bs_env%Sigma_c_exists(i_t, ispin)) THEN
    1782          120 :                CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_p_name, read_write_index)
    1783              :                CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_pos_tau(i_t, ispin)%matrix, &
    1784          120 :                                      keep_sparsity=.FALSE.)
    1785          120 :                CALL fm_read(bs_env%fm_work_mo(1), bs_env, bs_env%Sigma_n_name, read_write_index)
    1786              :                CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_neg_tau(i_t, ispin)%matrix, &
    1787          120 :                                      keep_sparsity=.FALSE.)
    1788          120 :                IF (bs_env%unit_nr > 0) THEN
    1789           60 :                   WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F7.1,A)') 'Read Σ^c(iτ,k=0) ', &
    1790           60 :                      'from file for time point  ', i_t, ' /', bs_env%num_time_freq_points, &
    1791          120 :                      ',    Execution time', m_walltime() - t1, ' s'
    1792              :                END IF
    1793              : 
    1794              :                CYCLE
    1795              : 
    1796              :             END IF
    1797              : 
    1798          244 :             tau = bs_env%imag_time_points(i_t)
    1799              : 
    1800          244 :             CALL G_occ_vir(bs_env, tau, bs_env%fm_Gocc, ispin, occ=.TRUE., vir=.FALSE.)
    1801          244 :             CALL G_occ_vir(bs_env, tau, bs_env%fm_Gvir, ispin, occ=.FALSE., vir=.TRUE.)
    1802              : 
    1803              :             ! fm G^occ, G^vir and W to local tensor
    1804              :             CALL fm_to_local_tensor(bs_env%fm_Gocc, bs_env%mat_ao_ao%matrix, &
    1805              :                                     bs_env%mat_ao_ao_tensor%matrix, t_2c_Gocc, bs_env, &
    1806          244 :                                     bs_env%atoms_i_t_group)
    1807              :             CALL fm_to_local_tensor(bs_env%fm_Gvir, bs_env%mat_ao_ao%matrix, &
    1808              :                                     bs_env%mat_ao_ao_tensor%matrix, t_2c_Gvir, bs_env, &
    1809          244 :                                     bs_env%atoms_i_t_group)
    1810              :             CALL fm_to_local_tensor(fm_W_MIC_time(i_t), bs_env%mat_RI_RI%matrix, &
    1811              :                                     bs_env%mat_RI_RI_tensor%matrix, t_2c_W, bs_env, &
    1812          244 :                                     bs_env%atoms_j_t_group)
    1813              : 
    1814              :             ! every group has its own range of i_atoms and j_atoms; only deal with a
    1815              :             ! limited number of i_atom-j_atom pairs simultaneously in a group to save memory
    1816          488 :             DO i_intval_idx = 1, bs_env%n_intervals_i
    1817          732 :                DO j_intval_idx = 1, bs_env%n_intervals_j
    1818          732 :                   i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
    1819          732 :                   j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
    1820              : 
    1821          244 :                   IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
    1822              :                       bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) CYCLE
    1823              : 
    1824              :                   ! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
    1825              :                   ! 2. tensor operation M_νσQ(iτ) = sum_P (νσ|P) W^MIC_PQ(iτ)
    1826          226 :                   CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
    1827              : 
    1828              :                   ! 3. Σ_λσ(iτ,k=0) = sum_νQ M_νσQ(iτ) sum_µ (λµ|Q) G^occ_µν(i|τ|) for τ < 0
    1829              :                   !    (recall M_νσQ(iτ) = M_νσQ(-iτ) because W^MIC_PQ(iτ) = W^MIC_PQ(-iτ) )
    1830              :                   CALL contract_to_Sigma(t_2c_Gocc, t_3c_x_W, t_2c_Sigma_neg_tau, i_atoms, j_atoms, &
    1831              :                                          qs_env, bs_env, occ=.TRUE., vir=.FALSE., clear_W=.FALSE., &
    1832          226 :                                          can_skip=bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx))
    1833              : 
    1834              :                   !    Σ_λσ(iτ,k=0) = sum_νQ M_νσQ(iτ) sum_µ (λµ|Q) G^vir_µν(i|τ|) for τ > 0
    1835              :                   CALL contract_to_Sigma(t_2c_Gvir, t_3c_x_W, t_2c_Sigma_pos_tau, i_atoms, j_atoms, &
    1836              :                                          qs_env, bs_env, occ=.FALSE., vir=.TRUE., clear_W=.TRUE., &
    1837          488 :                                          can_skip=bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx))
    1838              : 
    1839              :                END DO ! j_atoms
    1840              :             END DO ! i_atoms
    1841              : 
    1842              :             ! 4. communicate data tensor t_2c_Sigma (which is local in the subgroup)
    1843              :             !    to the global dbcsr matrix mat_Sigma_pos/neg_tau (which stores Σ for all iτ)
    1844              :             CALL local_dbt_to_global_mat(t_2c_Sigma_neg_tau, bs_env%mat_ao_ao_tensor%matrix, &
    1845          244 :                                          mat_Sigma_neg_tau(i_t, ispin)%matrix, bs_env%para_env)
    1846              :             CALL local_dbt_to_global_mat(t_2c_Sigma_pos_tau, bs_env%mat_ao_ao_tensor%matrix, &
    1847          244 :                                          mat_Sigma_pos_tau(i_t, ispin)%matrix, bs_env%para_env)
    1848              : 
    1849              :             CALL write_matrix(mat_Sigma_pos_tau(i_t, ispin)%matrix, read_write_index, &
    1850          244 :                               bs_env%Sigma_p_name, bs_env%fm_work_mo(1), qs_env)
    1851              :             CALL write_matrix(mat_Sigma_neg_tau(i_t, ispin)%matrix, read_write_index, &
    1852          244 :                               bs_env%Sigma_n_name, bs_env%fm_work_mo(1), qs_env)
    1853              : 
    1854          568 :             IF (bs_env%unit_nr > 0) THEN
    1855              :                WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
    1856          122 :                   'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
    1857          244 :                   ',    Execution time', m_walltime() - t1, ' s'
    1858              :             END IF
    1859              : 
    1860              :          END DO ! ispin
    1861              : 
    1862              :       END DO ! i_t
    1863              : 
    1864           22 :       IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
    1865              : 
    1866              :       CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
    1867           22 :                                       mat_Sigma_pos_tau, mat_Sigma_neg_tau)
    1868              : 
    1869           22 :       CALL print_skipping(bs_env)
    1870              : 
    1871              :       CALL destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
    1872              :                                t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
    1873           22 :                                mat_Sigma_neg_tau, mat_Sigma_pos_tau)
    1874              : 
    1875           22 :       CALL delete_unnecessary_files(bs_env)
    1876              : 
    1877           22 :       CALL timestop(handle)
    1878              : 
    1879           44 :    END SUBROUTINE get_Sigma_c
    1880              : 
    1881              : ! **************************************************************************************************
    1882              : !> \brief ...
    1883              : !> \param bs_env ...
    1884              : !> \param t_2c_Gocc ...
    1885              : !> \param t_2c_Gvir ...
    1886              : !> \param t_2c_W ...
    1887              : !> \param t_2c_Sigma_neg_tau ...
    1888              : !> \param t_2c_Sigma_pos_tau ...
    1889              : !> \param t_3c_x_W ...
    1890              : !> \param mat_Sigma_neg_tau ...
    1891              : !> \param mat_Sigma_pos_tau ...
    1892              : ! **************************************************************************************************
    1893           22 :    SUBROUTINE create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
    1894              :                                      t_2c_Sigma_pos_tau, t_3c_x_W, &
    1895              :                                      mat_Sigma_neg_tau, mat_Sigma_pos_tau)
    1896              : 
    1897              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1898              :       TYPE(dbt_type)                                     :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
    1899              :                                                             t_2c_Sigma_neg_tau, &
    1900              :                                                             t_2c_Sigma_pos_tau, t_3c_x_W
    1901              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
    1902              : 
    1903              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_mat_for_Sigma_c'
    1904              : 
    1905              :       INTEGER                                            :: handle, i_t, ispin
    1906              : 
    1907           22 :       CALL timeset(routineN, handle)
    1908              : 
    1909           22 :       CALL dbt_create(bs_env%t_G, t_2c_Gocc)
    1910           22 :       CALL dbt_create(bs_env%t_G, t_2c_Gvir)
    1911           22 :       CALL dbt_create(bs_env%t_W, t_2c_W)
    1912           22 :       CALL dbt_create(bs_env%t_G, t_2c_Sigma_neg_tau)
    1913           22 :       CALL dbt_create(bs_env%t_G, t_2c_Sigma_pos_tau)
    1914           22 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_W)
    1915              : 
    1916           22 :       NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
    1917          478 :       ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
    1918          478 :       ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
    1919              : 
    1920          346 :       DO i_t = 1, bs_env%num_time_freq_points
    1921          710 :          DO ispin = 1, bs_env%n_spin
    1922          364 :             ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
    1923          364 :             ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
    1924          364 :             CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
    1925          688 :             CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
    1926              :          END DO
    1927              :       END DO
    1928              : 
    1929           22 :       CALL timestop(handle)
    1930              : 
    1931           22 :    END SUBROUTINE create_mat_for_Sigma_c
    1932              : 
    1933              : ! **************************************************************************************************
    1934              : !> \brief ...
    1935              : !> \param qs_env ...
    1936              : !> \param bs_env ...
    1937              : !> \param i_atoms ...
    1938              : !> \param j_atoms ...
    1939              : !> \param t_3c_x_W ...
    1940              : !> \param t_2c_W ...
    1941              : ! **************************************************************************************************
    1942          244 :    SUBROUTINE compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
    1943              : 
    1944              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1945              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1946              :       INTEGER, DIMENSION(2)                              :: i_atoms, j_atoms
    1947              :       TYPE(dbt_type)                                     :: t_3c_x_W, t_2c_W
    1948              : 
    1949              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_and_contract_W'
    1950              : 
    1951              :       INTEGER                                            :: handle, RI_intval_idx
    1952              :       INTEGER, DIMENSION(2)                              :: bounds_j, RI_atoms
    1953         4148 :       TYPE(dbt_type)                                     :: t_3c_for_W, t_3c_x_W_tmp
    1954              : 
    1955          244 :       CALL timeset(routineN, handle)
    1956              : 
    1957          244 :       CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_W_tmp)
    1958          244 :       CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_for_W)
    1959              : 
    1960              :       bounds_j(1:2) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
    1961          732 :                        bs_env%i_RI_end_from_atom(j_atoms(2))]
    1962              : 
    1963          488 :       DO RI_intval_idx = 1, bs_env%n_intervals_inner_loop_atoms
    1964          732 :          RI_atoms = bs_env%inner_loop_atom_intervals(1:2, RI_intval_idx)
    1965              : 
    1966              :          ! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
    1967              :          CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_W, &
    1968          244 :                                    atoms_AO_1=i_atoms, atoms_RI=RI_atoms)
    1969              : 
    1970              :          ! 2. tensor operation M_νσQ(iτ) = sum_P (νσ|P) W^MIC_PQ(iτ)
    1971              :          CALL dbt_contract(alpha=1.0_dp, &
    1972              :                            tensor_1=t_2c_W, &
    1973              :                            tensor_2=t_3c_for_W, &
    1974              :                            beta=1.0_dp, &
    1975              :                            tensor_3=t_3c_x_W_tmp, &
    1976              :                            contract_1=[2], notcontract_1=[1], map_1=[1], &
    1977              :                            contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
    1978              :                            bounds_2=bounds_j, &
    1979          488 :                            filter_eps=bs_env%eps_filter)
    1980              : 
    1981              :       END DO ! RI_atoms
    1982              : 
    1983              :       ! 3. reorder tensor
    1984          244 :       CALL dbt_copy(t_3c_x_W_tmp, t_3c_x_W, order=[1, 2, 3], move_data=.TRUE.)
    1985              : 
    1986          244 :       CALL dbt_destroy(t_3c_x_W_tmp)
    1987          244 :       CALL dbt_destroy(t_3c_for_W)
    1988              : 
    1989          244 :       CALL timestop(handle)
    1990              : 
    1991          244 :    END SUBROUTINE compute_3c_and_contract_W
    1992              : 
    1993              : ! **************************************************************************************************
    1994              : !> \brief ...
    1995              : !> \param t_2c_G ...
    1996              : !> \param t_3c_x_W ...
    1997              : !> \param t_2c_Sigma ...
    1998              : !> \param i_atoms ...
    1999              : !> \param j_atoms ...
    2000              : !> \param qs_env ...
    2001              : !> \param bs_env ...
    2002              : !> \param occ ...
    2003              : !> \param vir ...
    2004              : !> \param clear_W ...
    2005              : !> \param can_skip ...
    2006              : ! **************************************************************************************************
    2007          470 :    SUBROUTINE contract_to_Sigma(t_2c_G, t_3c_x_W, t_2c_Sigma, i_atoms, j_atoms, qs_env, bs_env, &
    2008              :                                 occ, vir, clear_W, can_skip)
    2009              :       TYPE(dbt_type)                                     :: t_2c_G, t_3c_x_W, t_2c_Sigma
    2010              :       INTEGER, DIMENSION(2)                              :: i_atoms, j_atoms
    2011              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2012              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2013              :       LOGICAL                                            :: occ, vir, clear_W
    2014              :       LOGICAL, OPTIONAL                                  :: can_skip
    2015              : 
    2016              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'contract_to_Sigma'
    2017              : 
    2018              :       INTEGER :: handle, inner_loop_atoms_interval_index
    2019              :       INTEGER(KIND=int_8)                                :: flop
    2020              :       INTEGER, DIMENSION(2)                              :: bounds_i, IL_atoms
    2021              :       REAL(KIND=dp)                                      :: sign_Sigma
    2022        11750 :       TYPE(dbt_type)                                     :: t_3c_for_G, t_3c_x_G, t_3c_x_G_2
    2023              : 
    2024          470 :       CALL timeset(routineN, handle)
    2025              : 
    2026          470 :       CPASSERT(occ .EQV. (.NOT. vir))
    2027          470 :       IF (occ) sign_Sigma = -1.0_dp
    2028          470 :       IF (vir) sign_Sigma = 1.0_dp
    2029              : 
    2030          470 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_G)
    2031          470 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G)
    2032          470 :       CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G_2)
    2033              : 
    2034              :       bounds_i(1:2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
    2035         1410 :                        bs_env%i_ao_end_from_atom(i_atoms(2))]
    2036              : 
    2037          940 :       DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
    2038         1410 :          IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
    2039              : 
    2040              :          CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_G, &
    2041          470 :                                    atoms_RI=j_atoms, atoms_AO_2=IL_atoms)
    2042              : 
    2043              :          CALL dbt_contract(alpha=1.0_dp, &
    2044              :                            tensor_1=t_2c_G, &
    2045              :                            tensor_2=t_3c_for_G, &
    2046              :                            beta=1.0_dp, &
    2047              :                            tensor_3=t_3c_x_G, &
    2048              :                            contract_1=[2], notcontract_1=[1], map_1=[3], &
    2049              :                            contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
    2050              :                            bounds_2=bounds_i, &
    2051          940 :                            filter_eps=bs_env%eps_filter)
    2052              : 
    2053              :       END DO ! IL_atoms
    2054              : 
    2055          470 :       CALL dbt_copy(t_3c_x_G, t_3c_x_G_2, order=[1, 3, 2], move_data=.TRUE.)
    2056              : 
    2057              :       CALL dbt_contract(alpha=sign_Sigma, &
    2058              :                         tensor_1=t_3c_x_W, &
    2059              :                         tensor_2=t_3c_x_G_2, &
    2060              :                         beta=1.0_dp, &
    2061              :                         tensor_3=t_2c_Sigma, &
    2062              :                         contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
    2063              :                         contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
    2064          470 :                         filter_eps=bs_env%eps_filter, move_data=clear_W, flop=flop)
    2065              : 
    2066          470 :       IF (PRESENT(can_skip)) THEN
    2067          452 :          IF (flop == 0_int_8) can_skip = .TRUE.
    2068              :       END IF
    2069              : 
    2070          470 :       CALL dbt_destroy(t_3c_for_G)
    2071          470 :       CALL dbt_destroy(t_3c_x_G)
    2072          470 :       CALL dbt_destroy(t_3c_x_G_2)
    2073              : 
    2074          470 :       CALL timestop(handle)
    2075              : 
    2076          470 :    END SUBROUTINE contract_to_Sigma
    2077              : 
    2078              : ! **************************************************************************************************
    2079              : !> \brief ...
    2080              : !> \param fm_Sigma_c_Gamma_time ...
    2081              : !> \param bs_env ...
    2082              : !> \param mat_Sigma_pos_tau ...
    2083              : !> \param mat_Sigma_neg_tau ...
    2084              : ! **************************************************************************************************
    2085           22 :    SUBROUTINE fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
    2086              :                                          mat_Sigma_pos_tau, mat_Sigma_neg_tau)
    2087              : 
    2088              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :)  :: fm_Sigma_c_Gamma_time
    2089              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2090              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_Sigma_pos_tau, mat_Sigma_neg_tau
    2091              : 
    2092              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_fm_Sigma_c_Gamma_time'
    2093              : 
    2094              :       INTEGER                                            :: handle, i_t, ispin, pos_neg
    2095              : 
    2096           22 :       CALL timeset(routineN, handle)
    2097              : 
    2098          894 :       ALLOCATE (fm_Sigma_c_Gamma_time(bs_env%num_time_freq_points, 2, bs_env%n_spin))
    2099          346 :       DO i_t = 1, bs_env%num_time_freq_points
    2100          710 :          DO ispin = 1, bs_env%n_spin
    2101         1092 :             DO pos_neg = 1, 2
    2102              :                CALL cp_fm_create(fm_Sigma_c_Gamma_time(i_t, pos_neg, ispin), &
    2103         1092 :                                  bs_env%fm_s_Gamma%matrix_struct)
    2104              :             END DO
    2105              :             CALL copy_dbcsr_to_fm(mat_Sigma_pos_tau(i_t, ispin)%matrix, &
    2106          364 :                                   fm_Sigma_c_Gamma_time(i_t, 1, ispin))
    2107              :             CALL copy_dbcsr_to_fm(mat_Sigma_neg_tau(i_t, ispin)%matrix, &
    2108          688 :                                   fm_Sigma_c_Gamma_time(i_t, 2, ispin))
    2109              :          END DO
    2110              :       END DO
    2111              : 
    2112           22 :       CALL timestop(handle)
    2113              : 
    2114           22 :    END SUBROUTINE fill_fm_Sigma_c_Gamma_time
    2115              : 
    2116              : ! **************************************************************************************************
    2117              : !> \brief ...
    2118              : !> \param bs_env ...
    2119              : ! **************************************************************************************************
    2120           22 :    SUBROUTINE print_skipping(bs_env)
    2121              : 
    2122              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2123              : 
    2124              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'print_skipping'
    2125              : 
    2126              :       INTEGER                                            :: handle, i_intval_idx, j_intval_idx, &
    2127              :                                                             n_skip
    2128              : 
    2129           22 :       CALL timeset(routineN, handle)
    2130              : 
    2131           22 :       n_skip = 0
    2132              : 
    2133           44 :       DO i_intval_idx = 1, bs_env%n_intervals_i
    2134           66 :          DO j_intval_idx = 1, bs_env%n_intervals_j
    2135           22 :             IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
    2136           22 :                 bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
    2137            2 :                n_skip = n_skip + 1
    2138              :             END IF
    2139              :          END DO
    2140              :       END DO
    2141              : 
    2142           22 :       IF (bs_env%unit_nr > 0) THEN
    2143              :          WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
    2144           11 :             'Sparsity of Σ^c(iτ,k=0): Percentage of skipped atom pairs:', &
    2145           22 :             REAL(100*n_skip, KIND=dp)/REAL(i_intval_idx*j_intval_idx, KIND=dp), ' %'
    2146              :       END IF
    2147              : 
    2148           22 :       CALL timestop(handle)
    2149              : 
    2150           22 :    END SUBROUTINE print_skipping
    2151              : 
    2152              : ! **************************************************************************************************
    2153              : !> \brief ...
    2154              : !> \param t_2c_Gocc ...
    2155              : !> \param t_2c_Gvir ...
    2156              : !> \param t_2c_W ...
    2157              : !> \param t_2c_Sigma_neg_tau ...
    2158              : !> \param t_2c_Sigma_pos_tau ...
    2159              : !> \param t_3c_x_W ...
    2160              : !> \param fm_W_MIC_time ...
    2161              : !> \param mat_Sigma_neg_tau ...
    2162              : !> \param mat_Sigma_pos_tau ...
    2163              : ! **************************************************************************************************
    2164           22 :    SUBROUTINE destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg_tau, &
    2165              :                                   t_2c_Sigma_pos_tau, t_3c_x_W, fm_W_MIC_time, &
    2166              :                                   mat_Sigma_neg_tau, mat_Sigma_pos_tau)
    2167              : 
    2168              :       TYPE(dbt_type)                                     :: t_2c_Gocc, t_2c_Gvir, t_2c_W, &
    2169              :                                                             t_2c_Sigma_neg_tau, &
    2170              :                                                             t_2c_Sigma_pos_tau, t_3c_x_W
    2171              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_W_MIC_time
    2172              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
    2173              : 
    2174              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'destroy_mat_Sigma_c'
    2175              : 
    2176              :       INTEGER                                            :: handle
    2177              : 
    2178           22 :       CALL timeset(routineN, handle)
    2179              : 
    2180           22 :       CALL dbt_destroy(t_2c_Gocc)
    2181           22 :       CALL dbt_destroy(t_2c_Gvir)
    2182           22 :       CALL dbt_destroy(t_2c_W)
    2183           22 :       CALL dbt_destroy(t_2c_Sigma_neg_tau)
    2184           22 :       CALL dbt_destroy(t_2c_Sigma_pos_tau)
    2185           22 :       CALL dbt_destroy(t_3c_x_W)
    2186           22 :       CALL cp_fm_release(fm_W_MIC_time)
    2187           22 :       CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
    2188           22 :       CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
    2189              : 
    2190           22 :       CALL timestop(handle)
    2191              : 
    2192           22 :    END SUBROUTINE destroy_mat_Sigma_c
    2193              : 
    2194              : ! **************************************************************************************************
    2195              : !> \brief ...
    2196              : !> \param bs_env ...
    2197              : ! **************************************************************************************************
    2198           22 :    SUBROUTINE delete_unnecessary_files(bs_env)
    2199              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2200              : 
    2201              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'delete_unnecessary_files'
    2202              : 
    2203              :       CHARACTER(LEN=default_string_length)               :: f_chi, f_W_t, prefix
    2204              :       INTEGER                                            :: handle, i_t
    2205              : 
    2206           22 :       CALL timeset(routineN, handle)
    2207              : 
    2208           22 :       prefix = bs_env%prefix
    2209              : 
    2210          346 :       DO i_t = 1, bs_env%num_time_freq_points
    2211              : 
    2212          324 :          IF (i_t < 10) THEN
    2213          186 :             WRITE (f_chi, '(3A,I1,A)') TRIM(prefix), bs_env%chi_name, "_00", i_t, ".matrix"
    2214          186 :             WRITE (f_W_t, '(3A,I1,A)') TRIM(prefix), bs_env%W_time_name, "_00", i_t, ".matrix"
    2215          138 :          ELSE IF (i_t < 100) THEN
    2216          138 :             WRITE (f_chi, '(3A,I2,A)') TRIM(prefix), bs_env%chi_name, "_0", i_t, ".matrix"
    2217          138 :             WRITE (f_W_t, '(3A,I2,A)') TRIM(prefix), bs_env%W_time_name, "_0", i_t, ".matrix"
    2218              :          ELSE
    2219            0 :             CPABORT('Please implement more than 99 time/frequency points.')
    2220              :          END IF
    2221              : 
    2222          324 :          CALL safe_delete(f_chi, bs_env)
    2223          346 :          CALL safe_delete(f_W_t, bs_env)
    2224              : 
    2225              :       END DO
    2226              : 
    2227           22 :       CALL timestop(handle)
    2228              : 
    2229           22 :    END SUBROUTINE delete_unnecessary_files
    2230              : 
    2231              : ! **************************************************************************************************
    2232              : !> \brief ...
    2233              : !> \param filename ...
    2234              : !> \param bs_env ...
    2235              : ! **************************************************************************************************
    2236          648 :    SUBROUTINE safe_delete(filename, bs_env)
    2237              :       CHARACTER(LEN=*)                                   :: filename
    2238              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2239              : 
    2240              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'safe_delete'
    2241              : 
    2242              :       INTEGER                                            :: handle
    2243              :       LOGICAL                                            :: file_exists
    2244              : 
    2245          648 :       CALL timeset(routineN, handle)
    2246              : 
    2247          648 :       IF (bs_env%para_env%mepos == 0) THEN
    2248              : 
    2249          324 :          INQUIRE (file=TRIM(filename), exist=file_exists)
    2250          324 :          IF (file_exists) CALL mp_file_delete(TRIM(filename))
    2251              : 
    2252              :       END IF
    2253              : 
    2254          648 :       CALL timestop(handle)
    2255              : 
    2256          648 :    END SUBROUTINE safe_delete
    2257              : 
    2258              : ! **************************************************************************************************
    2259              : !> \brief ...
    2260              : !> \param bs_env ...
    2261              : !> \param qs_env ...
    2262              : !> \param fm_Sigma_x_Gamma ...
    2263              : !> \param fm_Sigma_c_Gamma_time ...
    2264              : ! **************************************************************************************************
    2265           22 :    SUBROUTINE compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
    2266              : 
    2267              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2268              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2269              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_Sigma_x_Gamma
    2270              :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :)  :: fm_Sigma_c_Gamma_time
    2271              : 
    2272              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
    2273              : 
    2274              :       INTEGER                                            :: handle, ikp, ispin, j_t
    2275              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Sigma_x_ikp_n, V_xc_ikp_n
    2276              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: Sigma_c_ikp_n_freq, Sigma_c_ikp_n_time
    2277              :       TYPE(cp_cfm_type)                                  :: cfm_ks_ikp, cfm_mos_ikp, cfm_s_ikp, &
    2278              :                                                             cfm_Sigma_x_ikp, cfm_work_ikp
    2279              : 
    2280           22 :       CALL timeset(routineN, handle)
    2281              : 
    2282           22 :       CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
    2283           22 :       CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_s_Gamma%matrix_struct)
    2284              :       ! JW TODO: fully distribute these arrays at given time; also eigenvalues in bs_env
    2285          110 :       ALLOCATE (V_xc_ikp_n(bs_env%n_ao), Sigma_x_ikp_n(bs_env%n_ao))
    2286          110 :       ALLOCATE (Sigma_c_ikp_n_time(bs_env%n_ao, bs_env%num_time_freq_points, 2))
    2287           88 :       ALLOCATE (Sigma_c_ikp_n_freq(bs_env%n_ao, bs_env%num_time_freq_points, 2))
    2288              : 
    2289           48 :       DO ispin = 1, bs_env%n_spin
    2290              : 
    2291           86 :          DO ikp = 1, bs_env%nkp_bs_and_DOS
    2292              : 
    2293              :             ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
    2294              :             CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
    2295           38 :                                        ikp, qs_env, bs_env%kpoints_DOS, "ORB")
    2296              : 
    2297              :             ! 2. get S_µν(k_i) from S_µν(k=0)
    2298              :             CALL cfm_ikp_from_fm_Gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
    2299           38 :                                        ikp, qs_env, bs_env%kpoints_DOS, "ORB")
    2300              : 
    2301              :             ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
    2302              :             CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp, cfm_mos_ikp, &
    2303           38 :                               bs_env%eigenval_scf(:, ikp, ispin), cfm_work_ikp)
    2304              : 
    2305              :             ! 4. V^xc_µν(k=0) -> V^xc_µν(k_i) -> V^xc_nn(k_i)
    2306              :             CALL to_ikp_and_mo(V_xc_ikp_n, bs_env%fm_V_xc_Gamma(ispin), &
    2307           38 :                                ikp, qs_env, bs_env, cfm_mos_ikp)
    2308              : 
    2309              :             ! 5. Σ^x_µν(k=0) -> Σ^x_µν(k_i) -> Σ^x_nn(k_i)
    2310              :             CALL to_ikp_and_mo(Sigma_x_ikp_n, fm_Sigma_x_Gamma(ispin), &
    2311           38 :                                ikp, qs_env, bs_env, cfm_mos_ikp)
    2312              : 
    2313              :             ! 6. Σ^c_µν(k=0,+/-i|τ_j|) -> Σ^c_µν(k_i,+/-i|τ_j|) -> Σ^c_nn(k_i,+/-i|τ_j|)
    2314          506 :             DO j_t = 1, bs_env%num_time_freq_points
    2315              :                CALL to_ikp_and_mo(Sigma_c_ikp_n_time(:, j_t, 1), &
    2316              :                                   fm_Sigma_c_Gamma_time(j_t, 1, ispin), &
    2317          468 :                                   ikp, qs_env, bs_env, cfm_mos_ikp)
    2318              :                CALL to_ikp_and_mo(Sigma_c_ikp_n_time(:, j_t, 2), &
    2319              :                                   fm_Sigma_c_Gamma_time(j_t, 2, ispin), &
    2320          506 :                                   ikp, qs_env, bs_env, cfm_mos_ikp)
    2321              :             END DO
    2322              : 
    2323              :             ! 7. Σ^c_nn(k_i,iτ) -> Σ^c_nn(k_i,iω)
    2324           38 :             CALL time_to_freq(bs_env, Sigma_c_ikp_n_time, Sigma_c_ikp_n_freq, ispin)
    2325              : 
    2326              :             ! 8. Analytic continuation Σ^c_nn(k_i,iω) -> Σ^c_nn(k_i,ϵ) and
    2327              :             !    ϵ_nk_i^GW = ϵ_nk_i^DFT + Σ^c_nn(k_i,ϵ) + Σ^x_nn(k_i) - v^xc_nn(k_i)
    2328              :             CALL analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_xc_ikp_n, &
    2329           64 :                                         bs_env%eigenval_scf(:, ikp, ispin), ikp, ispin)
    2330              : 
    2331              :          END DO ! ikp_DOS
    2332              : 
    2333              :       END DO ! ispin
    2334              : 
    2335           22 :       CALL get_all_VBM_CBM_bandgaps(bs_env)
    2336              : 
    2337           22 :       CALL cp_fm_release(fm_Sigma_x_Gamma)
    2338           22 :       CALL cp_fm_release(fm_Sigma_c_Gamma_time)
    2339           22 :       CALL cp_cfm_release(cfm_ks_ikp)
    2340           22 :       CALL cp_cfm_release(cfm_s_ikp)
    2341           22 :       CALL cp_cfm_release(cfm_mos_ikp)
    2342           22 :       CALL cp_cfm_release(cfm_work_ikp)
    2343           22 :       CALL cp_cfm_release(cfm_Sigma_x_ikp)
    2344              : 
    2345           22 :       CALL timestop(handle)
    2346              : 
    2347           44 :    END SUBROUTINE compute_QP_energies
    2348              : 
    2349              : ! **************************************************************************************************
    2350              : !> \brief ...
    2351              : !> \param array_ikp_n ...
    2352              : !> \param fm_Gamma ...
    2353              : !> \param ikp ...
    2354              : !> \param qs_env ...
    2355              : !> \param bs_env ...
    2356              : !> \param cfm_mos_ikp ...
    2357              : ! **************************************************************************************************
    2358         1012 :    SUBROUTINE to_ikp_and_mo(array_ikp_n, fm_Gamma, ikp, qs_env, bs_env, cfm_mos_ikp)
    2359              : 
    2360              :       REAL(KIND=dp), DIMENSION(:)                        :: array_ikp_n
    2361              :       TYPE(cp_fm_type)                                   :: fm_Gamma
    2362              :       INTEGER                                            :: ikp
    2363              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2364              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2365              :       TYPE(cp_cfm_type)                                  :: cfm_mos_ikp
    2366              : 
    2367              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'to_ikp_and_mo'
    2368              : 
    2369              :       INTEGER                                            :: handle
    2370              :       TYPE(cp_fm_type)                                   :: fm_ikp_mo_re
    2371              : 
    2372         1012 :       CALL timeset(routineN, handle)
    2373              : 
    2374         1012 :       CALL cp_fm_create(fm_ikp_mo_re, fm_Gamma%matrix_struct)
    2375              : 
    2376         1012 :       CALL fm_Gamma_ao_to_cfm_ikp_mo(fm_Gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
    2377              : 
    2378         1012 :       CALL cp_fm_get_diag(fm_ikp_mo_re, array_ikp_n)
    2379              : 
    2380         1012 :       CALL cp_fm_release(fm_ikp_mo_re)
    2381              : 
    2382         1012 :       CALL timestop(handle)
    2383              : 
    2384         1012 :    END SUBROUTINE to_ikp_and_mo
    2385              : 
    2386              : ! **************************************************************************************************
    2387              : !> \brief ...
    2388              : !> \param fm_Gamma ...
    2389              : !> \param fm_ikp_mo_re ...
    2390              : !> \param ikp ...
    2391              : !> \param qs_env ...
    2392              : !> \param bs_env ...
    2393              : !> \param cfm_mos_ikp ...
    2394              : ! **************************************************************************************************
    2395         4048 :    SUBROUTINE fm_Gamma_ao_to_cfm_ikp_mo(fm_Gamma, fm_ikp_mo_re, ikp, qs_env, bs_env, cfm_mos_ikp)
    2396              :       TYPE(cp_fm_type)                                   :: fm_Gamma, fm_ikp_mo_re
    2397              :       INTEGER                                            :: ikp
    2398              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2399              :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2400              :       TYPE(cp_cfm_type)                                  :: cfm_mos_ikp
    2401              : 
    2402              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_Gamma_ao_to_cfm_ikp_mo'
    2403              : 
    2404              :       INTEGER                                            :: handle, nmo
    2405              :       TYPE(cp_cfm_type)                                  :: cfm_ikp_ao, cfm_ikp_mo, cfm_tmp
    2406              : 
    2407         1012 :       CALL timeset(routineN, handle)
    2408              : 
    2409         1012 :       CALL cp_cfm_create(cfm_ikp_ao, fm_Gamma%matrix_struct)
    2410         1012 :       CALL cp_cfm_create(cfm_ikp_mo, fm_Gamma%matrix_struct)
    2411         1012 :       CALL cp_cfm_create(cfm_tmp, fm_Gamma%matrix_struct)
    2412              : 
    2413              :       ! get cfm_µν(k_i) from fm_µν(k=0)
    2414         1012 :       CALL cfm_ikp_from_fm_Gamma(cfm_ikp_ao, fm_Gamma, ikp, qs_env, bs_env%kpoints_DOS, "ORB")
    2415              : 
    2416         1012 :       nmo = bs_env%n_ao
    2417         1012 :       CALL parallel_gemm('N', 'N', nmo, nmo, nmo, z_one, cfm_ikp_ao, cfm_mos_ikp, z_zero, cfm_tmp)
    2418         1012 :       CALL parallel_gemm('C', 'N', nmo, nmo, nmo, z_one, cfm_mos_ikp, cfm_tmp, z_zero, cfm_ikp_mo)
    2419              : 
    2420         1012 :       CALL cp_cfm_to_fm(cfm_ikp_mo, fm_ikp_mo_re)
    2421              : 
    2422         1012 :       CALL cp_cfm_release(cfm_ikp_mo)
    2423         1012 :       CALL cp_cfm_release(cfm_ikp_ao)
    2424         1012 :       CALL cp_cfm_release(cfm_tmp)
    2425              : 
    2426         1012 :       CALL timestop(handle)
    2427              : 
    2428         1012 :    END SUBROUTINE fm_Gamma_ao_to_cfm_ikp_mo
    2429              : 
    2430              : END MODULE gw_large_cell_gamma
        

Generated by: LCOV version 2.0-1