gw_tensor_large_cell_gamma_static_cohsex.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Static COHSEX band structure in the tensor-based large-cell Gamma GW code:
!>        COHSEX DIAGONAL (diagonal energies) or COHSEX QP_UPDATE (canonical energies and MOs)
!> \author Ritaj Tyagi
!> \date Jan 2025
! **************************************************************************************************
MODULE gw_tensor_large_cell_gamma_static_cohsex
   USE atomic_kind_types,               ONLY: atomic_kind_type
   USE bibliography,                    ONLY: hedin1965,&
                                              hybertsen1986,&
                                              tyagi2024,&
                                              tyagi2026,&
                                              cite_reference
   USE cell_types,                      ONLY: cell_type
   USE cp_cfm_diag,                     ONLY: cp_cfm_geeig
   USE cp_cfm_types,                    ONLY: cp_cfm_create,&
                                              cp_cfm_release,&
                                              cp_cfm_to_fm,&
                                              cp_cfm_type
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
                                              dbcsr_p_type,&
                                              dbcsr_set
   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
                                              dbcsr_deallocate_matrix_set
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add,&
                                              cp_fm_transpose
   USE cp_fm_diag,                      ONLY: cp_fm_syevx
   USE cp_fm_types,                     ONLY: cp_fm_create,&
                                              cp_fm_get_diag,&
                                              cp_fm_get_info,&
                                              cp_fm_release,&
                                              cp_fm_set_all,&
                                              cp_fm_to_fm,&
                                              cp_fm_type,&
                                              cp_fm_write_formatted
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr
   USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
   USE dbt_api,                         ONLY: dbt_type
   USE gw_compute_v_and_w,              ONLY: apply_extrapol_factor,&
                                              compute_minvvsqrt_vsqrt,&
                                              compute_v_k_by_lattice_sum,&
                                              compute_fm_w_mic_freq_j_from_chi,&
                                              fm_read,&
                                              fm_write
   USE gw_tensor_large_cell_gamma,      ONLY: compute_sigma_c_occ_vir,&
                                              create_mat_for_sigma_c,&
                                              delete_unnecessary_files,&
                                              destroy_mat_sigma_c,&
                                              get_sigma_x,&
                                              get_mat_chi_gamma_tau,&
                                              print_skipping
   USE gw_utils,                        ONLY: de_init_bs_env
   USE gw_utils_fm,                     ONLY: cfm_contract_aba,&
                                              fm_contract_aba
   USE input_constants,                 ONLY: cohsex_qp_update
   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kinds,                           ONLY: default_path_length,&
                                              default_string_length,&
                                              dp
   USE machine,                         ONLY: m_walltime
   USE parallel_gemm_api,               ONLY: parallel_gemm
   USE particle_list_types,             ONLY: particle_list_type
   USE particle_types,                  ONLY: particle_type
   USE physcon,                         ONLY: evolt
   USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
   USE post_scf_bandstructure_utils,    ONLY: cfm_ikp_from_fm_gamma,&
                                              get_all_vbm_cbm_bandgaps
   USE pw_env_types,                    ONLY: pw_env_get,&
                                              pw_env_type
   USE pw_pool_types,                   ONLY: pw_pool_type
   USE pw_types,                        ONLY: pw_c1d_gs_type,&
                                              pw_r3d_rs_type
   USE qs_collocate_density,            ONLY: calculate_wavefunction
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_kind_types,                   ONLY: qs_kind_type
   USE qs_mo_io,                        ONLY: write_mo_set_to_output_unit
   USE qs_mo_types,                     ONLY: allocate_mo_set,&
                                              deallocate_mo_set,&
                                              init_mo_set,&
                                              mo_set_type
   USE qs_subsys_types,                 ONLY: qs_subsys_get,&
                                              qs_subsys_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_tensor_large_cell_Gamma_static_cohsex'
 
   PUBLIC :: gw_calc_tensor_large_cell_gamma_static_cohsex
 
CONTAINS
 
! **************************************************************************************************
!> \brief Perform static COHSEX band structure calculation
!> \param qs_env ...
!> \param bs_env ...
!> \par History
!>    * 12.2024 created [Ritaj Tyagi]
! **************************************************************************************************
   SUBROUTINE gw_calc_tensor_large_cell_gamma_static_cohsex(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'gw_calc_tensor_large_cell_Gamma_static_cohsex'
 
      INTEGER                                            :: handle
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_sigma_c_gamma, fm_sigma_x_gamma, &
                                                            fm_w_mic_freq_zero
 
      CALL timeset(routinen, handle)
 
      ! static COHSEX; with QP_UPDATE the quasiparticle update; scaled self-energy if requested
      CALL cite_reference(hedin1965)
      CALL cite_reference(hybertsen1986)
      IF (bs_env%cohsex_flavor == cohsex_qp_update) CALL cite_reference(tyagi2024)
      IF (bs_env%cohsex_scale_occ /= 1.0_dp .OR. bs_env%cohsex_scale_virt /= 1.0_dp) THEN
         CALL cite_reference(tyagi2026)
      END IF
 
      IF (bs_env%Sigma_c_freq_zero_exists) THEN
 
         ! restart: Σ^c_λσ(iω=0,k=0) from file; χ and W are not needed
         CALL read_sigma_c_freq_zero(bs_env, fm_sigma_c_gamma)
 
      ELSE
 
         ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
         ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
         ! χ_PQ(iτ,k=0) = sum_λν [sum_µ (µν|P) G^occ_µλ(i|τ|)] [sum_σ (σλ|Q) G^vir_σν(i|τ|)]
         CALL get_mat_chi_gamma_tau(bs_env, qs_env, bs_env%mat_chi_Gamma_tau)
 
         ! χ_PQ(iτ,k=0) -> χ_PQ(iω=0,k=0) -> χ_PQ(iω=0,k) -> ε_PQ(iω=0,k) -> W_PQ(iω=0,k)
         ! -> W^MIC_PQ(iω=0) -> M^-1*W^MIC(iω=0)*M^-1
         CALL get_w_mic_freq_zero(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_w_mic_freq_zero)
 
         ! Σ^c_λσ(iω=0,k=0) = 1/2 sum_n sgn(ϵ_n - ϵ_F) sum_νQ [sum_P (νσ|P) W^MIC_PQ(iω=0)]
         !                    [sum_µ (λµ|Q) C_µn C_νn]; written to the Σ^c(iω=0) restart file
         CALL get_sigma_c(bs_env, qs_env, fm_w_mic_freq_zero, fm_sigma_c_gamma)
 
      END IF
 
      ! 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>
      ! Σ^x_λσ(k=0) = sum_νQ [sum_P (νσ|P) V^trunc_PQ] [sum_µ (λµ|Q) D_µν)]
      CALL get_sigma_x(bs_env, qs_env, fm_sigma_x_gamma)
 
      ! Σ^c_λσ(k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
      CALL compute_qp_energies(bs_env, qs_env, fm_sigma_x_gamma, fm_sigma_c_gamma)
 
      CALL de_init_bs_env(qs_env, bs_env)
 
      CALL timestop(handle)
 
   END SUBROUTINE gw_calc_tensor_large_cell_gamma_static_cohsex
 
! **************************************************************************************************
!> \brief Static screened Coulomb interaction at iω = 0:
!>          χ_PQ(iω=0,k=0) = sum_i w^0_i χ_PQ(iτ_i,k=0)   (w^0: time weights at zero frequency)
!>          χ_PQ(iω=0,k=0) -> χ_PQ(iω=0,k) -> ε_PQ(iω=0,k) -> W_PQ(iω=0,k) -> W^MIC_PQ(iω=0)
!>        as for one frequency point of G0W0 (compute_fm_W_MIC_freq_j_from_chi); then the
!>        k-point extrapolation factor and M^-1(k=0) W^MIC(iω=0) M^-1(k=0)
!> \param bs_env ...
!> \param qs_env ...
!> \param mat_chi_Gamma_tau ...
!> \param fm_W_MIC_freq_zero ...
! **************************************************************************************************
   SUBROUTINE get_w_mic_freq_zero(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_freq_zero)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_chi_gamma_tau
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_w_mic_freq_zero
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'get_W_MIC_freq_zero'
 
      INTEGER                                            :: handle, i_t, ikp, ikp_batch, ikp_in_batch
      REAL(kind=dp)                                      :: t1
      TYPE(cp_cfm_type)                                  :: cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp
      TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:)       :: cfm_v_kp
 
      CALL timeset(routinen, handle)
 
      ! χ_PQ(iω=0,k=0) = sum_i w^0_i χ_PQ(iτ_i,k=0) -> bs_env%fm_chi_Gamma_freq
      CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
      DO i_t = 1, bs_env%num_time_freq_points
         CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_gamma_tau(i_t)%matrix, 1.0_dp, &
                        bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t))
      END DO
      CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, bs_env%fm_chi_Gamma_freq)
      CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
 
      ALLOCATE (fm_w_mic_freq_zero(1))
      CALL cp_fm_create(fm_w_mic_freq_zero(1), bs_env%fm_RI_RI%matrix_struct)
      CALL cp_fm_set_all(fm_w_mic_freq_zero(1), 0.0_dp)
 
      DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
 
         t1 = m_walltime()
 
         ! V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
         CALL compute_v_k_by_lattice_sum(bs_env, qs_env, cfm_v_kp, ikp_batch)
 
         DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
 
            ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
 
            IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) cycle
 
            CALL compute_minvvsqrt_vsqrt(bs_env, qs_env, cfm_v_kp, &
                                         cfm_v_sqrt_ikp, cfm_m_inv_v_sqrt_ikp, ikp)
 
            CALL bs_env%para_env%sync()
            CALL cp_cfm_release(cfm_v_kp(ikp))
 
            ! W^MIC contribution of k-point ikp at iω = 0 (j_w = 1: first frequency point for the
            ! approximate k-point extrapolation)
            CALL compute_fm_w_mic_freq_j_from_chi(bs_env, qs_env, bs_env%fm_W_MIC_freq, 1, ikp, &
                                                  cfm_m_inv_v_sqrt_ikp, cfm_v_sqrt_ikp)
 
            CALL cp_fm_scale_and_add(1.0_dp, fm_w_mic_freq_zero(1), 1.0_dp, bs_env%fm_W_MIC_freq)
 
         END DO ! ikp_in_batch
 
         DEALLOCATE (cfm_v_kp)
 
         IF (bs_env%unit_nr > 0) THEN
            WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F10.1,A)') &
               'Computed W(iω=0,k) for k-point batch', &
               ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
               ', Execution time', m_walltime() - t1, ' s'
         END IF
 
      END DO ! ikp_batch
 
      IF (bs_env%approx_kp_extrapol) THEN
         CALL apply_extrapol_factor(bs_env, fm_w_mic_freq_zero)
      END IF
 
      ! M^-1(k=0) W^MIC(iω=0) M^-1(k=0)
      CALL fm_contract_aba(bs_env%fm_Minv_Gamma, fm_w_mic_freq_zero(1))
 
      CALL cp_cfm_release(cfm_m_inv_v_sqrt_ikp)
      CALL cp_cfm_release(cfm_v_sqrt_ikp)
 
      IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
 
      CALL timestop(handle)
 
   END SUBROUTINE get_w_mic_freq_zero
 
! **************************************************************************************************
!> \brief Static COHSEX correlation self-energy at the Γ point in the AO basis,
!>          Σ^c_λσ = 1/2 sum_n sgn(ϵ_n - ϵ_F) sum_νQ [sum_P (νσ|P) W^MIC_PQ(iω=0)] [sum_µ (λµ|Q) C_µn C_νn]
!>                 = 1/2 (Σ^neg_λσ(iτ=0) + Σ^pos_λσ(iτ=0))   with W = W^MIC(iω=0),
!>        i.e. the occupied and virtual self-energy parts of G0W0 at τ = 0 (compute_Sigma_c_occ_vir),
!>        symmetrized
!> \param bs_env ...
!> \param qs_env ...
!> \param fm_W_MIC_freq_zero ...
!> \param fm_Sigma_c_Gamma ...
! **************************************************************************************************
   SUBROUTINE get_sigma_c(bs_env, qs_env, fm_W_MIC_freq_zero, fm_Sigma_c_Gamma)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_w_mic_freq_zero, fm_sigma_c_gamma
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'get_Sigma_c'
 
      INTEGER                                            :: handle, ispin
      REAL(kind=dp)                                      :: t1
      TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: mat_sigma_neg, mat_sigma_pos
      TYPE(dbt_type)                                     :: t_2c_gocc, t_2c_gvir, t_2c_sigma_neg, &
                                                            t_2c_sigma_pos, t_2c_w, t_3c_x_w
 
      CALL timeset(routinen, handle)
 
      CALL create_mat_for_sigma_c(bs_env, t_2c_gocc, t_2c_gvir, t_2c_w, t_2c_sigma_neg, &
                                  t_2c_sigma_pos, t_3c_x_w, mat_sigma_neg, mat_sigma_pos)
 
      ALLOCATE (fm_sigma_c_gamma(bs_env%n_spin))
      DO ispin = 1, bs_env%n_spin
 
         t1 = m_walltime()
 
         ! Σ^neg(iτ=0) = -G^occ(0) W^MIC(iω=0),  Σ^pos(iτ=0) = G^vir(0) W^MIC(iω=0)
         CALL compute_sigma_c_occ_vir(bs_env, qs_env, 0.0_dp, ispin, fm_w_mic_freq_zero(1), &
                                      t_2c_gocc, t_2c_gvir, t_2c_w, t_2c_sigma_neg, &
                                      t_2c_sigma_pos, t_3c_x_w, mat_sigma_neg(1, ispin)%matrix, &
                                      mat_sigma_pos(1, ispin)%matrix, count_skipped=.false.)
 
         ! Σ^c = 1/2 (Σ^neg + Σ^pos), symmetrized
         CALL cp_fm_create(fm_sigma_c_gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
         CALL copy_dbcsr_to_fm(mat_sigma_neg(1, ispin)%matrix, fm_sigma_c_gamma(ispin))
         CALL copy_dbcsr_to_fm(mat_sigma_pos(1, ispin)%matrix, bs_env%fm_work_mo(1))
         CALL cp_fm_scale_and_add(0.5_dp, fm_sigma_c_gamma(ispin), 0.5_dp, bs_env%fm_work_mo(1))
         CALL cp_fm_transpose(fm_sigma_c_gamma(ispin), bs_env%fm_work_mo(1))
         CALL cp_fm_scale_and_add(0.5_dp, fm_sigma_c_gamma(ispin), 0.5_dp, bs_env%fm_work_mo(1))
 
         ! restart file of Σ^c(iω=0,k=0), one per spin
         CALL fm_write(fm_sigma_c_gamma(ispin), ispin, bs_env%Sigma_c_freq_zero_name, qs_env)
 
         IF (bs_env%unit_nr > 0) THEN
            WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
               'Computed Σ^c(k=0),', ' Execution time', m_walltime() - t1, ' s'
         END IF
 
      END DO
 
      IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
 
      CALL print_skipping(bs_env)
 
      CALL destroy_mat_sigma_c(t_2c_gocc, t_2c_gvir, t_2c_w, t_2c_sigma_neg, t_2c_sigma_pos, &
                               t_3c_x_w, fm_w_mic_freq_zero, mat_sigma_neg, mat_sigma_pos)
 
      CALL delete_unnecessary_files(bs_env)
 
      CALL timestop(handle)
 
   END SUBROUTINE get_sigma_c
 
! **************************************************************************************************
!> \brief Read the static COHSEX restart file Σ^c_λσ(iω=0,k=0), one per spin
!> \param bs_env ...
!> \param fm_Sigma_c_Gamma Σ^c_λσ(iω=0,k=0) per spin, created here
! **************************************************************************************************
   SUBROUTINE read_sigma_c_freq_zero(bs_env, fm_Sigma_c_Gamma)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_sigma_c_gamma
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'read_Sigma_c_freq_zero'
 
      INTEGER                                            :: handle, ispin
      REAL(kind=dp)                                      :: t1
 
      CALL timeset(routinen, handle)
 
      ALLOCATE (fm_sigma_c_gamma(bs_env%n_spin))
      DO ispin = 1, bs_env%n_spin
         t1 = m_walltime()
         CALL cp_fm_create(fm_sigma_c_gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
         CALL fm_read(fm_sigma_c_gamma(ispin), bs_env, bs_env%Sigma_c_freq_zero_name, ispin)
         IF (bs_env%unit_nr > 0) THEN
            WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
               'Read Σ^c(k=0) from file,', ' Execution time', m_walltime() - t1, ' s'
         END IF
      END DO
 
      IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
 
      CALL timestop(handle)
 
   END SUBROUTINE read_sigma_c_freq_zero
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
!> \param qs_env ...
!> \param fm_Sigma_x_Gamma ...
!> \param fm_Sigma_c_Gamma ...
! **************************************************************************************************
   SUBROUTINE compute_qp_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma)
 
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_sigma_x_gamma, fm_sigma_c_gamma
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_QP_energies'
 
      INTEGER                                            :: handle, ikp, ispin
      TYPE(cp_cfm_type)                                  :: cfm_ks_ikp, cfm_mos_ikp, cfm_s_ikp, &
                                                            cfm_work_ikp
      TYPE(cp_fm_type)                                   :: fm_sigma_c_ikp_mo, fm_sigma_x_ikp_mo, &
                                                            fm_v_xc_ikp_mo
 
      CALL timeset(routinen, handle)
 
      CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
      CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_s_Gamma%matrix_struct)
 
      ! MO-basis operator matrices at k_i (n_ao x n_ao), reused for every k-point
      CALL cp_fm_create(fm_v_xc_ikp_mo, bs_env%fm_s_Gamma%matrix_struct)
      CALL cp_fm_create(fm_sigma_x_ikp_mo, bs_env%fm_s_Gamma%matrix_struct)
      CALL cp_fm_create(fm_sigma_c_ikp_mo, bs_env%fm_s_Gamma%matrix_struct)
 
      DO ispin = 1, bs_env%n_spin
 
         DO ikp = 1, bs_env%nkp_bs_and_DOS
 
            ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
            CALL cfm_ikp_from_fm_gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
                                       ikp, qs_env, bs_env%kpoints_DOS, "ORB")
 
            ! 2. get S_µν(k_i) from S_µν(k=0)
            CALL cfm_ikp_from_fm_gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
                                       ikp, qs_env, bs_env%kpoints_DOS, "ORB")
 
            ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
            CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp, cfm_mos_ikp, &
                              bs_env%eigenval_scf(:, ikp, ispin), cfm_work_ikp)
 
            ! Reuse the KS buffer for the full MO operators. Retain their real parts
            ! at the existing real COHSEX Hamiltonian boundary.
            ! 4. V^xc_µν(k=0) -> V^xc_µν(k_i) -> V^xc_mn(k_i)
            CALL cfm_ikp_from_fm_gamma(cfm_ks_ikp, bs_env%fm_V_xc_Gamma(ispin), &
                                       ikp, qs_env, bs_env%kpoints_DOS, "ORB")
            CALL cfm_contract_aba(cfm_mos_ikp, cfm_ks_ikp, workspace=cfm_work_ikp)
            CALL cp_cfm_to_fm(cfm_ks_ikp, fm_v_xc_ikp_mo)
 
            ! 5. Σ^x_µν(k=0) -> Σ^x_µν(k_i) -> Σ^x_mn(k_i)
            CALL cfm_ikp_from_fm_gamma(cfm_ks_ikp, fm_sigma_x_gamma(ispin), &
                                       ikp, qs_env, bs_env%kpoints_DOS, "ORB")
            CALL cfm_contract_aba(cfm_mos_ikp, cfm_ks_ikp, workspace=cfm_work_ikp)
            CALL cp_cfm_to_fm(cfm_ks_ikp, fm_sigma_x_ikp_mo)
 
            ! 6. Σ^c_µν(k=0) -> Σ^c_µν(k_i) -> Σ^c_mn(k_i)
            CALL cfm_ikp_from_fm_gamma(cfm_ks_ikp, fm_sigma_c_gamma(ispin), &
                                       ikp, qs_env, bs_env%kpoints_DOS, "ORB")
            CALL cfm_contract_aba(cfm_mos_ikp, cfm_ks_ikp, workspace=cfm_work_ikp)
            CALL cp_cfm_to_fm(cfm_ks_ikp, fm_sigma_c_ikp_mo)
 
            ! 7. build static COHSEX QP Hamiltonian, diagonalize, and print
            CALL diagonalize_and_print(bs_env, qs_env, fm_sigma_x_ikp_mo, fm_v_xc_ikp_mo, &
                                       fm_sigma_c_ikp_mo, ikp, ispin)
 
         END DO ! ikp_DOS
 
      END DO ! ispin
 
      CALL get_all_vbm_cbm_bandgaps(bs_env)
 
      CALL cp_fm_release(fm_sigma_x_gamma)
      CALL cp_fm_release(fm_sigma_c_gamma)
      CALL cp_fm_release(fm_v_xc_ikp_mo)
      CALL cp_fm_release(fm_sigma_x_ikp_mo)
      CALL cp_fm_release(fm_sigma_c_ikp_mo)
      CALL cp_cfm_release(cfm_ks_ikp)
      CALL cp_cfm_release(cfm_s_ikp)
      CALL cp_cfm_release(cfm_mos_ikp)
      CALL cp_cfm_release(cfm_work_ikp)
 
      CALL timestop(handle)
 
   END SUBROUTINE compute_qp_energies
 
! **************************************************************************************************
!> \brief Build the static COHSEX quasiparticle Hamiltonian H^QP at k-point ikp and print the
!>        diagonal energies ϵ^DG = H^QP_nn (COHSEX DIAGONAL); with COHSEX QP_UPDATE also
!>        diagonalize H^QP and print the canonical energies ϵ^CN and the canonical MOs
!> \param bs_env ...
!> \param qs_env ...
!> \param fm_Sigma_x_ikp_mo ...
!> \param fm_V_xc_ikp_mo ...
!> \param fm_Sigma_c_ikp_mo ...
!> \param ikp ...
!> \param ispin ...
! **************************************************************************************************
   SUBROUTINE diagonalize_and_print(bs_env, qs_env, fm_Sigma_x_ikp_mo, fm_V_xc_ikp_mo, &
                                    fm_Sigma_c_ikp_mo, ikp, ispin)
 
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_fm_type)                                   :: fm_sigma_x_ikp_mo, fm_v_xc_ikp_mo, &
                                                            fm_sigma_c_ikp_mo
      INTEGER                                            :: ikp, ispin
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'diagonalize_and_print'
 
      CHARACTER(len=3)                                   :: occ_vir
      CHARACTER(len=default_string_length)               :: fname
      INTEGER                                            :: handle, i_col, i_row, ikp_for_print, &
                                                            iunit, j_col, j_row, m, n_mo, &
                                                            ncol_local, nrow_local
      INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
      LOGICAL                                            :: is_bandstruc_kpoint, print_dos_kpoints, &
                                                            print_ikp, qp_update
      REAL(kind=dp)                                      :: val
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: diag_sigma_c, diag_sigma_x, diag_sp, &
                                                            diag_v_xc, sp_eigval
      TYPE(cp_fm_type)                                   :: fm_sp_mat, sp_eigvec, sp_eigvec_ao
 
      CALL timeset(routinen, handle)
 
      n_mo = bs_env%n_ao
      ! COHSEX QP_UPDATE: diagonalize H^QP (canonical energies and MOs); DIAGONAL: ϵ^DG only
      qp_update = (bs_env%cohsex_flavor == cohsex_qp_update)
 
      CALL cp_fm_create(fm_sp_mat, bs_env%fm_s_Gamma%matrix_struct)
 
      ALLOCATE (sp_eigval(n_mo))
      sp_eigval(:) = 0.0_dp
 
      ! Build the static COHSEX QP Hamiltonian directly in the distributed fm by looping
      ! over local elements only (global indices from row_indices/col_indices):
      !   occ-occ block : Σ^x - v^xc + α Σ^c   (α = COHSEX_SCALE_OCC)
      !   vir-vir block : Σ^x - v^xc + α Σ^c   (α = COHSEX_SCALE_VIRT)
      !   occ-vir / vir-occ : 0, i.e. occupied and virtual states are not mixed
      !   diagonal additionally gets + ϵ^DFT_nk
      CALL cp_fm_get_info(fm_sp_mat, nrow_local=nrow_local, ncol_local=ncol_local, &
                          row_indices=row_indices, col_indices=col_indices)
 
!$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE) &
!$OMP             SHARED(ncol_local, nrow_local, col_indices, row_indices, bs_env, ispin, ikp, &
!$OMP                    fm_Sigma_x_ikp_mo, fm_V_xc_ikp_mo, fm_Sigma_c_ikp_mo, fm_SP_mat) &
!$OMP             PRIVATE(i_col, i_row, j_col, j_row, val)
      DO i_col = 1, ncol_local
         DO i_row = 1, nrow_local
            j_col = col_indices(i_col)
            j_row = row_indices(i_row)
            IF (j_row <= bs_env%n_occ(ispin) .AND. j_col <= bs_env%n_occ(ispin)) THEN
               val = fm_sigma_x_ikp_mo%local_data(i_row, i_col) &
                     - fm_v_xc_ikp_mo%local_data(i_row, i_col) &
                     + bs_env%cohsex_scale_occ*fm_sigma_c_ikp_mo%local_data(i_row, i_col)
            ELSE IF (j_row > bs_env%n_occ(ispin) .AND. j_col > bs_env%n_occ(ispin)) THEN
               val = fm_sigma_x_ikp_mo%local_data(i_row, i_col) &
                     - fm_v_xc_ikp_mo%local_data(i_row, i_col) &
                     + bs_env%cohsex_scale_virt*fm_sigma_c_ikp_mo%local_data(i_row, i_col)
            ELSE
               val = 0.0_dp
            END IF
            IF (j_row == j_col) val = val + bs_env%eigenval_scf(j_row, ikp, ispin)
            fm_sp_mat%local_data(i_row, i_col) = val
         END DO
      END DO
!$OMP END PARALLEL DO
 
      ! diagonal elements (for the per-MO QP printout) before diagonalization
      ALLOCATE (diag_sp(n_mo), diag_sigma_c(n_mo), diag_sigma_x(n_mo), diag_v_xc(n_mo))
      CALL cp_fm_get_diag(fm_sp_mat, diag_sp)
      CALL cp_fm_get_diag(fm_sigma_c_ikp_mo, diag_sigma_c)
      CALL cp_fm_get_diag(fm_sigma_x_ikp_mo, diag_sigma_x)
      CALL cp_fm_get_diag(fm_v_xc_ikp_mo, diag_v_xc)
 
      IF (qp_update) THEN
         ! canonical static COHSEX: H^QP U = U ϵ^CN (the occ-virt blocks are zero)
         CALL cp_fm_create(sp_eigvec, bs_env%fm_s_Gamma%matrix_struct)
         CALL cp_fm_create(sp_eigvec_ao, bs_env%fm_s_Gamma%matrix_struct)
         CALL cp_fm_syevx(matrix=fm_sp_mat, eigenvectors=sp_eigvec, eigenvalues=sp_eigval)
         IF (ikp == bs_env%nkp_bs_and_DOS) THEN
            CALL fm_write_formatted(sp_eigvec, 0, "Singlepole_Eigenvector_Gamma", qs_env)
         END IF
 
         ! transform QP eigenvectors from the KS-MO basis to the AO basis
         CALL parallel_gemm(transa="N", transb="N", m=n_mo, n=n_mo, k=n_mo, alpha=1.0_dp, &
                            matrix_a=bs_env%fm_mo_coeff_Gamma(ispin), matrix_b=sp_eigvec, &
                            beta=0.0_dp, matrix_c=sp_eigvec_ao)
      ELSE
         ! diagonal static COHSEX: ϵ^DG = H^QP_nn, no diagonalization
         sp_eigval(:) = diag_sp(:)
      END IF
 
      bs_env%eigenval_GW(:, ikp, ispin) = sp_eigval(:)
      IF (ALLOCATED(bs_env%eigenval_G0W0)) THEN
         bs_env%eigenval_G0W0(:, ikp, ispin) = sp_eigval(:)
      END IF
      DO m = 1, n_mo
         bs_env%eigenval_HF(m, ikp, ispin) = bs_env%eigenval_scf(m, ikp, ispin) + &
                                             diag_sigma_x(m) - diag_v_xc(m)
      END DO
 
      ! only print eigenvalues of DOS k-points in case no bandstructure path has been given;
      ! in kpoints_DOS the last nkp_only_bs entries are the bandstructure-path k-points
      print_dos_kpoints = (bs_env%nkp_only_bs <= 0)
      is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
      print_ikp = print_dos_kpoints .OR. is_bandstruc_kpoint
 
      !! Printing Quasiparticle energies to a file
      IF (bs_env%para_env%is_source() .AND. print_ikp) THEN
 
         IF (print_dos_kpoints) THEN
            ikp_for_print = ikp
         ELSE
            ikp_for_print = ikp - bs_env%nkp_only_DOS
         END IF
 
         fname = "Quasiparticle_energies_singlepole"
         IF (ikp_for_print == 1 .AND. ispin == 1) THEN
            CALL open_file(trim(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
         ELSE
            CALL open_file(trim(fname), unit_number=iunit, file_status="OLD", file_action="WRITE", &
                           file_position="APPEND")
         END IF
 
         WRITE (iunit, "(A)") " "
         WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_print, "coordinate: ", &
            bs_env%kpoints_DOS%xkp(:, ikp)
         WRITE (iunit, "(A)") " "
         IF (qp_update) THEN
            WRITE (iunit, "(A5,A24,2A17,A16,A18,A18)") "n", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", &
               "Σ^x_nk (eV)", "v_nk^xc (eV)", "ϵ_nk^DG (eV)", "ϵ_nk^CN (eV)"
         ELSE
            WRITE (iunit, "(A5,A24,2A17,A16,A18)") "n", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", &
               "Σ^x_nk (eV)", "v_nk^xc (eV)", "ϵ_nk^DG (eV)"
         END IF
         WRITE (iunit, "(A)") " "
 
         DO m = 1, n_mo
            IF (m <= bs_env%n_occ(ispin)) occ_vir = 'occ'
            IF (m > bs_env%n_occ(ispin)) occ_vir = 'vir'
            IF (qp_update) THEN
               WRITE (iunit, "(I5,3A,4F16.3,F17.3,F17.3)") m, ' (', occ_vir, ') ', &
                  bs_env%eigenval_scf(m, ikp, ispin)*evolt, &
                  diag_sigma_c(m)*evolt, &
                  diag_sigma_x(m)*evolt, &
                  diag_v_xc(m)*evolt, &
                  diag_sp(m)*evolt, &
                  sp_eigval(m)*evolt
            ELSE
               WRITE (iunit, "(I5,3A,4F16.3,F17.3)") m, ' (', occ_vir, ') ', &
                  bs_env%eigenval_scf(m, ikp, ispin)*evolt, &
                  diag_sigma_c(m)*evolt, &
                  diag_sigma_x(m)*evolt, &
                  diag_v_xc(m)*evolt, &
                  diag_sp(m)*evolt
            END IF
         END DO
 
         WRITE (iunit, "(A)") " "
 
         CALL close_file(iunit)
      END IF
 
      IF (qp_update) THEN
         ! canonical MO energies/coefficients and MO cubes (DFT%PRINT%MO, DFT%PRINT%MO_CUBES),
         ! Gamma point only
         IF (ikp == 1) THEN
            CALL print_static_cohsex_mos(bs_env, qs_env, sp_eigvec_ao, sp_eigval, ispin)
         END IF
         CALL cp_fm_release(sp_eigvec)
         CALL cp_fm_release(sp_eigvec_ao)
      END IF
 
      CALL cp_fm_release(fm_sp_mat)
      DEALLOCATE (sp_eigval, diag_sp, diag_sigma_c, diag_sigma_x, diag_v_xc)
 
      CALL timestop(handle)
 
   END SUBROUTINE diagonalize_and_print
 
! **************************************************************************************************
!> \brief Print static COHSEX QP MO energies/occupations/coefficients and MO cubes at the Gamma point,
!>        reusing the standard DFT%PRINT%MO / DFT%PRINT%MO_CUBES print keys.
!> \param bs_env ...
!> \param qs_env ...
!> \param fm_qp_mo_ao ...
!> \param qp_eigval ...
!> \param ispin ...
! **************************************************************************************************
   SUBROUTINE print_static_cohsex_mos(bs_env, qs_env, fm_qp_mo_ao, qp_eigval, ispin)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_fm_type)                                   :: fm_qp_mo_ao
      REAL(kind=dp), DIMENSION(:)                        :: qp_eigval
      INTEGER                                            :: ispin
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'print_static_cohsex_mos'
 
      CHARACTER(LEN=5)                                   :: spin_label
      CHARACTER(LEN=default_path_length)                 :: fname, title
      INTEGER                                            :: first_mo, handle, ivector, last_mo, &
                                                            n_mo, nelectron, nhomo, nlumo, unit_nr
      INTEGER, DIMENSION(:), POINTER                     :: stride
      LOGICAL                                            :: print_coeff, print_cubes, &
                                                            print_energies, print_mo, print_occ
      REAL(kind=dp)                                      :: maxocc
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: cell
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(mo_set_type)                                  :: qp_mo_set
      TYPE(particle_list_type), POINTER                  :: particles
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(pw_c1d_gs_type)                               :: wf_g
      TYPE(pw_env_type), POINTER                         :: pw_env
      TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      TYPE(pw_r3d_rs_type)                               :: wf_r
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(section_vals_type), POINTER                   :: dft_section, input
 
      CALL timeset(routinen, handle)
 
      n_mo = bs_env%n_ao
 
      CALL get_qs_env(qs_env, input=input)
      dft_section => section_vals_get_subs_vals(input, "DFT")
      logger => cp_get_default_logger()
 
      ! ---------------------------------------------------------------------------------------
      ! MO energies / occupation numbers / coefficients
      ! ---------------------------------------------------------------------------------------
      CALL section_vals_val_get(dft_section, "PRINT%MO%ENERGIES", l_val=print_energies)
      CALL section_vals_val_get(dft_section, "PRINT%MO%COEFFICIENTS", l_val=print_coeff)
      CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occ)
      print_mo = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO") /= 0)
 
      IF (print_mo .AND. (print_energies .OR. print_occ .OR. print_coeff)) THEN
 
         ! occupation numbers of the QP states: the lowest n_occ are occupied
         IF (bs_env%n_spin == 1) THEN
            maxocc = 2.0_dp
         ELSE
            maxocc = 1.0_dp
         END IF
         nelectron = nint(maxocc*real(bs_env%n_occ(ispin), kind=dp))
 
         CALL allocate_mo_set(qp_mo_set, nao=n_mo, nmo=n_mo, nelectron=nelectron, &
                              n_el_f=real(nelectron, kind=dp), maxocc=maxocc, &
                              flexible_electron_count=0.0_dp)
         CALL init_mo_set(qp_mo_set, fm_struct=fm_qp_mo_ao%matrix_struct, name="static_cohsex")
         qp_mo_set%eigenvalues(:) = qp_eigval(:)
         CALL cp_fm_to_fm(fm_qp_mo_ao, qp_mo_set%mo_coeff)
 
         CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
 
         IF (bs_env%n_spin == 1) THEN
            CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
                                             4, 0, final_mos=.true.)
         ELSE
            IF (ispin == 1) THEN
               spin_label = "ALPHA"
            ELSE
               spin_label = "BETA"
            END IF
            CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
                                             4, 0, final_mos=.true., spin=trim(spin_label))
         END IF
 
         CALL deallocate_mo_set(qp_mo_set)
      END IF
 
      ! ---------------------------------------------------------------------------------------
      ! MO cubes
      ! ---------------------------------------------------------------------------------------
      print_cubes = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") /= 0)
 
      IF (print_cubes) THEN
         CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%NHOMO", i_val=nhomo)
         CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%NLUMO", i_val=nlumo)
         CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%STRIDE", i_vals=stride)
 
         IF (nhomo == -1) nhomo = bs_env%n_occ(ispin)
         IF (nlumo == -1) nlumo = n_mo - bs_env%n_occ(ispin)
         first_mo = max(1, bs_env%n_occ(ispin) - nhomo + 1)
         last_mo = min(n_mo, bs_env%n_occ(ispin) + nlumo)
 
         CALL get_qs_env(qs_env, &
                         atomic_kind_set=atomic_kind_set, &
                         qs_kind_set=qs_kind_set, &
                         cell=cell, &
                         dft_control=dft_control, &
                         particle_set=particle_set, &
                         pw_env=pw_env, &
                         subsys=subsys)
         CALL qs_subsys_get(subsys, particles=particles)
 
         CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
         CALL auxbas_pw_pool%create_pw(wf_r)
         CALL auxbas_pw_pool%create_pw(wf_g)
 
         DO ivector = first_mo, last_mo
            CALL calculate_wavefunction(fm_qp_mo_ao, ivector, wf_r, wf_g, atomic_kind_set, &
                                        qs_kind_set, cell, dft_control, particle_set, pw_env)
 
            WRITE (fname, '(A,I5.5,A,I1)') "static_cohsex_WFN_", ivector, "_", ispin
            unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO_CUBES", &
                                           extension=".cube", middle_name=trim(fname), &
                                           file_position="REWIND", log_filename=.false.)
            WRITE (title, *) "static COHSEX WAVEFUNCTION ", ivector, " spin ", ispin, &
               " i.e. HOMO ", ivector - bs_env%n_occ(ispin)
            CALL cp_pw_to_cube(wf_r, unit_nr, title, particles, stride=stride)
            CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%MO_CUBES")
         END DO
 
         CALL auxbas_pw_pool%give_back_pw(wf_r)
         CALL auxbas_pw_pool%give_back_pw(wf_g)
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE print_static_cohsex_mos
 
! **************************************************************************************************
!> \brief ...
!> \param fm ...
!> \param matrix_index ...
!> \param matrix_name ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE fm_write_formatted(fm, matrix_index, matrix_name, qs_env)
      TYPE(cp_fm_type)                                   :: fm
      INTEGER                                            :: matrix_index
      CHARACTER(LEN=*)                                   :: matrix_name
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
         routinen = 'fm_write_formatted'
 
      CHARACTER(LEN=default_string_length)               :: filename
      INTEGER                                            :: handle, unit_nr
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_vals_type), POINTER                   :: input
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, input=input)
 
      logger => cp_get_default_logger()
 
      IF (btest(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
 
         IF (matrix_index < 10) THEN
            WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
         ELSE IF (matrix_index < 100) THEN
            WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
         ELSE
            cpabort('Please implement more than 99 time/frequency points.')
         END IF
 
         unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
                                        file_form="FORMATTED", middle_name=trim(filename), &
                                        file_position="REWIND", file_action="WRITE")
 
         CALL cp_fm_write_formatted(fm, unit_nr)
         IF (unit_nr > 0) THEN
            CALL close_file(unit_nr)
         END IF
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE fm_write_formatted
 
END MODULE gw_tensor_large_cell_gamma_static_cohsex