LCOV - code coverage report
Current view: top level - src - gw_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 1145 1243 92.1 %
Date: 2025-02-18 08:24:35 Functions: 43 46 93.5 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief
      10             : !> \author Jan Wilhelm
      11             : !> \date 07.2023
      12             : ! **************************************************************************************************
      13             : MODULE gw_utils
      14             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15             :                                               get_atomic_kind_set
      16             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      17             :                                               gto_basis_set_type
      18             :    USE bibliography,                    ONLY: Graml2024,&
      19             :                                               cite_reference
      20             :    USE cell_types,                      ONLY: cell_type,&
      21             :                                               pbc
      22             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
      23             :                                               cp_blacs_env_release,&
      24             :                                               cp_blacs_env_type
      25             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      26             :                                               cp_cfm_release,&
      27             :                                               cp_cfm_to_cfm,&
      28             :                                               cp_cfm_to_fm,&
      29             :                                               cp_cfm_type
      30             :    USE cp_control_types,                ONLY: dft_control_type
      31             :    USE cp_dbcsr_api,                    ONLY: &
      32             :         dbcsr_create, dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, &
      33             :         dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
      34             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      35             :                                               copy_fm_to_dbcsr,&
      36             :                                               cp_dbcsr_dist2d_to_dist,&
      37             :                                               dbcsr_allocate_matrix_set,&
      38             :                                               dbcsr_deallocate_matrix_set
      39             :    USE cp_files,                        ONLY: close_file,&
      40             :                                               open_file
      41             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add
      42             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      43             :                                               cp_fm_struct_release,&
      44             :                                               cp_fm_struct_type
      45             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      46             :                                               cp_fm_get_diag,&
      47             :                                               cp_fm_release,&
      48             :                                               cp_fm_set_all,&
      49             :                                               cp_fm_type
      50             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      51             :                                               cp_logger_type
      52             :    USE cp_output_handling,              ONLY: cp_print_key_generate_filename
      53             :    USE dbt_api,                         ONLY: &
      54             :         dbt_clear, dbt_create, dbt_destroy, dbt_filter, dbt_iterator_blocks_left, &
      55             :         dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
      56             :         dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
      57             :    USE distribution_2d_types,           ONLY: distribution_2d_type
      58             :    USE gw_communication,                ONLY: fm_to_local_array
      59             :    USE gw_integrals,                    ONLY: build_3c_integral_block
      60             :    USE gw_kp_to_real_space_and_back,    ONLY: trafo_rs_to_ikp
      61             :    USE input_constants,                 ONLY: do_potential_truncated,&
      62             :                                               large_cell_Gamma,&
      63             :                                               ri_rpa_g0w0_crossing_newton,&
      64             :                                               rtp_method_bse,&
      65             :                                               small_cell_full_kp,&
      66             :                                               xc_none
      67             :    USE input_section_types,             ONLY: section_vals_get,&
      68             :                                               section_vals_get_subs_vals,&
      69             :                                               section_vals_type,&
      70             :                                               section_vals_val_get,&
      71             :                                               section_vals_val_set
      72             :    USE kinds,                           ONLY: default_string_length,&
      73             :                                               dp,&
      74             :                                               int_8
      75             :    USE kpoint_methods,                  ONLY: kpoint_init_cell_index
      76             :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      77             :                                               kpoint_create,&
      78             :                                               kpoint_type
      79             :    USE libint_2c_3c,                    ONLY: libint_potential_type
      80             :    USE libint_wrapper,                  ONLY: cp_libint_static_cleanup,&
      81             :                                               cp_libint_static_init
      82             :    USE machine,                         ONLY: m_memory,&
      83             :                                               m_walltime
      84             :    USE mathconstants,                   ONLY: gaussi,&
      85             :                                               z_one,&
      86             :                                               z_zero
      87             :    USE mathlib,                         ONLY: diag_complex,&
      88             :                                               gcd
      89             :    USE message_passing,                 ONLY: mp_cart_type,&
      90             :                                               mp_para_env_type
      91             :    USE minimax_exp,                     ONLY: get_exp_minimax_coeff
      92             :    USE minimax_exp_gw,                  ONLY: get_exp_minimax_coeff_gw
      93             :    USE minimax_rpa,                     ONLY: get_rpa_minimax_coeff,&
      94             :                                               get_rpa_minimax_coeff_larger_grid
      95             :    USE mp2_gpw,                         ONLY: create_mat_munu
      96             :    USE mp2_grids,                       ONLY: get_l_sq_wghts_cos_tf_t_to_w,&
      97             :                                               get_l_sq_wghts_cos_tf_w_to_t,&
      98             :                                               get_l_sq_wghts_sin_tf_t_to_w
      99             :    USE mp2_ri_2c,                       ONLY: trunc_coulomb_for_exchange
     100             :    USE parallel_gemm_api,               ONLY: parallel_gemm
     101             :    USE particle_methods,                ONLY: get_particle_set
     102             :    USE particle_types,                  ONLY: particle_type
     103             :    USE physcon,                         ONLY: angstrom,&
     104             :                                               evolt
     105             :    USE post_scf_bandstructure_types,    ONLY: post_scf_bandstructure_type
     106             :    USE post_scf_bandstructure_utils,    ONLY: rsmat_to_kp
     107             :    USE qs_energy_types,                 ONLY: qs_energy_type
     108             :    USE qs_environment_types,            ONLY: get_qs_env,&
     109             :                                               qs_env_part_release,&
     110             :                                               qs_environment_type
     111             :    USE qs_integral_utils,               ONLY: basis_set_list_setup
     112             :    USE qs_interactions,                 ONLY: init_interaction_radii_orb_basis
     113             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
     114             :                                               qs_kind_type
     115             :    USE qs_ks_methods,                   ONLY: qs_ks_build_kohn_sham_matrix
     116             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type,&
     117             :                                               release_neighbor_list_sets
     118             :    USE qs_tensors,                      ONLY: build_2c_integrals,&
     119             :                                               build_2c_neighbor_lists,&
     120             :                                               build_3c_integrals,&
     121             :                                               build_3c_neighbor_lists,&
     122             :                                               get_tensor_occupancy,&
     123             :                                               neighbor_list_3c_destroy
     124             :    USE qs_tensors_types,                ONLY: create_2c_tensor,&
     125             :                                               create_3c_tensor,&
     126             :                                               distribution_3d_create,&
     127             :                                               distribution_3d_type,&
     128             :                                               neighbor_list_3c_type
     129             :    USE rpa_gw,                          ONLY: continuation_pade
     130             : #include "base/base_uses.f90"
     131             : 
     132             :    IMPLICIT NONE
     133             : 
     134             :    PRIVATE
     135             : 
     136             :    PUBLIC :: create_and_init_bs_env_for_gw, de_init_bs_env, get_i_j_atoms, &
     137             :              kpoint_init_cell_index_simple, compute_xkp, time_to_freq, analyt_conti_and_print, &
     138             :              add_R, is_cell_in_index_to_cell, get_V_tr_R, power
     139             : 
     140             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_utils'
     141             : 
     142             : CONTAINS
     143             : 
     144             : ! **************************************************************************************************
     145             : !> \brief ...
     146             : !> \param qs_env ...
     147             : !> \param bs_env ...
     148             : !> \param bs_sec ...
     149             : ! **************************************************************************************************
     150          36 :    SUBROUTINE create_and_init_bs_env_for_gw(qs_env, bs_env, bs_sec)
     151             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     152             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     153             :       TYPE(section_vals_type), POINTER                   :: bs_sec
     154             : 
     155             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'create_and_init_bs_env_for_gw'
     156             : 
     157             :       INTEGER                                            :: handle
     158             : 
     159          36 :       CALL timeset(routineN, handle)
     160             : 
     161          36 :       CALL cite_reference(Graml2024)
     162             : 
     163          36 :       CALL read_gw_input_parameters(bs_env, bs_sec)
     164             : 
     165          36 :       CALL print_header_and_input_parameters(bs_env)
     166             : 
     167          36 :       CALL setup_AO_and_RI_basis_set(qs_env, bs_env)
     168             : 
     169          36 :       CALL get_RI_basis_and_basis_function_indices(qs_env, bs_env)
     170             : 
     171          36 :       CALL set_heuristic_parameters(bs_env, qs_env)
     172             : 
     173          36 :       CALL cp_libint_static_init()
     174             : 
     175          36 :       CALL setup_kpoints_chi_eps_W(bs_env, bs_env%kpoints_chi_eps_W)
     176             : 
     177          36 :       IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
     178           8 :          CALL setup_cells_3c(qs_env, bs_env)
     179             :       END IF
     180             : 
     181          36 :       CALL set_parallelization_parameters(qs_env, bs_env)
     182             : 
     183          36 :       CALL allocate_matrices(qs_env, bs_env)
     184             : 
     185          36 :       CALL compute_V_xc(qs_env, bs_env)
     186             : 
     187          36 :       CALL create_tensors(qs_env, bs_env)
     188             : 
     189          64 :       SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
     190             :       CASE (large_cell_Gamma)
     191             : 
     192          28 :          CALL allocate_GW_eigenvalues(bs_env)
     193             : 
     194          28 :          CALL check_sparsity_3c(qs_env, bs_env)
     195             : 
     196          28 :          CALL set_sparsity_parallelization_parameters(bs_env)
     197             : 
     198          28 :          CALL check_for_restart_files(qs_env, bs_env)
     199             : 
     200             :       CASE (small_cell_full_kp)
     201             : 
     202           8 :          CALL compute_3c_integrals(qs_env, bs_env)
     203             : 
     204           8 :          CALL setup_cells_Delta_R(bs_env)
     205             : 
     206           8 :          CALL setup_parallelization_Delta_R(bs_env)
     207             : 
     208           8 :          CALL allocate_matrices_small_cell_full_kp(qs_env, bs_env)
     209             : 
     210           8 :          CALL trafo_V_xc_R_to_kp(qs_env, bs_env)
     211             : 
     212          44 :          CALL heuristic_RI_regularization(qs_env, bs_env)
     213             : 
     214             :       END SELECT
     215             : 
     216          36 :       CALL setup_time_and_frequency_minimax_grid(bs_env)
     217             : 
     218             :       ! free memory in qs_env; only if one is not calculating the LDOS because
     219             :       ! we need real-space grid operations in pw_env, task_list for the LDOS
     220             :       ! Recommendation in case of memory issues: first perform GW calculation without calculating
     221             :       !                                          LDOS (to safe memor). Then, use GW restart files
     222             :       !                                          in a subsequent calculation to calculate the LDOS
     223             :       ! Marek : TODO - boolean that does not interfere with RTP init but sets this to correct value
     224             :       IF (.NOT. bs_env%do_ldos .AND. .FALSE.) THEN
     225             :          CALL qs_env_part_release(qs_env)
     226             :       END IF
     227             : 
     228          36 :       CALL timestop(handle)
     229             : 
     230          36 :    END SUBROUTINE create_and_init_bs_env_for_gw
     231             : 
     232             : ! **************************************************************************************************
     233             : !> \brief ...
     234             : !> \param bs_env ...
     235             : ! **************************************************************************************************
     236          36 :    SUBROUTINE de_init_bs_env(bs_env)
     237             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     238             : 
     239             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'de_init_bs_env'
     240             : 
     241             :       INTEGER                                            :: handle
     242             : 
     243          36 :       CALL timeset(routineN, handle)
     244             :       ! deallocate quantities here which:
     245             :       ! 1. cannot be deallocated in bs_env_release due to circular dependencies
     246             :       ! 2. consume a lot of memory and should not be kept until the quantity is
     247             :       !    deallocated in bs_env_release
     248             : 
     249          36 :       IF (ASSOCIATED(bs_env%nl_3c%ij_list) .AND. (bs_env%rtp_method == rtp_method_bse)) THEN
     250          12 :          IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, *) "Retaining nl_3c for RTBSE"
     251             :       ELSE
     252          24 :          CALL neighbor_list_3c_destroy(bs_env%nl_3c)
     253             :       END IF
     254             : 
     255          36 :       CALL cp_libint_static_cleanup()
     256             : 
     257          36 :       CALL timestop(handle)
     258             : 
     259          36 :    END SUBROUTINE de_init_bs_env
     260             : 
     261             : ! **************************************************************************************************
     262             : !> \brief ...
     263             : !> \param bs_env ...
     264             : !> \param bs_sec ...
     265             : ! **************************************************************************************************
     266          36 :    SUBROUTINE read_gw_input_parameters(bs_env, bs_sec)
     267             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     268             :       TYPE(section_vals_type), POINTER                   :: bs_sec
     269             : 
     270             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_gw_input_parameters'
     271             : 
     272             :       INTEGER                                            :: handle
     273             :       TYPE(section_vals_type), POINTER                   :: gw_sec
     274             : 
     275          36 :       CALL timeset(routineN, handle)
     276             : 
     277          36 :       NULLIFY (gw_sec)
     278          36 :       gw_sec => section_vals_get_subs_vals(bs_sec, "GW")
     279             : 
     280          36 :       CALL section_vals_val_get(gw_sec, "NUM_TIME_FREQ_POINTS", i_val=bs_env%num_time_freq_points)
     281          36 :       CALL section_vals_val_get(gw_sec, "EPS_FILTER", r_val=bs_env%eps_filter)
     282          36 :       CALL section_vals_val_get(gw_sec, "REGULARIZATION_RI", r_val=bs_env%input_regularization_RI)
     283          36 :       CALL section_vals_val_get(gw_sec, "CUTOFF_RADIUS_RI", r_val=bs_env%ri_metric%cutoff_radius)
     284          36 :       CALL section_vals_val_get(gw_sec, "MEMORY_PER_PROC", r_val=bs_env%input_memory_per_proc_GB)
     285          36 :       CALL section_vals_val_get(gw_sec, "APPROX_KP_EXTRAPOL", l_val=bs_env%approx_kp_extrapol)
     286          36 :       CALL section_vals_val_get(gw_sec, "SIZE_LATTICE_SUM", i_val=bs_env%size_lattice_sum_V)
     287          36 :       CALL section_vals_val_get(gw_sec, "KPOINTS_W", i_vals=bs_env%nkp_grid_chi_eps_W_input)
     288          36 :       CALL section_vals_val_get(gw_sec, "HEDIN_SHIFT", l_val=bs_env%do_hedin_shift)
     289          36 :       CALL section_vals_val_get(gw_sec, "FREQ_MAX_FIT", r_val=bs_env%freq_max_fit)
     290             : 
     291          36 :       CALL timestop(handle)
     292             : 
     293          36 :    END SUBROUTINE read_gw_input_parameters
     294             : 
     295             : ! **************************************************************************************************
     296             : !> \brief ...
     297             : !> \param qs_env ...
     298             : !> \param bs_env ...
     299             : ! **************************************************************************************************
     300          36 :    SUBROUTINE setup_AO_and_RI_basis_set(qs_env, bs_env)
     301             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     302             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     303             : 
     304             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_AO_and_RI_basis_set'
     305             : 
     306             :       INTEGER                                            :: handle, natom, nkind
     307          36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     308          36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     309             : 
     310          36 :       CALL timeset(routineN, handle)
     311             : 
     312             :       CALL get_qs_env(qs_env, &
     313             :                       qs_kind_set=qs_kind_set, &
     314             :                       particle_set=particle_set, &
     315          36 :                       natom=natom, nkind=nkind)
     316             : 
     317             :       ! set up basis
     318         180 :       ALLOCATE (bs_env%sizes_RI(natom), bs_env%sizes_AO(natom))
     319         300 :       ALLOCATE (bs_env%basis_set_RI(nkind), bs_env%basis_set_AO(nkind))
     320             : 
     321          36 :       CALL basis_set_list_setup(bs_env%basis_set_RI, "RI_AUX", qs_kind_set)
     322          36 :       CALL basis_set_list_setup(bs_env%basis_set_AO, "ORB", qs_kind_set)
     323             : 
     324             :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_RI, &
     325          36 :                             basis=bs_env%basis_set_RI)
     326             :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=bs_env%sizes_AO, &
     327          36 :                             basis=bs_env%basis_set_AO)
     328             : 
     329          36 :       CALL timestop(handle)
     330             : 
     331          36 :    END SUBROUTINE setup_AO_and_RI_basis_set
     332             : 
     333             : ! **************************************************************************************************
     334             : !> \brief ...
     335             : !> \param qs_env ...
     336             : !> \param bs_env ...
     337             : ! **************************************************************************************************
     338          36 :    SUBROUTINE get_RI_basis_and_basis_function_indices(qs_env, bs_env)
     339             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     340             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     341             : 
     342             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_RI_basis_and_basis_function_indices'
     343             : 
     344             :       INTEGER                                            :: handle, i_RI, iatom, ikind, iset, &
     345             :                                                             max_AO_bf_per_atom, n_ao_test, n_atom, &
     346             :                                                             n_kind, n_RI, nset, nsgf, u
     347          36 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     348          36 :       INTEGER, DIMENSION(:), POINTER                     :: l_max, l_min, nsgf_set
     349          36 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     350             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a
     351          36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     352             : 
     353          36 :       CALL timeset(routineN, handle)
     354             : 
     355             :       ! determine RI basis set size
     356          36 :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
     357             : 
     358          36 :       n_kind = SIZE(qs_kind_set)
     359          36 :       n_atom = bs_env%n_atom
     360             : 
     361          36 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     362             : 
     363          96 :       DO ikind = 1, n_kind
     364             :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, &
     365          60 :                           basis_type="RI_AUX")
     366          96 :          IF (.NOT. ASSOCIATED(basis_set_a)) THEN
     367             :             CALL cp_abort(__LOCATION__, &
     368           0 :                           "At least one RI_AUX basis set was not explicitly invoked in &KIND-section.")
     369             :          END IF
     370             :       END DO
     371             : 
     372         108 :       ALLOCATE (bs_env%i_RI_start_from_atom(n_atom))
     373          72 :       ALLOCATE (bs_env%i_RI_end_from_atom(n_atom))
     374          72 :       ALLOCATE (bs_env%i_ao_start_from_atom(n_atom))
     375          72 :       ALLOCATE (bs_env%i_ao_end_from_atom(n_atom))
     376             : 
     377          36 :       n_RI = 0
     378         116 :       DO iatom = 1, n_atom
     379          80 :          bs_env%i_RI_start_from_atom(iatom) = n_RI + 1
     380          80 :          ikind = kind_of(iatom)
     381          80 :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
     382          80 :          n_RI = n_RI + nsgf
     383         116 :          bs_env%i_RI_end_from_atom(iatom) = n_RI
     384             :       END DO
     385          36 :       bs_env%n_RI = n_RI
     386             : 
     387          36 :       max_AO_bf_per_atom = 0
     388          36 :       n_ao_test = 0
     389         116 :       DO iatom = 1, n_atom
     390          80 :          bs_env%i_ao_start_from_atom(iatom) = n_ao_test + 1
     391          80 :          ikind = kind_of(iatom)
     392          80 :          CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="ORB")
     393          80 :          n_ao_test = n_ao_test + nsgf
     394          80 :          bs_env%i_ao_end_from_atom(iatom) = n_ao_test
     395         116 :          max_AO_bf_per_atom = MAX(max_AO_bf_per_atom, nsgf)
     396             :       END DO
     397          36 :       CPASSERT(n_ao_test == bs_env%n_ao)
     398          36 :       bs_env%max_AO_bf_per_atom = max_AO_bf_per_atom
     399             : 
     400         108 :       ALLOCATE (bs_env%l_RI(n_RI))
     401          36 :       i_RI = 0
     402         116 :       DO iatom = 1, n_atom
     403          80 :          ikind = kind_of(iatom)
     404             : 
     405          80 :          nset = bs_env%basis_set_RI(ikind)%gto_basis_set%nset
     406          80 :          l_max => bs_env%basis_set_RI(ikind)%gto_basis_set%lmax
     407          80 :          l_min => bs_env%basis_set_RI(ikind)%gto_basis_set%lmin
     408          80 :          nsgf_set => bs_env%basis_set_RI(ikind)%gto_basis_set%nsgf_set
     409             : 
     410         316 :          DO iset = 1, nset
     411         200 :             CPASSERT(l_max(iset) == l_min(iset))
     412         600 :             bs_env%l_RI(i_RI + 1:i_RI + nsgf_set(iset)) = l_max(iset)
     413         280 :             i_RI = i_RI + nsgf_set(iset)
     414             :          END DO
     415             : 
     416             :       END DO
     417          36 :       CPASSERT(i_RI == n_RI)
     418             : 
     419          36 :       u = bs_env%unit_nr
     420             : 
     421          36 :       IF (u > 0) THEN
     422          18 :          WRITE (u, FMT="(T2,A)") " "
     423          18 :          WRITE (u, FMT="(T2,2A,T75,I8)") "Number of auxiliary Gaussian basis functions ", &
     424          36 :             "for χ, ε, W", n_RI
     425             :       END IF
     426             : 
     427          36 :       CALL timestop(handle)
     428             : 
     429          72 :    END SUBROUTINE get_RI_basis_and_basis_function_indices
     430             : 
     431             : ! **************************************************************************************************
     432             : !> \brief ...
     433             : !> \param bs_env ...
     434             : !> \param kpoints ...
     435             : ! **************************************************************************************************
     436          36 :    SUBROUTINE setup_kpoints_chi_eps_W(bs_env, kpoints)
     437             : 
     438             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     439             :       TYPE(kpoint_type), POINTER                         :: kpoints
     440             : 
     441             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_chi_eps_W'
     442             : 
     443             :       INTEGER                                            :: handle, i_dim, n_dim, nkp, nkp_extra, &
     444             :                                                             nkp_orig, u
     445             :       INTEGER, DIMENSION(3)                              :: nkp_grid, nkp_grid_extra, periodic
     446             :       REAL(KIND=dp)                                      :: exp_s_p, n_dim_inv
     447             : 
     448          36 :       CALL timeset(routineN, handle)
     449             : 
     450             :       ! routine adapted from mp2_integrals.F
     451          36 :       NULLIFY (kpoints)
     452          36 :       CALL kpoint_create(kpoints)
     453             : 
     454          36 :       kpoints%kp_scheme = "GENERAL"
     455             : 
     456         144 :       periodic(1:3) = bs_env%periodic(1:3)
     457             : 
     458          36 :       CPASSERT(SIZE(bs_env%nkp_grid_chi_eps_W_input) == 3)
     459             : 
     460             :       IF (bs_env%nkp_grid_chi_eps_W_input(1) > 0 .AND. &
     461          36 :           bs_env%nkp_grid_chi_eps_W_input(2) > 0 .AND. &
     462             :           bs_env%nkp_grid_chi_eps_W_input(3) > 0) THEN
     463             :          ! 1. k-point mesh for χ, ε, W from input
     464           0 :          DO i_dim = 1, 3
     465           0 :             SELECT CASE (periodic(i_dim))
     466             :             CASE (0)
     467           0 :                nkp_grid(i_dim) = 1
     468           0 :                nkp_grid_extra(i_dim) = 1
     469             :             CASE (1)
     470           0 :                nkp_grid(i_dim) = bs_env%nkp_grid_chi_eps_W_input(i_dim)
     471           0 :                nkp_grid_extra(i_dim) = nkp_grid(i_dim)*2
     472             :             CASE DEFAULT
     473           0 :                CPABORT("Error in periodicity.")
     474             :             END SELECT
     475             :          END DO
     476             : 
     477             :       ELSE IF (bs_env%nkp_grid_chi_eps_W_input(1) == -1 .AND. &
     478          36 :                bs_env%nkp_grid_chi_eps_W_input(2) == -1 .AND. &
     479             :                bs_env%nkp_grid_chi_eps_W_input(3) == -1) THEN
     480             :          ! 2. automatic k-point mesh for χ, ε, W
     481             : 
     482         144 :          DO i_dim = 1, 3
     483             : 
     484         108 :             CPASSERT(periodic(i_dim) == 0 .OR. periodic(i_dim) == 1)
     485             : 
     486          36 :             SELECT CASE (periodic(i_dim))
     487             :             CASE (0)
     488          60 :                nkp_grid(i_dim) = 1
     489          60 :                nkp_grid_extra(i_dim) = 1
     490             :             CASE (1)
     491          80 :                SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
     492             :                CASE (large_cell_Gamma)
     493          32 :                   nkp_grid(i_dim) = 4
     494          32 :                   nkp_grid_extra(i_dim) = 6
     495             :                CASE (small_cell_full_kp)
     496          16 :                   nkp_grid(i_dim) = bs_env%kpoints_scf_desymm%nkp_grid(i_dim)*4
     497          48 :                   nkp_grid_extra(i_dim) = bs_env%kpoints_scf_desymm%nkp_grid(i_dim)*8
     498             :                END SELECT
     499             :             CASE DEFAULT
     500         108 :                CPABORT("Error in periodicity.")
     501             :             END SELECT
     502             : 
     503             :          END DO
     504             : 
     505             :       ELSE
     506             : 
     507           0 :          CPABORT("An error occured when setting up the k-mesh for W.")
     508             : 
     509             :       END IF
     510             : 
     511          36 :       nkp_orig = MAX(nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2, 1)
     512             : 
     513          36 :       nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
     514             : 
     515          36 :       nkp = nkp_orig + nkp_extra
     516             : 
     517         144 :       kpoints%nkp_grid(1:3) = nkp_grid(1:3)
     518          36 :       kpoints%nkp = nkp
     519             : 
     520         144 :       bs_env%nkp_grid_chi_eps_W_orig(1:3) = nkp_grid(1:3)
     521         144 :       bs_env%nkp_grid_chi_eps_W_extra(1:3) = nkp_grid_extra(1:3)
     522          36 :       bs_env%nkp_chi_eps_W_orig = nkp_orig
     523          36 :       bs_env%nkp_chi_eps_W_extra = nkp_extra
     524          36 :       bs_env%nkp_chi_eps_W_orig_plus_extra = nkp
     525             : 
     526         180 :       ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
     527         180 :       ALLOCATE (bs_env%wkp_no_extra(nkp), bs_env%wkp_s_p(nkp))
     528             : 
     529          36 :       CALL compute_xkp(kpoints%xkp, 1, nkp_orig, nkp_grid)
     530          36 :       CALL compute_xkp(kpoints%xkp, nkp_orig + 1, nkp, nkp_grid_extra)
     531             : 
     532         144 :       n_dim = SUM(periodic)
     533          36 :       IF (n_dim == 0) THEN
     534             :          ! molecules
     535          12 :          kpoints%wkp(1) = 1.0_dp
     536          12 :          bs_env%wkp_s_p(1) = 1.0_dp
     537          12 :          bs_env%wkp_no_extra(1) = 1.0_dp
     538             :       ELSE
     539             : 
     540          24 :          n_dim_inv = 1.0_dp/REAL(n_dim, KIND=dp)
     541             : 
     542             :          ! k-point weights are chosen to automatically extrapolate the k-point mesh
     543          24 :          CALL compute_wkp(kpoints%wkp(1:nkp_orig), nkp_orig, nkp_extra, n_dim_inv)
     544          24 :          CALL compute_wkp(kpoints%wkp(nkp_orig + 1:nkp), nkp_extra, nkp_orig, n_dim_inv)
     545             : 
     546        1176 :          bs_env%wkp_no_extra(1:nkp_orig) = 0.0_dp
     547        4408 :          bs_env%wkp_no_extra(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp)
     548             : 
     549          24 :          IF (n_dim == 3) THEN
     550             :             ! W_PQ(k) for an s-function P and a p-function Q diverges as 1/k at k=0
     551             :             ! (instead of 1/k^2 for P and Q both being s-functions).
     552           0 :             exp_s_p = 2.0_dp*n_dim_inv
     553           0 :             CALL compute_wkp(bs_env%wkp_s_p(1:nkp_orig), nkp_orig, nkp_extra, exp_s_p)
     554           0 :             CALL compute_wkp(bs_env%wkp_s_p(nkp_orig + 1:nkp), nkp_extra, nkp_orig, exp_s_p)
     555             :          ELSE
     556        5560 :             bs_env%wkp_s_p(1:nkp) = bs_env%wkp_no_extra(1:nkp)
     557             :          END IF
     558             : 
     559             :       END IF
     560             : 
     561          36 :       IF (bs_env%approx_kp_extrapol) THEN
     562           2 :          bs_env%wkp_orig = 1.0_dp/REAL(nkp_orig, KIND=dp)
     563             :       END IF
     564             : 
     565             :       ! heuristic parameter: how many k-points for χ, ε, and W are used simultaneously
     566             :       ! (less simultaneous k-points: less memory, but more computational effort because of
     567             :       !  recomputation of V(k))
     568          36 :       bs_env%nkp_chi_eps_W_batch = 4
     569             : 
     570             :       bs_env%num_chi_eps_W_batches = (bs_env%nkp_chi_eps_W_orig_plus_extra - 1)/ &
     571          36 :                                      bs_env%nkp_chi_eps_W_batch + 1
     572             : 
     573          36 :       u = bs_env%unit_nr
     574             : 
     575          36 :       IF (u > 0) THEN
     576          18 :          WRITE (u, FMT="(T2,A)") " "
     577          18 :          WRITE (u, FMT="(T2,1A,T71,3I4)") "K-point mesh 1 for χ, ε, W", nkp_grid(1:3)
     578          18 :          WRITE (u, FMT="(T2,2A,T71,3I4)") "K-point mesh 2 for χ, ε, W ", &
     579          36 :             "(for k-point extrapolation of W)", nkp_grid_extra(1:3)
     580          18 :          WRITE (u, FMT="(T2,A,T80,L)") "Approximate the k-point extrapolation", &
     581          36 :             bs_env%approx_kp_extrapol
     582             :       END IF
     583             : 
     584          36 :       CALL timestop(handle)
     585             : 
     586          36 :    END SUBROUTINE setup_kpoints_chi_eps_W
     587             : 
     588             : ! **************************************************************************************************
     589             : !> \brief ...
     590             : !> \param kpoints ...
     591             : !> \param qs_env ...
     592             : ! **************************************************************************************************
     593           0 :    SUBROUTINE kpoint_init_cell_index_simple(kpoints, qs_env)
     594             : 
     595             :       TYPE(kpoint_type), POINTER                         :: kpoints
     596             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     597             : 
     598             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'kpoint_init_cell_index_simple'
     599             : 
     600             :       INTEGER                                            :: handle
     601             :       TYPE(dft_control_type), POINTER                    :: dft_control
     602             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     603             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     604           0 :          POINTER                                         :: sab_orb
     605             : 
     606           0 :       CALL timeset(routineN, handle)
     607             : 
     608           0 :       NULLIFY (dft_control, para_env, sab_orb)
     609           0 :       CALL get_qs_env(qs_env=qs_env, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb)
     610           0 :       CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control)
     611             : 
     612           0 :       CALL timestop(handle)
     613             : 
     614           0 :    END SUBROUTINE kpoint_init_cell_index_simple
     615             : 
     616             : ! **************************************************************************************************
     617             : !> \brief ...
     618             : !> \param xkp ...
     619             : !> \param ikp_start ...
     620             : !> \param ikp_end ...
     621             : !> \param grid ...
     622             : ! **************************************************************************************************
     623          72 :    SUBROUTINE compute_xkp(xkp, ikp_start, ikp_end, grid)
     624             : 
     625             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: xkp
     626             :       INTEGER                                            :: ikp_start, ikp_end
     627             :       INTEGER, DIMENSION(3)                              :: grid
     628             : 
     629             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_xkp'
     630             : 
     631             :       INTEGER                                            :: handle, i, ix, iy, iz
     632             : 
     633          72 :       CALL timeset(routineN, handle)
     634             : 
     635          72 :       i = ikp_start
     636         392 :       DO ix = 1, grid(1)
     637        6224 :          DO iy = 1, grid(2)
     638       17248 :             DO iz = 1, grid(3)
     639             : 
     640       11096 :                IF (i > ikp_end) CYCLE
     641             : 
     642        5548 :                xkp(1, i) = REAL(2*ix - grid(1) - 1, KIND=dp)/(2._dp*REAL(grid(1), KIND=dp))
     643        5548 :                xkp(2, i) = REAL(2*iy - grid(2) - 1, KIND=dp)/(2._dp*REAL(grid(2), KIND=dp))
     644        5548 :                xkp(3, i) = REAL(2*iz - grid(3) - 1, KIND=dp)/(2._dp*REAL(grid(3), KIND=dp))
     645       16928 :                i = i + 1
     646             : 
     647             :             END DO
     648             :          END DO
     649             :       END DO
     650             : 
     651          72 :       CALL timestop(handle)
     652             : 
     653          72 :    END SUBROUTINE compute_xkp
     654             : 
     655             : ! **************************************************************************************************
     656             : !> \brief ...
     657             : !> \param wkp ...
     658             : !> \param nkp_1 ...
     659             : !> \param nkp_2 ...
     660             : !> \param exponent ...
     661             : ! **************************************************************************************************
     662          48 :    SUBROUTINE compute_wkp(wkp, nkp_1, nkp_2, exponent)
     663             :       REAL(KIND=dp), DIMENSION(:)                        :: wkp
     664             :       INTEGER                                            :: nkp_1, nkp_2
     665             :       REAL(KIND=dp)                                      :: exponent
     666             : 
     667             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_wkp'
     668             : 
     669             :       INTEGER                                            :: handle
     670             :       REAL(KIND=dp)                                      :: nkp_ratio
     671             : 
     672          48 :       CALL timeset(routineN, handle)
     673             : 
     674          48 :       nkp_ratio = REAL(nkp_2, KIND=dp)/REAL(nkp_1, KIND=dp)
     675             : 
     676        5584 :       wkp(:) = 1.0_dp/REAL(nkp_1, KIND=dp)/(1.0_dp - nkp_ratio**exponent)
     677             : 
     678          48 :       CALL timestop(handle)
     679             : 
     680          48 :    END SUBROUTINE compute_wkp
     681             : 
     682             : ! **************************************************************************************************
     683             : !> \brief ...
     684             : !> \param qs_env ...
     685             : !> \param bs_env ...
     686             : ! **************************************************************************************************
     687          36 :    SUBROUTINE allocate_matrices(qs_env, bs_env)
     688             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     689             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     690             : 
     691             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'allocate_matrices'
     692             : 
     693             :       INTEGER                                            :: handle, i_t
     694             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env, blacs_env_tensor
     695             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_RI_global
     696             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     697             : 
     698          36 :       CALL timeset(routineN, handle)
     699             : 
     700          36 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
     701             : 
     702          36 :       fm_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
     703             : 
     704          36 :       CALL cp_fm_create(bs_env%fm_Gocc, fm_struct)
     705          36 :       CALL cp_fm_create(bs_env%fm_Gvir, fm_struct)
     706             : 
     707          36 :       NULLIFY (fm_struct_RI_global)
     708             :       CALL cp_fm_struct_create(fm_struct_RI_global, context=blacs_env, nrow_global=bs_env%n_RI, &
     709          36 :                                ncol_global=bs_env%n_RI, para_env=para_env)
     710          36 :       CALL cp_fm_create(bs_env%fm_RI_RI, fm_struct_RI_global)
     711          36 :       CALL cp_fm_create(bs_env%fm_chi_Gamma_freq, fm_struct_RI_global)
     712          36 :       CALL cp_fm_create(bs_env%fm_W_MIC_freq, fm_struct_RI_global)
     713          36 :       IF (bs_env%approx_kp_extrapol) THEN
     714           2 :          CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_extra, fm_struct_RI_global)
     715           2 :          CALL cp_fm_create(bs_env%fm_W_MIC_freq_1_no_extra, fm_struct_RI_global)
     716           2 :          CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_extra, 0.0_dp)
     717           2 :          CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_1_no_extra, 0.0_dp)
     718             :       END IF
     719          36 :       CALL cp_fm_struct_release(fm_struct_RI_global)
     720             : 
     721             :       ! create blacs_env for subgroups of tensor operations
     722          36 :       NULLIFY (blacs_env_tensor)
     723          36 :       CALL cp_blacs_env_create(blacs_env=blacs_env_tensor, para_env=bs_env%para_env_tensor)
     724             : 
     725             :       ! allocate dbcsr matrices in the tensor subgroup; actually, one only needs a small
     726             :       ! subset of blocks in the tensor subgroup, however, all atomic blocks are allocated.
     727             :       ! One might think of creating a dbcsr matrix with only the blocks that are needed
     728             :       ! in the tensor subgroup
     729             :       CALL create_mat_munu(bs_env%mat_ao_ao_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
     730          36 :                            blacs_env_tensor, do_ri_aux_basis=.FALSE.)
     731             : 
     732             :       CALL create_mat_munu(bs_env%mat_RI_RI_tensor, qs_env, bs_env%eps_atom_grid_2d_mat, &
     733          36 :                            blacs_env_tensor, do_ri_aux_basis=.TRUE.)
     734             : 
     735             :       CALL create_mat_munu(bs_env%mat_RI_RI, qs_env, bs_env%eps_atom_grid_2d_mat, &
     736          36 :                            blacs_env, do_ri_aux_basis=.TRUE.)
     737             : 
     738          36 :       CALL cp_blacs_env_release(blacs_env_tensor)
     739             : 
     740          36 :       NULLIFY (bs_env%mat_chi_Gamma_tau)
     741          36 :       CALL dbcsr_allocate_matrix_set(bs_env%mat_chi_Gamma_tau, bs_env%num_time_freq_points)
     742             : 
     743         484 :       DO i_t = 1, bs_env%num_time_freq_points
     744         448 :          ALLOCATE (bs_env%mat_chi_Gamma_tau(i_t)%matrix)
     745         484 :          CALL dbcsr_create(bs_env%mat_chi_Gamma_tau(i_t)%matrix, template=bs_env%mat_RI_RI%matrix)
     746             :       END DO
     747             : 
     748          36 :       CALL timestop(handle)
     749             : 
     750          36 :    END SUBROUTINE allocate_matrices
     751             : 
     752             : ! **************************************************************************************************
     753             : !> \brief ...
     754             : !> \param bs_env ...
     755             : ! **************************************************************************************************
     756          28 :    SUBROUTINE allocate_GW_eigenvalues(bs_env)
     757             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     758             : 
     759             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_GW_eigenvalues'
     760             : 
     761             :       INTEGER                                            :: handle
     762             : 
     763          28 :       CALL timeset(routineN, handle)
     764             : 
     765         140 :       ALLOCATE (bs_env%eigenval_G0W0(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
     766         140 :       ALLOCATE (bs_env%eigenval_HF(bs_env%n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
     767             : 
     768          28 :       CALL timestop(handle)
     769             : 
     770          28 :    END SUBROUTINE allocate_GW_eigenvalues
     771             : 
     772             : ! **************************************************************************************************
     773             : !> \brief ...
     774             : !> \param qs_env ...
     775             : !> \param bs_env ...
     776             : ! **************************************************************************************************
     777          36 :    SUBROUTINE create_tensors(qs_env, bs_env)
     778             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     779             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     780             : 
     781             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'create_tensors'
     782             : 
     783             :       INTEGER                                            :: handle
     784             : 
     785          36 :       CALL timeset(routineN, handle)
     786             : 
     787          36 :       CALL init_interaction_radii(bs_env)
     788             : 
     789             :       ! split blocks does not improve load balancing/efficienfy for tensor contraction, so we go
     790             :       ! with the standard atomic blocks
     791             :       CALL create_3c_t(bs_env%t_RI_AO__AO, bs_env%para_env_tensor, "(RI AO | AO)", [1, 2], [3], &
     792             :                        bs_env%sizes_RI, bs_env%sizes_AO, &
     793          36 :                        create_nl_3c=.TRUE., nl_3c=bs_env%nl_3c, qs_env=qs_env)
     794             :       CALL create_3c_t(bs_env%t_RI__AO_AO, bs_env%para_env_tensor, "(RI | AO AO)", [1], [2, 3], &
     795          36 :                        bs_env%sizes_RI, bs_env%sizes_AO)
     796             : 
     797          36 :       CALL create_2c_t(bs_env, bs_env%sizes_RI, bs_env%sizes_AO)
     798             : 
     799          36 :       CALL timestop(handle)
     800             : 
     801          36 :    END SUBROUTINE create_tensors
     802             : 
     803             : ! **************************************************************************************************
     804             : !> \brief ...
     805             : !> \param qs_env ...
     806             : !> \param bs_env ...
     807             : ! **************************************************************************************************
     808          28 :    SUBROUTINE check_sparsity_3c(qs_env, bs_env)
     809             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     810             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     811             : 
     812             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'check_sparsity_3c'
     813             : 
     814             :       INTEGER                                            :: handle, n_atom_step, RI_atom
     815             :       INTEGER(int_8)                                     :: mem, non_zero_elements_sum, nze
     816             :       REAL(dp)                                           :: max_dist_AO_atoms, occ, occupation_sum
     817             :       REAL(KIND=dp)                                      :: t1, t2
     818         196 :       TYPE(dbt_type)                                     :: t_3c_global
     819          28 :       TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :)       :: t_3c_global_array
     820             :       TYPE(neighbor_list_3c_type)                        :: nl_3c_global
     821             : 
     822          28 :       CALL timeset(routineN, handle)
     823             : 
     824             :       ! check the sparsity of 3c integral tensor (µν|P); calculate maximum distance between
     825             :       ! AO atoms µ, ν where at least a single integral (µν|P) is larger than the filter threshold
     826             :       CALL create_3c_t(t_3c_global, bs_env%para_env, "(RI AO | AO)", [1, 2], [3], &
     827             :                        bs_env%sizes_RI, bs_env%sizes_AO, &
     828          28 :                        create_nl_3c=.TRUE., nl_3c=nl_3c_global, qs_env=qs_env)
     829             : 
     830          28 :       CALL m_memory(mem)
     831          28 :       CALL bs_env%para_env%max(mem)
     832             : 
     833         252 :       ALLOCATE (t_3c_global_array(1, 1))
     834          28 :       CALL dbt_create(t_3c_global, t_3c_global_array(1, 1))
     835             : 
     836          28 :       CALL bs_env%para_env%sync()
     837          28 :       t1 = m_walltime()
     838             : 
     839          28 :       occupation_sum = 0.0_dp
     840          28 :       non_zero_elements_sum = 0
     841          28 :       max_dist_AO_atoms = 0.0_dp
     842          28 :       n_atom_step = INT(SQRT(REAL(bs_env%n_atom, KIND=dp)))
     843             :       ! do not compute full 3c integrals at once because it may cause out of memory
     844          88 :       DO RI_atom = 1, bs_env%n_atom, n_atom_step
     845             : 
     846             :          CALL build_3c_integrals(t_3c_global_array, &
     847             :                                  bs_env%eps_filter, &
     848             :                                  qs_env, &
     849             :                                  nl_3c_global, &
     850             :                                  int_eps=bs_env%eps_filter, &
     851             :                                  basis_i=bs_env%basis_set_RI, &
     852             :                                  basis_j=bs_env%basis_set_AO, &
     853             :                                  basis_k=bs_env%basis_set_AO, &
     854             :                                  bounds_i=[RI_atom, MIN(RI_atom + n_atom_step - 1, bs_env%n_atom)], &
     855             :                                  potential_parameter=bs_env%ri_metric, &
     856         180 :                                  desymmetrize=.FALSE.)
     857             : 
     858          60 :          CALL dbt_filter(t_3c_global_array(1, 1), bs_env%eps_filter)
     859             : 
     860          60 :          CALL bs_env%para_env%sync()
     861             : 
     862          60 :          CALL get_tensor_occupancy(t_3c_global_array(1, 1), nze, occ)
     863          60 :          non_zero_elements_sum = non_zero_elements_sum + nze
     864          60 :          occupation_sum = occupation_sum + occ
     865             : 
     866          60 :          CALL get_max_dist_AO_atoms(t_3c_global_array(1, 1), max_dist_AO_atoms, qs_env)
     867             : 
     868         148 :          CALL dbt_clear(t_3c_global_array(1, 1))
     869             : 
     870             :       END DO
     871             : 
     872          28 :       t2 = m_walltime()
     873             : 
     874          28 :       bs_env%occupation_3c_int = occupation_sum
     875          28 :       bs_env%max_dist_AO_atoms = max_dist_AO_atoms
     876             : 
     877          28 :       CALL dbt_destroy(t_3c_global)
     878          28 :       CALL dbt_destroy(t_3c_global_array(1, 1))
     879          56 :       DEALLOCATE (t_3c_global_array)
     880             : 
     881          28 :       CALL neighbor_list_3c_destroy(nl_3c_global)
     882             : 
     883          28 :       IF (bs_env%unit_nr > 0) THEN
     884          14 :          WRITE (bs_env%unit_nr, '(T2,A)') ''
     885             :          WRITE (bs_env%unit_nr, '(T2,A,F27.1,A)') &
     886          14 :             'Computed 3-center integrals (µν|P), execution time', t2 - t1, ' s'
     887          14 :          WRITE (bs_env%unit_nr, '(T2,A,F48.3,A)') 'Percentage of non-zero (µν|P)', &
     888          28 :             occupation_sum*100, ' %'
     889          14 :          WRITE (bs_env%unit_nr, '(T2,A,F33.1,A)') 'Max. distance between µ,ν in non-zero (µν|P)', &
     890          28 :             max_dist_AO_atoms*angstrom, ' A'
     891          14 :          WRITE (bs_env%unit_nr, '(T2,2A,I20,A)') 'Required memory if storing all 3-center ', &
     892          28 :             'integrals (µν|P)', INT(REAL(non_zero_elements_sum, KIND=dp)*8.0E-9_dp), ' GB'
     893             :       END IF
     894             : 
     895          28 :       CALL timestop(handle)
     896             : 
     897         112 :    END SUBROUTINE check_sparsity_3c
     898             : 
     899             : ! **************************************************************************************************
     900             : !> \brief ...
     901             : !> \param bs_env ...
     902             : !> \param sizes_RI ...
     903             : !> \param sizes_AO ...
     904             : ! **************************************************************************************************
     905          36 :    SUBROUTINE create_2c_t(bs_env, sizes_RI, sizes_AO)
     906             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
     907             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sizes_RI, sizes_AO
     908             : 
     909             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'create_2c_t'
     910             : 
     911             :       INTEGER                                            :: handle
     912          36 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: dist_1, dist_2
     913             :       INTEGER, DIMENSION(2)                              :: pdims_2d
     914         108 :       TYPE(dbt_pgrid_type)                               :: pgrid_2d
     915             : 
     916          36 :       CALL timeset(routineN, handle)
     917             : 
     918             :       ! inspired from rpa_im_time.F / hfx_types.F
     919             : 
     920          36 :       pdims_2d = 0
     921          36 :       CALL dbt_pgrid_create(bs_env%para_env_tensor, pdims_2d, pgrid_2d)
     922             : 
     923             :       CALL create_2c_tensor(bs_env%t_G, dist_1, dist_2, pgrid_2d, sizes_AO, sizes_AO, &
     924          36 :                             name="(AO | AO)")
     925          36 :       DEALLOCATE (dist_1, dist_2)
     926             :       CALL create_2c_tensor(bs_env%t_chi, dist_1, dist_2, pgrid_2d, sizes_RI, sizes_RI, &
     927          36 :                             name="(RI | RI)")
     928          36 :       DEALLOCATE (dist_1, dist_2)
     929             :       CALL create_2c_tensor(bs_env%t_W, dist_1, dist_2, pgrid_2d, sizes_RI, sizes_RI, &
     930          36 :                             name="(RI | RI)")
     931          36 :       DEALLOCATE (dist_1, dist_2)
     932          36 :       CALL dbt_pgrid_destroy(pgrid_2d)
     933             : 
     934          36 :       CALL timestop(handle)
     935             : 
     936          36 :    END SUBROUTINE create_2c_t
     937             : 
     938             : ! **************************************************************************************************
     939             : !> \brief ...
     940             : !> \param tensor ...
     941             : !> \param para_env ...
     942             : !> \param tensor_name ...
     943             : !> \param map1 ...
     944             : !> \param map2 ...
     945             : !> \param sizes_RI ...
     946             : !> \param sizes_AO ...
     947             : !> \param create_nl_3c ...
     948             : !> \param nl_3c ...
     949             : !> \param qs_env ...
     950             : ! **************************************************************************************************
     951         100 :    SUBROUTINE create_3c_t(tensor, para_env, tensor_name, map1, map2, sizes_RI, sizes_AO, &
     952             :                           create_nl_3c, nl_3c, qs_env)
     953             :       TYPE(dbt_type)                                     :: tensor
     954             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     955             :       CHARACTER(LEN=12)                                  :: tensor_name
     956             :       INTEGER, DIMENSION(:)                              :: map1, map2
     957             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sizes_RI, sizes_AO
     958             :       LOGICAL, OPTIONAL                                  :: create_nl_3c
     959             :       TYPE(neighbor_list_3c_type), OPTIONAL              :: nl_3c
     960             :       TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
     961             : 
     962             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'create_3c_t'
     963             : 
     964             :       INTEGER                                            :: handle, nkind
     965         100 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: dist_AO_1, dist_AO_2, dist_RI
     966             :       INTEGER, DIMENSION(3)                              :: pcoord, pdims, pdims_3d
     967             :       LOGICAL                                            :: my_create_nl_3c
     968         300 :       TYPE(dbt_pgrid_type)                               :: pgrid_3d
     969             :       TYPE(distribution_3d_type)                         :: dist_3d
     970         100 :       TYPE(mp_cart_type)                                 :: mp_comm_t3c_2
     971         100 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     972             : 
     973         100 :       CALL timeset(routineN, handle)
     974             : 
     975         100 :       pdims_3d = 0
     976         100 :       CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d)
     977             :       CALL create_3c_tensor(tensor, dist_RI, dist_AO_1, dist_AO_2, &
     978             :                             pgrid_3d, sizes_RI, sizes_AO, sizes_AO, &
     979         100 :                             map1=map1, map2=map2, name=tensor_name)
     980             : 
     981         100 :       IF (PRESENT(create_nl_3c)) THEN
     982          64 :          my_create_nl_3c = create_nl_3c
     983             :       ELSE
     984             :          my_create_nl_3c = .FALSE.
     985             :       END IF
     986             : 
     987          64 :       IF (my_create_nl_3c) THEN
     988          64 :          CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
     989          64 :          CALL dbt_mp_environ_pgrid(pgrid_3d, pdims, pcoord)
     990          64 :          CALL mp_comm_t3c_2%create(pgrid_3d%mp_comm_2d, 3, pdims)
     991             :          CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, &
     992          64 :                                      nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.)
     993             : 
     994             :          CALL build_3c_neighbor_lists(nl_3c, &
     995             :                                       qs_env%bs_env%basis_set_RI, &
     996             :                                       qs_env%bs_env%basis_set_AO, &
     997             :                                       qs_env%bs_env%basis_set_AO, &
     998             :                                       dist_3d, qs_env%bs_env%ri_metric, &
     999          64 :                                       "GW_3c_nl", qs_env, own_dist=.TRUE.)
    1000             :       END IF
    1001             : 
    1002         100 :       DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
    1003         100 :       CALL dbt_pgrid_destroy(pgrid_3d)
    1004             : 
    1005         100 :       CALL timestop(handle)
    1006             : 
    1007         200 :    END SUBROUTINE create_3c_t
    1008             : 
    1009             : ! **************************************************************************************************
    1010             : !> \brief ...
    1011             : !> \param bs_env ...
    1012             : ! **************************************************************************************************
    1013          36 :    SUBROUTINE init_interaction_radii(bs_env)
    1014             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1015             : 
    1016             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'init_interaction_radii'
    1017             : 
    1018             :       INTEGER                                            :: handle, ibasis
    1019             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis, ri_basis
    1020             : 
    1021          36 :       CALL timeset(routineN, handle)
    1022             : 
    1023          96 :       DO ibasis = 1, SIZE(bs_env%basis_set_AO)
    1024             : 
    1025          60 :          orb_basis => bs_env%basis_set_AO(ibasis)%gto_basis_set
    1026          60 :          CALL init_interaction_radii_orb_basis(orb_basis, bs_env%eps_filter)
    1027             : 
    1028          60 :          ri_basis => bs_env%basis_set_RI(ibasis)%gto_basis_set
    1029          96 :          CALL init_interaction_radii_orb_basis(ri_basis, bs_env%eps_filter)
    1030             : 
    1031             :       END DO
    1032             : 
    1033          36 :       CALL timestop(handle)
    1034             : 
    1035          36 :    END SUBROUTINE init_interaction_radii
    1036             : 
    1037             : ! **************************************************************************************************
    1038             : !> \brief ...
    1039             : !> \param t_3c_int ...
    1040             : !> \param max_dist_AO_atoms ...
    1041             : !> \param qs_env ...
    1042             : ! **************************************************************************************************
    1043          60 :    SUBROUTINE get_max_dist_AO_atoms(t_3c_int, max_dist_AO_atoms, qs_env)
    1044             :       TYPE(dbt_type)                                     :: t_3c_int
    1045             :       REAL(KIND=dp)                                      :: max_dist_AO_atoms
    1046             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1047             : 
    1048             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_max_dist_AO_atoms'
    1049             : 
    1050             :       INTEGER                                            :: atom_1, atom_2, handle, num_cells
    1051             :       INTEGER, DIMENSION(3)                              :: atom_ind
    1052          60 :       INTEGER, DIMENSION(:, :), POINTER                  :: index_to_cell
    1053             :       REAL(KIND=dp)                                      :: abs_rab
    1054             :       REAL(KIND=dp), DIMENSION(3)                        :: rab
    1055             :       TYPE(cell_type), POINTER                           :: cell
    1056             :       TYPE(dbt_iterator_type)                            :: iter
    1057             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1058          60 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1059             : 
    1060          60 :       CALL timeset(routineN, handle)
    1061             : 
    1062          60 :       NULLIFY (cell, particle_set, para_env)
    1063          60 :       CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, para_env=para_env)
    1064             : 
    1065             : !$OMP PARALLEL DEFAULT(NONE) &
    1066             : !$OMP SHARED(t_3c_int, max_dist_AO_atoms, num_cells, index_to_cell, particle_set, cell) &
    1067          60 : !$OMP PRIVATE(iter,atom_ind,rab, abs_rab, atom_1, atom_2)
    1068             :       CALL dbt_iterator_start(iter, t_3c_int)
    1069             :       DO WHILE (dbt_iterator_blocks_left(iter))
    1070             :          CALL dbt_iterator_next_block(iter, atom_ind)
    1071             : 
    1072             :          atom_1 = atom_ind(2)
    1073             :          atom_2 = atom_ind(3)
    1074             : 
    1075             :          rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
    1076             : 
    1077             :          abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
    1078             : 
    1079             :          max_dist_AO_atoms = MAX(max_dist_AO_atoms, abs_rab)
    1080             : 
    1081             :       END DO
    1082             :       CALL dbt_iterator_stop(iter)
    1083             : !$OMP END PARALLEL
    1084             : 
    1085          60 :       CALL para_env%max(max_dist_AO_atoms)
    1086             : 
    1087          60 :       CALL timestop(handle)
    1088             : 
    1089          60 :    END SUBROUTINE get_max_dist_AO_atoms
    1090             : 
    1091             : ! **************************************************************************************************
    1092             : !> \brief ...
    1093             : !> \param bs_env ...
    1094             : ! **************************************************************************************************
    1095          28 :    SUBROUTINE set_sparsity_parallelization_parameters(bs_env)
    1096             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1097             : 
    1098             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'set_sparsity_parallelization_parameters'
    1099             : 
    1100             :       INTEGER :: handle, i_ivl, IL_ivl, j_ivl, n_atom_per_IL_ivl, n_atom_per_ivl, n_intervals_i, &
    1101             :          n_intervals_inner_loop_atoms, n_intervals_j, u
    1102             :       INTEGER(KIND=int_8)                                :: input_memory_per_proc
    1103             : 
    1104          28 :       CALL timeset(routineN, handle)
    1105             : 
    1106             :       ! heuristic parameter to prevent out of memory
    1107          28 :       bs_env%safety_factor_memory = 0.10_dp
    1108             : 
    1109          28 :       input_memory_per_proc = INT(bs_env%input_memory_per_proc_GB*1.0E9_dp, KIND=int_8)
    1110             : 
    1111             :       ! choose atomic range for λ ("i_atom"), ν ("j_atom") in
    1112             :       ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
    1113             :       ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
    1114             :       ! such that M and N fit into the memory
    1115             :       n_atom_per_ivl = INT(SQRT(bs_env%safety_factor_memory*input_memory_per_proc &
    1116             :                                 *bs_env%group_size_tensor/24/bs_env%n_RI &
    1117          28 :                                 /SQRT(bs_env%occupation_3c_int)))/bs_env%max_AO_bf_per_atom
    1118             : 
    1119          28 :       n_intervals_i = (bs_env%n_atom_i - 1)/n_atom_per_ivl + 1
    1120          28 :       n_intervals_j = (bs_env%n_atom_j - 1)/n_atom_per_ivl + 1
    1121             : 
    1122          28 :       bs_env%n_atom_per_interval_ij = n_atom_per_ivl
    1123          28 :       bs_env%n_intervals_i = n_intervals_i
    1124          28 :       bs_env%n_intervals_j = n_intervals_j
    1125             : 
    1126          84 :       ALLOCATE (bs_env%i_atom_intervals(2, n_intervals_i))
    1127          84 :       ALLOCATE (bs_env%j_atom_intervals(2, n_intervals_j))
    1128             : 
    1129          56 :       DO i_ivl = 1, n_intervals_i
    1130          28 :          bs_env%i_atom_intervals(1, i_ivl) = (i_ivl - 1)*n_atom_per_ivl + bs_env%atoms_i(1)
    1131             :          bs_env%i_atom_intervals(2, i_ivl) = MIN(i_ivl*n_atom_per_ivl + bs_env%atoms_i(1) - 1, &
    1132          56 :                                                  bs_env%atoms_i(2))
    1133             :       END DO
    1134             : 
    1135          56 :       DO j_ivl = 1, n_intervals_j
    1136          28 :          bs_env%j_atom_intervals(1, j_ivl) = (j_ivl - 1)*n_atom_per_ivl + bs_env%atoms_j(1)
    1137             :          bs_env%j_atom_intervals(2, j_ivl) = MIN(j_ivl*n_atom_per_ivl + bs_env%atoms_j(1) - 1, &
    1138          56 :                                                  bs_env%atoms_j(2))
    1139             :       END DO
    1140             : 
    1141         112 :       ALLOCATE (bs_env%skip_Sigma_occ(n_intervals_i, n_intervals_j))
    1142          84 :       ALLOCATE (bs_env%skip_Sigma_vir(n_intervals_i, n_intervals_j))
    1143          84 :       bs_env%skip_Sigma_occ(:, :) = .FALSE.
    1144          84 :       bs_env%skip_Sigma_vir(:, :) = .FALSE.
    1145             : 
    1146             :       ! choose atomic range for µ and σ ("inner loop (IL) atom") in
    1147             :       ! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
    1148             :       ! N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
    1149             :       n_atom_per_IL_ivl = MIN(INT(bs_env%safety_factor_memory*input_memory_per_proc &
    1150             :                                   *bs_env%group_size_tensor/n_atom_per_ivl &
    1151             :                                   /bs_env%max_AO_bf_per_atom &
    1152             :                                   /bs_env%n_RI/8/SQRT(bs_env%occupation_3c_int) &
    1153          28 :                                   /bs_env%max_AO_bf_per_atom), bs_env%n_atom)
    1154             : 
    1155          28 :       n_intervals_inner_loop_atoms = (bs_env%n_atom - 1)/n_atom_per_IL_ivl + 1
    1156             : 
    1157          28 :       bs_env%n_atom_per_IL_interval = n_atom_per_IL_ivl
    1158          28 :       bs_env%n_intervals_inner_loop_atoms = n_intervals_inner_loop_atoms
    1159             : 
    1160          84 :       ALLOCATE (bs_env%inner_loop_atom_intervals(2, n_intervals_inner_loop_atoms))
    1161          56 :       DO IL_ivl = 1, n_intervals_inner_loop_atoms
    1162          28 :          bs_env%inner_loop_atom_intervals(1, IL_ivl) = (IL_ivl - 1)*n_atom_per_IL_ivl + 1
    1163          56 :          bs_env%inner_loop_atom_intervals(2, IL_ivl) = MIN(IL_ivl*n_atom_per_IL_ivl, bs_env%n_atom)
    1164             :       END DO
    1165             : 
    1166          28 :       u = bs_env%unit_nr
    1167          28 :       IF (u > 0) THEN
    1168          14 :          WRITE (u, '(T2,A)') ''
    1169          14 :          WRITE (u, '(T2,A,I33)') 'Number of i and j atoms in M_λνP(τ), N_νλQ(τ):', n_atom_per_ivl
    1170          14 :          WRITE (u, '(T2,A,I18)') 'Number of inner loop atoms for µ in M_λνP = sum_µ (µν|P) G_µλ', &
    1171          28 :             n_atom_per_IL_ivl
    1172             :       END IF
    1173             : 
    1174          28 :       CALL timestop(handle)
    1175             : 
    1176          28 :    END SUBROUTINE set_sparsity_parallelization_parameters
    1177             : 
    1178             : ! **************************************************************************************************
    1179             : !> \brief ...
    1180             : !> \param qs_env ...
    1181             : !> \param bs_env ...
    1182             : ! **************************************************************************************************
    1183          28 :    SUBROUTINE check_for_restart_files(qs_env, bs_env)
    1184             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1185             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1186             : 
    1187             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'check_for_restart_files'
    1188             : 
    1189             :       CHARACTER(LEN=9)                                   :: frmt
    1190             :       CHARACTER(LEN=default_string_length)               :: f_chi, f_S_n, f_S_p, f_S_x, f_W_t, &
    1191             :                                                             prefix, project_name
    1192             :       INTEGER                                            :: handle, i_spin, i_t_or_w, ind, n_spin, &
    1193             :                                                             num_time_freq_points
    1194             :       LOGICAL                                            :: chi_exists, Sigma_neg_time_exists, &
    1195             :                                                             Sigma_pos_time_exists, &
    1196             :                                                             Sigma_x_spin_exists, W_time_exists
    1197             :       TYPE(cp_logger_type), POINTER                      :: logger
    1198             :       TYPE(section_vals_type), POINTER                   :: input, print_key
    1199             : 
    1200          28 :       CALL timeset(routineN, handle)
    1201             : 
    1202          28 :       num_time_freq_points = bs_env%num_time_freq_points
    1203          28 :       n_spin = bs_env%n_spin
    1204             : 
    1205          84 :       ALLOCATE (bs_env%read_chi(num_time_freq_points))
    1206          56 :       ALLOCATE (bs_env%calc_chi(num_time_freq_points))
    1207         112 :       ALLOCATE (bs_env%Sigma_c_exists(num_time_freq_points, n_spin))
    1208             : 
    1209          28 :       CALL get_qs_env(qs_env, input=input)
    1210             : 
    1211          28 :       logger => cp_get_default_logger()
    1212          28 :       print_key => section_vals_get_subs_vals(input, 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART')
    1213             :       project_name = cp_print_key_generate_filename(logger, print_key, extension="", &
    1214          28 :                                                     my_local=.FALSE.)
    1215          28 :       WRITE (prefix, '(2A)') TRIM(project_name), "-RESTART_"
    1216          28 :       bs_env%prefix = prefix
    1217             : 
    1218          28 :       bs_env%all_W_exist = .TRUE.
    1219             : 
    1220         412 :       DO i_t_or_w = 1, num_time_freq_points
    1221             : 
    1222         384 :          IF (i_t_or_w < 10) THEN
    1223         240 :             WRITE (frmt, '(A)') '(3A,I1,A)'
    1224         240 :             WRITE (f_chi, frmt) TRIM(prefix), bs_env%chi_name, "_0", i_t_or_w, ".matrix"
    1225         240 :             WRITE (f_W_t, frmt) TRIM(prefix), bs_env%W_time_name, "_0", i_t_or_w, ".matrix"
    1226         144 :          ELSE IF (i_t_or_w < 100) THEN
    1227         144 :             WRITE (frmt, '(A)') '(3A,I2,A)'
    1228         144 :             WRITE (f_chi, frmt) TRIM(prefix), bs_env%chi_name, "_", i_t_or_w, ".matrix"
    1229         144 :             WRITE (f_W_t, frmt) TRIM(prefix), bs_env%W_time_name, "_", i_t_or_w, ".matrix"
    1230             :          ELSE
    1231           0 :             CPABORT('Please implement more than 99 time/frequency points.')
    1232             :          END IF
    1233             : 
    1234         384 :          INQUIRE (file=TRIM(f_chi), exist=chi_exists)
    1235         384 :          INQUIRE (file=TRIM(f_W_t), exist=W_time_exists)
    1236             : 
    1237         384 :          bs_env%read_chi(i_t_or_w) = chi_exists
    1238         384 :          bs_env%calc_chi(i_t_or_w) = .NOT. chi_exists
    1239             : 
    1240         384 :          bs_env%all_W_exist = bs_env%all_W_exist .AND. W_time_exists
    1241             : 
    1242             :          ! the self-energy is spin-dependent
    1243         876 :          DO i_spin = 1, n_spin
    1244             : 
    1245         464 :             ind = i_t_or_w + (i_spin - 1)*num_time_freq_points
    1246             : 
    1247         464 :             IF (ind < 10) THEN
    1248         240 :                WRITE (frmt, '(A)') '(3A,I1,A)'
    1249         240 :                WRITE (f_S_p, frmt) TRIM(prefix), bs_env%Sigma_p_name, "_0", ind, ".matrix"
    1250         240 :                WRITE (f_S_n, frmt) TRIM(prefix), bs_env%Sigma_n_name, "_0", ind, ".matrix"
    1251         224 :             ELSE IF (i_t_or_w < 100) THEN
    1252         224 :                WRITE (frmt, '(A)') '(3A,I2,A)'
    1253         224 :                WRITE (f_S_p, frmt) TRIM(prefix), bs_env%Sigma_p_name, "_", ind, ".matrix"
    1254         224 :                WRITE (f_S_n, frmt) TRIM(prefix), bs_env%Sigma_n_name, "_", ind, ".matrix"
    1255             :             END IF
    1256             : 
    1257         464 :             INQUIRE (file=TRIM(f_S_p), exist=Sigma_pos_time_exists)
    1258         464 :             INQUIRE (file=TRIM(f_S_n), exist=Sigma_neg_time_exists)
    1259             : 
    1260             :             bs_env%Sigma_c_exists(i_t_or_w, i_spin) = Sigma_pos_time_exists .AND. &
    1261        1072 :                                                       Sigma_neg_time_exists
    1262             : 
    1263             :          END DO
    1264             : 
    1265             :       END DO
    1266             : 
    1267             :       ! Marek : In the RTBSE run, check also for zero frequency W
    1268          28 :       IF (bs_env%rtp_method == rtp_method_bse) THEN
    1269          12 :          WRITE (f_W_t, '(3A,I1,A)') TRIM(prefix), "W_freq_rtp", "_0", 0, ".matrix"
    1270          12 :          INQUIRE (file=TRIM(f_W_t), exist=W_time_exists)
    1271          18 :          bs_env%all_W_exist = bs_env%all_W_exist .AND. W_time_exists
    1272             :       END IF
    1273             : 
    1274          28 :       IF (bs_env%all_W_exist) THEN
    1275         192 :          bs_env%read_chi(:) = .FALSE.
    1276         192 :          bs_env%calc_chi(:) = .FALSE.
    1277             :       END IF
    1278             : 
    1279          28 :       bs_env%Sigma_x_exists = .TRUE.
    1280          64 :       DO i_spin = 1, n_spin
    1281          36 :          WRITE (f_S_x, '(3A,I1,A)') TRIM(prefix), bs_env%Sigma_x_name, "_0", i_spin, ".matrix"
    1282          36 :          INQUIRE (file=TRIM(f_S_x), exist=Sigma_x_spin_exists)
    1283          82 :          bs_env%Sigma_x_exists = bs_env%Sigma_x_exists .AND. Sigma_x_spin_exists
    1284             :       END DO
    1285             : 
    1286          28 :       CALL timestop(handle)
    1287             : 
    1288          28 :    END SUBROUTINE check_for_restart_files
    1289             : 
    1290             : ! **************************************************************************************************
    1291             : !> \brief ...
    1292             : !> \param qs_env ...
    1293             : !> \param bs_env ...
    1294             : ! **************************************************************************************************
    1295          36 :    SUBROUTINE set_parallelization_parameters(qs_env, bs_env)
    1296             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1297             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1298             : 
    1299             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'set_parallelization_parameters'
    1300             : 
    1301             :       INTEGER                                            :: color_sub, dummy_1, dummy_2, handle, &
    1302             :                                                             num_pe, num_t_groups, u
    1303             :       INTEGER(KIND=int_8)                                :: mem
    1304             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1305             : 
    1306          36 :       CALL timeset(routineN, handle)
    1307             : 
    1308          36 :       CALL get_qs_env(qs_env, para_env=para_env)
    1309             : 
    1310          36 :       num_pe = para_env%num_pe
    1311             :       ! if not already set, use all processors for the group (for large-cell GW, performance
    1312             :       ! seems to be best for a single group with all MPI processes per group)
    1313          36 :       IF (bs_env%group_size_tensor < 0 .OR. bs_env%group_size_tensor > num_pe) &
    1314          28 :          bs_env%group_size_tensor = num_pe
    1315             : 
    1316             :       ! group_size_tensor must divide num_pe without rest; otherwise everything will be complicated
    1317          36 :       IF (MODULO(num_pe, bs_env%group_size_tensor) .NE. 0) THEN
    1318           0 :          CALL find_good_group_size(num_pe, bs_env%group_size_tensor)
    1319             :       END IF
    1320             : 
    1321             :       ! para_env_tensor for tensor subgroups
    1322          36 :       color_sub = para_env%mepos/bs_env%group_size_tensor
    1323          36 :       bs_env%tensor_group_color = color_sub
    1324             : 
    1325          36 :       ALLOCATE (bs_env%para_env_tensor)
    1326          36 :       CALL bs_env%para_env_tensor%from_split(para_env, color_sub)
    1327             : 
    1328          36 :       num_t_groups = para_env%num_pe/bs_env%group_size_tensor
    1329          36 :       bs_env%num_tensor_groups = num_t_groups
    1330             : 
    1331             :       CALL get_i_j_atoms(bs_env%atoms_i, bs_env%atoms_j, bs_env%n_atom_i, bs_env%n_atom_j, &
    1332          36 :                          color_sub, bs_env)
    1333             : 
    1334         108 :       ALLOCATE (bs_env%atoms_i_t_group(2, num_t_groups))
    1335         108 :       ALLOCATE (bs_env%atoms_j_t_group(2, num_t_groups))
    1336          80 :       DO color_sub = 0, num_t_groups - 1
    1337             :          CALL get_i_j_atoms(bs_env%atoms_i_t_group(1:2, color_sub + 1), &
    1338             :                             bs_env%atoms_j_t_group(1:2, color_sub + 1), &
    1339          80 :                             dummy_1, dummy_2, color_sub, bs_env)
    1340             :       END DO
    1341             : 
    1342          36 :       CALL m_memory(mem)
    1343          36 :       CALL bs_env%para_env%max(mem)
    1344             : 
    1345          36 :       u = bs_env%unit_nr
    1346          36 :       IF (u > 0) THEN
    1347          18 :          WRITE (u, '(T2,A,I47)') 'Group size for tensor operations', bs_env%group_size_tensor
    1348          18 :          IF (bs_env%group_size_tensor > 1 .AND. bs_env%n_atom < 5) THEN
    1349          14 :             WRITE (u, '(T2,A)') 'The requested group size is > 1 which can lead to bad performance.'
    1350          14 :             WRITE (u, '(T2,A)') 'Using more memory per MPI process might improve performance.'
    1351          14 :             WRITE (u, '(T2,A)') '(Also increase MEMORY_PER_PROC when using more memory per process.)'
    1352             :          END IF
    1353             :       END IF
    1354             : 
    1355          36 :       CALL timestop(handle)
    1356             : 
    1357          72 :    END SUBROUTINE set_parallelization_parameters
    1358             : 
    1359             : ! **************************************************************************************************
    1360             : !> \brief ...
    1361             : !> \param num_pe ...
    1362             : !> \param group_size ...
    1363             : ! **************************************************************************************************
    1364           0 :    SUBROUTINE find_good_group_size(num_pe, group_size)
    1365             : 
    1366             :       INTEGER                                            :: num_pe, group_size
    1367             : 
    1368             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'find_good_group_size'
    1369             : 
    1370             :       INTEGER                                            :: group_size_minus, group_size_orig, &
    1371             :                                                             group_size_plus, handle, i_diff
    1372             : 
    1373           0 :       CALL timeset(routineN, handle)
    1374             : 
    1375           0 :       group_size_orig = group_size
    1376             : 
    1377           0 :       DO i_diff = 1, num_pe
    1378             : 
    1379           0 :          group_size_minus = group_size - i_diff
    1380             : 
    1381           0 :          IF (MODULO(num_pe, group_size_minus) == 0 .AND. group_size_minus > 0) THEN
    1382           0 :             group_size = group_size_minus
    1383           0 :             EXIT
    1384             :          END IF
    1385             : 
    1386           0 :          group_size_plus = group_size + i_diff
    1387             : 
    1388           0 :          IF (MODULO(num_pe, group_size_plus) == 0 .AND. group_size_plus <= num_pe) THEN
    1389           0 :             group_size = group_size_plus
    1390           0 :             EXIT
    1391             :          END IF
    1392             : 
    1393             :       END DO
    1394             : 
    1395           0 :       IF (group_size_orig == group_size) CPABORT("Group size error")
    1396             : 
    1397           0 :       CALL timestop(handle)
    1398             : 
    1399           0 :    END SUBROUTINE find_good_group_size
    1400             : 
    1401             : ! **************************************************************************************************
    1402             : !> \brief ...
    1403             : !> \param atoms_i ...
    1404             : !> \param atoms_j ...
    1405             : !> \param n_atom_i ...
    1406             : !> \param n_atom_j ...
    1407             : !> \param color_sub ...
    1408             : !> \param bs_env ...
    1409             : ! **************************************************************************************************
    1410          80 :    SUBROUTINE get_i_j_atoms(atoms_i, atoms_j, n_atom_i, n_atom_j, color_sub, bs_env)
    1411             : 
    1412             :       INTEGER, DIMENSION(2)                              :: atoms_i, atoms_j
    1413             :       INTEGER                                            :: n_atom_i, n_atom_j, color_sub
    1414             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1415             : 
    1416             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_i_j_atoms'
    1417             : 
    1418             :       INTEGER                                            :: handle, i_atoms_per_group, i_group, &
    1419             :                                                             ipcol, ipcol_loop, iprow, iprow_loop, &
    1420             :                                                             j_atoms_per_group, npcol, nprow
    1421             : 
    1422          80 :       CALL timeset(routineN, handle)
    1423             : 
    1424             :       ! create a square mesh of tensor groups for iatom and jatom; code from blacs_env_create
    1425          80 :       CALL square_mesh(nprow, npcol, bs_env%num_tensor_groups)
    1426             : 
    1427          80 :       i_group = 0
    1428         160 :       DO ipcol_loop = 0, npcol - 1
    1429         264 :          DO iprow_loop = 0, nprow - 1
    1430         104 :             IF (i_group == color_sub) THEN
    1431          80 :                iprow = iprow_loop
    1432          80 :                ipcol = ipcol_loop
    1433             :             END IF
    1434         184 :             i_group = i_group + 1
    1435             :          END DO
    1436             :       END DO
    1437             : 
    1438          80 :       IF (MODULO(bs_env%n_atom, nprow) == 0) THEN
    1439          68 :          i_atoms_per_group = bs_env%n_atom/nprow
    1440             :       ELSE
    1441          12 :          i_atoms_per_group = bs_env%n_atom/nprow + 1
    1442             :       END IF
    1443             : 
    1444          80 :       IF (MODULO(bs_env%n_atom, npcol) == 0) THEN
    1445          80 :          j_atoms_per_group = bs_env%n_atom/npcol
    1446             :       ELSE
    1447           0 :          j_atoms_per_group = bs_env%n_atom/npcol + 1
    1448             :       END IF
    1449             : 
    1450          80 :       atoms_i(1) = iprow*i_atoms_per_group + 1
    1451          80 :       atoms_i(2) = MIN((iprow + 1)*i_atoms_per_group, bs_env%n_atom)
    1452          80 :       n_atom_i = atoms_i(2) - atoms_i(1) + 1
    1453             : 
    1454          80 :       atoms_j(1) = ipcol*j_atoms_per_group + 1
    1455          80 :       atoms_j(2) = MIN((ipcol + 1)*j_atoms_per_group, bs_env%n_atom)
    1456          80 :       n_atom_j = atoms_j(2) - atoms_j(1) + 1
    1457             : 
    1458          80 :       CALL timestop(handle)
    1459             : 
    1460          80 :    END SUBROUTINE get_i_j_atoms
    1461             : 
    1462             : ! **************************************************************************************************
    1463             : !> \brief ...
    1464             : !> \param nprow ...
    1465             : !> \param npcol ...
    1466             : !> \param nproc ...
    1467             : ! **************************************************************************************************
    1468          80 :    SUBROUTINE square_mesh(nprow, npcol, nproc)
    1469             :       INTEGER                                            :: nprow, npcol, nproc
    1470             : 
    1471             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'square_mesh'
    1472             : 
    1473             :       INTEGER                                            :: gcd_max, handle, ipe, jpe
    1474             : 
    1475          80 :       CALL timeset(routineN, handle)
    1476             : 
    1477          80 :       gcd_max = -1
    1478         184 :       DO ipe = 1, CEILING(SQRT(REAL(nproc, dp)))
    1479         104 :          jpe = nproc/ipe
    1480         104 :          IF (ipe*jpe .NE. nproc) CYCLE
    1481         184 :          IF (gcd(ipe, jpe) >= gcd_max) THEN
    1482         104 :             nprow = ipe
    1483         104 :             npcol = jpe
    1484         104 :             gcd_max = gcd(ipe, jpe)
    1485             :          END IF
    1486             :       END DO
    1487             : 
    1488          80 :       CALL timestop(handle)
    1489             : 
    1490          80 :    END SUBROUTINE square_mesh
    1491             : 
    1492             : ! **************************************************************************************************
    1493             : !> \brief ...
    1494             : !> \param bs_env ...
    1495             : !> \param qs_env ...
    1496             : ! **************************************************************************************************
    1497          36 :    SUBROUTINE set_heuristic_parameters(bs_env, qs_env)
    1498             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1499             :       TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
    1500             : 
    1501             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'set_heuristic_parameters'
    1502             : 
    1503             :       INTEGER                                            :: handle, u
    1504             :       LOGICAL                                            :: do_BvK_cell
    1505             : 
    1506          36 :       CALL timeset(routineN, handle)
    1507             : 
    1508             :       ! for generating numerically stable minimax Fourier integration weights
    1509          36 :       bs_env%num_points_per_magnitude = 200
    1510             : 
    1511             :       ! for periodic systems and for 20 minimax points, we use a regularized minimax mesh
    1512             :       ! (from experience: regularized minimax meshes converges faster for periodic systems
    1513             :       !  and for 20 pts)
    1514         144 :       IF (SUM(bs_env%periodic) .NE. 0 .OR. bs_env%num_time_freq_points >= 20) THEN
    1515          36 :          bs_env%regularization_minimax = 1.0E-6_dp
    1516             :       ELSE
    1517           0 :          bs_env%regularization_minimax = 0.0_dp
    1518             :       END IF
    1519             : 
    1520          36 :       bs_env%stabilize_exp = 70.0_dp
    1521          36 :       bs_env%eps_atom_grid_2d_mat = 1.0E-50_dp
    1522             : 
    1523             :       ! use a 16-parameter Padé fit
    1524          36 :       bs_env%nparam_pade = 16
    1525             : 
    1526             :       ! resolution of the identity with the truncated Coulomb metric, cutoff radius 3 Angström
    1527          36 :       bs_env%ri_metric%potential_type = do_potential_truncated
    1528          36 :       bs_env%ri_metric%omega = 0.0_dp
    1529             :       ! cutoff radius is specified in the input
    1530          36 :       bs_env%ri_metric%filename = "t_c_g.dat"
    1531             : 
    1532          36 :       bs_env%eps_eigval_mat_RI = 0.0_dp
    1533             : 
    1534          36 :       IF (bs_env%input_regularization_RI > -1.0E-12_dp) THEN
    1535           0 :          bs_env%regularization_RI = bs_env%input_regularization_RI
    1536             :       ELSE
    1537             :          ! default case:
    1538             : 
    1539             :          ! 1. for periodic systems, we use the regularized resolution of the identity per default
    1540          36 :          bs_env%regularization_RI = 1.0E-2_dp
    1541             : 
    1542             :          ! 2. for molecules, no regularization is necessary
    1543         144 :          IF (SUM(bs_env%periodic) == 0) bs_env%regularization_RI = 0.0_dp
    1544             : 
    1545             :       END IF
    1546             : 
    1547             :       ! truncated Coulomb operator for exchange self-energy
    1548             :       ! (see details in Guidon, VandeVondele, Hutter, JCTC 5, 3010 (2009) and references therein)
    1549          36 :       do_BvK_cell = bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp
    1550             :       CALL trunc_coulomb_for_exchange(qs_env, bs_env%trunc_coulomb, &
    1551             :                                       rel_cutoff_trunc_coulomb_ri_x=0.5_dp, &
    1552             :                                       cell_grid=bs_env%cell_grid_scf_desymm, &
    1553          36 :                                       do_BvK_cell=do_BvK_cell)
    1554             : 
    1555             :       ! for small-cell GW, we need more cells than normally used by the filter bs_env%eps_filter
    1556             :       ! (in particular for computing the self-energy because of higher number of cells needed)
    1557          36 :       bs_env%heuristic_filter_factor = 1.0E-4
    1558             : 
    1559          36 :       u = bs_env%unit_nr
    1560          36 :       IF (u > 0) THEN
    1561          18 :          WRITE (u, FMT="(T2,2A,F21.1,A)") "Cutoff radius for the truncated Coulomb ", &
    1562          36 :             "operator in Σ^x:", bs_env%trunc_coulomb%cutoff_radius*angstrom, " Å"
    1563          18 :          WRITE (u, FMT="(T2,2A,F15.1,A)") "Cutoff radius for the truncated Coulomb ", &
    1564          36 :             "operator in RI metric:", bs_env%ri_metric%cutoff_radius*angstrom, " Å"
    1565          18 :          WRITE (u, FMT="(T2,A,ES48.1)") "Regularization parameter of RI ", bs_env%regularization_RI
    1566          18 :          WRITE (u, FMT="(T2,A,I53)") "Lattice sum size for V(k):", bs_env%size_lattice_sum_V
    1567             :       END IF
    1568             : 
    1569          36 :       CALL timestop(handle)
    1570             : 
    1571          36 :    END SUBROUTINE set_heuristic_parameters
    1572             : 
    1573             : ! **************************************************************************************************
    1574             : !> \brief ...
    1575             : !> \param bs_env ...
    1576             : ! **************************************************************************************************
    1577          36 :    SUBROUTINE print_header_and_input_parameters(bs_env)
    1578             : 
    1579             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1580             : 
    1581             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_header_and_input_parameters'
    1582             : 
    1583             :       INTEGER                                            :: handle, u
    1584             : 
    1585          36 :       CALL timeset(routineN, handle)
    1586             : 
    1587          36 :       u = bs_env%unit_nr
    1588             : 
    1589          36 :       IF (u > 0) THEN
    1590          18 :          WRITE (u, '(T2,A)') ' '
    1591          18 :          WRITE (u, '(T2,A)') REPEAT('-', 79)
    1592          18 :          WRITE (u, '(T2,A,A78)') '-', '-'
    1593          18 :          WRITE (u, '(T2,A,A46,A32)') '-', 'GW CALCULATION', '-'
    1594          18 :          WRITE (u, '(T2,A,A78)') '-', '-'
    1595          18 :          WRITE (u, '(T2,A)') REPEAT('-', 79)
    1596          18 :          WRITE (u, '(T2,A)') ' '
    1597          18 :          WRITE (u, '(T2,A,I45)') 'Input: Number of time/freq. points', bs_env%num_time_freq_points
    1598          18 :          WRITE (u, "(T2,A,F44.1,A)") 'Input: ω_max for fitting Σ(iω) (eV)', bs_env%freq_max_fit*evolt
    1599          18 :          WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
    1600          36 :             bs_env%eps_filter
    1601          18 :          WRITE (u, "(T2,A,L55)") 'Input: Apply Hedin shift', bs_env%do_hedin_shift
    1602             :       END IF
    1603             : 
    1604          36 :       CALL timestop(handle)
    1605             : 
    1606          36 :    END SUBROUTINE print_header_and_input_parameters
    1607             : 
    1608             : ! **************************************************************************************************
    1609             : !> \brief ...
    1610             : !> \param qs_env ...
    1611             : !> \param bs_env ...
    1612             : ! **************************************************************************************************
    1613          72 :    SUBROUTINE compute_V_xc(qs_env, bs_env)
    1614             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1615             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1616             : 
    1617             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'compute_V_xc'
    1618             : 
    1619             :       INTEGER                                            :: handle, img, ispin, myfun, nimages
    1620             :       LOGICAL                                            :: hf_present
    1621             :       REAL(KIND=dp)                                      :: energy_ex, energy_exc, energy_total, &
    1622             :                                                             myfraction
    1623          36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mat_ks_without_v_xc
    1624          36 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks_kp
    1625             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1626             :       TYPE(qs_energy_type), POINTER                      :: energy
    1627             :       TYPE(section_vals_type), POINTER                   :: hf_section, input, xc_section
    1628             : 
    1629          36 :       CALL timeset(routineN, handle)
    1630             : 
    1631          36 :       CALL get_qs_env(qs_env, input=input, energy=energy, dft_control=dft_control)
    1632             : 
    1633             :       ! previously, dft_control%nimages set to # neighbor cells, revert for Γ-only KS matrix
    1634          36 :       nimages = dft_control%nimages
    1635          36 :       dft_control%nimages = bs_env%nimages_scf
    1636             : 
    1637             :       ! we need to reset XC functional, therefore, get XC input
    1638          36 :       xc_section => section_vals_get_subs_vals(input, "DFT%XC")
    1639          36 :       CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
    1640          36 :       CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=xc_none)
    1641             :       ! IF (ASSOCIATED(section_vals_get_subs_vals(xc_section, "HF", can_return_null=.TRUE.))) THEN
    1642          36 :       hf_section => section_vals_get_subs_vals(input, "DFT%XC%HF", can_return_null=.TRUE.)
    1643          36 :       hf_present = .FALSE.
    1644          36 :       IF (ASSOCIATED(hf_section)) THEN
    1645          36 :          CALL section_vals_get(hf_section, explicit=hf_present)
    1646             :       END IF
    1647          36 :       IF (hf_present) THEN
    1648             :          ! Special case for handling hfx
    1649           0 :          CALL section_vals_val_get(xc_section, "HF%FRACTION", r_val=myfraction)
    1650           0 :          CALL section_vals_val_set(xc_section, "HF%FRACTION", r_val=0.0_dp)
    1651             :       END IF
    1652             : 
    1653             :       ! save the energy before the energy gets updated
    1654          36 :       energy_total = energy%total
    1655          36 :       energy_exc = energy%exc
    1656          36 :       energy_ex = energy%ex
    1657             : 
    1658          64 :       SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
    1659             :       CASE (large_cell_Gamma)
    1660             : 
    1661          28 :          NULLIFY (mat_ks_without_v_xc)
    1662          28 :          CALL dbcsr_allocate_matrix_set(mat_ks_without_v_xc, bs_env%n_spin)
    1663             : 
    1664          64 :          DO ispin = 1, bs_env%n_spin
    1665          36 :             ALLOCATE (mat_ks_without_v_xc(ispin)%matrix)
    1666          64 :             IF (hf_present) THEN
    1667             :                CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix, &
    1668           0 :                                  matrix_type=dbcsr_type_symmetric)
    1669             :             ELSE
    1670          36 :                CALL dbcsr_create(mat_ks_without_v_xc(ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
    1671             :             END IF
    1672             :          END DO
    1673             : 
    1674             :          ! calculate KS-matrix without XC
    1675             :          CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
    1676          28 :                                            ext_ks_matrix=mat_ks_without_v_xc)
    1677             : 
    1678          64 :          DO ispin = 1, bs_env%n_spin
    1679             :             ! transfer dbcsr matrix to fm
    1680          36 :             CALL cp_fm_create(bs_env%fm_V_xc_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
    1681          36 :             CALL copy_dbcsr_to_fm(mat_ks_without_v_xc(ispin)%matrix, bs_env%fm_V_xc_Gamma(ispin))
    1682             : 
    1683             :             ! v_xc = h_ks - h_ks(v_xc = 0)
    1684             :             CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_Gamma(ispin), &
    1685          64 :                                      beta=1.0_dp, matrix_b=bs_env%fm_ks_Gamma(ispin))
    1686             :          END DO
    1687             : 
    1688          28 :          CALL dbcsr_deallocate_matrix_set(mat_ks_without_v_xc)
    1689             : 
    1690             :       CASE (small_cell_full_kp)
    1691             : 
    1692             :          ! calculate KS-matrix without XC
    1693           8 :          CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
    1694           8 :          CALL get_qs_env(qs_env=qs_env, matrix_ks_kp=matrix_ks_kp)
    1695             : 
    1696         256 :          ALLOCATE (bs_env%fm_V_xc_R(dft_control%nimages, bs_env%n_spin))
    1697          52 :          DO ispin = 1, bs_env%n_spin
    1698         232 :             DO img = 1, dft_control%nimages
    1699             :                ! safe fm_V_xc_R in fm_matrix because saving in dbcsr matrix caused trouble...
    1700         216 :                CALL copy_dbcsr_to_fm(matrix_ks_kp(ispin, img)%matrix, bs_env%fm_work_mo(1))
    1701         216 :                CALL cp_fm_create(bs_env%fm_V_xc_R(img, ispin), bs_env%fm_work_mo(1)%matrix_struct)
    1702             :                ! store h_ks(v_xc = 0) in fm_V_xc_R
    1703             :                CALL cp_fm_scale_and_add(alpha=1.0_dp, matrix_a=bs_env%fm_V_xc_R(img, ispin), &
    1704         224 :                                         beta=1.0_dp, matrix_b=bs_env%fm_work_mo(1))
    1705             :             END DO
    1706             :          END DO
    1707             : 
    1708             :       END SELECT
    1709             : 
    1710             :       ! set back the energy
    1711          36 :       energy%total = energy_total
    1712          36 :       energy%exc = energy_exc
    1713          36 :       energy%ex = energy_ex
    1714             : 
    1715             :       ! set back nimages
    1716          36 :       dft_control%nimages = nimages
    1717             : 
    1718             :       ! set the DFT functional and HF fraction back
    1719             :       CALL section_vals_val_set(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
    1720          36 :                                 i_val=myfun)
    1721          36 :       IF (hf_present) THEN
    1722             :          CALL section_vals_val_set(xc_section, "HF%FRACTION", &
    1723           0 :                                    r_val=myfraction)
    1724             :       END IF
    1725             : 
    1726          36 :       IF (bs_env%small_cell_full_kp_or_large_cell_Gamma == small_cell_full_kp) THEN
    1727             :          ! calculate KS-matrix again with XC
    1728           8 :          CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
    1729          16 :          DO ispin = 1, bs_env%n_spin
    1730         232 :             DO img = 1, dft_control%nimages
    1731             :                ! store h_ks in fm_work_mo
    1732         216 :                CALL copy_dbcsr_to_fm(matrix_ks_kp(ispin, img)%matrix, bs_env%fm_work_mo(1))
    1733             :                ! v_xc = h_ks - h_ks(v_xc = 0)
    1734             :                CALL cp_fm_scale_and_add(alpha=-1.0_dp, matrix_a=bs_env%fm_V_xc_R(img, ispin), &
    1735         224 :                                         beta=1.0_dp, matrix_b=bs_env%fm_work_mo(1))
    1736             :             END DO
    1737             :          END DO
    1738             :       END IF
    1739             : 
    1740          36 :       CALL timestop(handle)
    1741             : 
    1742          36 :    END SUBROUTINE compute_V_xc
    1743             : 
    1744             : ! **************************************************************************************************
    1745             : !> \brief ...
    1746             : !> \param bs_env ...
    1747             : ! **************************************************************************************************
    1748          36 :    SUBROUTINE setup_time_and_frequency_minimax_grid(bs_env)
    1749             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1750             : 
    1751             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_time_and_frequency_minimax_grid'
    1752             : 
    1753             :       INTEGER                                            :: handle, homo, i_w, ierr, ispin, j_w, &
    1754             :                                                             n_mo, num_time_freq_points, u
    1755             :       REAL(KIND=dp)                                      :: E_max, E_max_ispin, E_min, E_min_ispin, &
    1756             :                                                             E_range, max_error_min
    1757          36 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: points_and_weights
    1758             : 
    1759          36 :       CALL timeset(routineN, handle)
    1760             : 
    1761          36 :       n_mo = bs_env%n_ao
    1762          36 :       num_time_freq_points = bs_env%num_time_freq_points
    1763             : 
    1764         108 :       ALLOCATE (bs_env%imag_freq_points(num_time_freq_points))
    1765         108 :       ALLOCATE (bs_env%imag_time_points(num_time_freq_points))
    1766         108 :       ALLOCATE (bs_env%imag_time_weights_freq_zero(num_time_freq_points))
    1767         144 :       ALLOCATE (bs_env%weights_cos_t_to_w(num_time_freq_points, num_time_freq_points))
    1768         144 :       ALLOCATE (bs_env%weights_cos_w_to_t(num_time_freq_points, num_time_freq_points))
    1769         144 :       ALLOCATE (bs_env%weights_sin_t_to_w(num_time_freq_points, num_time_freq_points))
    1770             : 
    1771             :       ! minimum and maximum difference between eigenvalues of unoccupied and an occupied MOs
    1772          36 :       E_min = 1000.0_dp
    1773          36 :       E_max = -1000.0_dp
    1774          80 :       DO ispin = 1, bs_env%n_spin
    1775          44 :          homo = bs_env%n_occ(ispin)
    1776          80 :          SELECT CASE (bs_env%small_cell_full_kp_or_large_cell_Gamma)
    1777             :          CASE (large_cell_Gamma)
    1778             :             E_min_ispin = bs_env%eigenval_scf_Gamma(homo + 1, ispin) - &
    1779          36 :                           bs_env%eigenval_scf_Gamma(homo, ispin)
    1780             :             E_max_ispin = bs_env%eigenval_scf_Gamma(n_mo, ispin) - &
    1781          36 :                           bs_env%eigenval_scf_Gamma(1, ispin)
    1782             :          CASE (small_cell_full_kp)
    1783             :             E_min_ispin = MINVAL(bs_env%eigenval_scf(homo + 1, :, ispin)) - &
    1784         528 :                           MAXVAL(bs_env%eigenval_scf(homo, :, ispin))
    1785             :             E_max_ispin = MAXVAL(bs_env%eigenval_scf(n_mo, :, ispin)) - &
    1786         572 :                           MINVAL(bs_env%eigenval_scf(1, :, ispin))
    1787             :          END SELECT
    1788          44 :          E_min = MIN(E_min, E_min_ispin)
    1789          80 :          E_max = MAX(E_max, E_max_ispin)
    1790             :       END DO
    1791             : 
    1792          36 :       E_range = E_max/E_min
    1793             : 
    1794         108 :       ALLOCATE (points_and_weights(2*num_time_freq_points))
    1795             : 
    1796             :       ! frequency points
    1797          36 :       IF (num_time_freq_points .LE. 20) THEN
    1798          36 :          CALL get_rpa_minimax_coeff(num_time_freq_points, E_range, points_and_weights, ierr, .FALSE.)
    1799             :       ELSE
    1800           0 :          CALL get_rpa_minimax_coeff_larger_grid(num_time_freq_points, E_range, points_and_weights)
    1801             :       END IF
    1802             : 
    1803             :       ! one needs to scale the minimax grids, see Azizi, Wilhelm, Golze, Panades-Barrueta,
    1804             :       ! Giantomassi, Rinke, Draxl, Gonze et al., 2 publications
    1805         484 :       bs_env%imag_freq_points(:) = points_and_weights(1:num_time_freq_points)*E_min
    1806             : 
    1807             :       ! determine number of fit points in the interval [0,ω_max] for virt, or [-ω_max,0] for occ
    1808          36 :       bs_env%num_freq_points_fit = 0
    1809         484 :       DO i_w = 1, num_time_freq_points
    1810         484 :          IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
    1811         184 :             bs_env%num_freq_points_fit = bs_env%num_freq_points_fit + 1
    1812             :          END IF
    1813             :       END DO
    1814             : 
    1815             :       ! iω values for the analytic continuation Σ^c_n(iω,k) -> Σ^c_n(ϵ,k)
    1816         108 :       ALLOCATE (bs_env%imag_freq_points_fit(bs_env%num_freq_points_fit))
    1817          36 :       j_w = 0
    1818         484 :       DO i_w = 1, num_time_freq_points
    1819         484 :          IF (bs_env%imag_freq_points(i_w) < bs_env%freq_max_fit) THEN
    1820         184 :             j_w = j_w + 1
    1821         184 :             bs_env%imag_freq_points_fit(j_w) = bs_env%imag_freq_points(i_w)
    1822             :          END IF
    1823             :       END DO
    1824             : 
    1825             :       ! reset the number of Padé parameters if smaller than the number of
    1826             :       ! imaginary-frequency points for the fit
    1827          36 :       IF (bs_env%num_freq_points_fit < bs_env%nparam_pade) THEN
    1828          36 :          bs_env%nparam_pade = bs_env%num_freq_points_fit
    1829             :       END IF
    1830             : 
    1831             :       ! time points
    1832          36 :       IF (num_time_freq_points .LE. 20) THEN
    1833          36 :          CALL get_exp_minimax_coeff(num_time_freq_points, E_range, points_and_weights)
    1834             :       ELSE
    1835           0 :          CALL get_exp_minimax_coeff_gw(num_time_freq_points, E_range, points_and_weights)
    1836             :       END IF
    1837             : 
    1838         484 :       bs_env%imag_time_points(:) = points_and_weights(1:num_time_freq_points)/(2.0_dp*E_min)
    1839         484 :       bs_env%imag_time_weights_freq_zero(:) = points_and_weights(num_time_freq_points + 1:)/(E_min)
    1840             : 
    1841          36 :       DEALLOCATE (points_and_weights)
    1842             : 
    1843          36 :       u = bs_env%unit_nr
    1844          36 :       IF (u > 0) THEN
    1845          18 :          WRITE (u, '(T2,A)') ''
    1846          18 :          WRITE (u, '(T2,A,F55.2)') 'SCF direct band gap (eV)', E_min*evolt
    1847          18 :          WRITE (u, '(T2,A,F53.2)') 'Max. SCF eigval diff. (eV)', E_max*evolt
    1848          18 :          WRITE (u, '(T2,A,F55.2)') 'E-Range for minimax grid', E_range
    1849          18 :          WRITE (u, '(T2,A,I27)') 'Number of Padé parameters for analytic continuation:', &
    1850          36 :             bs_env%nparam_pade
    1851          18 :          WRITE (u, '(T2,A)') ''
    1852             :       END IF
    1853             : 
    1854             :       ! in minimax grids, Fourier transforms t -> w and w -> t are split using
    1855             :       ! e^(iwt) = cos(wt) + i sin(wt); we thus calculate weights for trafos with a cos and
    1856             :       ! sine prefactor; details in Azizi, Wilhelm, Golze, Giantomassi, Panades-Barrueta,
    1857             :       ! Rinke, Draxl, Gonze et al., 2 publications
    1858             : 
    1859             :       ! cosine transform weights imaginary time to imaginary frequency
    1860             :       CALL get_l_sq_wghts_cos_tf_t_to_w(num_time_freq_points, &
    1861             :                                         bs_env%imag_time_points, &
    1862             :                                         bs_env%weights_cos_t_to_w, &
    1863             :                                         bs_env%imag_freq_points, &
    1864             :                                         E_min, E_max, max_error_min, &
    1865             :                                         bs_env%num_points_per_magnitude, &
    1866          36 :                                         bs_env%regularization_minimax)
    1867             : 
    1868             :       ! cosine transform weights imaginary frequency to imaginary time
    1869             :       CALL get_l_sq_wghts_cos_tf_w_to_t(num_time_freq_points, &
    1870             :                                         bs_env%imag_time_points, &
    1871             :                                         bs_env%weights_cos_w_to_t, &
    1872             :                                         bs_env%imag_freq_points, &
    1873             :                                         E_min, E_max, max_error_min, &
    1874             :                                         bs_env%num_points_per_magnitude, &
    1875          36 :                                         bs_env%regularization_minimax)
    1876             : 
    1877             :       ! sine transform weights imaginary time to imaginary frequency
    1878             :       CALL get_l_sq_wghts_sin_tf_t_to_w(num_time_freq_points, &
    1879             :                                         bs_env%imag_time_points, &
    1880             :                                         bs_env%weights_sin_t_to_w, &
    1881             :                                         bs_env%imag_freq_points, &
    1882             :                                         E_min, E_max, max_error_min, &
    1883             :                                         bs_env%num_points_per_magnitude, &
    1884          36 :                                         bs_env%regularization_minimax)
    1885             : 
    1886          36 :       CALL timestop(handle)
    1887             : 
    1888          72 :    END SUBROUTINE setup_time_and_frequency_minimax_grid
    1889             : 
    1890             : ! **************************************************************************************************
    1891             : !> \brief ...
    1892             : !> \param qs_env ...
    1893             : !> \param bs_env ...
    1894             : ! **************************************************************************************************
    1895           8 :    SUBROUTINE setup_cells_3c(qs_env, bs_env)
    1896             : 
    1897             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1898             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    1899             : 
    1900             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'setup_cells_3c'
    1901             : 
    1902             :       INTEGER :: atom_i, atom_j, atom_k, cell_pair_count, handle, i, i_cell_x, i_cell_x_max, &
    1903             :          i_cell_x_min, i_size, ikind, img, j, j_cell, j_cell_max, j_cell_y, j_cell_y_max, &
    1904             :          j_cell_y_min, j_size, k_cell, k_cell_max, k_cell_z, k_cell_z_max, k_cell_z_min, k_size, &
    1905             :          nimage_pairs_3c, nimages_3c, nimages_3c_max, nkind, u
    1906             :       INTEGER(KIND=int_8)                                :: mem_occ_per_proc
    1907           8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: n_other_3c_images_max
    1908           8 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: index_to_cell_3c_max, nblocks_3c_max
    1909             :       INTEGER, DIMENSION(3)                              :: cell_index, n_max
    1910             :       REAL(KIND=dp) :: avail_mem_per_proc_GB, cell_dist, cell_radius_3c, eps, exp_min_ao, &
    1911             :          exp_min_RI, frobenius_norm, mem_3c_GB, mem_occ_per_proc_GB, radius_ao, radius_ao_product, &
    1912             :          radius_RI
    1913           8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: int_3c
    1914           8 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: exp_ao, exp_RI
    1915             : 
    1916           8 :       CALL timeset(routineN, handle)
    1917             : 
    1918           8 :       CALL get_qs_env(qs_env, nkind=nkind)
    1919             : 
    1920           8 :       exp_min_RI = 10.0_dp
    1921           8 :       exp_min_ao = 10.0_dp
    1922             : 
    1923          24 :       DO ikind = 1, nkind
    1924             : 
    1925          16 :          CALL get_gto_basis_set(bs_env%basis_set_RI(ikind)%gto_basis_set, zet=exp_RI)
    1926          16 :          CALL get_gto_basis_set(bs_env%basis_set_ao(ikind)%gto_basis_set, zet=exp_ao)
    1927             : 
    1928             :          ! we need to remove all exponents lower than a lower bound, e.g. 1E-3, because
    1929             :          ! for contracted basis sets, there might be exponents = 0 in zet
    1930          32 :          DO i = 1, SIZE(exp_RI, 1)
    1931          56 :             DO j = 1, SIZE(exp_RI, 2)
    1932          40 :                IF (exp_RI(i, j) < exp_min_RI .AND. exp_RI(i, j) > 1E-3_dp) exp_min_RI = exp_RI(i, j)
    1933             :             END DO
    1934             :          END DO
    1935          88 :          DO i = 1, SIZE(exp_ao, 1)
    1936         192 :             DO j = 1, SIZE(exp_ao, 2)
    1937         176 :                IF (exp_ao(i, j) < exp_min_ao .AND. exp_ao(i, j) > 1E-3_dp) exp_min_ao = exp_ao(i, j)
    1938             :             END DO
    1939             :          END DO
    1940             : 
    1941             :       END DO
    1942             : 
    1943           8 :       eps = bs_env%eps_filter*bs_env%heuristic_filter_factor
    1944             : 
    1945           8 :       radius_ao = SQRT(-LOG(eps)/exp_min_ao)
    1946           8 :       radius_ao_product = SQRT(-LOG(eps)/(2.0_dp*exp_min_ao))
    1947           8 :       radius_RI = SQRT(-LOG(eps)/exp_min_RI)
    1948             : 
    1949             :       ! For a 3c integral (μR υS | P0) we have that cell R and cell S need to be within radius_3c
    1950           8 :       cell_radius_3c = radius_ao_product + radius_RI + bs_env%ri_metric%cutoff_radius
    1951             : 
    1952          32 :       n_max(1:3) = bs_env%periodic(1:3)*30
    1953             : 
    1954           8 :       nimages_3c_max = 0
    1955             : 
    1956           8 :       i_cell_x_min = 0
    1957           8 :       i_cell_x_max = 0
    1958           8 :       j_cell_y_min = 0
    1959           8 :       j_cell_y_max = 0
    1960           8 :       k_cell_z_min = 0
    1961           8 :       k_cell_z_max = 0
    1962             : 
    1963         256 :       DO i_cell_x = -n_max(1), n_max(1)
    1964       15384 :          DO j_cell_y = -n_max(2), n_max(2)
    1965       45144 :             DO k_cell_z = -n_max(3), n_max(3)
    1966             : 
    1967      119072 :                cell_index(1:3) = (/i_cell_x, j_cell_y, k_cell_z/)
    1968             : 
    1969       29768 :                CALL get_cell_dist(cell_index, bs_env%hmat, cell_dist)
    1970             : 
    1971       44896 :                IF (cell_dist < cell_radius_3c) THEN
    1972         184 :                   nimages_3c_max = nimages_3c_max + 1
    1973         184 :                   i_cell_x_min = MIN(i_cell_x_min, i_cell_x)
    1974         184 :                   i_cell_x_max = MAX(i_cell_x_max, i_cell_x)
    1975         184 :                   j_cell_y_min = MIN(j_cell_y_min, j_cell_y)
    1976         184 :                   j_cell_y_max = MAX(j_cell_y_max, j_cell_y)
    1977         184 :                   k_cell_z_min = MIN(k_cell_z_min, k_cell_z)
    1978         184 :                   k_cell_z_max = MAX(k_cell_z_max, k_cell_z)
    1979             :                END IF
    1980             : 
    1981             :             END DO
    1982             :          END DO
    1983             :       END DO
    1984             : 
    1985             :       ! get index_to_cell_3c_max for the maximum possible cell range;
    1986             :       ! compute 3c integrals later in this routine and check really which cell is needed
    1987          24 :       ALLOCATE (index_to_cell_3c_max(nimages_3c_max, 3))
    1988             : 
    1989           8 :       img = 0
    1990         256 :       DO i_cell_x = -n_max(1), n_max(1)
    1991       15384 :          DO j_cell_y = -n_max(2), n_max(2)
    1992       45144 :             DO k_cell_z = -n_max(3), n_max(3)
    1993             : 
    1994      119072 :                cell_index(1:3) = (/i_cell_x, j_cell_y, k_cell_z/)
    1995             : 
    1996       29768 :                CALL get_cell_dist(cell_index, bs_env%hmat, cell_dist)
    1997             : 
    1998       44896 :                IF (cell_dist < cell_radius_3c) THEN
    1999         184 :                   img = img + 1
    2000         736 :                   index_to_cell_3c_max(img, 1:3) = cell_index(1:3)
    2001             :                END IF
    2002             : 
    2003             :             END DO
    2004             :          END DO
    2005             :       END DO
    2006             : 
    2007             :       ! get pairs of R and S which have non-zero 3c integral (μR υS | P0)
    2008          32 :       ALLOCATE (nblocks_3c_max(nimages_3c_max, nimages_3c_max))
    2009        4456 :       nblocks_3c_max(:, :) = 0
    2010             : 
    2011             :       cell_pair_count = 0
    2012         192 :       DO j_cell = 1, nimages_3c_max
    2013        4456 :          DO k_cell = 1, nimages_3c_max
    2014             : 
    2015        4264 :             cell_pair_count = cell_pair_count + 1
    2016             : 
    2017             :             ! trivial parallelization over cell pairs
    2018        4264 :             IF (MODULO(cell_pair_count, bs_env%para_env%num_pe) .NE. bs_env%para_env%mepos) CYCLE
    2019             : 
    2020        7830 :             DO atom_j = 1, bs_env%n_atom
    2021       24556 :             DO atom_k = 1, bs_env%n_atom
    2022       61098 :             DO atom_i = 1, bs_env%n_atom
    2023             : 
    2024       40806 :                j_size = bs_env%i_ao_end_from_atom(atom_j) - bs_env%i_ao_start_from_atom(atom_j) + 1
    2025       40806 :                k_size = bs_env%i_ao_end_from_atom(atom_k) - bs_env%i_ao_start_from_atom(atom_k) + 1
    2026       40806 :                i_size = bs_env%i_RI_end_from_atom(atom_i) - bs_env%i_RI_start_from_atom(atom_i) + 1
    2027             : 
    2028      204030 :                ALLOCATE (int_3c(j_size, k_size, i_size))
    2029             : 
    2030             :                ! compute 3-c int. ( μ(atom j) R , ν (atom k) S | P (atom i) 0 )
    2031             :                ! ("|": truncated Coulomb operator), inside build_3c_integrals: (j k | i)
    2032             :                CALL build_3c_integral_block(int_3c, qs_env, bs_env%ri_metric, &
    2033             :                                             basis_j=bs_env%basis_set_AO, &
    2034             :                                             basis_k=bs_env%basis_set_AO, &
    2035             :                                             basis_i=bs_env%basis_set_RI, &
    2036             :                                             cell_j=index_to_cell_3c_max(j_cell, 1:3), &
    2037             :                                             cell_k=index_to_cell_3c_max(k_cell, 1:3), &
    2038      285642 :                                             atom_k=atom_k, atom_j=atom_j, atom_i=atom_i)
    2039             : 
    2040     2347536 :                frobenius_norm = SQRT(SUM(int_3c(:, :, :)**2))
    2041             : 
    2042       40806 :                DEALLOCATE (int_3c)
    2043             : 
    2044             :                ! we use a higher threshold here to safe memory when storing the 3c integrals
    2045             :                ! in every tensor group
    2046       55584 :                IF (frobenius_norm > eps) THEN
    2047         892 :                   nblocks_3c_max(j_cell, k_cell) = nblocks_3c_max(j_cell, k_cell) + 1
    2048             :                END IF
    2049             : 
    2050             :             END DO
    2051             :             END DO
    2052             :             END DO
    2053             : 
    2054             :          END DO
    2055             :       END DO
    2056             : 
    2057           8 :       CALL bs_env%para_env%sum(nblocks_3c_max)
    2058             : 
    2059          24 :       ALLOCATE (n_other_3c_images_max(nimages_3c_max))
    2060         192 :       n_other_3c_images_max(:) = 0
    2061             : 
    2062           8 :       nimages_3c = 0
    2063           8 :       nimage_pairs_3c = 0
    2064             : 
    2065         192 :       DO j_cell = 1, nimages_3c_max
    2066        4448 :          DO k_cell = 1, nimages_3c_max
    2067        4448 :             IF (nblocks_3c_max(j_cell, k_cell) > 0) THEN
    2068         312 :                n_other_3c_images_max(j_cell) = n_other_3c_images_max(j_cell) + 1
    2069         312 :                nimage_pairs_3c = nimage_pairs_3c + 1
    2070             :             END IF
    2071             :          END DO
    2072             : 
    2073         192 :          IF (n_other_3c_images_max(j_cell) > 0) nimages_3c = nimages_3c + 1
    2074             : 
    2075             :       END DO
    2076             : 
    2077           8 :       bs_env%nimages_3c = nimages_3c
    2078          24 :       ALLOCATE (bs_env%index_to_cell_3c(nimages_3c, 3))
    2079             :       ALLOCATE (bs_env%cell_to_index_3c(i_cell_x_min:i_cell_x_max, &
    2080             :                                         j_cell_y_min:j_cell_y_max, &
    2081          40 :                                         k_cell_z_min:k_cell_z_max))
    2082         352 :       bs_env%cell_to_index_3c(:, :, :) = -1
    2083             : 
    2084          32 :       ALLOCATE (bs_env%nblocks_3c(nimages_3c, nimages_3c))
    2085           8 :       bs_env%nblocks_3c(nimages_3c, nimages_3c) = 0
    2086             : 
    2087           8 :       j_cell = 0
    2088         192 :       DO j_cell_max = 1, nimages_3c_max
    2089         184 :          IF (n_other_3c_images_max(j_cell_max) == 0) CYCLE
    2090          68 :          j_cell = j_cell + 1
    2091         272 :          cell_index(1:3) = index_to_cell_3c_max(j_cell_max, 1:3)
    2092         272 :          bs_env%index_to_cell_3c(j_cell, 1:3) = cell_index(1:3)
    2093          68 :          bs_env%cell_to_index_3c(cell_index(1), cell_index(2), cell_index(3)) = j_cell
    2094             : 
    2095          68 :          k_cell = 0
    2096        1696 :          DO k_cell_max = 1, nimages_3c_max
    2097        1620 :             IF (n_other_3c_images_max(k_cell_max) == 0) CYCLE
    2098         676 :             k_cell = k_cell + 1
    2099             : 
    2100        1804 :             bs_env%nblocks_3c(j_cell, k_cell) = nblocks_3c_max(j_cell_max, k_cell_max)
    2101             :          END DO
    2102             : 
    2103             :       END DO
    2104             : 
    2105             :       ! we use: 8*10^-9 GB / double precision number
    2106             :       mem_3c_GB = REAL(bs_env%n_RI, KIND=dp)*REAL(bs_env%n_ao, KIND=dp)**2 &
    2107           8 :                   *REAL(nimage_pairs_3c, KIND=dp)*8E-9_dp
    2108             : 
    2109           8 :       CALL m_memory(mem_occ_per_proc)
    2110           8 :       CALL bs_env%para_env%max(mem_occ_per_proc)
    2111             : 
    2112           8 :       mem_occ_per_proc_GB = REAL(mem_occ_per_proc, KIND=dp)/1.0E9_dp
    2113             : 
    2114             :       ! number of processors per group that entirely stores the 3c integrals and does tensor ops
    2115           8 :       avail_mem_per_proc_GB = bs_env%input_memory_per_proc_GB - mem_occ_per_proc_GB
    2116             : 
    2117             :       ! careful: downconvering real to integer, 1.9 -> 1; thus add 1.0 for upconversion, 1.9 -> 2
    2118           8 :       bs_env%group_size_tensor = MAX(INT(mem_3c_GB/avail_mem_per_proc_GB + 1.0_dp), 1)
    2119             : 
    2120           8 :       u = bs_env%unit_nr
    2121             : 
    2122           8 :       IF (u > 0) THEN
    2123           4 :          WRITE (u, FMT="(T2,A,F52.1,A)") "Radius of atomic orbitals", radius_ao*angstrom, " Å"
    2124           4 :          WRITE (u, FMT="(T2,A,F55.1,A)") "Radius of RI functions", radius_RI*angstrom, " Å"
    2125           4 :          WRITE (u, FMT="(T2,A,I47)") "Number of cells for 3c integrals", nimages_3c
    2126           4 :          WRITE (u, FMT="(T2,A,I42)") "Number of cell pairs for 3c integrals", nimage_pairs_3c
    2127           4 :          WRITE (u, '(T2,A)') ''
    2128           4 :          WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
    2129           8 :             bs_env%input_memory_per_proc_GB, ' GB'
    2130           4 :          WRITE (u, '(T2,A,F35.1,A)') 'Used memory per MPI process before GW run', &
    2131           8 :             mem_occ_per_proc_GB, ' GB'
    2132           4 :          WRITE (u, '(T2,A,F44.1,A)') 'Memory of three-center integrals', mem_3c_GB, ' GB'
    2133             :       END IF
    2134             : 
    2135           8 :       CALL timestop(handle)
    2136             : 
    2137          24 :    END SUBROUTINE setup_cells_3c
    2138             : 
    2139             : ! **************************************************************************************************
    2140             : !> \brief ...
    2141             : !> \param index_to_cell_1 ...
    2142             : !> \param index_to_cell_2 ...
    2143             : !> \param nimages_1 ...
    2144             : !> \param nimages_2 ...
    2145             : !> \param index_to_cell ...
    2146             : !> \param cell_to_index ...
    2147             : !> \param nimages ...
    2148             : ! **************************************************************************************************
    2149           8 :    SUBROUTINE sum_two_R_grids(index_to_cell_1, index_to_cell_2, nimages_1, nimages_2, &
    2150             :                               index_to_cell, cell_to_index, nimages)
    2151             : 
    2152             :       INTEGER, DIMENSION(:, :)                           :: index_to_cell_1, index_to_cell_2
    2153             :       INTEGER                                            :: nimages_1, nimages_2
    2154             :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: index_to_cell
    2155             :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    2156             :       INTEGER                                            :: nimages
    2157             : 
    2158             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'sum_two_R_grids'
    2159             : 
    2160             :       INTEGER                                            :: handle, i_dim, img_1, img_2, nimages_max
    2161           8 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: index_to_cell_tmp
    2162             :       INTEGER, DIMENSION(3)                              :: cell_1, cell_2, R, R_max, R_min
    2163             : 
    2164           8 :       CALL timeset(routineN, handle)
    2165             : 
    2166          32 :       DO i_dim = 1, 3
    2167         432 :          R_min(i_dim) = MINVAL(index_to_cell_1(:, i_dim)) + MINVAL(index_to_cell_2(:, i_dim))
    2168         464 :          R_max(i_dim) = MAXVAL(index_to_cell_1(:, i_dim)) + MAXVAL(index_to_cell_2(:, i_dim))
    2169             :       END DO
    2170             : 
    2171           8 :       nimages_max = (R_max(1) - R_min(1) + 1)*(R_max(2) - R_min(2) + 1)*(R_max(3) - R_min(3) + 1)
    2172             : 
    2173          24 :       ALLOCATE (index_to_cell_tmp(nimages_max, 3))
    2174         752 :       index_to_cell_tmp(:, :) = -1
    2175             : 
    2176          40 :       ALLOCATE (cell_to_index(R_min(1):R_max(1), R_min(2):R_max(2), R_min(3):R_max(3)))
    2177         440 :       cell_to_index(:, :, :) = -1
    2178             : 
    2179           8 :       nimages = 0
    2180             : 
    2181          76 :       DO img_1 = 1, nimages_1
    2182             : 
    2183         752 :          DO img_2 = 1, nimages_2
    2184             : 
    2185        2704 :             cell_1(1:3) = index_to_cell_1(img_1, 1:3)
    2186        2704 :             cell_2(1:3) = index_to_cell_2(img_2, 1:3)
    2187             : 
    2188        2704 :             R(1:3) = cell_1(1:3) + cell_2(1:3)
    2189             : 
    2190             :             ! check whether we have found a new cell
    2191         744 :             IF (cell_to_index(R(1), R(2), R(3)) == -1) THEN
    2192             : 
    2193         192 :                nimages = nimages + 1
    2194         192 :                cell_to_index(R(1), R(2), R(3)) = nimages
    2195         768 :                index_to_cell_tmp(nimages, 1:3) = R(1:3)
    2196             : 
    2197             :             END IF
    2198             : 
    2199             :          END DO
    2200             : 
    2201             :       END DO
    2202             : 
    2203          24 :       ALLOCATE (index_to_cell(nimages, 3))
    2204         608 :       index_to_cell(:, :) = index_to_cell_tmp(1:nimages, 1:3)
    2205             : 
    2206           8 :       CALL timestop(handle)
    2207             : 
    2208          16 :    END SUBROUTINE sum_two_R_grids
    2209             : 
    2210             : ! **************************************************************************************************
    2211             : !> \brief ...
    2212             : !> \param qs_env ...
    2213             : !> \param bs_env ...
    2214             : ! **************************************************************************************************
    2215           8 :    SUBROUTINE compute_3c_integrals(qs_env, bs_env)
    2216             : 
    2217             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2218             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2219             : 
    2220             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_integrals'
    2221             : 
    2222             :       INTEGER                                            :: handle, j_cell, k_cell, nimages_3c
    2223             : 
    2224           8 :       CALL timeset(routineN, handle)
    2225             : 
    2226           8 :       nimages_3c = bs_env%nimages_3c
    2227         840 :       ALLOCATE (bs_env%t_3c_int(nimages_3c, nimages_3c))
    2228          76 :       DO j_cell = 1, nimages_3c
    2229         752 :          DO k_cell = 1, nimages_3c
    2230         744 :             CALL dbt_create(bs_env%t_RI_AO__AO, bs_env%t_3c_int(j_cell, k_cell))
    2231             :          END DO
    2232             :       END DO
    2233             : 
    2234             :       CALL build_3c_integrals(bs_env%t_3c_int, &
    2235             :                               bs_env%eps_filter, &
    2236             :                               qs_env, &
    2237             :                               bs_env%nl_3c, &
    2238             :                               int_eps=bs_env%eps_filter*0.05_dp, &
    2239             :                               basis_i=bs_env%basis_set_RI, &
    2240             :                               basis_j=bs_env%basis_set_AO, &
    2241             :                               basis_k=bs_env%basis_set_AO, &
    2242             :                               potential_parameter=bs_env%ri_metric, &
    2243             :                               desymmetrize=.FALSE., do_kpoints=.TRUE., cell_sym=.TRUE., &
    2244           8 :                               cell_to_index_ext=bs_env%cell_to_index_3c)
    2245             : 
    2246           8 :       CALL bs_env%para_env%sync()
    2247             : 
    2248           8 :       CALL timestop(handle)
    2249             : 
    2250           8 :    END SUBROUTINE compute_3c_integrals
    2251             : 
    2252             : ! **************************************************************************************************
    2253             : !> \brief ...
    2254             : !> \param cell_index ...
    2255             : !> \param hmat ...
    2256             : !> \param cell_dist ...
    2257             : ! **************************************************************************************************
    2258       59536 :    SUBROUTINE get_cell_dist(cell_index, hmat, cell_dist)
    2259             : 
    2260             :       INTEGER, DIMENSION(3)                              :: cell_index
    2261             :       REAL(KIND=dp)                                      :: hmat(3, 3), cell_dist
    2262             : 
    2263             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_cell_dist'
    2264             : 
    2265             :       INTEGER                                            :: handle, i_dim
    2266             :       INTEGER, DIMENSION(3)                              :: cell_index_adj
    2267             :       REAL(KIND=dp)                                      :: cell_dist_3(3)
    2268             : 
    2269       59536 :       CALL timeset(routineN, handle)
    2270             : 
    2271             :       ! the distance of cells needs to be taken to adjacent neighbors, not
    2272             :       ! between the center of the cells. We thus need to rescale the cell index
    2273      238144 :       DO i_dim = 1, 3
    2274      178608 :          IF (cell_index(i_dim) > 0) cell_index_adj(i_dim) = cell_index(i_dim) - 1
    2275      178608 :          IF (cell_index(i_dim) < 0) cell_index_adj(i_dim) = cell_index(i_dim) + 1
    2276      238144 :          IF (cell_index(i_dim) == 0) cell_index_adj(i_dim) = cell_index(i_dim)
    2277             :       END DO
    2278             : 
    2279      952576 :       cell_dist_3(1:3) = MATMUL(hmat, REAL(cell_index_adj, KIND=dp))
    2280             : 
    2281      238144 :       cell_dist = SQRT(ABS(SUM(cell_dist_3(1:3)**2)))
    2282             : 
    2283       59536 :       CALL timestop(handle)
    2284             : 
    2285       59536 :    END SUBROUTINE get_cell_dist
    2286             : 
    2287             : ! **************************************************************************************************
    2288             : !> \brief ...
    2289             : !> \param qs_env ...
    2290             : !> \param bs_env ...
    2291             : !> \param kpoints ...
    2292             : !> \param do_print ...
    2293             : ! **************************************************************************************************
    2294           0 :    SUBROUTINE setup_kpoints_scf_desymm(qs_env, bs_env, kpoints, do_print)
    2295             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2296             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2297             :       TYPE(kpoint_type), POINTER                         :: kpoints
    2298             : 
    2299             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_kpoints_scf_desymm'
    2300             : 
    2301             :       INTEGER                                            :: handle, i_cell_x, i_dim, img, j_cell_y, &
    2302             :                                                             k_cell_z, nimages, nkp, u
    2303             :       INTEGER, DIMENSION(3)                              :: cell_grid, cixd, nkp_grid
    2304             :       TYPE(kpoint_type), POINTER                         :: kpoints_scf
    2305             : 
    2306             :       LOGICAL:: do_print
    2307             : 
    2308           0 :       CALL timeset(routineN, handle)
    2309             : 
    2310           0 :       NULLIFY (kpoints)
    2311           0 :       CALL kpoint_create(kpoints)
    2312             : 
    2313           0 :       CALL get_qs_env(qs_env=qs_env, kpoints=kpoints_scf)
    2314             : 
    2315           0 :       nkp_grid(1:3) = kpoints_scf%nkp_grid(1:3)
    2316           0 :       nkp = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)
    2317             : 
    2318             :       ! we need in periodic directions at least 2 k-points in the SCF
    2319           0 :       DO i_dim = 1, 3
    2320           0 :          IF (bs_env%periodic(i_dim) == 1) THEN
    2321           0 :             CPASSERT(nkp_grid(i_dim) > 1)
    2322             :          END IF
    2323             :       END DO
    2324             : 
    2325           0 :       kpoints%kp_scheme = "GENERAL"
    2326           0 :       kpoints%nkp_grid(1:3) = nkp_grid(1:3)
    2327           0 :       kpoints%nkp = nkp
    2328           0 :       bs_env%nkp_scf_desymm = nkp
    2329             : 
    2330           0 :       ALLOCATE (kpoints%xkp(1:3, nkp))
    2331           0 :       CALL compute_xkp(kpoints%xkp, 1, nkp, nkp_grid)
    2332             : 
    2333           0 :       ALLOCATE (kpoints%wkp(nkp))
    2334           0 :       kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp)
    2335             : 
    2336             :       ! for example 4x3x6 kpoint grid -> 3x3x5 cell grid because we need the same number of
    2337             :       ! neighbor cells on both sides of the unit cell
    2338           0 :       cell_grid(1:3) = nkp_grid(1:3) - MODULO(nkp_grid(1:3) + 1, 2)
    2339             :       ! cell index: for example for x: from -n_x/2 to +n_x/2, n_x: number of cells in x direction
    2340           0 :       cixd(1:3) = cell_grid(1:3)/2
    2341             : 
    2342           0 :       nimages = cell_grid(1)*cell_grid(2)*cell_grid(3)
    2343             : 
    2344           0 :       bs_env%nimages_scf_desymm = nimages
    2345             : 
    2346           0 :       ALLOCATE (kpoints%cell_to_index(-cixd(1):cixd(1), -cixd(2):cixd(2), -cixd(3):cixd(3)))
    2347           0 :       ALLOCATE (kpoints%index_to_cell(nimages, 3))
    2348             : 
    2349           0 :       img = 0
    2350           0 :       DO i_cell_x = -cixd(1), cixd(1)
    2351           0 :          DO j_cell_y = -cixd(2), cixd(2)
    2352           0 :             DO k_cell_z = -cixd(3), cixd(3)
    2353           0 :                img = img + 1
    2354           0 :                kpoints%cell_to_index(i_cell_x, j_cell_y, k_cell_z) = img
    2355           0 :                kpoints%index_to_cell(img, 1:3) = (/i_cell_x, j_cell_y, k_cell_z/)
    2356             :             END DO
    2357             :          END DO
    2358             :       END DO
    2359             : 
    2360           0 :       u = bs_env%unit_nr
    2361           0 :       IF (u > 0 .AND. do_print) THEN
    2362           0 :          WRITE (u, FMT="(T2,A,I49)") "Number of cells for G, χ, W, Σ", nimages
    2363             :       END IF
    2364             : 
    2365           0 :       CALL timestop(handle)
    2366             : 
    2367           0 :    END SUBROUTINE setup_kpoints_scf_desymm
    2368             : 
    2369             : ! **************************************************************************************************
    2370             : !> \brief ...
    2371             : !> \param bs_env ...
    2372             : ! **************************************************************************************************
    2373           8 :    SUBROUTINE setup_cells_Delta_R(bs_env)
    2374             : 
    2375             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2376             : 
    2377             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_cells_Delta_R'
    2378             : 
    2379             :       INTEGER                                            :: handle
    2380             : 
    2381           8 :       CALL timeset(routineN, handle)
    2382             : 
    2383             :       ! cell sums batch wise for fixed ΔR = S_1 - R_1; for example:
    2384             :       ! Σ_λσ^R = sum_PR1νS1 M^G_λ0,νS1,PR1 M^W_σR,νS1,PR1
    2385             : 
    2386             :       CALL sum_two_R_grids(bs_env%index_to_cell_3c, &
    2387             :                            bs_env%index_to_cell_3c, &
    2388             :                            bs_env%nimages_3c, bs_env%nimages_3c, &
    2389             :                            bs_env%index_to_cell_Delta_R, &
    2390             :                            bs_env%cell_to_index_Delta_R, &
    2391           8 :                            bs_env%nimages_Delta_R)
    2392             : 
    2393           8 :       IF (bs_env%unit_nr > 0) THEN
    2394           4 :          WRITE (bs_env%unit_nr, FMT="(T2,A,I61)") "Number of cells ΔR", bs_env%nimages_Delta_R
    2395             :       END IF
    2396             : 
    2397           8 :       CALL timestop(handle)
    2398             : 
    2399           8 :    END SUBROUTINE setup_cells_Delta_R
    2400             : 
    2401             : ! **************************************************************************************************
    2402             : !> \brief ...
    2403             : !> \param bs_env ...
    2404             : ! **************************************************************************************************
    2405           8 :    SUBROUTINE setup_parallelization_Delta_R(bs_env)
    2406             : 
    2407             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2408             : 
    2409             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_parallelization_Delta_R'
    2410             : 
    2411             :       INTEGER                                            :: handle, i_cell_Delta_R, i_task_local, &
    2412             :                                                             n_tasks_local
    2413           8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: i_cell_Delta_R_group, &
    2414           8 :                                                             n_tensor_ops_Delta_R
    2415             : 
    2416           8 :       CALL timeset(routineN, handle)
    2417             : 
    2418           8 :       CALL compute_n_tensor_ops_Delta_R(bs_env, n_tensor_ops_Delta_R)
    2419             : 
    2420           8 :       CALL compute_Delta_R_dist(bs_env, n_tensor_ops_Delta_R, i_cell_Delta_R_group, n_tasks_local)
    2421             : 
    2422           8 :       bs_env%n_tasks_Delta_R_local = n_tasks_local
    2423             : 
    2424          24 :       ALLOCATE (bs_env%task_Delta_R(n_tasks_local))
    2425             : 
    2426           8 :       i_task_local = 0
    2427         200 :       DO i_cell_Delta_R = 1, bs_env%nimages_Delta_R
    2428             : 
    2429         192 :          IF (i_cell_Delta_R_group(i_cell_Delta_R) /= bs_env%tensor_group_color) CYCLE
    2430             : 
    2431          86 :          i_task_local = i_task_local + 1
    2432             : 
    2433         200 :          bs_env%task_Delta_R(i_task_local) = i_cell_Delta_R
    2434             : 
    2435             :       END DO
    2436             : 
    2437           8 :       CALL timestop(handle)
    2438             : 
    2439          16 :    END SUBROUTINE setup_parallelization_Delta_R
    2440             : 
    2441             : ! **************************************************************************************************
    2442             : !> \brief ...
    2443             : !> \param bs_env ...
    2444             : !> \param n_tensor_ops_Delta_R ...
    2445             : !> \param i_cell_Delta_R_group ...
    2446             : !> \param n_tasks_local ...
    2447             : ! **************************************************************************************************
    2448           8 :    SUBROUTINE compute_Delta_R_dist(bs_env, n_tensor_ops_Delta_R, i_cell_Delta_R_group, n_tasks_local)
    2449             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2450             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: n_tensor_ops_Delta_R, &
    2451             :                                                             i_cell_Delta_R_group
    2452             :       INTEGER                                            :: n_tasks_local
    2453             : 
    2454             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Delta_R_dist'
    2455             : 
    2456             :       INTEGER                                            :: handle, i_Delta_R_max_op, i_group_min, &
    2457             :                                                             nimages_Delta_R, u
    2458           8 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: n_tensor_ops_Delta_R_in_group
    2459             : 
    2460           8 :       CALL timeset(routineN, handle)
    2461             : 
    2462           8 :       nimages_Delta_R = bs_env%nimages_Delta_R
    2463             : 
    2464           8 :       u = bs_env%unit_nr
    2465             : 
    2466           8 :       IF (u > 0 .AND. nimages_Delta_R < bs_env%num_tensor_groups) THEN
    2467           0 :          WRITE (u, FMT="(T2,A,I5,A,I5,A)") "There are only ", nimages_Delta_R, &
    2468           0 :             " tasks to work on but there are ", bs_env%num_tensor_groups, " groups."
    2469           0 :          WRITE (u, FMT="(T2,A)") "Please reduce the number of MPI processes."
    2470           0 :          WRITE (u, '(T2,A)') ''
    2471             :       END IF
    2472             : 
    2473          24 :       ALLOCATE (n_tensor_ops_Delta_R_in_group(bs_env%num_tensor_groups))
    2474          24 :       n_tensor_ops_Delta_R_in_group(:) = 0
    2475          24 :       ALLOCATE (i_cell_Delta_R_group(nimages_Delta_R))
    2476         200 :       i_cell_Delta_R_group(:) = -1
    2477             : 
    2478           8 :       n_tasks_local = 0
    2479             : 
    2480         732 :       DO WHILE (ANY(n_tensor_ops_Delta_R(:) .NE. 0))
    2481             : 
    2482             :          ! get largest element of n_tensor_ops_Delta_R
    2483        5048 :          i_Delta_R_max_op = MAXLOC(n_tensor_ops_Delta_R, 1)
    2484             : 
    2485             :          ! distribute i_Delta_R_max_op to tensor group which has currently the smallest load
    2486         688 :          i_group_min = MINLOC(n_tensor_ops_Delta_R_in_group, 1)
    2487             : 
    2488             :          ! the tensor groups are 0-index based; but i_group_min is 1-index based
    2489         172 :          i_cell_Delta_R_group(i_Delta_R_max_op) = i_group_min - 1
    2490             :          n_tensor_ops_Delta_R_in_group(i_group_min) = n_tensor_ops_Delta_R_in_group(i_group_min) + &
    2491         172 :                                                       n_tensor_ops_Delta_R(i_Delta_R_max_op)
    2492             : 
    2493             :          ! remove i_Delta_R_max_op from n_tensor_ops_Delta_R
    2494         172 :          n_tensor_ops_Delta_R(i_Delta_R_max_op) = 0
    2495             : 
    2496         180 :          IF (bs_env%tensor_group_color == i_group_min - 1) n_tasks_local = n_tasks_local + 1
    2497             : 
    2498             :       END DO
    2499             : 
    2500           8 :       CALL timestop(handle)
    2501             : 
    2502          16 :    END SUBROUTINE compute_Delta_R_dist
    2503             : 
    2504             : ! **************************************************************************************************
    2505             : !> \brief ...
    2506             : !> \param bs_env ...
    2507             : !> \param n_tensor_ops_Delta_R ...
    2508             : ! **************************************************************************************************
    2509           8 :    SUBROUTINE compute_n_tensor_ops_Delta_R(bs_env, n_tensor_ops_Delta_R)
    2510             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2511             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: n_tensor_ops_Delta_R
    2512             : 
    2513             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_n_tensor_ops_Delta_R'
    2514             : 
    2515             :       INTEGER :: handle, i_cell_Delta_R, i_cell_R, i_cell_R1, i_cell_R1_minus_R, i_cell_R2, &
    2516             :          i_cell_R2_m_R1, i_cell_S1, i_cell_S1_m_R1_p_R2, i_cell_S1_minus_R, i_cell_S2, &
    2517             :          nimages_Delta_R
    2518             :       INTEGER, DIMENSION(3) :: cell_DR, cell_m_R1, cell_R, cell_R1, cell_R1_minus_R, cell_R2, &
    2519             :          cell_R2_m_R1, cell_S1, cell_S1_m_R2_p_R1, cell_S1_minus_R, cell_S1_p_S2_m_R1, cell_S2
    2520             :       LOGICAL                                            :: cell_found
    2521             : 
    2522           8 :       CALL timeset(routineN, handle)
    2523             : 
    2524           8 :       nimages_Delta_R = bs_env%nimages_Delta_R
    2525             : 
    2526          24 :       ALLOCATE (n_tensor_ops_Delta_R(nimages_Delta_R))
    2527         200 :       n_tensor_ops_Delta_R(:) = 0
    2528             : 
    2529             :       ! compute number of tensor operations for specific Delta_R
    2530         200 :       DO i_cell_Delta_R = 1, nimages_Delta_R
    2531             : 
    2532         192 :          IF (MODULO(i_cell_Delta_R, bs_env%num_tensor_groups) /= bs_env%tensor_group_color) CYCLE
    2533             : 
    2534        1074 :          DO i_cell_R1 = 1, bs_env%nimages_3c
    2535             : 
    2536        3880 :             cell_R1(1:3) = bs_env%index_to_cell_3c(i_cell_R1, 1:3)
    2537        3880 :             cell_DR(1:3) = bs_env%index_to_cell_Delta_R(i_cell_Delta_R, 1:3)
    2538             : 
    2539             :             ! S_1 = R_1 + ΔR (from ΔR = S_1 - R_1)
    2540             :             CALL add_R(cell_R1, cell_DR, bs_env%index_to_cell_3c, cell_S1, &
    2541         970 :                        cell_found, bs_env%cell_to_index_3c, i_cell_S1)
    2542         970 :             IF (.NOT. cell_found) CYCLE
    2543             : 
    2544        3200 :             DO i_cell_R2 = 1, bs_env%nimages_scf_desymm
    2545             : 
    2546       11520 :                cell_R2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(i_cell_R2, 1:3)
    2547             : 
    2548             :                ! R_2 - R_1
    2549             :                CALL add_R(cell_R2, -cell_R1, bs_env%index_to_cell_3c, cell_R2_m_R1, &
    2550       11520 :                           cell_found, bs_env%cell_to_index_3c, i_cell_R2_m_R1)
    2551        2880 :                IF (.NOT. cell_found) CYCLE
    2552             : 
    2553             :                ! S_1 - R_1 + R_2
    2554             :                CALL add_R(cell_S1, cell_R2_m_R1, bs_env%index_to_cell_3c, cell_S1_m_R2_p_R1, &
    2555        1680 :                           cell_found, bs_env%cell_to_index_3c, i_cell_S1_m_R1_p_R2)
    2556        1680 :                IF (.NOT. cell_found) CYCLE
    2557             : 
    2558        4300 :                n_tensor_ops_Delta_R(i_cell_Delta_R) = n_tensor_ops_Delta_R(i_cell_Delta_R) + 1
    2559             : 
    2560             :             END DO ! i_cell_R2
    2561             : 
    2562        3200 :             DO i_cell_S2 = 1, bs_env%nimages_scf_desymm
    2563             : 
    2564       11520 :                cell_S2(1:3) = bs_env%kpoints_scf_desymm%index_to_cell(i_cell_S2, 1:3)
    2565       11520 :                cell_m_R1(1:3) = -cell_R1(1:3)
    2566       11520 :                cell_S1_p_S2_m_R1(1:3) = cell_S1(1:3) + cell_S2(1:3) - cell_R1(1:3)
    2567             : 
    2568        2880 :                CALL is_cell_in_index_to_cell(cell_m_R1, bs_env%index_to_cell_3c, cell_found)
    2569        2880 :                IF (.NOT. cell_found) CYCLE
    2570             : 
    2571        2394 :                CALL is_cell_in_index_to_cell(cell_S1_p_S2_m_R1, bs_env%index_to_cell_3c, cell_found)
    2572         320 :                IF (.NOT. cell_found) CYCLE
    2573             : 
    2574             :             END DO ! i_cell_S2
    2575             : 
    2576        4362 :             DO i_cell_R = 1, bs_env%nimages_scf_desymm
    2577             : 
    2578       11520 :                cell_R = bs_env%kpoints_scf_desymm%index_to_cell(i_cell_R, 1:3)
    2579             : 
    2580             :                ! R_1 - R
    2581             :                CALL add_R(cell_R1, -cell_R, bs_env%index_to_cell_3c, cell_R1_minus_R, &
    2582       11520 :                           cell_found, bs_env%cell_to_index_3c, i_cell_R1_minus_R)
    2583        2880 :                IF (.NOT. cell_found) CYCLE
    2584             : 
    2585             :                ! S_1 - R
    2586             :                CALL add_R(cell_S1, -cell_R, bs_env%index_to_cell_3c, cell_S1_minus_R, &
    2587        7224 :                           cell_found, bs_env%cell_to_index_3c, i_cell_S1_minus_R)
    2588         970 :                IF (.NOT. cell_found) CYCLE
    2589             : 
    2590             :             END DO ! i_cell_R
    2591             : 
    2592             :          END DO ! i_cell_R1
    2593             : 
    2594             :       END DO ! i_cell_Delta_R
    2595             : 
    2596           8 :       CALL bs_env%para_env%sum(n_tensor_ops_Delta_R)
    2597             : 
    2598           8 :       CALL timestop(handle)
    2599             : 
    2600           8 :    END SUBROUTINE compute_n_tensor_ops_Delta_R
    2601             : 
    2602             : ! **************************************************************************************************
    2603             : !> \brief ...
    2604             : !> \param cell_1 ...
    2605             : !> \param cell_2 ...
    2606             : !> \param index_to_cell ...
    2607             : !> \param cell_1_plus_2 ...
    2608             : !> \param cell_found ...
    2609             : !> \param cell_to_index ...
    2610             : !> \param i_cell_1_plus_2 ...
    2611             : ! **************************************************************************************************
    2612      216490 :    SUBROUTINE add_R(cell_1, cell_2, index_to_cell, cell_1_plus_2, cell_found, &
    2613             :                     cell_to_index, i_cell_1_plus_2)
    2614             : 
    2615             :       INTEGER, DIMENSION(3)                              :: cell_1, cell_2
    2616             :       INTEGER, DIMENSION(:, :)                           :: index_to_cell
    2617             :       INTEGER, DIMENSION(3)                              :: cell_1_plus_2
    2618             :       LOGICAL                                            :: cell_found
    2619             :       INTEGER, DIMENSION(:, :, :), INTENT(IN), &
    2620             :          OPTIONAL, POINTER                               :: cell_to_index
    2621             :       INTEGER, INTENT(OUT), OPTIONAL                     :: i_cell_1_plus_2
    2622             : 
    2623             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'add_R'
    2624             : 
    2625             :       INTEGER                                            :: handle
    2626             : 
    2627      216490 :       CALL timeset(routineN, handle)
    2628             : 
    2629      865960 :       cell_1_plus_2(1:3) = cell_1(1:3) + cell_2(1:3)
    2630             : 
    2631      216490 :       CALL is_cell_in_index_to_cell(cell_1_plus_2, index_to_cell, cell_found)
    2632             : 
    2633      216490 :       IF (PRESENT(i_cell_1_plus_2)) THEN
    2634      216490 :          IF (cell_found) THEN
    2635      116784 :             CPASSERT(PRESENT(cell_to_index))
    2636      116784 :             i_cell_1_plus_2 = cell_to_index(cell_1_plus_2(1), cell_1_plus_2(2), cell_1_plus_2(3))
    2637             :          ELSE
    2638       99706 :             i_cell_1_plus_2 = -1000
    2639             :          END IF
    2640             :       END IF
    2641             : 
    2642      216490 :       CALL timestop(handle)
    2643             : 
    2644      216490 :    END SUBROUTINE add_R
    2645             : 
    2646             : ! **************************************************************************************************
    2647             : !> \brief ...
    2648             : !> \param cell ...
    2649             : !> \param index_to_cell ...
    2650             : !> \param cell_found ...
    2651             : ! **************************************************************************************************
    2652      335686 :    SUBROUTINE is_cell_in_index_to_cell(cell, index_to_cell, cell_found)
    2653             :       INTEGER, DIMENSION(3)                              :: cell
    2654             :       INTEGER, DIMENSION(:, :)                           :: index_to_cell
    2655             :       LOGICAL                                            :: cell_found
    2656             : 
    2657             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'is_cell_in_index_to_cell'
    2658             : 
    2659             :       INTEGER                                            :: handle, i_cell, nimg
    2660             :       INTEGER, DIMENSION(3)                              :: cell_i
    2661             : 
    2662      335686 :       CALL timeset(routineN, handle)
    2663             : 
    2664      335686 :       nimg = SIZE(index_to_cell, 1)
    2665             : 
    2666      335686 :       cell_found = .FALSE.
    2667             : 
    2668     3838008 :       DO i_cell = 1, nimg
    2669             : 
    2670    14009288 :          cell_i(1:3) = index_to_cell(i_cell, 1:3)
    2671             : 
    2672     3838008 :          IF (cell_i(1) == cell(1) .AND. cell_i(2) == cell(2) .AND. cell_i(3) == cell(3)) THEN
    2673      195888 :             cell_found = .TRUE.
    2674             :          END IF
    2675             : 
    2676             :       END DO
    2677             : 
    2678      335686 :       CALL timestop(handle)
    2679             : 
    2680      335686 :    END SUBROUTINE is_cell_in_index_to_cell
    2681             : 
    2682             : ! **************************************************************************************************
    2683             : !> \brief ...
    2684             : !> \param qs_env ...
    2685             : !> \param bs_env ...
    2686             : ! **************************************************************************************************
    2687           8 :    SUBROUTINE allocate_matrices_small_cell_full_kp(qs_env, bs_env)
    2688             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2689             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2690             : 
    2691             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_matrices_small_cell_full_kp'
    2692             : 
    2693             :       INTEGER                                            :: handle, i_spin, i_t, img, n_spin, &
    2694             :                                                             nimages_scf, num_time_freq_points
    2695             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    2696             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2697             : 
    2698           8 :       CALL timeset(routineN, handle)
    2699             : 
    2700           8 :       nimages_scf = bs_env%nimages_scf_desymm
    2701           8 :       num_time_freq_points = bs_env%num_time_freq_points
    2702           8 :       n_spin = bs_env%n_spin
    2703             : 
    2704           8 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
    2705             : 
    2706          96 :       ALLOCATE (bs_env%fm_G_S(nimages_scf))
    2707          96 :       ALLOCATE (bs_env%fm_Sigma_x_R(nimages_scf))
    2708         672 :       ALLOCATE (bs_env%fm_chi_R_t(nimages_scf, num_time_freq_points))
    2709         672 :       ALLOCATE (bs_env%fm_MWM_R_t(nimages_scf, num_time_freq_points))
    2710         688 :       ALLOCATE (bs_env%fm_Sigma_c_R_neg_tau(nimages_scf, num_time_freq_points, n_spin))
    2711         688 :       ALLOCATE (bs_env%fm_Sigma_c_R_pos_tau(nimages_scf, num_time_freq_points, n_spin))
    2712          80 :       DO img = 1, nimages_scf
    2713          72 :          CALL cp_fm_create(bs_env%fm_G_S(img), bs_env%fm_work_mo(1)%matrix_struct)
    2714          72 :          CALL cp_fm_create(bs_env%fm_Sigma_x_R(img), bs_env%fm_work_mo(1)%matrix_struct)
    2715         656 :          DO i_t = 1, num_time_freq_points
    2716         576 :             CALL cp_fm_create(bs_env%fm_chi_R_t(img, i_t), bs_env%fm_RI_RI%matrix_struct)
    2717         576 :             CALL cp_fm_create(bs_env%fm_MWM_R_t(img, i_t), bs_env%fm_RI_RI%matrix_struct)
    2718         576 :             CALL cp_fm_set_all(bs_env%fm_MWM_R_t(img, i_t), 0.0_dp)
    2719        1224 :             DO i_spin = 1, n_spin
    2720             :                CALL cp_fm_create(bs_env%fm_Sigma_c_R_neg_tau(img, i_t, i_spin), &
    2721         576 :                                  bs_env%fm_work_mo(1)%matrix_struct)
    2722             :                CALL cp_fm_create(bs_env%fm_Sigma_c_R_pos_tau(img, i_t, i_spin), &
    2723         576 :                                  bs_env%fm_work_mo(1)%matrix_struct)
    2724         576 :                CALL cp_fm_set_all(bs_env%fm_Sigma_c_R_neg_tau(img, i_t, i_spin), 0.0_dp)
    2725        1152 :                CALL cp_fm_set_all(bs_env%fm_Sigma_c_R_pos_tau(img, i_t, i_spin), 0.0_dp)
    2726             :             END DO
    2727             :          END DO
    2728             :       END DO
    2729             : 
    2730           8 :       CALL timestop(handle)
    2731             : 
    2732           8 :    END SUBROUTINE allocate_matrices_small_cell_full_kp
    2733             : 
    2734             : ! **************************************************************************************************
    2735             : !> \brief ...
    2736             : !> \param qs_env ...
    2737             : !> \param bs_env ...
    2738             : ! **************************************************************************************************
    2739           8 :    SUBROUTINE trafo_V_xc_R_to_kp(qs_env, bs_env)
    2740             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2741             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2742             : 
    2743             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'trafo_V_xc_R_to_kp'
    2744             : 
    2745             :       INTEGER                                            :: handle, ikp, img, ispin, n_ao
    2746           8 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index_scf
    2747             :       TYPE(cp_cfm_type)                                  :: cfm_mo_coeff, cfm_tmp, cfm_V_xc
    2748             :       TYPE(cp_fm_type)                                   :: fm_V_xc_re
    2749           8 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_ks
    2750             :       TYPE(kpoint_type), POINTER                         :: kpoints_scf
    2751             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2752           8 :          POINTER                                         :: sab_nl
    2753             : 
    2754           8 :       CALL timeset(routineN, handle)
    2755             : 
    2756           8 :       n_ao = bs_env%n_ao
    2757             : 
    2758           8 :       CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, kpoints=kpoints_scf)
    2759             : 
    2760           8 :       NULLIFY (sab_nl)
    2761           8 :       CALL get_kpoint_info(kpoints_scf, sab_nl=sab_nl, cell_to_index=cell_to_index_scf)
    2762             : 
    2763           8 :       CALL cp_cfm_create(cfm_V_xc, bs_env%cfm_work_mo%matrix_struct)
    2764           8 :       CALL cp_cfm_create(cfm_mo_coeff, bs_env%cfm_work_mo%matrix_struct)
    2765           8 :       CALL cp_cfm_create(cfm_tmp, bs_env%cfm_work_mo%matrix_struct)
    2766           8 :       CALL cp_fm_create(fm_V_xc_re, bs_env%cfm_work_mo%matrix_struct)
    2767             : 
    2768         224 :       DO img = 1, bs_env%nimages_scf
    2769         440 :          DO ispin = 1, bs_env%n_spin
    2770             :             ! JW kind of hack because the format of matrix_ks remains dubious...
    2771         216 :             CALL dbcsr_set(matrix_ks(ispin, img)%matrix, 0.0_dp)
    2772         432 :             CALL copy_fm_to_dbcsr(bs_env%fm_V_xc_R(img, ispin), matrix_ks(ispin, img)%matrix)
    2773             :          END DO
    2774             :       END DO
    2775             : 
    2776          40 :       ALLOCATE (bs_env%v_xc_n(n_ao, bs_env%nkp_bs_and_DOS, bs_env%n_spin))
    2777             : 
    2778          16 :       DO ispin = 1, bs_env%n_spin
    2779         268 :          DO ikp = 1, bs_env%nkp_bs_and_DOS
    2780             : 
    2781             :             ! v^xc^R -> v^xc(k)  (matrix_ks stores v^xc^R, see SUBROUTINE compute_V_xc)
    2782             :             CALL rsmat_to_kp(matrix_ks, ispin, bs_env%kpoints_DOS%xkp(1:3, ikp), &
    2783         252 :                              cell_to_index_scf, sab_nl, bs_env, cfm_V_xc)
    2784             : 
    2785             :             ! get C_µn(k)
    2786         252 :             CALL cp_cfm_to_cfm(bs_env%cfm_mo_coeff_kp(ikp, ispin), cfm_mo_coeff)
    2787             : 
    2788             :             ! v^xc_nm(k_i) = sum_µν C^*_µn(k_i) v^xc_µν(k_i) C_νn(k_i)
    2789             :             CALL parallel_gemm('N', 'N', n_ao, n_ao, n_ao, z_one, cfm_V_xc, cfm_mo_coeff, &
    2790         252 :                                z_zero, cfm_tmp)
    2791             :             CALL parallel_gemm('C', 'N', n_ao, n_ao, n_ao, z_one, cfm_mo_coeff, cfm_tmp, &
    2792         252 :                                z_zero, cfm_V_xc)
    2793             : 
    2794             :             ! get v^xc_nn(k_i) which is a real quantity as v^xc is Hermitian
    2795         252 :             CALL cp_cfm_to_fm(cfm_V_xc, fm_V_xc_re)
    2796         260 :             CALL cp_fm_get_diag(fm_V_xc_re, bs_env%v_xc_n(:, ikp, ispin))
    2797             : 
    2798             :          END DO
    2799             : 
    2800             :       END DO
    2801             : 
    2802             :       ! just rebuild the overwritten KS matrix again
    2803           8 :       CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
    2804             : 
    2805           8 :       CALL cp_cfm_release(cfm_V_xc)
    2806           8 :       CALL cp_cfm_release(cfm_mo_coeff)
    2807           8 :       CALL cp_cfm_release(cfm_tmp)
    2808           8 :       CALL cp_fm_release(fm_V_xc_re)
    2809             : 
    2810           8 :       CALL timestop(handle)
    2811             : 
    2812          16 :    END SUBROUTINE trafo_V_xc_R_to_kp
    2813             : 
    2814             : ! **************************************************************************************************
    2815             : !> \brief ...
    2816             : !> \param qs_env ...
    2817             : !> \param bs_env ...
    2818             : ! **************************************************************************************************
    2819           8 :    SUBROUTINE heuristic_RI_regularization(qs_env, bs_env)
    2820             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2821             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2822             : 
    2823             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'heuristic_RI_regularization'
    2824             : 
    2825           8 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)  :: M
    2826             :       INTEGER                                            :: handle, ikp, ikp_local, n_RI, nkp, &
    2827             :                                                             nkp_local, u
    2828             :       REAL(KIND=dp)                                      :: cond_nr, cond_nr_max, max_ev, &
    2829             :                                                             max_ev_ikp, min_ev, min_ev_ikp
    2830           8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: M_R
    2831             : 
    2832           8 :       CALL timeset(routineN, handle)
    2833             : 
    2834             :       ! compute M^R_PQ = <phi_P,0|V^tr(rc)|phi_Q,R> for RI metric
    2835           8 :       CALL get_V_tr_R(M_R, bs_env%ri_metric, 0.0_dp, bs_env, qs_env)
    2836             : 
    2837           8 :       nkp = bs_env%nkp_chi_eps_W_orig_plus_extra
    2838           8 :       n_RI = bs_env%n_RI
    2839             : 
    2840           8 :       nkp_local = 0
    2841        5128 :       DO ikp = 1, nkp
    2842             :          ! trivial parallelization over k-points
    2843        5120 :          IF (MODULO(ikp, bs_env%para_env%num_pe) .NE. bs_env%para_env%mepos) CYCLE
    2844        5128 :          nkp_local = nkp_local + 1
    2845             :       END DO
    2846             : 
    2847          40 :       ALLOCATE (M(n_RI, n_RI, nkp_local))
    2848             : 
    2849           8 :       ikp_local = 0
    2850           8 :       cond_nr_max = 0.0_dp
    2851           8 :       min_ev = 1000.0_dp
    2852           8 :       max_ev = -1000.0_dp
    2853             : 
    2854        5128 :       DO ikp = 1, nkp
    2855             : 
    2856             :          ! trivial parallelization
    2857        5120 :          IF (MODULO(ikp, bs_env%para_env%num_pe) .NE. bs_env%para_env%mepos) CYCLE
    2858             : 
    2859        2560 :          ikp_local = ikp_local + 1
    2860             : 
    2861             :          ! M(k) = sum_R e^ikR M^R
    2862             :          CALL trafo_rs_to_ikp(M_R, M(:, :, ikp_local), &
    2863             :                               bs_env%kpoints_scf_desymm%index_to_cell, &
    2864        2560 :                               bs_env%kpoints_chi_eps_W%xkp(1:3, ikp))
    2865             : 
    2866             :          ! compute condition number of M_PQ(k)
    2867        2560 :          CALL power(M(:, :, ikp_local), 1.0_dp, 0.0_dp, cond_nr, min_ev_ikp, max_ev_ikp)
    2868             : 
    2869        2560 :          IF (cond_nr > cond_nr_max) cond_nr_max = cond_nr
    2870        2560 :          IF (max_ev_ikp > max_ev) max_ev = max_ev_ikp
    2871        2568 :          IF (min_ev_ikp < min_ev) min_ev = min_ev_ikp
    2872             : 
    2873             :       END DO ! ikp
    2874             : 
    2875           8 :       CALL bs_env%para_env%max(cond_nr_max)
    2876             : 
    2877           8 :       u = bs_env%unit_nr
    2878           8 :       IF (u > 0) THEN
    2879           4 :          WRITE (u, FMT="(T2,A,ES34.1)") "Min. abs. eigenvalue of RI metric matrix M(k)", min_ev
    2880           4 :          WRITE (u, FMT="(T2,A,ES34.1)") "Max. abs. eigenvalue of RI metric matrix M(k)", max_ev
    2881           4 :          WRITE (u, FMT="(T2,A,ES50.1)") "Max. condition number of M(k)", cond_nr_max
    2882             :       END IF
    2883             : 
    2884           8 :       CALL timestop(handle)
    2885             : 
    2886          16 :    END SUBROUTINE heuristic_RI_regularization
    2887             : 
    2888             : ! **************************************************************************************************
    2889             : !> \brief ...
    2890             : !> \param V_tr_R ...
    2891             : !> \param pot_type ...
    2892             : !> \param regularization_RI ...
    2893             : !> \param bs_env ...
    2894             : !> \param qs_env ...
    2895             : ! **************************************************************************************************
    2896          96 :    SUBROUTINE get_V_tr_R(V_tr_R, pot_type, regularization_RI, bs_env, qs_env)
    2897             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: V_tr_R
    2898             :       TYPE(libint_potential_type)                        :: pot_type
    2899             :       REAL(KIND=dp)                                      :: regularization_RI
    2900             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    2901             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2902             : 
    2903             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_V_tr_R'
    2904             : 
    2905             :       INTEGER                                            :: handle, img, nimages_scf_desymm
    2906             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sizes_RI
    2907          96 :       INTEGER, DIMENSION(:), POINTER                     :: col_bsize, row_bsize
    2908             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    2909          96 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_V_tr_R
    2910             :       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
    2911          96 :       TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: mat_V_tr_R
    2912             :       TYPE(distribution_2d_type), POINTER                :: dist_2d
    2913             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    2914          96 :          POINTER                                         :: sab_RI
    2915          96 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    2916          96 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2917             : 
    2918          96 :       CALL timeset(routineN, handle)
    2919             : 
    2920          96 :       NULLIFY (sab_RI, dist_2d)
    2921             : 
    2922             :       CALL get_qs_env(qs_env=qs_env, &
    2923             :                       blacs_env=blacs_env, &
    2924             :                       distribution_2d=dist_2d, &
    2925             :                       qs_kind_set=qs_kind_set, &
    2926          96 :                       particle_set=particle_set)
    2927             : 
    2928         288 :       ALLOCATE (sizes_RI(bs_env%n_atom))
    2929          96 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=bs_env%basis_set_RI)
    2930             :       CALL build_2c_neighbor_lists(sab_RI, bs_env%basis_set_RI, bs_env%basis_set_RI, &
    2931             :                                    pot_type, "2c_nl_RI", qs_env, sym_ij=.FALSE., &
    2932          96 :                                    dist_2d=dist_2d)
    2933          96 :       CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
    2934         288 :       ALLOCATE (row_bsize(SIZE(sizes_RI)))
    2935         288 :       ALLOCATE (col_bsize(SIZE(sizes_RI)))
    2936         328 :       row_bsize(:) = sizes_RI
    2937         328 :       col_bsize(:) = sizes_RI
    2938             : 
    2939          96 :       nimages_scf_desymm = bs_env%nimages_scf_desymm
    2940        1152 :       ALLOCATE (mat_V_tr_R(nimages_scf_desymm))
    2941             :       CALL dbcsr_create(mat_V_tr_R(1), "(RI|RI)", dbcsr_dist, dbcsr_type_no_symmetry, &
    2942          96 :                         row_bsize, col_bsize)
    2943          96 :       DEALLOCATE (row_bsize, col_bsize)
    2944             : 
    2945         864 :       DO img = 2, nimages_scf_desymm
    2946         864 :          CALL dbcsr_create(mat_V_tr_R(img), template=mat_V_tr_R(1))
    2947             :       END DO
    2948             : 
    2949             :       CALL build_2c_integrals(mat_V_tr_R, 0.0_dp, qs_env, sab_RI, bs_env%basis_set_RI, &
    2950             :                               bs_env%basis_set_RI, pot_type, do_kpoints=.TRUE., &
    2951             :                               ext_kpoints=bs_env%kpoints_scf_desymm, &
    2952          96 :                               regularization_RI=regularization_RI)
    2953             : 
    2954        1152 :       ALLOCATE (fm_V_tr_R(nimages_scf_desymm))
    2955         960 :       DO img = 1, nimages_scf_desymm
    2956         864 :          CALL cp_fm_create(fm_V_tr_R(img), bs_env%fm_RI_RI%matrix_struct)
    2957         864 :          CALL copy_dbcsr_to_fm(mat_V_tr_R(img), fm_V_tr_R(img))
    2958         960 :          CALL dbcsr_release(mat_V_tr_R(img))
    2959             :       END DO
    2960             : 
    2961          96 :       IF (.NOT. ALLOCATED(V_tr_R)) THEN
    2962         480 :          ALLOCATE (V_tr_R(bs_env%n_RI, bs_env%n_RI, nimages_scf_desymm))
    2963             :       END IF
    2964             : 
    2965          96 :       CALL fm_to_local_array(fm_V_tr_R, V_tr_R)
    2966             : 
    2967          96 :       CALL cp_fm_release(fm_V_tr_R)
    2968          96 :       CALL dbcsr_distribution_release(dbcsr_dist)
    2969          96 :       CALL release_neighbor_list_sets(sab_RI)
    2970             : 
    2971          96 :       CALL timestop(handle)
    2972             : 
    2973         288 :    END SUBROUTINE get_V_tr_R
    2974             : 
    2975             : ! **************************************************************************************************
    2976             : !> \brief ...
    2977             : !> \param matrix ...
    2978             : !> \param exponent ...
    2979             : !> \param eps ...
    2980             : !> \param cond_nr ...
    2981             : !> \param min_ev ...
    2982             : !> \param max_ev ...
    2983             : ! **************************************************************************************************
    2984       49216 :    SUBROUTINE power(matrix, exponent, eps, cond_nr, min_ev, max_ev)
    2985             :       COMPLEX(KIND=dp), DIMENSION(:, :)                  :: matrix
    2986             :       REAL(KIND=dp)                                      :: exponent, eps
    2987             :       REAL(KIND=dp), OPTIONAL                            :: cond_nr, min_ev, max_ev
    2988             : 
    2989             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'power'
    2990             : 
    2991       49216 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :)     :: eigenvectors
    2992             :       INTEGER                                            :: handle, i, n
    2993             :       REAL(KIND=dp)                                      :: pos_eval
    2994       49216 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
    2995             : 
    2996       49216 :       CALL timeset(routineN, handle)
    2997             : 
    2998             :       ! make matrix perfectly Hermitian
    2999     3468352 :       matrix(:, :) = 0.5_dp*(matrix(:, :) + CONJG(TRANSPOSE(matrix(:, :))))
    3000             : 
    3001       49216 :       n = SIZE(matrix, 1)
    3002      295296 :       ALLOCATE (eigenvalues(n), eigenvectors(n, n))
    3003       49216 :       CALL diag_complex(matrix, eigenvectors, eigenvalues)
    3004             : 
    3005       82496 :       IF (PRESENT(cond_nr)) cond_nr = MAXVAL(ABS(eigenvalues))/MINVAL(ABS(eigenvalues))
    3006       65856 :       IF (PRESENT(min_ev)) min_ev = MINVAL(ABS(eigenvalues))
    3007       65856 :       IF (PRESENT(max_ev)) max_ev = MAXVAL(ABS(eigenvalues))
    3008             : 
    3009      314720 :       DO i = 1, n
    3010      265504 :          IF (eps < eigenvalues(i)) THEN
    3011      265504 :             pos_eval = (eigenvalues(i))**(0.5_dp*exponent)
    3012             :          ELSE
    3013             :             pos_eval = 0.0_dp
    3014             :          END IF
    3015     1758784 :          eigenvectors(:, i) = eigenvectors(:, i)*pos_eval
    3016             :       END DO
    3017             : 
    3018       49216 :       CALL ZGEMM("N", "C", n, n, n, z_one, eigenvectors, n, eigenvectors, n, z_zero, matrix, n)
    3019             : 
    3020       49216 :       DEALLOCATE (eigenvalues, eigenvectors)
    3021             : 
    3022       49216 :       CALL timestop(handle)
    3023             : 
    3024       49216 :    END SUBROUTINE power
    3025             : 
    3026             : ! **************************************************************************************************
    3027             : !> \brief ...
    3028             : !> \param bs_env ...
    3029             : !> \param Sigma_c_n_time ...
    3030             : !> \param Sigma_c_n_freq ...
    3031             : !> \param ispin ...
    3032             : ! **************************************************************************************************
    3033         308 :    SUBROUTINE time_to_freq(bs_env, Sigma_c_n_time, Sigma_c_n_freq, ispin)
    3034             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    3035             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: Sigma_c_n_time, Sigma_c_n_freq
    3036             :       INTEGER                                            :: ispin
    3037             : 
    3038             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'time_to_freq'
    3039             : 
    3040             :       INTEGER                                            :: handle, i_t, j_w, n_occ
    3041             :       REAL(KIND=dp)                                      :: freq_j, time_i, w_cos_ij, w_sin_ij
    3042         308 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: Sigma_c_n_cos_time, Sigma_c_n_sin_time
    3043             : 
    3044         308 :       CALL timeset(routineN, handle)
    3045             : 
    3046        1232 :       ALLOCATE (Sigma_c_n_cos_time(bs_env%n_ao, bs_env%num_time_freq_points))
    3047         924 :       ALLOCATE (Sigma_c_n_sin_time(bs_env%n_ao, bs_env%num_time_freq_points))
    3048             : 
    3049       30340 :       Sigma_c_n_cos_time(:, :) = 0.5_dp*(Sigma_c_n_time(:, :, 1) + Sigma_c_n_time(:, :, 2))
    3050       30340 :       Sigma_c_n_sin_time(:, :) = 0.5_dp*(Sigma_c_n_time(:, :, 1) - Sigma_c_n_time(:, :, 2))
    3051             : 
    3052       60988 :       Sigma_c_n_freq(:, :, :) = 0.0_dp
    3053             : 
    3054        3220 :       DO i_t = 1, bs_env%num_time_freq_points
    3055             : 
    3056       33012 :          DO j_w = 1, bs_env%num_time_freq_points
    3057             : 
    3058       29792 :             freq_j = bs_env%imag_freq_points(j_w)
    3059       29792 :             time_i = bs_env%imag_time_points(i_t)
    3060             :             ! integration weights for cosine and sine transform
    3061       29792 :             w_cos_ij = bs_env%weights_cos_t_to_w(j_w, i_t)*COS(freq_j*time_i)
    3062       29792 :             w_sin_ij = bs_env%weights_sin_t_to_w(j_w, i_t)*SIN(freq_j*time_i)
    3063             : 
    3064             :             ! 1. Re(Σ^c_nn(k_i,iω)) from cosine transform
    3065             :             Sigma_c_n_freq(:, j_w, 1) = Sigma_c_n_freq(:, j_w, 1) + &
    3066      279872 :                                         w_cos_ij*Sigma_c_n_cos_time(:, i_t)
    3067             : 
    3068             :             ! 2. Im(Σ^c_nn(k_i,iω)) from sine transform
    3069             :             Sigma_c_n_freq(:, j_w, 2) = Sigma_c_n_freq(:, j_w, 2) + &
    3070      282784 :                                         w_sin_ij*Sigma_c_n_sin_time(:, i_t)
    3071             : 
    3072             :          END DO
    3073             : 
    3074             :       END DO
    3075             : 
    3076             :       ! for occupied levels, we need the correlation self-energy for negative omega.
    3077             :       ! Therefore, weight_sin should be computed with -omega, which results in an
    3078             :       ! additional minus for the imaginary part:
    3079         308 :       n_occ = bs_env%n_occ(ispin)
    3080       13988 :       Sigma_c_n_freq(1:n_occ, :, 2) = -Sigma_c_n_freq(1:n_occ, :, 2)
    3081             : 
    3082         308 :       CALL timestop(handle)
    3083             : 
    3084         616 :    END SUBROUTINE time_to_freq
    3085             : 
    3086             : ! **************************************************************************************************
    3087             : !> \brief ...
    3088             : !> \param bs_env ...
    3089             : !> \param Sigma_c_ikp_n_freq ...
    3090             : !> \param Sigma_x_ikp_n ...
    3091             : !> \param V_xc_ikp_n ...
    3092             : !> \param eigenval_scf ...
    3093             : !> \param ikp ...
    3094             : !> \param ispin ...
    3095             : ! **************************************************************************************************
    3096         308 :    SUBROUTINE analyt_conti_and_print(bs_env, Sigma_c_ikp_n_freq, Sigma_x_ikp_n, V_xc_ikp_n, &
    3097         308 :                                      eigenval_scf, ikp, ispin)
    3098             : 
    3099             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    3100             :       REAL(KIND=dp), DIMENSION(:, :, :)                  :: Sigma_c_ikp_n_freq
    3101             :       REAL(KIND=dp), DIMENSION(:)                        :: Sigma_x_ikp_n, V_xc_ikp_n, eigenval_scf
    3102             :       INTEGER                                            :: ikp, ispin
    3103             : 
    3104             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'analyt_conti_and_print'
    3105             : 
    3106             :       CHARACTER(len=3)                                   :: occ_vir
    3107             :       CHARACTER(len=default_string_length)               :: fname
    3108             :       INTEGER                                            :: handle, i_mo, ikp_for_print, iunit, &
    3109             :                                                             n_mo, nkp
    3110             :       LOGICAL                                            :: is_bandstruc_kpoint, print_DOS_kpoints, &
    3111             :                                                             print_ikp
    3112             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dummy, Sigma_c_ikp_n_qp
    3113             : 
    3114         308 :       CALL timeset(routineN, handle)
    3115             : 
    3116         308 :       n_mo = bs_env%n_ao
    3117        1232 :       ALLOCATE (dummy(n_mo), Sigma_c_ikp_n_qp(n_mo))
    3118        3428 :       Sigma_c_ikp_n_qp(:) = 0.0_dp
    3119             : 
    3120        3428 :       DO i_mo = 1, n_mo
    3121             : 
    3122             :          ! parallelization
    3123        3120 :          IF (MODULO(i_mo, bs_env%para_env%num_pe) /= bs_env%para_env%mepos) CYCLE
    3124             : 
    3125             :          CALL continuation_pade(Sigma_c_ikp_n_qp, &
    3126             :                                 bs_env%imag_freq_points_fit, dummy, dummy, &
    3127             :                                 Sigma_c_ikp_n_freq(:, 1:bs_env%num_freq_points_fit, 1)*z_one + &
    3128             :                                 Sigma_c_ikp_n_freq(:, 1:bs_env%num_freq_points_fit, 2)*gaussi, &
    3129             :                                 Sigma_x_ikp_n(:) - V_xc_ikp_n(:), &
    3130             :                                 eigenval_scf(:), eigenval_scf(:), &
    3131             :                                 bs_env%do_hedin_shift, &
    3132             :                                 i_mo, bs_env%n_occ(ispin), bs_env%n_vir(ispin), &
    3133             :                                 bs_env%nparam_pade, bs_env%num_freq_points_fit, &
    3134             :                                 ri_rpa_g0w0_crossing_newton, bs_env%n_occ(ispin), &
    3135       97748 :                                 0.0_dp, .TRUE., .FALSE., 1, e_fermi_ext=bs_env%e_fermi(ispin))
    3136             :       END DO
    3137             : 
    3138         308 :       CALL bs_env%para_env%sum(Sigma_c_ikp_n_qp)
    3139             : 
    3140         308 :       CALL correct_obvious_fitting_fails(Sigma_c_ikp_n_qp, ispin, bs_env)
    3141             : 
    3142             :       bs_env%eigenval_G0W0(:, ikp, ispin) = eigenval_scf(:) + &
    3143             :                                             Sigma_c_ikp_n_qp(:) + &
    3144             :                                             Sigma_x_ikp_n(:) - &
    3145        3428 :                                             V_xc_ikp_n(:)
    3146             : 
    3147        3428 :       bs_env%eigenval_HF(:, ikp, ispin) = eigenval_scf(:) + Sigma_x_ikp_n(:) - V_xc_ikp_n(:)
    3148             : 
    3149             :       ! only print eigenvalues of DOS k-points in case no bandstructure path has been given
    3150         308 :       print_DOS_kpoints = (bs_env%nkp_only_bs .LE. 0)
    3151             :       ! in kpoints_DOS, the last nkp_only_bs are bandstructure k-points
    3152         308 :       is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
    3153         308 :       print_ikp = print_DOS_kpoints .OR. is_bandstruc_kpoint
    3154             : 
    3155         308 :       IF (bs_env%para_env%is_source() .AND. print_ikp) THEN
    3156             : 
    3157         122 :          IF (print_DOS_kpoints) THEN
    3158          60 :             nkp = bs_env%nkp_only_DOS
    3159          60 :             ikp_for_print = ikp
    3160             :          ELSE
    3161          62 :             nkp = bs_env%nkp_only_bs
    3162          62 :             ikp_for_print = ikp - bs_env%nkp_only_DOS
    3163             :          END IF
    3164             : 
    3165         122 :          fname = "bandstructure_SCF_and_G0W0"
    3166             : 
    3167         122 :          IF (ikp_for_print == 1) THEN
    3168             :             CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", &
    3169          22 :                            file_action="WRITE")
    3170             :          ELSE
    3171             :             CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", &
    3172         100 :                            file_action="WRITE", file_position="APPEND")
    3173             :          END IF
    3174             : 
    3175         122 :          WRITE (iunit, "(A)") " "
    3176         122 :          WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_print, "coordinate: ", &
    3177         244 :             bs_env%kpoints_DOS%xkp(:, ikp)
    3178         122 :          WRITE (iunit, "(A)") " "
    3179         122 :          WRITE (iunit, "(A5,A12,3A17,A16,A18)") "n", "k", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", &
    3180         244 :             "Σ^x_nk (eV)", "v_nk^xc (eV)", "ϵ_nk^G0W0 (eV)"
    3181         122 :          WRITE (iunit, "(A)") " "
    3182             : 
    3183        1394 :          DO i_mo = 1, n_mo
    3184        1272 :             IF (i_mo .LE. bs_env%n_occ(ispin)) occ_vir = 'occ'
    3185        1272 :             IF (i_mo > bs_env%n_occ(ispin)) occ_vir = 'vir'
    3186        1272 :             WRITE (iunit, "(I5,3A,I5,4F16.3,F17.3)") i_mo, ' (', occ_vir, ') ', ikp_for_print, &
    3187        1272 :                eigenval_scf(i_mo)*evolt, &
    3188        1272 :                Sigma_c_ikp_n_qp(i_mo)*evolt, &
    3189        1272 :                Sigma_x_ikp_n(i_mo)*evolt, &
    3190        1272 :                V_xc_ikp_n(i_mo)*evolt, &
    3191        2666 :                bs_env%eigenval_G0W0(i_mo, ikp, ispin)*evolt
    3192             :          END DO
    3193             : 
    3194         122 :          WRITE (iunit, "(A)") " "
    3195             : 
    3196         122 :          CALL close_file(iunit)
    3197             : 
    3198             :       END IF
    3199             : 
    3200         308 :       CALL timestop(handle)
    3201             : 
    3202         616 :    END SUBROUTINE analyt_conti_and_print
    3203             : 
    3204             : ! **************************************************************************************************
    3205             : !> \brief ...
    3206             : !> \param Sigma_c_ikp_n_qp ...
    3207             : !> \param ispin ...
    3208             : !> \param bs_env ...
    3209             : ! **************************************************************************************************
    3210         308 :    SUBROUTINE correct_obvious_fitting_fails(Sigma_c_ikp_n_qp, ispin, bs_env)
    3211             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Sigma_c_ikp_n_qp
    3212             :       INTEGER                                            :: ispin
    3213             :       TYPE(post_scf_bandstructure_type), POINTER         :: bs_env
    3214             : 
    3215             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'correct_obvious_fitting_fails'
    3216             : 
    3217             :       INTEGER                                            :: handle, homo, i_mo, j_mo, &
    3218             :                                                             n_levels_scissor, n_mo
    3219             :       LOGICAL                                            :: is_occ, is_vir
    3220             :       REAL(KIND=dp)                                      :: sum_Sigma_c
    3221             : 
    3222         308 :       CALL timeset(routineN, handle)
    3223             : 
    3224         308 :       n_mo = bs_env%n_ao
    3225         308 :       homo = bs_env%n_occ(ispin)
    3226             : 
    3227        3428 :       DO i_mo = 1, n_mo
    3228             : 
    3229             :          ! if |𝚺^c| > 13 eV, we use a scissors shift
    3230        3428 :          IF (ABS(Sigma_c_ikp_n_qp(i_mo)) > 13.0_dp/evolt) THEN
    3231             : 
    3232           0 :             is_occ = (i_mo .LE. homo)
    3233           0 :             is_vir = (i_mo > homo)
    3234             : 
    3235           0 :             n_levels_scissor = 0
    3236           0 :             sum_Sigma_c = 0.0_dp
    3237             : 
    3238             :             ! compute scissor
    3239           0 :             DO j_mo = 1, n_mo
    3240             : 
    3241             :                ! only compute scissor from other GW levels close in energy
    3242           0 :                IF (is_occ .AND. j_mo > homo) CYCLE
    3243           0 :                IF (is_vir .AND. j_mo .LE. homo) CYCLE
    3244           0 :                IF (ABS(i_mo - j_mo) > 10) CYCLE
    3245           0 :                IF (i_mo == j_mo) CYCLE
    3246             : 
    3247           0 :                n_levels_scissor = n_levels_scissor + 1
    3248           0 :                sum_Sigma_c = sum_Sigma_c + Sigma_c_ikp_n_qp(j_mo)
    3249             : 
    3250             :             END DO
    3251             : 
    3252             :             ! overwrite the self-energy with scissor shift
    3253           0 :             Sigma_c_ikp_n_qp(i_mo) = sum_Sigma_c/REAL(n_levels_scissor, KIND=dp)
    3254             : 
    3255             :          END IF
    3256             : 
    3257             :       END DO ! i_mo
    3258             : 
    3259         308 :       CALL timestop(handle)
    3260             : 
    3261         308 :    END SUBROUTINE correct_obvious_fitting_fails
    3262             : 
    3263             : END MODULE gw_utils

Generated by: LCOV version 1.15