LCOV - code coverage report
Current view: top level - src - gw_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:3158929) Lines: 1194 1293 92.3 %
Date: 2025-07-22 22:41:15 Functions: 44 47 93.6 %

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

Generated by: LCOV version 1.15