gw_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief
!> \author Jan Wilhelm
!> \date 07.2023
! **************************************************************************************************
MODULE gw_utils
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind_set
   USE basis_set_types,                 ONLY: gto_basis_set_type
   USE bibliography,                    ONLY: graml2024,&
                                              cite_reference
   USE cell_types,                      ONLY: cell_type,&
                                              pbc
   USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
                                              cp_blacs_env_release,&
                                              cp_blacs_env_type
   USE cp_control_types,                ONLY: dft_control_type
   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
                                              dbcsr_allocate_matrix_set,&
                                              dbcsr_deallocate_matrix_set
   USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
   USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
                                              cp_fm_struct_release,&
                                              cp_fm_struct_type
   USE cp_fm_types,                     ONLY: cp_fm_create,&
                                              cp_fm_set_all
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_print_key_generate_filename
   USE dbcsr_api,                       ONLY: dbcsr_create,&
                                              dbcsr_p_type
   USE dbt_api,                         ONLY: &
        dbt_clear, dbt_create, dbt_destroy, dbt_filter, dbt_iterator_blocks_left, &
        dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
        dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
   USE input_constants,                 ONLY: do_potential_truncated,&
                                              xc_none
   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get,&
                                              section_vals_val_set
   USE kinds,                           ONLY: default_string_length,&
                                              dp,&
                                              int_8
   USE kpoint_methods,                  ONLY: kpoint_init_cell_index
   USE kpoint_types,                    ONLY: kpoint_create,&
                                              kpoint_type
   USE libint_wrapper,                  ONLY: cp_libint_static_cleanup,&
                                              cp_libint_static_init
   USE machine,                         ONLY: m_memory,&
                                              m_walltime
   USE mathlib,                         ONLY: gcd
   USE message_passing,                 ONLY: mp_cart_type,&
                                              mp_para_env_type
   USE minimax_exp,                     ONLY: get_exp_minimax_coeff
   USE minimax_exp_gw,                  ONLY: get_exp_minimax_coeff_gw
   USE minimax_rpa,                     ONLY: get_rpa_minimax_coeff,&
                                              get_rpa_minimax_coeff_larger_grid
   USE mp2_gpw,                         ONLY: create_mat_munu
   USE mp2_grids,                       ONLY: get_l_sq_wghts_cos_tf_t_to_w,&
                                              get_l_sq_wghts_cos_tf_w_to_t,&
                                              get_l_sq_wghts_sin_tf_t_to_w
   USE mp2_ri_2c,                       ONLY: setup_trunc_coulomb_pot_for_exchange_self_energy
   USE particle_methods,                ONLY: get_particle_set
   USE particle_types,                  ONLY: particle_type
   USE physcon,                         ONLY: angstrom,&
                                              evolt
   USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
   USE qs_energy_types,                 ONLY: qs_energy_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_env_part_release,&
                                              qs_environment_type
   USE qs_integral_utils,               ONLY: basis_set_list_setup
   USE qs_interactions,                 ONLY: init_interaction_radii_orb_basis
   USE qs_kind_types,                   ONLY: get_qs_kind,&
                                              qs_kind_type
   USE qs_ks_methods,                   ONLY: qs_ks_build_kohn_sham_matrix
   USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
   USE qs_tensors,                      ONLY: build_3c_integrals,&
                                              build_3c_neighbor_lists,&
                                              get_tensor_occupancy,&
                                              neighbor_list_3c_destroy
   USE qs_tensors_types,                ONLY: create_2c_tensor,&
                                              create_3c_tensor,&
                                              distribution_3d_create,&
                                              distribution_3d_type,&
                                              neighbor_list_3c_type
#include "base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   PUBLIC :: create_and_init_bs_env_for_gw, de_init_bs_env, get_i_j_atoms, &
             kpoint_init_cell_index_simple, compute_xkp
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils'
 
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param bs_sec ...
! **************************************************************************************************
   SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(section_vals_type), POINTER                   :: bs_sec
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'create_and_init_bs_env_for_gw'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
 
      CALL cite_reference(graml2024)
 
      CALL read_gw_input_parameters(bs_env, bs_sec)
 
      CALL print_header_and_input_parameters(bs_env)
 
      CALL setup_ao_and_ri_basis_set(qs_env, bs_env)
 
      CALL get_ri_basis_and_basis_function_indices(qs_env, bs_env)
 
      CALL setup_kpoints_chi_eps_w(qs_env, bs_env, bs_env%kpoints_chi_eps_W)
 
      CALL set_heuristic_parameters(bs_env, qs_env)
 
      CALL set_parallelization_parameters(qs_env, bs_env)
 
      CALL allocate_and_fill_matrices_and_arrays(qs_env, bs_env)
 
      CALL cp_libint_static_init()
 
      CALL create_tensors(qs_env, bs_env)
 
      CALL set_sparsity_parallelization_parameters(bs_env)
 
      CALL check_for_restart_files(qs_env, bs_env)
 
      CALL compute_fm_v_xc_gamma(qs_env, bs_env)
 
      CALL setup_time_and_frequency_minimax_grid(bs_env)
 
      ! free memory in qs_env; only if one is not calculating the LDOS because
      ! we need real-space grid operations in pw_env, task_list for the LDOS
      ! Recommendation in case of memory issues: first perform GW calculation without calculating
      !                                          LDOS (to safe memor). Then, use GW restart files
      !                                          in a subsequent calculation to calculate the LDOS
      IF (.NOT. bs_env%do_ldos) THEN
         CALL qs_env_part_release(qs_env)
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE create_and_init_bs_env_for_gw
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE de_init_bs_env(bs_env)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'de_init_bs_env'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
      ! deallocate quantities here which:
      ! 1. cannot be deallocated in bs_env_release due to circular dependencies
      ! 2. consume a lot of memory and should not be kept until the quantity is
      !    deallocated in bs_env_release
 
      IF (ASSOCIATED(bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(bs_env%nl_3c)
 
      CALL cp_libint_static_cleanup()
 
      CALL timestop(handle)
 
   END SUBROUTINE de_init_bs_env
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
!> \param bs_sec ...
! **************************************************************************************************
   SUBROUTINE read_gw_input_parameters(bs_env, bs_sec)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(section_vals_type), POINTER                   :: bs_sec
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'read_gw_input_parameters'
 
      INTEGER                                            :: handle
      TYPE(section_vals_type), POINTER                   :: gw_sec
 
      CALL timeset(routinen, handle)
 
      NULLIFY (gw_sec)
      gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
 
      CALL section_vals_val_get(gw_sec, "NUM_TIME_FREQ_POINTS", i_val=bs_env%num_time_freq_points)
      CALL section_vals_val_get(gw_sec, "EPS_FILTER", r_val=bs_env%eps_filter)
      CALL section_vals_val_get(gw_sec, "MEMORY_PER_PROC", r_val=bs_env%input_memory_per_proc_GB)
      CALL section_vals_val_get(gw_sec, "APPROX_KP_EXTRAPOL", l_val=bs_env%approx_kp_extrapol)
 
      CALL timestop(handle)
 
   END SUBROUTINE read_gw_input_parameters
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE setup_ao_and_ri_basis_set(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_AO_and_RI_basis_set'
 
      INTEGER                                            :: handle, natom, nkind
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, &
                      qs_kind_set=qs_kind_set, &
                      particle_set=particle_set, &
                      natom=natom, nkind=nkind)
 
      ! set up basis
      ALLOCATE (bs_env%sizes_RI(natom), bs_env%sizes_AO(natom))
      ALLOCATE (bs_env%basis_set_RI(nkind), bs_env%basis_set_AO(nkind))
 
      CALL basis_set_list_setup(bs_env%basis_set_RI, "RI_AUX", qs_kind_set)
      CALL basis_set_list_setup(bs_env%basis_set_AO, "ORB", qs_kind_set)
 
      CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_RI, &
                            basis=bs_env%basis_set_RI)
      CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_AO, &
                            basis=bs_env%basis_set_AO)
 
      CALL timestop(handle)
 
   END SUBROUTINE setup_ao_and_ri_basis_set
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE get_ri_basis_and_basis_function_indices(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'get_RI_basis_and_basis_function_indices'
 
      INTEGER                                            :: handle, i_ri, iatom, ikind, iset, &
                                                            max_ao_bf_per_atom, n_ao_test, n_atom, &
                                                            n_kind, n_ri, nset, nsgf, u
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
      INTEGER, DIMENSION(:), POINTER                     :: l_max, l_min, nsgf_set
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(gto_basis_set_type), POINTER                  :: basis_set_a
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      CALL timeset(routinen, handle)
 
      ! determine RI basis set size
      CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
 
      n_kind = SIZE(qs_kind_set)
      n_atom = bs_env%n_atom
 
      CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
 
      DO ikind = 1, n_kind
         CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
                          basis_type="RI_AUX")
         cpassert(ASSOCIATED(basis_set_a))
      END DO
 
      ALLOCATE (bs_env%i_RI_start_from_atom(n_atom))
      ALLOCATE (bs_env%i_RI_end_from_atom(n_atom))
      ALLOCATE (bs_env%i_ao_start_from_atom(n_atom))
      ALLOCATE (bs_env%i_ao_end_from_atom(n_atom))
 
      n_ri = 0
      DO iatom = 1, n_atom
         bs_env%i_RI_start_from_atom(iatom) = n_ri + 1
         ikind = kind_of(iatom)
         CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
         n_ri = n_ri + nsgf
         bs_env%i_RI_end_from_atom(iatom) = n_ri
      END DO
      bs_env%n_RI = n_ri
 
      max_ao_bf_per_atom = 0
      n_ao_test = 0
      DO iatom = 1, n_atom
         bs_env%i_ao_start_from_atom(iatom) = n_ao_test + 1
         ikind = kind_of(iatom)
         CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
         n_ao_test = n_ao_test + nsgf
         bs_env%i_ao_end_from_atom(iatom) = n_ao_test
         max_ao_bf_per_atom = max(max_ao_bf_per_atom, nsgf)
      END DO
      cpassert(n_ao_test == bs_env%n_ao)
      bs_env%max_AO_bf_per_atom = max_ao_bf_per_atom
 
      ALLOCATE (bs_env%l_RI(n_ri))
      i_ri = 0
      DO iatom = 1, n_atom
         ikind = kind_of(iatom)
 
         nset = bs_env%basis_set_RI(ikind)%gto_basis_set%nset
         l_max => bs_env%basis_set_RI(ikind)%gto_basis_set%lmax
         l_min => bs_env%basis_set_RI(ikind)%gto_basis_set%lmin
         nsgf_set => bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf_set
 
         DO iset = 1, nset
            cpassert(l_max(iset) == l_min(iset))
            bs_env%l_RI(i_ri + 1:i_ri + nsgf_set(iset)) = l_max(iset)
            i_ri = i_ri + nsgf_set(iset)
         END DO
 
      END DO
      cpassert(i_ri == n_ri)
 
      u = bs_env%unit_nr
 
      IF (u > 0) THEN
         WRITE (unit=u, fmt="(T2,A)") " "
         WRITE (unit=u, fmt="(T2,2A,T75,I8)") "Number of auxiliary Gaussian basis functions ", &
            χε"for , , W", n_ri
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE get_ri_basis_and_basis_function_indices
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
!> \param kpoints ...
! **************************************************************************************************
   SUBROUTINE setup_kpoints_chi_eps_w(qs_env, bs_env, kpoints)
 
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(kpoint_type), POINTER                         :: kpoints
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_kpoints_chi_eps_W'
 
      INTEGER                                            :: handle, i_dim, n_dim, nkp, nkp_extra, &
                                                            nkp_orig, u
      INTEGER, DIMENSION(3)                              :: nkp_grid, nkp_grid_extra, periodic
      REAL(kind=dp)                                      :: exp_s_p, n_dim_inv
 
      CALL timeset(routinen, handle)
 
      ! routine adapted from mp2_integrals.F
      NULLIFY (kpoints)
      CALL kpoint_create(kpoints)
 
      kpoints%kp_scheme = "GENERAL"
 
      periodic(1:3) = bs_env%periodic(1:3)
 
      DO i_dim = 1, 3
 
         cpassert(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
 
         IF (periodic(i_dim) == 1) nkp_grid(i_dim) = 4
         IF (periodic(i_dim) == 0) nkp_grid(i_dim) = 1
 
         IF (periodic(i_dim) == 1) THEN
            ! only even k-meshes in periodic direction implemented
            cpassert(modulo(nkp_grid(i_dim), 2) == 0)
         END IF
         IF (periodic(i_dim) == 0) THEN
            ! single k-kpoint in non-periodic direction needed
            cpassert(nkp_grid(i_dim) == 1)
         END IF
 
         IF (periodic(i_dim) == 1) nkp_grid_extra(i_dim) = nkp_grid(i_dim) + 2
         IF (periodic(i_dim) == 0) nkp_grid_extra(i_dim) = 1
 
      END DO
 
      nkp_orig = max(nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2, 1)
 
      nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
 
      nkp = nkp_orig + nkp_extra
 
      kpoints%nkp_grid(1:3) = nkp_grid(1:3)
      kpoints%nkp = nkp
 
      bs_env%nkp_grid_chi_eps_W_orig(1:3) = nkp_grid(1:3)
      bs_env%nkp_grid_chi_eps_W_extra(1:3) = nkp_grid_extra(1:3)
      bs_env%nkp_chi_eps_W_orig = nkp_orig
      bs_env%nkp_chi_eps_W_extra = nkp_extra
      bs_env%nkp_chi_eps_W_orig_plus_extra = nkp
 
      ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
      ALLOCATE (bs_env%wkp_no_extra(nkp), bs_env%wkp_s_p(nkp))
 
      CALL compute_xkp(kpoints%xkp, 1, nkp_orig, nkp_grid)
      CALL compute_xkp(kpoints%xkp, nkp_orig + 1, nkp, nkp_grid_extra)
 
      n_dim = sum(periodic)
      IF (n_dim == 0) THEN
         ! molecules
         kpoints%wkp(1) = 1.0_dp
         bs_env%wkp_s_p(1) = 1.0_dp
         bs_env%wkp_no_extra(1) = 1.0_dp
      ELSE
 
         n_dim_inv = 1.0_dp/real(n_dim, kind=dp)
 
         ! k-point weights are chosen to automatically extrapolate the k-point mesh
         CALL compute_wkp(kpoints%wkp(1:nkp_orig), nkp_orig, nkp_extra, n_dim_inv)
         CALL compute_wkp(kpoints%wkp(nkp_orig + 1:nkp), nkp_extra, nkp_orig, n_dim_inv)
 
         bs_env%wkp_no_extra(1:nkp_orig) = 0.0_dp
         bs_env%wkp_no_extra(nkp_orig + 1:nkp) = 1.0_dp/real(nkp_extra, kind=dp)
 
         IF (n_dim == 3) THEN
            ! W_PQ(k) for an s-function P and a p-function Q diverges as 1/k at k=0
            ! (instead of 1/k^2 for P and Q both being s-functions).
            exp_s_p = 2.0_dp*n_dim_inv
            CALL compute_wkp(bs_env%wkp_s_p(1:nkp_orig), nkp_orig, nkp_extra, exp_s_p)
            CALL compute_wkp(bs_env%wkp_s_p(nkp_orig + 1:nkp), nkp_extra, nkp_orig, exp_s_p)
         ELSE
            bs_env%wkp_s_p(1:nkp) = bs_env%wkp_no_extra(1:nkp)
         END IF
 
      END IF
 
      IF (bs_env%approx_kp_extrapol) THEN
         bs_env%wkp_orig = 1.0_dp/real(nkp_orig, kind=dp)
      END IF
 
      CALL kpoint_init_cell_index_simple(kpoints, qs_env)
 
      ! heuristic parameter: how many k-points for χ, ε, and W are used simultaneously
      ! (less simultaneous k-points: less memory, but more computational effort because of
      !  recomputation of V(k))
      bs_env%nkp_chi_eps_W_batch = 4
 
      bs_env%num_chi_eps_W_batches = (bs_env%nkp_chi_eps_W_orig_plus_extra - 1)/ &
                                     bs_env%nkp_chi_eps_W_batch + 1
 
      u = bs_env%unit_nr
 
      IF (u > 0) THEN
         WRITE (unit=u, fmt="(T2,A)") " "
         WRITE (unit=u, fmt="(T2,1A,T71,3I4)") χε"K-point mesh 1 for , , W", nkp_grid(1:3)
         WRITE (unit=u, fmt="(T2,2A,T71,3I4)") χε"K-point mesh 2 for , , W ", &
            "(for k-point extrapolation of W)", nkp_grid_extra(1:3)
         WRITE (unit=u, fmt="(T2,A,T80,L)") "Approximate the k-point extrapolation", &
            bs_env%approx_kp_extrapol
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE setup_kpoints_chi_eps_w
 
! **************************************************************************************************
!> \brief ...
!> \param kpoints ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE kpoint_init_cell_index_simple(kpoints, qs_env)
 
      TYPE(kpoint_type), POINTER                         :: kpoints
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'kpoint_init_cell_index_simple'
 
      INTEGER                                            :: handle
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
         POINTER                                         :: sab_orb
 
      CALL timeset(routinen, handle)
 
      NULLIFY (dft_control, para_env, sab_orb)
      CALL get_qs_env(qs_env=qs_env, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb)
      CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control)
 
      CALL timestop(handle)
 
   END SUBROUTINE kpoint_init_cell_index_simple
 
! **************************************************************************************************
!> \brief ...
!> \param xkp ...
!> \param ikp_start ...
!> \param ikp_end ...
!> \param grid ...
! **************************************************************************************************
   SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
 
      REAL(kind=dp), DIMENSION(:, :), POINTER            :: xkp
      INTEGER                                            :: ikp_start, ikp_end
      INTEGER, DIMENSION(3)                              :: grid
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'compute_xkp'
 
      INTEGER                                            :: handle, i, ix, iy, iz
 
      CALL timeset(routinen, handle)
 
      i = ikp_start
      DO ix = 1, grid(1)
         DO iy = 1, grid(2)
            DO iz = 1, grid(3)
 
               IF (i > ikp_end) cycle
 
               xkp(1, i) = real(2*ix - grid(1) - 1, kind=dp)/(2._dp*real(grid(1), kind=dp))
               xkp(2, i) = real(2*iy - grid(2) - 1, kind=dp)/(2._dp*real(grid(2), kind=dp))
               xkp(3, i) = real(2*iz - grid(3) - 1, kind=dp)/(2._dp*real(grid(3), kind=dp))
               i = i + 1
 
            END DO
         END DO
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE compute_xkp
 
! **************************************************************************************************
!> \brief ...
!> \param wkp ...
!> \param nkp_1 ...
!> \param nkp_2 ...
!> \param exponent ...
! **************************************************************************************************
   SUBROUTINE compute_wkp(wkp, nkp_1, nkp_2, exponent)
      REAL(kind=dp), DIMENSION(:)                        :: wkp
      INTEGER                                            :: nkp_1, nkp_2
      REAL(kind=dp)                                      :: exponent
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'compute_wkp'
 
      INTEGER                                            :: handle
      REAL(kind=dp)                                      :: nkp_ratio
 
      CALL timeset(routinen, handle)
 
      nkp_ratio = real(nkp_2, kind=dp)/real(nkp_1, kind=dp)
 
      wkp(:) = 1.0_dp/real(nkp_1, kind=dp)/(1.0_dp - nkp_ratio**exponent)
 
      CALL timestop(handle)
 
   END SUBROUTINE
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE allocate_and_fill_matrices_and_arrays(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'allocate_and_fill_matrices_and_arrays'
 
      INTEGER                                            :: handle, i_t, num_time_freq_points
      TYPE(cp_blacs_env_type), POINTER                   :: blacs_env, blacs_env_tensor
      TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_ri_global
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      TYPE(mp_para_env_type), POINTER                    :: para_env
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env, matrix_s=matrix_s)
 
      fm_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
 
      CALL cp_fm_create(bs_env%fm_Gocc, fm_struct)
      CALL cp_fm_create(bs_env%fm_Gvir, fm_struct)
      CALL cp_fm_create(bs_env%fm_h_G0W0_Gamma, fm_struct)
 
      NULLIFY (fm_struct_ri_global)
      CALL cp_fm_struct_create(fm_struct_ri_global, context=blacs_env, nrow_global=bs_env%n_RI, &
                               ncol_global=bs_env%n_RI, para_env=para_env)
      CALL cp_fm_create(bs_env%fm_RI_RI, fm_struct_ri_global)
      CALL cp_fm_create(bs_env%fm_chi_Gamma_freq, fm_struct_ri_global)
      CALL cp_fm_create(bs_env%fm_W_MIC_freq, fm_struct_ri_global)
      IF (bs_env%approx_kp_extrapol) THEN
         CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_extra, fm_struct_ri_global)
         CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_no_extra, fm_struct_ri_global)
         CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_extra, 0.0_dp)
         CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_no_extra, 0.0_dp)
      END IF
      CALL cp_fm_struct_release(fm_struct_ri_global)
 
      ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, bs_env%nkp_DOS, bs_env%n_spin))
 
      num_time_freq_points = bs_env%num_time_freq_points
 
      ALLOCATE (bs_env%imag_freq_points(num_time_freq_points))
      ALLOCATE (bs_env%imag_time_points(num_time_freq_points))
      ALLOCATE (bs_env%weights_cos_t_to_w(num_time_freq_points, num_time_freq_points))
      ALLOCATE (bs_env%weights_cos_w_to_t(num_time_freq_points, num_time_freq_points))
      ALLOCATE (bs_env%weights_sin_t_to_w(num_time_freq_points, num_time_freq_points))
 
      ! create blacs_env for subgroups of tensor operations
      NULLIFY (blacs_env_tensor)
      CALL cp_blacs_env_create(blacs_env=blacs_env_tensor, para_env=bs_env%para_env_tensor)
 
      ! allocate dbcsr matrices in the tensor subgroup; actually, one only needs a small
      ! subset of blocks in the tensor subgroup, however, all atomic blocks are allocated.
      ! One might think of creating a dbcsr matrix with only the blocks that are needed
      ! in the tensor subgroup
      CALL create_mat_munu(bs_env%mat_ao_ao_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
                           blacs_env_tensor, do_ri_aux_basis=.false.)
 
      CALL create_mat_munu(bs_env%mat_RI_RI_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
                           blacs_env_tensor, do_ri_aux_basis=.true.)
 
      CALL create_mat_munu(bs_env%mat_RI_RI, qs_env, bs_env%eps_atom_grid_2d_mat, &
                           blacs_env, do_ri_aux_basis=.true.)
 
      CALL cp_blacs_env_release(blacs_env_tensor)
 
      NULLIFY (bs_env%mat_chi_Gamma_tau)
      CALL dbcsr_allocate_matrix_set(bs_env%mat_chi_Gamma_tau, bs_env%num_time_freq_points)
 
      DO i_t = 1, bs_env%num_time_freq_points
         ALLOCATE (bs_env%mat_chi_Gamma_tau(i_t)%matrix)
         CALL dbcsr_create(bs_env%mat_chi_Gamma_tau(i_t)%matrix, template=bs_env%mat_RI_RI%matrix)
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE allocate_and_fill_matrices_and_arrays
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE create_tensors(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'create_tensors'
 
      INTEGER                                            :: handle, n_atom_step, ri_atom
      INTEGER(int_8)                                     :: mem, non_zero_elements_sum, nze
      REAL(dp)                                           :: max_dist_ao_atoms, occ, occupation_sum
      TYPE(dbt_type)                                     :: t_3c_global
      TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_global_array
      TYPE(neighbor_list_3c_type)                        :: nl_3c_global
 
      CALL timeset(routinen, handle)
 
      ! be careful: routine needs bs_env%eps_3c_int which is set in set_heuristic_parameters
      CALL init_interaction_radii(bs_env)
 
      ! split blocks does not improve load balancing/efficienfy for tensor contraction, so we go
      ! with the standard atomic blocks
      CALL create_3c_t(bs_env%t_RI_AO__AO, bs_env%para_env_tensor, "(RI AO | AO)", [1, 2], [3], &
                       bs_env%sizes_RI, bs_env%sizes_AO, &
                       create_nl_3c=.true., nl_3c=bs_env%nl_3c, qs_env=qs_env)
      CALL create_3c_t(bs_env%t_RI__AO_AO, bs_env%para_env_tensor, "(RI | AO AO)", [1], [2, 3], &
                       bs_env%sizes_RI, bs_env%sizes_AO)
 
      CALL create_2c_t(bs_env, bs_env%sizes_RI, bs_env%sizes_AO)
 
      ! check the sparsity of 3c integral tensor (µν|P); calculate maximum distance between
      ! AO atoms µ, ν where at least a single integral (µν|P) is larger than the filter threshold
      CALL create_3c_t(t_3c_global, bs_env%para_env, "(RI AO | AO)", [1, 2], [3], &
                       bs_env%sizes_RI, bs_env%sizes_AO, &
                       create_nl_3c=.true., nl_3c=nl_3c_global, qs_env=qs_env)
 
      CALL m_memory(mem)
      CALL bs_env%para_env%max(mem)
 
      ALLOCATE (t_3c_global_array(1, 1))
      CALL dbt_create(t_3c_global, t_3c_global_array(1, 1))
 
      CALL bs_env%para_env%sync()
      bs_env%t1 = m_walltime()
 
      occupation_sum = 0.0_dp
      non_zero_elements_sum = 0
      max_dist_ao_atoms = 0.0_dp
      n_atom_step = int(sqrt(real(bs_env%n_atom, kind=dp)))
      ! do not compute full 3c integrals at once because it may cause out of memory
      DO ri_atom = 1, bs_env%n_atom, n_atom_step
 
         CALL build_3c_integrals(t_3c_global_array, &
                                 bs_env%eps_filter, &
                                 qs_env, &
                                 nl_3c_global, &
                                 int_eps=bs_env%eps_3c_int, &
                                 basis_i=bs_env%basis_set_RI, &
                                 basis_j=bs_env%basis_set_AO, &
                                 basis_k=bs_env%basis_set_AO, &
                                 bounds_i=[ri_atom, min(ri_atom + n_atom_step - 1, bs_env%n_atom)], &
                                 potential_parameter=bs_env%ri_metric, &
                                 desymmetrize=.false.)
 
         CALL dbt_filter(t_3c_global_array(1, 1), bs_env%eps_filter)
 
         CALL bs_env%para_env%sync()
 
         CALL get_tensor_occupancy(t_3c_global_array(1, 1), nze, occ)
         non_zero_elements_sum = non_zero_elements_sum + nze
         occupation_sum = occupation_sum + occ
 
         CALL get_max_dist_ao_atoms(t_3c_global_array(1, 1), max_dist_ao_atoms, qs_env)
 
         CALL dbt_clear(t_3c_global_array(1, 1))
 
      END DO
 
      bs_env%t2 = m_walltime()
 
      bs_env%occupation_3c_int = occupation_sum
      bs_env%max_dist_AO_atoms = max_dist_ao_atoms
 
      CALL dbt_destroy(t_3c_global)
      CALL dbt_destroy(t_3c_global_array(1, 1))
      DEALLOCATE (t_3c_global_array)
 
      CALL neighbor_list_3c_destroy(nl_3c_global)
 
      IF (bs_env%unit_nr > 0) THEN
         WRITE (bs_env%unit_nr, '(T2,A)') ''
         WRITE (bs_env%unit_nr, '(T2,A,F27.1,A)') &
            µν'Computed 3-center integrals (|P), execution time', bs_env%t2 - bs_env%t1, ' s'
         WRITE (bs_env%unit_nr, '(T2,A,F48.3,A)') µν'Percentage of non-zero (|P)', &
            occupation_sum*100, ' %'
         WRITE (bs_env%unit_nr, '(T2,A,F33.1,A)') µνµν'Max. distance between , in non-zero (|P)', &
            max_dist_ao_atoms*angstrom, ' A'
         WRITE (bs_env%unit_nr, '(T2,2A,I20,A)') 'Required memory if storing all 3-center ', &
            µν'integrals (|P)', int(real(non_zero_elements_sum, kind=dp)*8.0e-9_dp), ' GB'
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE create_tensors
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
!> \param sizes_RI ...
!> \param sizes_AO ...
! **************************************************************************************************
   SUBROUTINE create_2c_t(bs_env, sizes_RI, sizes_AO)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sizes_ri, sizes_ao
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'create_2c_t'
 
      INTEGER                                            :: handle
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: dist_1, dist_2
      INTEGER, DIMENSION(2)                              :: pdims_2d
      TYPE(dbt_pgrid_type)                               :: pgrid_2d
 
      CALL timeset(routinen, handle)
 
      ! inspired from rpa_im_time.F / hfx_types.F
 
      pdims_2d = 0
      CALL dbt_pgrid_create(bs_env%para_env_tensor, pdims_2d, pgrid_2d)
 
      CALL create_2c_tensor(bs_env%t_G, dist_1, dist_2, pgrid_2d, sizes_ao, sizes_ao, &
                            name="(AO | AO)")
      DEALLOCATE (dist_1, dist_2)
      CALL create_2c_tensor(bs_env%t_chi, dist_1, dist_2, pgrid_2d, sizes_ri, sizes_ri, &
                            name="(RI | RI)")
      DEALLOCATE (dist_1, dist_2)
      CALL create_2c_tensor(bs_env%t_W, dist_1, dist_2, pgrid_2d, sizes_ri, sizes_ri, &
                            name="(RI | RI)")
      DEALLOCATE (dist_1, dist_2)
      CALL dbt_pgrid_destroy(pgrid_2d)
 
      CALL timestop(handle)
 
   END SUBROUTINE create_2c_t
 
! **************************************************************************************************
!> \brief ...
!> \param tensor ...
!> \param para_env ...
!> \param tensor_name ...
!> \param map1 ...
!> \param map2 ...
!> \param sizes_RI ...
!> \param sizes_AO ...
!> \param create_nl_3c ...
!> \param nl_3c ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE create_3c_t(tensor, para_env, tensor_name, map1, map2, sizes_RI, sizes_AO, &
                          create_nl_3c, nl_3c, qs_env)
      TYPE(dbt_type)                                     :: tensor
      TYPE(mp_para_env_type), POINTER                    :: para_env
      CHARACTER(LEN=12)                                  :: tensor_name
      INTEGER, DIMENSION(:)                              :: map1, map2
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sizes_ri, sizes_ao
      LOGICAL, OPTIONAL                                  :: create_nl_3c
      TYPE(neighbor_list_3c_type), OPTIONAL              :: nl_3c
      TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'create_3c_t'
 
      INTEGER                                            :: handle, nkind
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: dist_ao_1, dist_ao_2, dist_ri
      INTEGER, DIMENSION(3)                              :: pcoord, pdims, pdims_3d
      LOGICAL                                            :: my_create_nl_3c
      TYPE(dbt_pgrid_type)                               :: pgrid_3d
      TYPE(distribution_3d_type)                         :: dist_3d
      TYPE(mp_cart_type)                                 :: mp_comm_t3c_2
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      CALL timeset(routinen, handle)
 
      pdims_3d = 0
      CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d)
      CALL create_3c_tensor(tensor, dist_ri, dist_ao_1, dist_ao_2, &
                            pgrid_3d, sizes_ri, sizes_ao, sizes_ao, &
                            map1=map1, map2=map2, name=tensor_name)
 
      IF (PRESENT(create_nl_3c)) THEN
         my_create_nl_3c = create_nl_3c
      ELSE
         my_create_nl_3c = .false.
      END IF
 
      IF (my_create_nl_3c) THEN
         CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
         CALL dbt_mp_environ_pgrid(pgrid_3d, pdims, pcoord)
         CALL mp_comm_t3c_2%create(pgrid_3d%mp_comm_2d, 3, pdims)
         CALL distribution_3d_create(dist_3d, dist_ri, dist_ao_1, dist_ao_2, &
                                     nkind, particle_set, mp_comm_t3c_2, own_comm=.true.)
 
         CALL build_3c_neighbor_lists(nl_3c, &
                                      qs_env%bs_env%basis_set_RI, &
                                      qs_env%bs_env%basis_set_AO, &
                                      qs_env%bs_env%basis_set_AO, &
                                      dist_3d, qs_env%bs_env%ri_metric, &
                                      "GW_3c_nl", qs_env, own_dist=.true.)
      END IF
 
      DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
      CALL dbt_pgrid_destroy(pgrid_3d)
 
      CALL timestop(handle)
 
   END SUBROUTINE create_3c_t
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE init_interaction_radii(bs_env)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'init_interaction_radii'
 
      INTEGER                                            :: handle, ibasis
      TYPE(gto_basis_set_type), POINTER                  :: orb_basis, ri_basis
 
      CALL timeset(routinen, handle)
 
      DO ibasis = 1, SIZE(bs_env%basis_set_AO)
 
         orb_basis => bs_env%basis_set_AO(ibasis)%gto_basis_set
         CALL init_interaction_radii_orb_basis(orb_basis, bs_env%eps_3c_int)
 
         ri_basis => bs_env%basis_set_RI(ibasis)%gto_basis_set
         CALL init_interaction_radii_orb_basis(ri_basis, bs_env%eps_3c_int)
 
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE init_interaction_radii
 
! **************************************************************************************************
!> \brief ...
!> \param t_3c_int ...
!> \param max_dist_AO_atoms ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE get_max_dist_ao_atoms(t_3c_int, max_dist_AO_atoms, qs_env)
      TYPE(dbt_type)                                     :: t_3c_int
      REAL(kind=dp)                                      :: max_dist_ao_atoms
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'get_max_dist_AO_atoms'
 
      INTEGER                                            :: atom_1, atom_2, handle, num_cells
      INTEGER, DIMENSION(3)                              :: atom_ind
      INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
      REAL(kind=dp)                                      :: abs_rab
      REAL(kind=dp), DIMENSION(3)                        :: rab
      REAL(kind=dp), DIMENSION(3, 3)                     :: hmat
      TYPE(cell_type), POINTER                           :: cell
      TYPE(dbt_iterator_type)                            :: iter
      TYPE(mp_para_env_type), POINTER                    :: para_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
 
      CALL timeset(routinen, handle)
 
      NULLIFY (cell, particle_set, para_env)
      CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, para_env=para_env)
 
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(t_3c_int, max_dist_AO_atoms, num_cells, index_to_cell, hmat, particle_set, cell) &
!$OMP PRIVATE(iter,atom_ind,rab, abs_rab, atom_1, atom_2)
      CALL dbt_iterator_start(iter, t_3c_int)
      DO WHILE (dbt_iterator_blocks_left(iter))
         CALL dbt_iterator_next_block(iter, atom_ind)
 
         atom_1 = atom_ind(2)
         atom_2 = atom_ind(3)
 
         rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
 
         abs_rab = sqrt(rab(1)**2 + rab(2)**2 + rab(3)**2)
 
         max_dist_ao_atoms = max(max_dist_ao_atoms, abs_rab)
 
      END DO
      CALL dbt_iterator_stop(iter)
!$OMP END PARALLEL
 
      CALL para_env%max(max_dist_ao_atoms)
 
      CALL timestop(handle)
 
   END SUBROUTINE get_max_dist_ao_atoms
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE set_sparsity_parallelization_parameters(bs_env)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'set_sparsity_parallelization_parameters'
 
      INTEGER :: handle, i_ivl, il_ivl, j_ivl, n_atom_per_il_ivl, n_atom_per_ivl, n_intervals_i, &
         n_intervals_inner_loop_atoms, n_intervals_j, u
 
      CALL timeset(routinen, handle)
 
      ! heuristic parameter to prevent out of memory
      bs_env%safety_factor_memory = 0.10_dp
 
      ! choose atomic range for λ ("i_atom"), ν ("j_atom") in
      ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
      ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
      ! such that M and N fit into the memory
      n_atom_per_ivl = int(sqrt(bs_env%safety_factor_memory*bs_env%input_memory_per_proc &
                                *bs_env%group_size_tensor/24/bs_env%n_RI &
                                /sqrt(bs_env%occupation_3c_int)))/bs_env%max_AO_bf_per_atom
 
      n_intervals_i = (bs_env%n_atom_i - 1)/n_atom_per_ivl + 1
      n_intervals_j = (bs_env%n_atom_j - 1)/n_atom_per_ivl + 1
 
      bs_env%n_atom_per_interval_ij = n_atom_per_ivl
      bs_env%n_intervals_i = n_intervals_i
      bs_env%n_intervals_j = n_intervals_j
 
      ALLOCATE (bs_env%i_atom_intervals(2, n_intervals_i))
      ALLOCATE (bs_env%j_atom_intervals(2, n_intervals_j))
 
      DO i_ivl = 1, n_intervals_i
         bs_env%i_atom_intervals(1, i_ivl) = (i_ivl - 1)*n_atom_per_ivl + bs_env%atoms_i(1)
         bs_env%i_atom_intervals(2, i_ivl) = min(i_ivl*n_atom_per_ivl + bs_env%atoms_i(1) - 1, &
                                                 bs_env%atoms_i(2))
      END DO
 
      DO j_ivl = 1, n_intervals_j
         bs_env%j_atom_intervals(1, j_ivl) = (j_ivl - 1)*n_atom_per_ivl + bs_env%atoms_j(1)
         bs_env%j_atom_intervals(2, j_ivl) = min(j_ivl*n_atom_per_ivl + bs_env%atoms_j(1) - 1, &
                                                 bs_env%atoms_j(2))
      END DO
 
      ALLOCATE (bs_env%skip_Sigma_occ(n_intervals_i, n_intervals_j))
      ALLOCATE (bs_env%skip_Sigma_vir(n_intervals_i, n_intervals_j))
      bs_env%skip_Sigma_occ(:, :) = .false.
      bs_env%skip_Sigma_vir(:, :) = .false.
 
      ! choose atomic range for µ and σ ("inner loop (IL) atom") in
      ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
      ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
      n_atom_per_il_ivl = min(int(bs_env%safety_factor_memory*bs_env%input_memory_per_proc &
                                  *bs_env%group_size_tensor/n_atom_per_ivl &
                                  /bs_env%max_AO_bf_per_atom &
                                  /bs_env%n_RI/8/sqrt(bs_env%occupation_3c_int) &
                                  /bs_env%max_AO_bf_per_atom), bs_env%n_atom)
 
      n_intervals_inner_loop_atoms = (bs_env%n_atom - 1)/n_atom_per_il_ivl + 1
 
      bs_env%n_atom_per_IL_interval = n_atom_per_il_ivl
      bs_env%n_intervals_inner_loop_atoms = n_intervals_inner_loop_atoms
 
      ALLOCATE (bs_env%inner_loop_atom_intervals(2, n_intervals_inner_loop_atoms))
      DO il_ivl = 1, n_intervals_inner_loop_atoms
         bs_env%inner_loop_atom_intervals(1, il_ivl) = (il_ivl - 1)*n_atom_per_il_ivl + 1
         bs_env%inner_loop_atom_intervals(2, il_ivl) = min(il_ivl*n_atom_per_il_ivl, bs_env%n_atom)
      END DO
 
      u = bs_env%unit_nr
      IF (u > 0) THEN
         WRITE (u, '(T2,A)') ''
         WRITE (u, '(T2,A,I33)') λντνλτ'Number of i and j atoms in M_P(), N_Q():', n_atom_per_ivl
         WRITE (u, '(T2,A,I18)') µλνµµνµλ'Number of inner loop atoms for  in M_P = sum_ (|P) G_', &
            n_atom_per_il_ivl
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE set_sparsity_parallelization_parameters
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE check_for_restart_files(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'check_for_restart_files'
 
      CHARACTER(LEN=9)                                   :: frmt
      CHARACTER(LEN=default_string_length)               :: f_chi, f_s_n, f_s_p, f_s_x, f_w_t, &
                                                            prefix, project_name
      INTEGER                                            :: handle, i_spin, i_t_or_w, ind, n_spin, &
                                                            num_time_freq_points
      LOGICAL                                            :: chi_exists, sigma_neg_time_exists, &
                                                            sigma_pos_time_exists, &
                                                            sigma_x_spin_exists, w_time_exists
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(section_vals_type), POINTER                   :: input, print_key
 
      CALL timeset(routinen, handle)
 
      num_time_freq_points = bs_env%num_time_freq_points
      n_spin = bs_env%n_spin
 
      ALLOCATE (bs_env%read_chi(num_time_freq_points))
      ALLOCATE (bs_env%calc_chi(num_time_freq_points))
      ALLOCATE (bs_env%Sigma_c_exists(num_time_freq_points, n_spin))
 
      CALL get_qs_env(qs_env, input=input)
 
      logger => cp_get_default_logger()
      print_key => section_vals_get_subs_vals(input, 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART')
      project_name = cp_print_key_generate_filename(logger, print_key, extension="", &
                                                    my_local=.false.)
      WRITE (prefix, '(2A)') trim(project_name), "-RESTART_"
      bs_env%prefix = prefix
 
      bs_env%all_W_exist = .true.
 
      DO i_t_or_w = 1, num_time_freq_points
 
         IF (i_t_or_w < 10) THEN
            WRITE (frmt, '(A)') '(3A,I1,A)'
            WRITE (f_chi, frmt) trim(prefix), bs_env%chi_name, "_0", i_t_or_w, ".matrix"
            WRITE (f_w_t, frmt) trim(prefix), bs_env%W_time_name, "_0", i_t_or_w, ".matrix"
         ELSE IF (i_t_or_w < 100) THEN
            WRITE (frmt, '(A)') '(3A,I2,A)'
            WRITE (f_chi, frmt) trim(prefix), bs_env%chi_name, "_", i_t_or_w, ".matrix"
            WRITE (f_w_t, frmt) trim(prefix), bs_env%W_time_name, "_", i_t_or_w, ".matrix"
         ELSE
            cpabort('Please implement more than 99 time/frequency points.')
         END IF
 
         INQUIRE (file=trim(f_chi), exist=chi_exists)
         INQUIRE (file=trim(f_w_t), exist=w_time_exists)
 
         bs_env%read_chi(i_t_or_w) = chi_exists
         bs_env%calc_chi(i_t_or_w) = .NOT. chi_exists
 
         bs_env%all_W_exist = bs_env%all_W_exist .AND. w_time_exists
 
         ! the self-energy is spin-dependent
         DO i_spin = 1, n_spin
 
            ind = i_t_or_w + (i_spin - 1)*num_time_freq_points
 
            IF (ind < 10) THEN
               WRITE (frmt, '(A)') '(3A,I1,A)'
               WRITE (f_s_p, frmt) trim(prefix), bs_env%Sigma_p_name, "_0", ind, ".matrix"
               WRITE (f_s_n, frmt) trim(prefix), bs_env%Sigma_n_name, "_0", ind, ".matrix"
            ELSE IF (i_t_or_w < 100) THEN
               WRITE (frmt, '(A)') '(3A,I2,A)'
               WRITE (f_s_p, frmt) trim(prefix), bs_env%Sigma_p_name, "_", ind, ".matrix"
               WRITE (f_s_n, frmt) trim(prefix), bs_env%Sigma_n_name, "_", ind, ".matrix"
            END IF
 
            INQUIRE (file=trim(f_s_p), exist=sigma_pos_time_exists)
            INQUIRE (file=trim(f_s_n), exist=sigma_neg_time_exists)
 
            bs_env%Sigma_c_exists(i_t_or_w, i_spin) = sigma_pos_time_exists .AND. &
                                                      sigma_neg_time_exists
 
         END DO
 
      END DO
 
      IF (bs_env%all_W_exist) THEN
         bs_env%read_chi(:) = .false.
         bs_env%calc_chi(:) = .false.
      END IF
 
      bs_env%Sigma_x_exists = .true.
      DO i_spin = 1, n_spin
         WRITE (f_s_x, '(3A,I1,A)') trim(prefix), bs_env%Sigma_x_name, "_0", i_spin, ".matrix"
         INQUIRE (file=trim(f_s_x), exist=sigma_x_spin_exists)
         bs_env%Sigma_x_exists = bs_env%Sigma_x_exists .AND. sigma_x_spin_exists
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE check_for_restart_files
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE set_parallelization_parameters(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'set_parallelization_parameters'
 
      INTEGER                                            :: color_sub, dummy_1, dummy_2, handle, &
                                                            num_pe, num_t_groups, u
      INTEGER(KIND=int_8)                                :: mem
      TYPE(mp_para_env_type), POINTER                    :: para_env
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, para_env=para_env)
 
      num_pe = para_env%num_pe
      ! use all processors for the group (in principle, number could be changed, but performance
      ! seems to be best for a single group with all MPI processes per group)
      bs_env%group_size_tensor = num_pe
 
      ! group_size_tensor must divide num_pe without rest; otherwise everything will be complicated
      IF (modulo(num_pe, bs_env%group_size_tensor) .NE. 0) THEN
         CALL find_good_group_size(num_pe, bs_env%group_size_tensor)
      END IF
 
      ! para_env_tensor for tensor subgroups
      color_sub = para_env%mepos/bs_env%group_size_tensor
      bs_env%tensor_group_color = color_sub
 
      ALLOCATE (bs_env%para_env_tensor)
      CALL bs_env%para_env_tensor%from_split(para_env, color_sub)
 
      num_t_groups = para_env%num_pe/bs_env%group_size_tensor
      bs_env%num_tensor_groups = num_t_groups
 
      CALL get_i_j_atoms(bs_env%atoms_i, bs_env%atoms_j, bs_env%n_atom_i, bs_env%n_atom_j, &
                         color_sub, bs_env)
 
      ALLOCATE (bs_env%atoms_i_t_group(2, num_t_groups))
      ALLOCATE (bs_env%atoms_j_t_group(2, num_t_groups))
      DO color_sub = 0, num_t_groups - 1
         CALL get_i_j_atoms(bs_env%atoms_i_t_group(1:2, color_sub + 1), &
                            bs_env%atoms_j_t_group(1:2, color_sub + 1), &
                            dummy_1, dummy_2, color_sub, bs_env)
      END DO
 
      CALL m_memory(mem)
      CALL bs_env%para_env%max(mem)
 
      bs_env%input_memory_per_proc = int(bs_env%input_memory_per_proc_GB*1.0e9_dp, kind=int_8)
 
      u = bs_env%unit_nr
      IF (u > 0) THEN
         WRITE (u, '(T2,A)') ''
         WRITE (u, '(T2,A,I47)') 'Group size for tensor operations', bs_env%group_size_tensor
         WRITE (u, '(T2,A)') ''
         WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
            bs_env%input_memory_per_proc_GB, ' GB'
         WRITE (u, '(T2,A,F35.1,A)') 'Used memory per MPI process before GW run', &
            REAL(mem, kind=dp)/1.e9_dp, ' GB'
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE set_parallelization_parameters
 
! **************************************************************************************************
!> \brief ...
!> \param num_pe ...
!> \param group_size ...
! **************************************************************************************************
   SUBROUTINE find_good_group_size(num_pe, group_size)
 
      INTEGER                                            :: num_pe, group_size
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'find_good_group_size'
 
      INTEGER                                            :: group_size_minus, group_size_orig, &
                                                            group_size_plus, handle, i_diff
 
      CALL timeset(routinen, handle)
 
      group_size_orig = group_size
 
      DO i_diff = 1, num_pe
 
         group_size_minus = group_size - i_diff
 
         IF (modulo(num_pe, group_size_minus) == 0 .AND. group_size_minus > 0) THEN
            group_size = group_size_minus
            EXIT
         END IF
 
         group_size_plus = group_size + i_diff
 
         IF (modulo(num_pe, group_size_plus) == 0 .AND. group_size_plus <= num_pe) THEN
            group_size = group_size_plus
            EXIT
         END IF
 
      END DO
 
      IF (group_size_orig == group_size) cpabort("Group size error")
 
      CALL timestop(handle)
 
   END SUBROUTINE find_good_group_size
 
! **************************************************************************************************
!> \brief ...
!> \param atoms_i ...
!> \param atoms_j ...
!> \param n_atom_i ...
!> \param n_atom_j ...
!> \param color_sub ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE get_i_j_atoms(atoms_i, atoms_j, n_atom_i, n_atom_j, color_sub, bs_env)
 
      INTEGER, DIMENSION(2)                              :: atoms_i, atoms_j
      INTEGER                                            :: n_atom_i, n_atom_j, color_sub
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'get_i_j_atoms'
 
      INTEGER                                            :: handle, i_atoms_per_group, i_group, &
                                                            ipcol, ipcol_loop, iprow, iprow_loop, &
                                                            j_atoms_per_group, npcol, nprow
 
      CALL timeset(routinen, handle)
 
      ! create a square mesh of tensor groups for iatom and jatom; code from blacs_env_create
      CALL square_mesh(nprow, npcol, bs_env%num_tensor_groups)
 
      i_group = 0
      DO ipcol_loop = 0, npcol - 1
         DO iprow_loop = 0, nprow - 1
            IF (i_group == color_sub) THEN
               iprow = iprow_loop
               ipcol = ipcol_loop
            END IF
            i_group = i_group + 1
         END DO
      END DO
 
      IF (modulo(bs_env%n_atom, nprow) == 0) THEN
         i_atoms_per_group = bs_env%n_atom/nprow
      ELSE
         i_atoms_per_group = bs_env%n_atom/nprow + 1
      END IF
 
      IF (modulo(bs_env%n_atom, npcol) == 0) THEN
         j_atoms_per_group = bs_env%n_atom/npcol
      ELSE
         j_atoms_per_group = bs_env%n_atom/npcol + 1
      END IF
 
      atoms_i(1) = iprow*i_atoms_per_group + 1
      atoms_i(2) = min((iprow + 1)*i_atoms_per_group, bs_env%n_atom)
      n_atom_i = atoms_i(2) - atoms_i(1) + 1
 
      atoms_j(1) = ipcol*j_atoms_per_group + 1
      atoms_j(2) = min((ipcol + 1)*j_atoms_per_group, bs_env%n_atom)
      n_atom_j = atoms_j(2) - atoms_j(1) + 1
 
      CALL timestop(handle)
 
   END SUBROUTINE get_i_j_atoms
 
! **************************************************************************************************
!> \brief ...
!> \param nprow ...
!> \param npcol ...
!> \param nproc ...
! **************************************************************************************************
   SUBROUTINE square_mesh(nprow, npcol, nproc)
      INTEGER                                            :: nprow, npcol, nproc
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'square_mesh'
 
      INTEGER                                            :: gcd_max, handle, ipe, jpe
 
      CALL timeset(routinen, handle)
 
      gcd_max = -1
      DO ipe = 1, ceiling(sqrt(real(nproc, dp)))
         jpe = nproc/ipe
         IF (ipe*jpe .NE. nproc) cycle
         IF (gcd(ipe, jpe) >= gcd_max) THEN
            nprow = ipe
            npcol = jpe
            gcd_max = gcd(ipe, jpe)
         END IF
      END DO
 
      CALL timestop(handle)
 
   END SUBROUTINE square_mesh
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
      TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'set_heuristic_parameters'
 
      INTEGER                                            :: handle
 
      CALL timeset(routinen, handle)
 
      ! use the same threshold for computing 3-center integrals (µν|P) as for filtering
      ! tensor operations
      bs_env%eps_3c_int = bs_env%eps_filter
 
      ! Determines number of cells used for summing the cells R in the Coulomb matrix,
      ! V_PQ(k) = \sum_R <P,cell=0 | 1/r | Q,cell=R>. SIZE_LATTICE_SUM_V 3 gives
      ! good convergence
      bs_env%size_lattice_sum_V = 3
 
      ! for generating numerically stable minimax Fourier integration weights
      bs_env%num_points_per_magnitude = 200
 
      ! for periodic systems and for 20 minimax points, we use a regularized minimax mesh
      ! (from experience: regularized minimax meshes converges faster for periodic systems
      !  and for 20 pts)
      IF (sum(bs_env%periodic) .NE. 0 .OR. bs_env%num_time_freq_points == 20) THEN
         bs_env%regularization_minimax = 1.0e-6_dp
      ELSE
         bs_env%regularization_minimax = 0.0_dp
      END IF
 
      bs_env%stabilize_exp = 70.0_dp
      bs_env%eps_atom_grid_2d_mat = 1.0e-50_dp
 
      ! only use interval ω in [0, 27.211 eV] (1 Hartree = 27.211 eV) for virt, and ω in
      ! [-27.211 eV, 0] for occ for use in analytic continuation of
      ! self-energy Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
      bs_env%freq_max_fit = 1.0_dp
 
      ! use a 16-parameter Padé fit
      bs_env%nparam_pade = 16
 
      ! minimum block size for tensor operations, taken from MP2/RPA input
      bs_env%min_block_size = 5
 
      ! resolution of the identity with the truncated Coulomb metric, cutoff radius 3 Angström
      bs_env%ri_metric%potential_type = do_potential_truncated
      bs_env%ri_metric%omega = 0.0_dp
      ! cutoff radius = 3 Angström
      bs_env%ri_metric%cutoff_radius = 3.0_dp/angstrom
      bs_env%ri_metric%filename = "t_c_g.dat"
 
      bs_env%eps_eigval_mat_RI = 0.0_dp
 
      ! for periodic systems, we use the regularized resolution of the identity
      IF (sum(bs_env%periodic) == 0) THEN
         bs_env%regularization_RI = 0.0_dp
      ELSE
         bs_env%regularization_RI = 1.0e-2_dp
      END IF
 
      ! truncated Coulomb operator for exchange self-energy
      ! (see details in Guidon, VandeVondele, Hutter, JCTC 5, 3010 (2009) and references therein)
      CALL setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, bs_env%trunc_coulomb, &
                                                            rel_cutoff_trunc_coulomb_ri_x=0.5_dp)
 
      CALL timestop(handle)
 
   END SUBROUTINE set_heuristic_parameters
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE print_header_and_input_parameters(bs_env)
 
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'print_header_and_input_parameters'
 
      INTEGER                                            :: handle, u
 
      CALL timeset(routinen, handle)
 
      u = bs_env%unit_nr
 
      IF (u > 0) THEN
         WRITE (u, *) ' '
         WRITE (u, '(T2,2A)') '------------------------------------------------', &
            '-------------------------------'
         WRITE (u, '(T2,2A)') '-                                               ', &
            '                              -'
         WRITE (u, '(T2,2A)') '-                              GW CALCULATION   ', &
            '                              -'
         WRITE (u, '(T2,2A)') '-                                               ', &
            '                              -'
         WRITE (u, '(T2,2A)') '------------------------------------------------', &
            '-------------------------------'
         WRITE (u, '(T2,A)') ' '
         WRITE (u, '(T2,A,I45)') 'Input: Number of time/freq. points', bs_env%num_time_freq_points
         WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
            bs_env%eps_filter
      END IF
 
      CALL timestop(handle)
 
   END SUBROUTINE print_header_and_input_parameters
 
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE compute_fm_v_xc_gamma(qs_env, bs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_fm_V_xc_Gamma'
 
      INTEGER                                            :: handle, ispin, myfun, nimages
      REAL(kind=dp)                                      :: energy_ex, energy_exc, energy_total
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_ks_without_v_xc
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(qs_energy_type), POINTER                      :: energy
      TYPE(section_vals_type), POINTER                   :: input, xc_section
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, input=input, energy=energy, dft_control=dft_control)
 
      ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
      nimages = dft_control%nimages
      dft_control%nimages = 1
 
      ! we need to reset XC functional, therefore, get XC input
      xc_section => section_vals_get_subs_vals(input, "DFT%XC")
      CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
      CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=xc_none)
 
      ! save the energy before the energy gets updated
      energy_total = energy%total
      energy_exc = energy%exc
      energy_ex = energy%ex
 
      NULLIFY (mat_ks_without_v_xc)
      CALL dbcsr_allocate_matrix_set(mat_ks_without_v_xc, bs_env%n_spin)
 
      DO ispin = 1, bs_env%n_spin
         ALLOCATE (mat_ks_without_v_xc(ispin)%matrix)
         CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
      END DO
 
      ! calculate KS-matrix without XC
      CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.false., just_energy=.false., &
                                        ext_ks_matrix=mat_ks_without_v_xc)
 
      DO ispin = 1, bs_env%n_spin
         ! transfer dbcsr matrix to fm
         CALL cp_fm_create(bs_env%fm_V_xc_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
         CALL copy_dbcsr_to_fm(mat_ks_without_v_xc(ispin)%matrix, bs_env%fm_V_xc_Gamma(ispin))
 
         ! finally compute the xc potential in the AO basis
         CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_Gamma(ispin), &
                                  beta=1.0_dp, matrix_b=bs_env%fm_ks_Gamma(ispin))
      END DO
 
      ! set back the energy
      energy%total = energy_total
      energy%exc = energy_exc
      energy%ex = energy_ex
 
      ! set back nimages
      dft_control%nimages = nimages
 
      ! set the DFT functional and HF fraction back
      CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
                                i_val=myfun)
 
      CALL dbcsr_deallocate_matrix_set(mat_ks_without_v_xc)
 
      CALL timestop(handle)
 
   END SUBROUTINE compute_fm_v_xc_gamma
 
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
! **************************************************************************************************
   SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env)
      TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_time_and_frequency_minimax_grid'
 
      INTEGER                                            :: handle, homo, i_w, ierr, j_w, n_mo, &
                                                            num_time_freq_points, u
      REAL(kind=dp)                                      :: e_max, e_min, e_range, max_error_min
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: points_and_weights
 
      CALL timeset(routinen, handle)
 
      homo = bs_env%n_occ(1)
      n_mo = bs_env%n_ao
      num_time_freq_points = bs_env%num_time_freq_points
 
      ! minimum and maximum difference between eigenvalues of unoccupied and an occupied MOs
      e_min = minval(bs_env%eigenval_scf_Gamma(homo + 1, :)) - &
              maxval(bs_env%eigenval_scf_Gamma(homo, :))
      e_max = maxval(bs_env%eigenval_scf_Gamma(n_mo, :)) - &
              minval(bs_env%eigenval_scf_Gamma(1, :))
 
      e_range = e_max/e_min
 
      ALLOCATE (points_and_weights(2*num_time_freq_points))
 
      ! frequency points
      IF (num_time_freq_points .LE. 20) THEN
         CALL get_rpa_minimax_coeff(num_time_freq_points, e_range, points_and_weights, ierr, .false.)
      ELSE
         CALL get_rpa_minimax_coeff_larger_grid(num_time_freq_points, e_range, points_and_weights)
      END IF
 
      ! one needs to scale the minimax grids, see Azizi, Wilhelm, Golze, Panades-Barrueta,
      ! Giantomassi, Rinke, Draxl, Gonze et al., 2 publications
      bs_env%imag_freq_points(:) = points_and_weights(1:num_time_freq_points)*e_min
 
      ! determine number of fit points in the interval [0,ω_max] for virt, or [-ω_max,0] for occ
      bs_env%num_freq_points_fit = 0
      DO i_w = 1, bs_env%num_time_freq_points
         IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
            bs_env%num_freq_points_fit = bs_env%num_freq_points_fit + 1
         END IF
      END DO
 
      ! iω values for the analytic continuation Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
      ALLOCATE (bs_env%imag_freq_points_fit(bs_env%num_freq_points_fit))
      j_w = 0
      DO i_w = 1, bs_env%num_time_freq_points
         IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
            j_w = j_w + 1
            bs_env%imag_freq_points_fit(j_w) = bs_env%imag_freq_points(i_w)
         END IF
      END DO
 
      ! reset the number of Padé parameters if smaller than the number of
      ! imaginary-frequency points for the fit
      IF (bs_env%num_freq_points_fit < bs_env%nparam_pade) THEN
         bs_env%nparam_pade = bs_env%num_freq_points_fit
      END IF
 
      ! time points
      IF (num_time_freq_points .LE. 20) THEN
         CALL get_exp_minimax_coeff(num_time_freq_points, e_range, points_and_weights)
      ELSE
         CALL get_exp_minimax_coeff_gw(num_time_freq_points, e_range, points_and_weights)
      END IF
 
      bs_env%imag_time_points(:) = points_and_weights(1:num_time_freq_points)/(2.0_dp*e_min)
 
      DEALLOCATE (points_and_weights)
 
      ! in minimax grids, Fourier transforms t -> w and w -> t are split using
      ! e^(iwt) = cos(wt) + i sin(wt); we thus calculate weights for trafos with a cos and
      ! sine prefactor; details in Azizi, Wilhelm, Golze, Giantomassi, Panades-Barrueta,
      ! Rinke, Draxl, Gonze et al., 2 publications
 
      ! cosine transform weights imaginary time to imaginary frequency
      CALL get_l_sq_wghts_cos_tf_t_to_w(num_time_freq_points, &
                                        bs_env%imag_time_points, &
                                        bs_env%weights_cos_t_to_w, &
                                        bs_env%imag_freq_points, &
                                        e_min, e_max, max_error_min, &
                                        bs_env%num_points_per_magnitude, &
                                        bs_env%regularization_minimax)
 
      ! cosine transform weights imaginary frequency to imaginary time
      CALL get_l_sq_wghts_cos_tf_w_to_t(num_time_freq_points, &
                                        bs_env%imag_time_points, &
                                        bs_env%weights_cos_w_to_t, &
                                        bs_env%imag_freq_points, &
                                        e_min, e_max, max_error_min, &
                                        bs_env%num_points_per_magnitude, &
                                        bs_env%regularization_minimax)
 
      ! sine transform weights imaginary time to imaginary frequency
      CALL get_l_sq_wghts_sin_tf_t_to_w(num_time_freq_points, &
                                        bs_env%imag_time_points, &
                                        bs_env%weights_sin_t_to_w, &
                                        bs_env%imag_freq_points, &
                                        e_min, e_max, max_error_min, &
                                        bs_env%num_points_per_magnitude, &
                                        bs_env%regularization_minimax)
 
      u = bs_env%unit_nr
      IF (u > 0) THEN
         WRITE (u, '(T2,A)') ''
         WRITE (u, '(T2,A,F44.2)') Γ'SCF direct band gap at -point (eV)', e_min*evolt
         WRITE (u, '(T2,A,F42.2)') Γ'Max. SCF eigval diff. at -point (eV)', e_max*evolt
         WRITE (u, '(T2,A,F55.2)') 'E-Range for minimax grid', e_range
         WRITE (u, '(T2,A,I27)') é'Number of Pad parameters for analytic continuation:', &
            bs_env%nparam_pade
         WRITE (u, '(T2,A)') ''
      END IF
      CALL timestop(handle)
 
   END SUBROUTINE setup_time_and_frequency_minimax_grid
 
END MODULE gw_utils