From 6a308192463da362cbf2e6b0309d41648d69dd4d Mon Sep 17 00:00:00 2001 From: Max Graml <78024843+MGraml@users.noreply.github.com> Date: Fri, 8 Mar 2024 20:30:58 +0100 Subject: [PATCH] Bethe Salpeter equation for molecules (#3308) --- src/CMakeLists.txt | 5 +- src/bse.F | 646 ----------- src/bse_full_diag.F | 712 ++++++++++++ src/bse_iterative.F | 1587 ++++++++++++++++++++++++++ src/bse_main.F | 221 ++++ src/bse_util.F | 1351 ++++++++++++++++++++++ src/input_constants.F | 17 + src/input_cp2k_mp2.F | 288 +++-- src/mp2_integrals.F | 10 +- src/mp2_setup.F | 50 +- src/mp2_types.F | 24 +- src/rpa_gw.F | 11 +- src/rpa_main.F | 66 +- src/rpa_util.F | 10 +- tests/QS/regtest-bse/BSE_H2O_PBE.inp | 80 ++ tests/QS/regtest-bse/H2O_gas.xyz | 5 + tests/QS/regtest-bse/TDA_H2O_PBE.inp | 80 ++ tests/QS/regtest-bse/TEST_FILES | 3 + tests/TEST_DIRS | 1 + tests/TEST_TYPES | 4 +- 20 files changed, 4398 insertions(+), 773 deletions(-) delete mode 100644 src/bse.F create mode 100644 src/bse_full_diag.F create mode 100644 src/bse_iterative.F create mode 100644 src/bse_main.F create mode 100644 src/bse_util.F create mode 100644 tests/QS/regtest-bse/BSE_H2O_PBE.inp create mode 100644 tests/QS/regtest-bse/H2O_gas.xyz create mode 100644 tests/QS/regtest-bse/TDA_H2O_PBE.inp create mode 100644 tests/QS/regtest-bse/TEST_FILES diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4f0a0f910a..828f2b7335 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -69,7 +69,10 @@ list( basis_set_output.F beta_gamma_psi.F block_p_types.F - bse.F + bse_main.F + bse_full_diag.F + bse_iterative.F + bse_util.F bsse.F cell_methods.F colvar_methods.F diff --git a/src/bse.F b/src/bse.F deleted file mode 100644 index 225725ef20..0000000000 --- a/src/bse.F +++ /dev/null @@ -1,646 +0,0 @@ -!--------------------------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright 2000-2024 CP2K developers group ! -! ! -! SPDX-License-Identifier: GPL-2.0-or-later ! -!--------------------------------------------------------------------------------------------------! - -! ************************************************************************************************** -!> \brief Routines for GW + Bethe-Salpeter for computing electronic excitations -!> \par History -!> 04.2017 created [Jan Wilhelm] -! ************************************************************************************************** -MODULE bse - USE cp_fm_basic_linalg, ONLY: cp_fm_upper_to_full - USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& - cp_fm_cholesky_invert - USE cp_fm_types, ONLY: cp_fm_create,& - cp_fm_get_info,& - cp_fm_release,& - cp_fm_set_all,& - cp_fm_to_fm,& - cp_fm_type - USE group_dist_types, ONLY: get_group_dist,& - group_dist_d1_type - USE kinds, ONLY: dp - USE message_passing, ONLY: mp_para_env_type,& - mp_request_type - USE mp2_types, ONLY: integ_mat_buffer_type - USE parallel_gemm_api, ONLY: parallel_gemm - USE rpa_communication, ONLY: communicate_buffer -#include "./base/base_uses.f90" - - IMPLICIT NONE - - PRIVATE - - CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse' - - PUBLIC :: mult_B_with_W_and_fill_local_3c_arrays, do_subspace_iterations - -CONTAINS - -! ************************************************************************************************** -!> \brief ... -!> \param B_bar_ijQ_bse_local ... -!> \param B_abQ_bse_local ... -!> \param B_bar_iaQ_bse_local ... -!> \param B_iaQ_bse_local ... -!> \param homo ... -!> \param virtual ... -!> \param num_Z_vectors ... -!> \param max_iter ... -!> \param threshold_min_trans ... -!> \param Eigenval ... -!> \param para_env ... -! ************************************************************************************************** - SUBROUTINE do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & - B_iaQ_bse_local, homo, virtual, num_Z_vectors, & - max_iter, threshold_min_trans, Eigenval, para_env) - - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_bar_ijQ_bse_local, B_abQ_bse_local, & - B_bar_iaQ_bse_local, B_iaQ_bse_local - INTEGER :: homo, virtual, num_Z_vectors, max_iter - REAL(KIND=dp) :: threshold_min_trans - REAL(KIND=dp), DIMENSION(:) :: Eigenval - TYPE(mp_para_env_type), INTENT(IN) :: para_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'do_subspace_iterations' - - INTEGER :: handle, i_iter, local_RI_size - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_ia_tmp, M_ji_tmp, RI_vector - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, BZ, Z_vectors - - CALL timeset(routineN, handle) - - ! JW hack 2del - threshold_min_trans = 0.01_dp - - ALLOCATE (Z_vectors(homo, virtual, num_Z_vectors)) - Z_vectors = 0.0_dp - - ALLOCATE (AZ(homo, virtual, num_Z_vectors)) - AZ = 0.0_dp - - ALLOCATE (BZ(homo, virtual, num_Z_vectors)) - BZ = 0.0_dp - - local_RI_size = SIZE(B_iaQ_bse_local, 3) - - ALLOCATE (M_ia_tmp(homo, virtual)) - M_ia_tmp = 0.0_dp - - ALLOCATE (M_ji_tmp(homo, homo)) - M_ji_tmp = 0.0_dp - - ALLOCATE (RI_vector(local_RI_size, num_Z_vectors)) - RI_vector = 0.0_dp - - CALL initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual) - - DO i_iter = 1, max_iter - - CALL compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, & - M_ia_tmp, RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, & - para_env) - - CALL compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, & - M_ji_tmp, RI_vector, homo, virtual, num_Z_vectors, local_RI_size, & - para_env) - - END DO - - DEALLOCATE (AZ, BZ, Z_vectors, M_ia_tmp, M_ji_tmp, RI_vector) - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param BZ ... -!> \param Z_vectors ... -!> \param B_iaQ_bse_local ... -!> \param B_bar_iaQ_bse_local ... -!> \param M_ji_tmp ... -!> \param RI_vector ... -!> \param homo ... -!> \param virtual ... -!> \param num_Z_vectors ... -!> \param local_RI_size ... -!> \param para_env ... -! ************************************************************************************************** - SUBROUTINE compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, & - M_ji_tmp, RI_vector, homo, virtual, num_Z_vectors, local_RI_size, para_env) - - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BZ, Z_vectors, B_iaQ_bse_local, & - B_bar_iaQ_bse_local - REAL(KIND=dp), DIMENSION(:, :) :: M_ji_tmp, RI_vector - INTEGER :: homo, virtual, num_Z_vectors, & - local_RI_size - TYPE(mp_para_env_type), INTENT(IN) :: para_env - - INTEGER :: i_Z_vector, LLL - - BZ(:, :, :) = 0.0_dp - - CALL compute_v_ia_jb_part(BZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, & - num_Z_vectors, homo, virtual) - - DO i_Z_vector = 1, num_Z_vectors - - DO LLL = 1, local_RI_size - - ! M_ji^P = sum_b Z_jb*B_bi^P - CALL DGEMM("N", "T", homo, homo, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, & - B_iaQ_bse_local(:, :, LLL), homo, 0.0_dp, M_ji_tmp, homo) - ! (BZ)_ia = sum_jP M_ij^P*B^bar_ja^P - CALL DGEMM("T", "N", homo, virtual, homo, 1.0_dp, M_ji_tmp, homo, & - B_bar_iaQ_bse_local, homo, 1.0_dp, BZ(:, :, i_Z_vector), homo) - - END DO - - END DO - - ! we make the sum over all RI basis functions - CALL para_env%sum(BZ) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param AZ ... -!> \param Z_vectors ... -!> \param B_iaQ_bse_local ... -!> \param B_bar_ijQ_bse_local ... -!> \param B_abQ_bse_local ... -!> \param M_ia_tmp ... -!> \param RI_vector ... -!> \param Eigenval ... -!> \param homo ... -!> \param virtual ... -!> \param num_Z_vectors ... -!> \param local_RI_size ... -!> \param para_env ... -! ************************************************************************************************** - SUBROUTINE compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, M_ia_tmp, & - RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, para_env) - - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, Z_vectors, B_iaQ_bse_local, & - B_bar_ijQ_bse_local, B_abQ_bse_local - REAL(KIND=dp), DIMENSION(:, :) :: M_ia_tmp, RI_vector - REAL(KIND=dp), DIMENSION(:) :: Eigenval - INTEGER :: homo, virtual, num_Z_vectors, & - local_RI_size - TYPE(mp_para_env_type), INTENT(IN) :: para_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_AZ' - - INTEGER :: a_virt, handle, i_occ, i_Z_vector, LLL - REAL(KIND=dp) :: eigen_diff - - CALL timeset(routineN, handle) - - AZ(:, :, :) = 0.0_dp - - CALL compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, & - num_Z_vectors, homo, virtual) - - DO i_Z_vector = 1, num_Z_vectors - - ! JW TO DO: OMP PARALLELIZATION - DO LLL = 1, local_RI_size - - ! M_ja^P = sum_j Z_jb*B_ba^P - CALL DGEMM("N", "N", homo, virtual, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, & - B_abQ_bse_local(:, :, LLL), virtual, 0.0_dp, M_ia_tmp, homo) - - ! (AZ)_ia = sum_jP B_bar_ij^P*M_ja^P - CALL DGEMM("N", "N", homo, virtual, homo, 1.0_dp, B_bar_ijQ_bse_local(:, :, LLL), homo, & - M_ia_tmp, homo, 1.0_dp, AZ(:, :, i_Z_vector), homo) - - END DO - - END DO - - ! we make the sum over all RI basis functions - CALL para_env%sum(AZ) - - ! add (e_a-e_i)*Z_ia - DO i_occ = 1, homo - DO a_virt = 1, virtual - - eigen_diff = Eigenval(a_virt + homo) - Eigenval(i_occ) - - AZ(i_occ, a_virt, :) = AZ(i_occ, a_virt, :) + Z_vectors(i_occ, a_virt, :)*eigen_diff - - END DO - END DO - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param AZ ... -!> \param Z_vectors ... -!> \param B_iaQ_bse_local ... -!> \param RI_vector ... -!> \param local_RI_size ... -!> \param num_Z_vectors ... -!> \param homo ... -!> \param virtual ... -! ************************************************************************************************** - SUBROUTINE compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, & - num_Z_vectors, homo, virtual) - - REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: AZ, Z_vectors, B_iaQ_bse_local - REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: RI_vector - INTEGER, INTENT(IN) :: local_RI_size, num_Z_vectors, homo, & - virtual - - CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_ia_jb_part' - - INTEGER :: a_virt, handle, i_occ, i_Z_vector, LLL - - CALL timeset(routineN, handle) - - RI_vector = 0.0_dp - - ! v_P = sum_jb B_jb^P Z_jb - DO LLL = 1, local_RI_size - DO i_Z_vector = 1, num_Z_vectors - DO i_occ = 1, homo - DO a_virt = 1, virtual - - RI_vector(LLL, i_Z_vector) = RI_vector(LLL, i_Z_vector) + & - Z_vectors(i_occ, a_virt, i_Z_vector)* & - B_iaQ_bse_local(i_occ, a_virt, LLL) - - END DO - END DO - END DO - END DO - - ! AZ = sum_P B_ia^P*v_P + ... - DO LLL = 1, local_RI_size - DO i_Z_vector = 1, num_Z_vectors - DO i_occ = 1, homo - DO a_virt = 1, virtual - - AZ(i_occ, a_virt, i_Z_vector) = AZ(i_occ, a_virt, i_Z_vector) + & - RI_vector(LLL, i_Z_vector)* & - B_iaQ_bse_local(i_occ, a_virt, LLL) - - END DO - END DO - END DO - END DO - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param Z_vectors ... -!> \param Eigenval ... -!> \param num_Z_vectors ... -!> \param homo ... -!> \param virtual ... -! ************************************************************************************************** - SUBROUTINE initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual) - - REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: Z_vectors - REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval - INTEGER, INTENT(IN) :: num_Z_vectors, homo, virtual - - CHARACTER(LEN=*), PARAMETER :: routineN = 'initial_guess_Z_vectors' - - INTEGER :: a_virt, handle, i_occ, i_Z_vector, & - min_loc(2) - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigen_diff_ia - - CALL timeset(routineN, handle) - - ALLOCATE (eigen_diff_ia(homo, virtual)) - - DO i_occ = 1, homo - DO a_virt = 1, virtual - eigen_diff_ia(i_occ, a_virt) = Eigenval(a_virt + homo) - Eigenval(i_occ) - END DO - END DO - - DO i_Z_vector = 1, num_Z_vectors - - min_loc = MINLOC(eigen_diff_ia) - - Z_vectors(min_loc(1), min_loc(2), i_Z_vector) = 1.0_dp - - eigen_diff_ia(min_loc(1), min_loc(2)) = 1.0E20_dp - - END DO - - DEALLOCATE (eigen_diff_ia) - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param fm_mat_S_ij_bse ... -!> \param fm_mat_S_ab_bse ... -!> \param fm_mat_S ... -!> \param fm_mat_Q_static_bse ... -!> \param fm_mat_Q_static_bse_gemm ... -!> \param B_bar_ijQ_bse_local ... -!> \param B_abQ_bse_local ... -!> \param B_bar_iaQ_bse_local ... -!> \param B_iaQ_bse_local ... -!> \param dimen_RI ... -!> \param homo ... -!> \param virtual ... -!> \param dimen_ia ... -!> \param gd_array ... -!> \param color_sub ... -!> \param para_env ... -! ************************************************************************************************** - SUBROUTINE mult_B_with_W_and_fill_local_3c_arrays(fm_mat_S_ij_bse, fm_mat_S_ab_bse, fm_mat_S, fm_mat_Q_static_bse, & - fm_mat_Q_static_bse_gemm, & - B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & - B_iaQ_bse_local, dimen_RI, homo, virtual, dimen_ia, & - gd_array, color_sub, para_env) - - TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ij_bse, fm_mat_S_ab_bse, & - fm_mat_S, fm_mat_Q_static_bse, & - fm_mat_Q_static_bse_gemm - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & - INTENT(OUT) :: B_bar_ijQ_bse_local, B_abQ_bse_local, & - B_bar_iaQ_bse_local, B_iaQ_bse_local - INTEGER, INTENT(IN) :: dimen_RI, homo, virtual, dimen_ia - TYPE(group_dist_d1_type), INTENT(IN) :: gd_array - INTEGER, INTENT(IN) :: color_sub - TYPE(mp_para_env_type), INTENT(IN) :: para_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'mult_B_with_W_and_fill_local_3c_arrays' - - INTEGER :: handle, i_global, iiB, info_chol, & - j_global, jjB, ncol_local, nrow_local - INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices - TYPE(cp_fm_type) :: fm_mat_S_bar_ia_bse, & - fm_mat_S_bar_ij_bse, fm_mat_work - - CALL timeset(routineN, handle) - - CALL cp_fm_create(fm_mat_S_bar_ia_bse, fm_mat_S%matrix_struct) - CALL cp_fm_to_fm(fm_mat_S, fm_mat_S_bar_ia_bse) - CALL cp_fm_set_all(fm_mat_S_bar_ia_bse, 0.0_dp) - - CALL cp_fm_create(fm_mat_S_bar_ij_bse, fm_mat_S_ij_bse%matrix_struct) - CALL cp_fm_to_fm(fm_mat_S_ij_bse, fm_mat_S_bar_ij_bse) - CALL cp_fm_set_all(fm_mat_S_bar_ij_bse, 0.0_dp) - - CALL cp_fm_create(fm_mat_work, fm_mat_Q_static_bse_gemm%matrix_struct) - CALL cp_fm_to_fm(fm_mat_Q_static_bse_gemm, fm_mat_work) - CALL cp_fm_set_all(fm_mat_work, 0.0_dp) - - ! get info of fm_mat_Q_static_bse and compute ((1+Q(0))^-1-1) - CALL cp_fm_get_info(matrix=fm_mat_Q_static_bse_gemm, & - nrow_local=nrow_local, & - ncol_local=ncol_local, & - row_indices=row_indices, & - col_indices=col_indices) - - DO jjB = 1, ncol_local - j_global = col_indices(jjB) - DO iiB = 1, nrow_local - i_global = row_indices(iiB) - IF (j_global == i_global .AND. i_global <= dimen_RI) THEN - fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) = fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) + 1.0_dp - END IF - END DO - END DO - - ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition - CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q_static_bse_gemm, n=dimen_RI, info_out=info_chol) - CPASSERT(info_chol == 0) - - ! calculate [1+Q(i0)]^-1 - CALL cp_fm_cholesky_invert(fm_mat_Q_static_bse_gemm) - ! symmetrize the result - CALL cp_fm_upper_to_full(fm_mat_Q_static_bse_gemm, fm_mat_work) - - CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=homo**2, k=dimen_RI, alpha=1.0_dp, & - matrix_a=fm_mat_Q_static_bse, matrix_b=fm_mat_S_ij_bse, beta=0.0_dp, & - matrix_c=fm_mat_S_bar_ij_bse) - - ! fm_mat_S_bar_ia_bse has a different blacs_env as fm_mat_S_ij_bse since we take - ! fm_mat_S from RPA. Therefore, we also need a different fm_mat_Q_static_bse_gemm - CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_ia, k=dimen_RI, alpha=1.0_dp, & - matrix_a=fm_mat_Q_static_bse_gemm, matrix_b=fm_mat_S, beta=0.0_dp, & - matrix_c=fm_mat_S_bar_ia_bse) - - CALL allocate_and_fill_local_array(B_iaQ_bse_local, fm_mat_S, gd_array, color_sub, homo, virtual, dimen_RI, para_env) - - CALL allocate_and_fill_local_array(B_bar_iaQ_bse_local, fm_mat_S_bar_ia_bse, gd_array, color_sub, homo, virtual, & - dimen_RI, para_env) - - CALL allocate_and_fill_local_array(B_bar_ijQ_bse_local, fm_mat_S_bar_ij_bse, gd_array, color_sub, homo, homo, & - dimen_RI, para_env) - - CALL allocate_and_fill_local_array(B_abQ_bse_local, fm_mat_S_ab_bse, gd_array, color_sub, virtual, virtual, & - dimen_RI, para_env) - - CALL cp_fm_release(fm_mat_S_bar_ia_bse) - CALL cp_fm_release(fm_mat_S_bar_ij_bse) - CALL cp_fm_release(fm_mat_work) - - CALL timestop(handle) - - END SUBROUTINE - -! ************************************************************************************************** -!> \brief ... -!> \param B_local ... -!> \param fm_mat_S ... -!> \param gd_array ... -!> \param color_sub ... -!> \param small_size ... -!> \param big_size ... -!> \param dimen_RI ... -!> \param para_env ... -! ************************************************************************************************** - SUBROUTINE allocate_and_fill_local_array(B_local, fm_mat_S, gd_array, & - color_sub, small_size, big_size, dimen_RI, para_env) - - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & - INTENT(OUT) :: B_local - TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S - TYPE(group_dist_d1_type), INTENT(IN) :: gd_array - INTEGER, INTENT(IN) :: color_sub, small_size, big_size, dimen_RI - TYPE(mp_para_env_type), INTENT(IN) :: para_env - - CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_and_fill_local_array' - - INTEGER :: combi_index, end_RI, handle, handle1, i_comm, i_entry, iiB, imepos, jjB, & - level_big_size, level_small_size, ncol_local, nrow_local, num_comm_cycles, RI_index, & - size_RI, start_RI - INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, mepos_from_RI_index, & - num_entries_rec, num_entries_send - INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices - REAL(KIND=dp) :: matrix_el - TYPE(integ_mat_buffer_type), ALLOCATABLE, & - DIMENSION(:) :: buffer_rec, buffer_send - TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array - - CALL timeset(routineN, handle) - - ALLOCATE (mepos_from_RI_index(dimen_RI)) - mepos_from_RI_index = 0 - - DO imepos = 0, para_env%num_pe - 1 - - CALL get_group_dist(gd_array, pos=imepos, starts=start_RI, ends=end_RI) - - mepos_from_RI_index(start_RI:end_RI) = imepos - - END DO - - ! color_sub is automatically the number of the process since every subgroup has only one MPI rank - CALL get_group_dist(gd_array, color_sub, start_RI, end_RI, size_RI) - - ALLOCATE (B_local(small_size, big_size, 1:size_RI)) - - ALLOCATE (num_entries_send(0:para_env%num_pe - 1)) - ALLOCATE (num_entries_rec(0:para_env%num_pe - 1)) - - ALLOCATE (req_array(1:para_env%num_pe, 4)) - - ALLOCATE (entry_counter(0:para_env%num_pe - 1)) - - CALL cp_fm_get_info(matrix=fm_mat_S, & - nrow_local=nrow_local, & - ncol_local=ncol_local, & - row_indices=row_indices, & - col_indices=col_indices) - - num_comm_cycles = 10 - - ! communicate not all due to huge memory overhead, since for every number in fm_mat_S, we store - ! three additional ones (RI index, first MO index, second MO index!!) - DO i_comm = 0, num_comm_cycles - 1 - - num_entries_send = 0 - num_entries_rec = 0 - - ! loop over RI index to get the number of sent entries - DO jjB = 1, nrow_local - - RI_index = row_indices(jjB) - - IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE - - imepos = mepos_from_RI_index(RI_index) - - num_entries_send(imepos) = num_entries_send(imepos) + ncol_local - - END DO - - CALL para_env%alltoall(num_entries_send, num_entries_rec, 1) - - ALLOCATE (buffer_rec(0:para_env%num_pe - 1)) - ALLOCATE (buffer_send(0:para_env%num_pe - 1)) - - ! allocate data message and corresponding indices - DO imepos = 0, para_env%num_pe - 1 - - ALLOCATE (buffer_rec(imepos)%msg(num_entries_rec(imepos))) - buffer_rec(imepos)%msg = 0.0_dp - - ALLOCATE (buffer_send(imepos)%msg(num_entries_send(imepos))) - buffer_send(imepos)%msg = 0.0_dp - - ALLOCATE (buffer_rec(imepos)%indx(num_entries_rec(imepos), 3)) - buffer_rec(imepos)%indx = 0 - - ALLOCATE (buffer_send(imepos)%indx(num_entries_send(imepos), 3)) - buffer_send(imepos)%indx = 0 - - END DO - - entry_counter(:) = 0 - - ! loop over RI index for filling the send-buffer - DO jjB = 1, nrow_local - - RI_index = row_indices(jjB) - - IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE - - imepos = mepos_from_RI_index(RI_index) - - DO iiB = 1, ncol_local - - combi_index = col_indices(iiB) - level_small_size = MAX(1, combi_index - 1)/big_size + 1 - level_big_size = combi_index - (level_small_size - 1)*big_size - - entry_counter(imepos) = entry_counter(imepos) + 1 - - buffer_send(imepos)%msg(entry_counter(imepos)) = fm_mat_S%local_data(jjB, iiB) - - buffer_send(imepos)%indx(entry_counter(imepos), 1) = RI_index - buffer_send(imepos)%indx(entry_counter(imepos), 2) = level_small_size - buffer_send(imepos)%indx(entry_counter(imepos), 3) = level_big_size - - END DO - - END DO - - CALL timeset("BSE_comm_data", handle1) - - CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array) - - CALL timestop(handle1) - - ! fill B_local - DO imepos = 0, para_env%num_pe - 1 - - DO i_entry = 1, num_entries_rec(imepos) - - RI_index = buffer_rec(imepos)%indx(i_entry, 1) - start_RI + 1 - level_small_size = buffer_rec(imepos)%indx(i_entry, 2) - level_big_size = buffer_rec(imepos)%indx(i_entry, 3) - - matrix_el = buffer_rec(imepos)%msg(i_entry) - - B_local(level_small_size, level_big_size, RI_index) = matrix_el - - END DO - - END DO - - DO imepos = 0, para_env%num_pe - 1 - DEALLOCATE (buffer_send(imepos)%msg) - DEALLOCATE (buffer_send(imepos)%indx) - DEALLOCATE (buffer_rec(imepos)%msg) - DEALLOCATE (buffer_rec(imepos)%indx) - END DO - - DEALLOCATE (buffer_rec, buffer_send) - - END DO - - DEALLOCATE (num_entries_send, num_entries_rec) - - DEALLOCATE (mepos_from_RI_index) - - DEALLOCATE (entry_counter, req_array) - - CALL timestop(handle) - - END SUBROUTINE - -END MODULE bse diff --git a/src/bse_full_diag.F b/src/bse_full_diag.F new file mode 100644 index 0000000000..e4c2d26447 --- /dev/null +++ b/src/bse_full_diag.F @@ -0,0 +1,712 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Routines for the full diagonalization of GW + Bethe-Salpeter for computing +!> electronic excitations +!> \par History +!> 10.2023 created [Maximilian Graml] +! ************************************************************************************************** +MODULE bse_full_diag + + USE bse_util, ONLY: comp_eigvec_coeff_BSE,& + filter_eigvec_contrib,& + fm_general_add_bse + USE cp_blacs_env, ONLY: cp_blacs_env_create,& + cp_blacs_env_release,& + cp_blacs_env_type + USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add + USE cp_fm_diag, ONLY: choose_eigv_solver,& + cp_fm_power + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type + USE cp_fm_types, ONLY: cp_fm_create,& + cp_fm_get_info,& + cp_fm_release,& + cp_fm_set_all,& + cp_fm_to_fm,& + cp_fm_type + USE input_constants, ONLY: bse_singlet,& + bse_triplet + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type + USE mp2_types, ONLY: mp2_type + USE parallel_gemm_api, ONLY: parallel_gemm + USE physcon, ONLY: evolt +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_full_diag' + + PUBLIC :: create_A, diagonalize_A, create_B, create_hermitian_form_of_ABBA, & + diagonalize_C + +CONTAINS + +! ************************************************************************************************** +!> \brief Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W) +!> A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab +!> ε_a, ε_i are GW singleparticle energies from Eigenval_reduced +!> α is a spin-dependent factor +!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction) +!> W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction) +!> \param fm_mat_S_ia_bse ... +!> \param fm_mat_S_bar_ij_bse ... +!> \param fm_mat_S_ab_bse ... +!> \param fm_A ... +!> \param Eigenval ... +!> \param unit_nr ... +!> \param homo ... +!> \param virtual ... +!> \param dimen_RI ... +!> \param mp2_env ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, & + fm_A, Eigenval, unit_nr, & + homo, virtual, dimen_RI, mp2_env, & + para_env) + + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, & + fm_mat_S_ab_bse + TYPE(cp_fm_type), INTENT(INOUT) :: fm_A + REAL(KIND=dp), DIMENSION(:) :: Eigenval + INTEGER, INTENT(IN) :: unit_nr, homo, virtual, dimen_RI + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + TYPE(mp_para_env_type), INTENT(INOUT) :: para_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'create_A' + + INTEGER :: a_virt_row, handle, i_occ_row, & + i_row_global, ii, j_col_global, jj, & + ncol_local_A, nrow_local_A + INTEGER, DIMENSION(4) :: reordering + INTEGER, DIMENSION(:), POINTER :: col_indices_A, row_indices_A + REAL(KIND=dp) :: alpha, eigen_diff + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_fm_struct_type), POINTER :: fm_struct_A, fm_struct_W + TYPE(cp_fm_type) :: fm_W + + CALL timeset(routineN, handle) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A' + END IF + + !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F) + SELECT CASE (mp2_env%ri_g0w0%bse_spin_config) + CASE (bse_singlet) + alpha = 2.0_dp + CASE (bse_triplet) + alpha = 0.0_dp + END SELECT + + ! create the blacs env for ij matrices (NOT fm_mat_S_ia_bse%matrix_struct related parallel_gemms!) + NULLIFY (blacs_env) + CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env) + + ! We have to use the same blacs_env for A as for the matrices fm_mat_S_ia_bse from RPA + ! Logic: A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab + ! We create v_ia,jb and W_ij,ab, then we communicate entries from local W_ij,ab + ! to the full matrix v_ia,jb. By adding these and the energy diffenences: v_ia,jb -> A_ia,jb + ! We use the A matrix already from the start instead of v + CALL cp_fm_struct_create(fm_struct_A, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, & + ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env) + CALL cp_fm_create(fm_A, fm_struct_A, name="fm_A_iajb") + CALL cp_fm_set_all(fm_A, 0.0_dp) + + CALL cp_fm_struct_create(fm_struct_W, context=fm_mat_S_ab_bse%matrix_struct%context, nrow_global=homo**2, & + ncol_global=virtual**2, para_env=fm_mat_S_ab_bse%matrix_struct%para_env) + CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ijab") + CALL cp_fm_set_all(fm_W, 0.0_dp) + + ! Create A matrix from GW Energies, v_ia,jb and W_ij,ab (different blacs_env!) + ! v_ia,jb, which is directly initialized in A (with a factor of alpha) + ! v_ia,jb = \sum_P B^P_ia B^P_jb + CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, & + matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, & + matrix_c=fm_A) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb' + END IF + + !W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab + CALL parallel_gemm(transa="T", transb="N", m=homo**2, n=virtual**2, k=dimen_RI, alpha=1.0_dp, & + matrix_a=fm_mat_S_bar_ij_bse, matrix_b=fm_mat_S_ab_bse, beta=0.0_dp, & + matrix_c=fm_W) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab' + END IF + + ! We start by moving data from local parts of W_ij,ab to the full matrix A_ia,jb using buffers + CALL cp_fm_get_info(matrix=fm_A, & + nrow_local=nrow_local_A, & + ncol_local=ncol_local_A, & + row_indices=row_indices_A, & + col_indices=col_indices_A) + ! Writing -1.0_dp * W_ij,ab to A_ia,jb, i.e. beta = -1.0_dp, + ! W_ij,ab: nrow_secidx_in = homo, ncol_secidx_in = virtual + ! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual + reordering = (/1, 3, 2, 4/) + CALL fm_general_add_bse(fm_A, fm_W, -1.0_dp, homo, virtual, & + virtual, virtual, unit_nr, reordering, mp2_env) + !full matrix W is not needed anymore, release it to save memory + CALL cp_fm_struct_release(fm_struct_W) + CALL cp_fm_release(fm_W) + CALL cp_fm_struct_release(fm_struct_A) + + !Now add the energy differences (ε_a-ε_i) on the diagonal (i.e. δ_ij δ_ab) of A_ia,jb + DO ii = 1, nrow_local_A + + i_row_global = row_indices_A(ii) + + DO jj = 1, ncol_local_A + + j_col_global = col_indices_A(jj) + + IF (i_row_global == j_col_global) THEN + i_occ_row = (i_row_global - 1)/virtual + 1 + a_virt_row = MOD(i_row_global - 1, virtual) + 1 + eigen_diff = Eigenval(a_virt_row + homo) - Eigenval(i_occ_row) + fm_A%local_data(ii, jj) = fm_A%local_data(ii, jj) + eigen_diff + + END IF + END DO + END DO + + CALL cp_fm_struct_release(fm_struct_A) + CALL cp_fm_struct_release(fm_struct_W) + + CALL cp_blacs_env_release(blacs_env) + + CALL timestop(handle) + + END SUBROUTINE create_A + +! ************************************************************************************************** +!> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W) +!> B_ia,jb = α * v_ia,jb - W_ib,aj +!> α is a spin-dependent factor +!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction) +!> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction) +!> \param fm_mat_S_ia_bse ... +!> \param fm_mat_S_bar_ia_bse ... +!> \param fm_B ... +!> \param homo ... +!> \param virtual ... +!> \param dimen_RI ... +!> \param unit_nr ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, & + homo, virtual, dimen_RI, unit_nr, mp2_env) + + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse + TYPE(cp_fm_type), INTENT(INOUT) :: fm_B + INTEGER, INTENT(IN) :: homo, virtual, dimen_RI, unit_nr + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'create_B' + + INTEGER :: handle + INTEGER, DIMENSION(4) :: reordering + REAL(KIND=dp) :: alpha + TYPE(cp_fm_struct_type), POINTER :: fm_struct_v + TYPE(cp_fm_type) :: fm_W + + CALL timeset(routineN, handle) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B' + END IF + + ! Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F) + SELECT CASE (mp2_env%ri_g0w0%bse_spin_config) + CASE (bse_singlet) + alpha = 2.0_dp + CASE (bse_triplet) + alpha = 0.0_dp + END SELECT + + CALL cp_fm_struct_create(fm_struct_v, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, & + ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env) + CALL cp_fm_create(fm_B, fm_struct_v, name="fm_B_iajb") + CALL cp_fm_set_all(fm_B, 0.0_dp) + + CALL cp_fm_create(fm_W, fm_struct_v, name="fm_W_ibaj") + CALL cp_fm_set_all(fm_W, 0.0_dp) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb' + END IF + ! v_ia,jb = \sum_P B^P_ia B^P_jb + CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, & + matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, & + matrix_c=fm_B) + + ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj + CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=1.0_dp, & + matrix_a=fm_mat_S_bar_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, & + matrix_c=fm_W) + ! from W_ib,ja to A_ia,jb (formally: W_ib,aj, but our internal indexorder is different) + ! Writing -1.0_dp * W_ib,ja to A_ia,jb, i.e. beta = -1.0_dp, + ! W_ib,ja: nrow_secidx_in = virtual, ncol_secidx_in = virtual + ! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual + reordering = (/1, 4, 3, 2/) + CALL fm_general_add_bse(fm_B, fm_W, -1.0_dp, virtual, virtual, & + virtual, virtual, unit_nr, reordering, mp2_env) + + CALL cp_fm_release(fm_W) + CALL cp_fm_struct_release(fm_struct_v) + CALL timestop(handle) + + END SUBROUTINE create_B + + ! ************************************************************************************************** +!> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem +!> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)). +!> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions +!> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001). +!> \param fm_A ... +!> \param fm_B ... +!> \param fm_C ... +!> \param fm_sqrt_A_minus_B ... +!> \param fm_inv_sqrt_A_minus_B ... +!> \param homo ... +!> \param virtual ... +!> \param unit_nr ... +!> \param mp2_env ... +!> \param diag_est ... +! ************************************************************************************************** + SUBROUTINE create_hermitian_form_of_ABBA(fm_A, fm_B, fm_C, & + fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, & + homo, virtual, unit_nr, mp2_env, diag_est) + + TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B + TYPE(cp_fm_type), INTENT(INOUT) :: fm_C, fm_sqrt_A_minus_B, & + fm_inv_sqrt_A_minus_B + INTEGER, INTENT(IN) :: homo, virtual, unit_nr + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + REAL(KIND=dp), INTENT(IN) :: diag_est + + CHARACTER(LEN=*), PARAMETER :: routineN = 'create_hermitian_form_of_ABBA' + + INTEGER :: dim_mat, handle, n_dependent + REAL(KIND=dp), DIMENSION(2) :: eigvals_AB_diff + TYPE(cp_fm_type) :: fm_A_minus_B, fm_A_plus_B, fm_dummy, & + fm_work_product + + CALL timeset(routineN, handle) + + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', & + ' with size of A. This will take around ', diag_est, " s." + END IF + + ! Create work matrices, which will hold A+B and A-B and their powers + ! C is created afterwards to save memory + ! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I + ! \_______/ \___/ \______/ + ! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B + ! (EQ.Ia) (EQ.Ib) (EQ.Ia) + ! Intermediate work matrices: + ! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II + ! fm_A_minus_B: (A-B) EQ.III + ! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV + CALL cp_fm_create(fm_A_plus_B, fm_A%matrix_struct) + CALL cp_fm_to_fm(fm_A, fm_A_plus_B) + CALL cp_fm_create(fm_A_minus_B, fm_A%matrix_struct) + CALL cp_fm_to_fm(fm_A, fm_A_minus_B) + CALL cp_fm_create(fm_sqrt_A_minus_B, fm_A%matrix_struct) + CALL cp_fm_set_all(fm_sqrt_A_minus_B, 0.0_dp) + CALL cp_fm_create(fm_inv_sqrt_A_minus_B, fm_A%matrix_struct) + CALL cp_fm_set_all(fm_inv_sqrt_A_minus_B, 0.0_dp) + + CALL cp_fm_create(fm_work_product, fm_A%matrix_struct) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays' + END IF + + ! Add/Substract B (cf. EQs. Ib and III) + CALL cp_fm_scale_and_add(1.0_dp, fm_A_plus_B, 1.0_dp, fm_B) + CALL cp_fm_scale_and_add(1.0_dp, fm_A_minus_B, -1.0_dp, fm_B) + + ! cp_fm_power will overwrite matrix, therefore we create copies + CALL cp_fm_to_fm(fm_A_minus_B, fm_inv_sqrt_A_minus_B) + + ! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia) + ! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards. + + ! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp + CALL cp_fm_create(fm_dummy, fm_A%matrix_struct) + ! Create (A-B)^-0.5 (cf. EQ.II) + CALL cp_fm_power(fm_inv_sqrt_A_minus_B, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_AB_diff) + CALL cp_fm_release(fm_dummy) + ! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues) + ! In this case, the procedure for hermitian form of ABBA is not applicable + IF (eigvals_AB_diff(1) < 0) THEN + CALL cp_abort(__LOCATION__, & + "Matrix (A-B) is not positive definite. "// & + "Hermitian diagonalization of full BSE matrix is ill-defined.") + END IF + + ! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA + ! We further create (A-B)^0.5 for the singleparticle transitions of ABBA + ! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia) + dim_mat = homo*virtual + CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_inv_sqrt_A_minus_B, fm_A_minus_B, 0.0_dp, & + fm_sqrt_A_minus_B) + + ! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV) + CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_sqrt_A_minus_B, fm_A_plus_B, 0.0_dp, & + fm_work_product) + + ! Release to save memory + CALL cp_fm_release(fm_A_plus_B) + CALL cp_fm_release(fm_A_minus_B) + + ! Now create full + CALL cp_fm_create(fm_C, fm_A%matrix_struct) + CALL cp_fm_set_all(fm_C, 0.0_dp) + ! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I) + CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_work_product, fm_sqrt_A_minus_B, 0.0_dp, & + fm_C) + CALL cp_fm_release(fm_work_product) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5' + END IF + + CALL timestop(handle) + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n . +!> Here, the eigenvectors Z^n relate to X^n via +!> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001). +!> \param fm_C ... +!> \param homo ... +!> \param virtual ... +!> \param homo_irred ... +!> \param fm_sqrt_A_minus_B ... +!> \param fm_inv_sqrt_A_minus_B ... +!> \param unit_nr ... +!> \param diag_est ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, & + fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, & + unit_nr, diag_est, mp2_env) + + TYPE(cp_fm_type), INTENT(INOUT) :: fm_C + INTEGER, INTENT(IN) :: homo, virtual, homo_irred + TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B + INTEGER, INTENT(IN) :: unit_nr + REAL(KIND=dp), INTENT(IN) :: diag_est + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_C' + + INTEGER :: diag_info, handle + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens + TYPE(cp_fm_type) :: fm_eigvec, fm_mat_eigvec_transform, & + fm_mat_eigvec_transform_neg + + CALL timeset(routineN, handle) + + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', & + 'This will take around ', diag_est, ' s.' + END IF + + !We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 + !Now: Diagonalize it + CALL cp_fm_create(fm_eigvec, fm_C%matrix_struct) + + ALLOCATE (Exc_ens(homo*virtual)) + + CALL choose_eigv_solver(fm_C, fm_eigvec, Exc_ens, diag_info) + CPASSERT(diag_info == 0) + Exc_ens = SQRT(Exc_ens) + + ! Prepare eigenvector for interpretation of singleparticle transitions + ! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001) + ! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result + + ! Following Furche, we basically use Eqs. (A10): First, we multiply + ! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues + ! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore! + + ! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^0.5 (A-B)^0.5 T_n + CALL cp_fm_create(fm_mat_eigvec_transform, fm_C%matrix_struct) + CALL cp_fm_set_all(fm_mat_eigvec_transform, 0.0_dp) + CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, & + matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec, beta=0.0_dp, & + matrix_c=fm_mat_eigvec_transform) + CALL cp_fm_release(fm_sqrt_A_minus_B) + CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.) + + ! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n + CALL cp_fm_create(fm_mat_eigvec_transform_neg, fm_C%matrix_struct) + CALL cp_fm_set_all(fm_mat_eigvec_transform_neg, 0.0_dp) + CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, & + matrix_a=fm_inv_sqrt_A_minus_B, matrix_b=fm_eigvec, beta=0.0_dp, & + matrix_c=fm_mat_eigvec_transform_neg) + CALL cp_fm_release(fm_inv_sqrt_A_minus_B) + CALL cp_fm_release(fm_eigvec) + CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_neg, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.) + + ! Now, we add the two equations to obtain X_n + CALL cp_fm_scale_and_add(1.0_dp, fm_mat_eigvec_transform, 1.0_dp, fm_mat_eigvec_transform_neg) + + !Cleanup + CALL cp_fm_release(fm_mat_eigvec_transform_neg) + + CALL success_message(Exc_ens, fm_mat_eigvec_transform, mp2_env, & + homo, virtual, homo_irred, unit_nr, .FALSE.) + + DEALLOCATE (Exc_ens) + CALL cp_fm_release(fm_mat_eigvec_transform) + + CALL timestop(handle) + + END SUBROUTINE diagonalize_C + +! ************************************************************************************************** +!> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n +!> \param fm_A ... +!> \param homo ... +!> \param virtual ... +!> \param homo_irred ... +!> \param unit_nr ... +!> \param diag_est ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, & + unit_nr, diag_est, mp2_env) + + TYPE(cp_fm_type), INTENT(INOUT) :: fm_A + INTEGER, INTENT(IN) :: homo, virtual, homo_irred, unit_nr + REAL(KIND=dp), INTENT(IN) :: diag_est + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_A' + + INTEGER :: diag_info, handle + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens + TYPE(cp_fm_type) :: fm_eigvec + + CALL timeset(routineN, handle) + + !Continue with formatting of subroutine create_A + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', & + 'This will take around ', diag_est, ' s.' + END IF + + !We have now the full matrix A_iajb, distributed over all ranks + !Now: Diagonalize it + CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct) + + ALLOCATE (Exc_ens(homo*virtual)) + + CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info) + CPASSERT(diag_info == 0) + CALL success_message(Exc_ens, fm_eigvec, mp2_env, & + homo, virtual, homo_irred, unit_nr, .TRUE.) + + CALL cp_fm_release(fm_eigvec) + DEALLOCATE (Exc_ens) + + CALL timestop(handle) + + END SUBROUTINE diagonalize_A + +! ************************************************************************************************** +!> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag) +!> \param Exc_ens ... +!> \param fm_eigvec ... +!> \param mp2_env ... +!> \param homo ... +!> \param virtual ... +!> \param homo_irred ... +!> \param unit_nr ... +!> \param flag_TDA ... +! ************************************************************************************************** + SUBROUTINE success_message(Exc_ens, fm_eigvec, mp2_env, & + homo, virtual, homo_irred, unit_nr, flag_TDA) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens + TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + INTEGER :: homo, virtual, homo_irred, unit_nr + LOGICAL, OPTIONAL :: flag_TDA + + CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message' + + CHARACTER(LEN=10) :: info_approximation, multiplet + INTEGER :: handle, i_exc, k, num_entries + INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt + REAL(KIND=dp) :: alpha + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries + + CALL timeset(routineN, handle) + + !Prepare variables for printing + IF (mp2_env%ri_g0w0%bse_spin_config == 0) THEN + multiplet = "Singlet" + alpha = 2.0_dp + ELSE + multiplet = "Triplet" + alpha = 0.0_dp + END IF + IF (.NOT. PRESENT(flag_TDA)) THEN + flag_TDA = .FALSE. + END IF + IF (flag_TDA) THEN + info_approximation = " -TDA- " + ELSE + info_approximation = "-full-" + END IF + + !Get information about eigenvector + + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + IF (flag_TDA) THEN + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************' + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA) *' + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', & + 'the BSE within the TDA:' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb X_jb^n ) = Ω^n X_ia^n' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0' + WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .' + ELSE + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************' + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Full Bethe Salpeter equation (BSE) (i.e. without TDA) *' + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', & + 'the BSE without the TDA:' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n| |1 0| |X^n|' + WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', ' sum_jb ( A_ia,jb X_jb^n + B_ia,jb Y_jb^n ) = Ω^n X_ia^n' + WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', '- sum_jb ( B_ia,jb X_jb^n + A_ia,jb Y_jb^n ) = Ω^n Y_ia^n' + END IF + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', & + 'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']' + WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', & + 'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']' + WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab' + IF (.NOT. flag_TDA) THEN + WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)' + WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .' + END IF + IF (.NOT. flag_TDA) THEN + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', 'The BSE is solved for Ω^n and X_ia^n as a hermitian problem, e.g. Eq.(42)' + WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .' + ! WRITE (unit_nr, '(T2,A4,T7,A69)') 'BSE|', 'C_ia,jb = sum_kc,ld ((A-B)^0.5)_ia,kc (A+B)_kc,ld ((A-B)^0.5)_ld,jb' + ! WRITE (unit_nr, '(T2,A4)') 'BSE|' + ! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X+Y)_ia,n = sum_jb,m ((A-B)^0.5)ia,jb Z_jb,m (Ω_m)^-0.5' + ! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X-Y)_ia,n = sum_jb,m ((A-B)^-0.5)ia,jb Z_jb,m (Ω_m)^0.5' + END IF + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A7,T19,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy' + WRITE (unit_nr, '(T2,A4,T7,A7,T19,A15)') 'BSE|', 'δ_...:', 'Kronecker delta' + WRITE (unit_nr, '(T2,A4,T7,A3,T19,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)' + WRITE (unit_nr, '(T2,A4,T7,A6,T18,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction' + WRITE (unit_nr, '(T2,A4,T7,A6,T18,A27)') 'BSE|', 'W_...:', 'Screened direct interaction' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1))') 'BSE|', & + 'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha + WRITE (unit_nr, '(T2,A4)') 'BSE|' + IF (flag_TDA) THEN + WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:' + ELSE + WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:' + END IF + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T13,A10,T27,A13,T42,A12,T59,A22)') 'BSE|', & + 'Excitation', "Multiplet", 'TDA/full BSE', 'Excitation energy (eV)' + END IF + !prints actual energies values + IF (unit_nr > 0) THEN + DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc) + WRITE (unit_nr, '(T2,A4,T7,I16,T27,A7,A6,T48,A6,T59,F22.4)') & + 'BSE|', i_exc, multiplet, " State", info_approximation, Exc_ens(i_exc)*evolt + END DO + END IF + !prints single particle transitions + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4)') 'BSE|' + + WRITE (unit_nr, '(T2,A4,T7,A70)') & + 'BSE|', "Excitations are built up by the following single-particle transitions," + WRITE (unit_nr, '(T2,A4,T7,A42,F5.2,A2)') & + 'BSE|', "neglecting contributions where |X_ia^n| < ", mp2_env%ri_g0w0%eps_x, " :" + + WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', & + homo_irred, ' and LUMO a =', homo_irred + 1, " --" + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A9,T20,A1,T22,A2,T29,A1,T42,A12,T71,A10)') & + "BSE|", "n-th exc.", "i", "=>", "a", 'TDA/full BSE', "|X_ia^n|" + END IF + DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc) + !Iterate through eigenvector and print values above threshold + CALL filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, & + i_exc, virtual, num_entries, mp2_env) + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4)') 'BSE|' + DO k = 1, num_entries + WRITE (unit_nr, '(T2,A4,T11,I5,T16,I5,T22,A2,T25,I5,T48,A6,T59,F22.4)') & + "BSE|", i_exc, homo_irred - homo + idx_homo(k), "=>", & + homo_irred + idx_virt(k), info_approximation, ABS(eigvec_entries(k)) + END DO + END IF + + DEALLOCATE (idx_homo) + DEALLOCATE (idx_virt) + DEALLOCATE (eigvec_entries) + END DO + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!' + END IF + + CALL timestop(handle) + END SUBROUTINE success_message + +END MODULE bse_full_diag diff --git a/src/bse_iterative.F b/src/bse_iterative.F new file mode 100644 index 0000000000..658d9ca6f6 --- /dev/null +++ b/src/bse_iterative.F @@ -0,0 +1,1587 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Iterative routines for GW + Bethe-Salpeter for computing electronic excitations +!> \par History +!> 04.2017 created [Jan Wilhelm] +!> 11.2023 Davidson solver implemented [Maximilian Graml] +! ************************************************************************************************** +MODULE bse_iterative + USE cp_fm_types, ONLY: cp_fm_get_info,& + cp_fm_type + USE group_dist_types, ONLY: get_group_dist,& + group_dist_d1_type + USE input_constants, ONLY: bse_singlet,& + bse_triplet + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type,& + mp_request_type + USE mp2_types, ONLY: integ_mat_buffer_type,& + mp2_type + USE physcon, ONLY: evolt + USE rpa_communication, ONLY: communicate_buffer +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_iterative' + + PUBLIC :: fill_local_3c_arrays, do_subspace_iterations + +CONTAINS + +! ************************************************************************************************** +!> \brief ... +!> \param B_bar_ijQ_bse_local ... +!> \param B_abQ_bse_local ... +!> \param B_bar_iaQ_bse_local ... +!> \param B_iaQ_bse_local ... +!> \param homo ... +!> \param virtual ... +!> \param bse_spin_config ... +!> \param unit_nr ... +!> \param Eigenval ... +!> \param para_env ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & + B_iaQ_bse_local, homo, virtual, bse_spin_config, unit_nr, & + Eigenval, para_env, mp2_env) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_bar_ijQ_bse_local, B_abQ_bse_local, & + B_bar_iaQ_bse_local, B_iaQ_bse_local + INTEGER :: homo, virtual, bse_spin_config, unit_nr + REAL(KIND=dp), DIMENSION(:) :: Eigenval + TYPE(mp_para_env_type), INTENT(IN) :: para_env + TYPE(mp2_type) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'do_subspace_iterations' + + CHARACTER(LEN=10) :: bse_davidson_abort_cond_string, & + success_abort_string + INTEGER :: bse_davidson_abort_cond, davidson_converged, fac_max_z_space, handle, i_iter, & + j_print, local_RI_size, num_add_start_z_space, num_davidson_iter, & + num_en_unconverged = -1, num_exact_en_unconverged = -1, num_exc_en, num_max_z_space, & + num_new_t, num_res_unconverged = -1, num_Z_vectors, num_Z_vectors_init + LOGICAL :: bse_full_diag_debug + REAL(kind=dp) :: eps_exc_en, eps_res, max_en_diff, & + max_res_norm, z_space_energy_cutoff + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: En_diffs, En_diffs_exact, Full_exc_spectrum, & + Res_norms, Subspace_full_eigenval, Subspace_new_eigenval, Subspace_prev_eigenval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: AZ_reshaped, M_ia_tmp, M_ji_tmp, RI_vector, & + Subspace_new_eigenvec, Subspace_residuals_reshaped, Z_vectors_reshaped + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, BZ, Subspace_add_dir, & + Subspace_ritzvec, W_vectors, Z_vectors + + CALL timeset(routineN, handle) + + !MG to del + !Debug flag for exact diagonalization (only using lapack!!!) + bse_full_diag_debug = .TRUE. + + bse_davidson_abort_cond = mp2_env%ri_g0w0%davidson_abort_cond + num_exc_en = mp2_env%ri_g0w0%num_exc_en + num_add_start_z_space = mp2_env%ri_g0w0%num_add_start_z_space + fac_max_z_space = mp2_env%ri_g0w0%fac_max_z_space + num_new_t = mp2_env%ri_g0w0%num_new_t + num_davidson_iter = mp2_env%ri_g0w0%num_davidson_iter + eps_res = mp2_env%ri_g0w0%eps_res + eps_exc_en = mp2_env%ri_g0w0%eps_exc_en + z_space_energy_cutoff = mp2_env%ri_g0w0%z_space_energy_cutoff + + num_Z_vectors_init = num_exc_en + num_add_start_z_space + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "bse_spin_config", bse_spin_config + WRITE (unit_nr, *) "num_exc_en", num_exc_en + WRITE (unit_nr, *) "num_add_start_z_space", num_add_start_z_space + WRITE (unit_nr, *) "num_Z_vectors_init", num_Z_vectors_init + WRITE (unit_nr, *) "fac_max_z_space", fac_max_z_space + WRITE (unit_nr, *) "num_new_t", num_new_t + WRITE (unit_nr, *) "eps_res", eps_res + WRITE (unit_nr, *) "num_davidson_iter", num_davidson_iter + WRITE (unit_nr, *) "eps_exc_en", eps_exc_en + WRITE (unit_nr, *) "bse_davidson_abort_cond", bse_davidson_abort_cond + WRITE (unit_nr, *) "z_space_energy_cutoff", z_space_energy_cutoff + WRITE (unit_nr, *) "Printing B_bar_iaQ_bse_local of shape", SHAPE(B_bar_iaQ_bse_local) + END IF + + local_RI_size = SIZE(B_iaQ_bse_local, 3) + + num_Z_vectors = num_Z_vectors_init + num_max_z_space = num_Z_vectors_init*fac_max_z_space + + !Check input parameters and correct them if necessary + IF (num_new_t > num_Z_vectors_init) THEN + num_new_t = num_Z_vectors_init + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, "Number of added directions has to be smaller/equals than "// & + "initial dimension. Corrected num_new_t accordingly.") + END IF + END IF + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Between BSE correction Warnings" + END IF + !If initial number is too big, already the first iteration causes trouble in LAPACK diagonal. (DORGQR) + IF (2*num_Z_vectors_init > homo*virtual) THEN + CALL cp_abort(__LOCATION__, "Initial dimension was too large and could not be corrected. "// & + "Choose another num_exc_en and num_add_start_z_space or adapt your basis set.") + END IF + IF (num_max_z_space .GE. homo*virtual) THEN + fac_max_z_space = homo*virtual/num_Z_vectors_init + num_max_z_space = num_Z_vectors_init*fac_max_z_space + + IF (fac_max_z_space == 0) THEN + CALL cp_abort(__LOCATION__, "Maximal dimension was too large and could not be corrected. "// & + "Choose another fac_max_z_space and num_Z_vectors_init or adapt your basis set.") + ELSE + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, "Maximal dimension of Z space has to be smaller than homo*virtual. "// & + "Corrected fac_max_z_space accordingly.") + END IF + END IF + END IF + + DO i_iter = 1, num_davidson_iter + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Allocating Z_vec,AZ,BZ with dimensions (homo,virt,num_Z)", homo, virtual, num_Z_vectors + WRITE (unit_nr, *) 'ProcNr', para_env%mepos, 'you really enter here for i_iter', i_iter + END IF + ALLOCATE (Z_vectors(homo, virtual, num_Z_vectors)) + Z_vectors = 0.0_dp + + !Dellocation procedures are a bit intricate, W_/Z_vectors and eigenvalues are needed for the next iteration, + ! therefore we have to deallocate them separately from the other quantities + IF (i_iter == 1) THEN + CALL initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual) + ALLOCATE (Subspace_prev_eigenval(num_exc_en)) + Subspace_prev_eigenval = 0.0_dp + ELSE + Z_vectors(:, :, :) = W_vectors(:, :, :) + DEALLOCATE (W_vectors) + END IF + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Allocated/rewritten Z arrays" + END IF + + CALL create_bse_work_arrays(AZ, Z_vectors_reshaped, AZ_reshaped, BZ, M_ia_tmp, M_ji_tmp, & + RI_vector, Subspace_new_eigenval, Subspace_full_eigenval, Subspace_new_eigenvec, & + Subspace_residuals_reshaped, Subspace_ritzvec, Subspace_add_dir, W_vectors, & + homo, virtual, num_Z_vectors, local_RI_size, num_new_t) + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Allocated Work arrays" + END IF + + CALL compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, & + M_ia_tmp, RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, & + para_env, bse_spin_config, z_space_energy_cutoff, i_iter, bse_full_diag_debug, & + Full_exc_spectrum, unit_nr) + + !MG: functionality of BZ not checked (issue with fm_mat_Q_static_bse_gemm in rpa_util needs to be checked!) + !CALL compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, & + ! M_ji_tmp, homo, virtual, num_Z_vectors, local_RI_size, & + ! para_env) + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Computed AZ" + END IF + + !MG to check: Reshaping correct? + AZ_reshaped(:, :) = RESHAPE(AZ, (/homo*virtual, num_Z_vectors/)) + Z_vectors_reshaped(:, :) = RESHAPE(Z_vectors, (/homo*virtual, num_Z_vectors/)) + + ! Diagonalize M and extract smallest eigenvalues/corresponding eigenvector + CALL compute_diagonalize_ZAZ(AZ_reshaped, Z_vectors_reshaped, num_Z_vectors, Subspace_new_eigenval, & + Subspace_new_eigenvec, num_new_t, Subspace_full_eigenval, para_env, unit_nr) + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Eigenval (eV) in iter=", i_iter, " is:", Subspace_new_eigenval(:6)*evolt + END IF + + ! Threshold in energies + CALL check_en_convergence(Subspace_full_eigenval, Subspace_prev_eigenval, eps_exc_en, num_en_unconverged, & + num_exc_en, max_en_diff, En_diffs) + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Largest change of desired exc ens =", max_en_diff + END IF + ! Compute residuals + CALL compute_residuals(AZ_reshaped, Z_vectors_reshaped, Subspace_new_eigenval, Subspace_new_eigenvec, & + Subspace_residuals_reshaped, homo, virtual, num_new_t, num_Z_vectors, Subspace_ritzvec) + + !Abort, if residuals are small enough w.r.t threshold + CALL check_res_convergence(Subspace_residuals_reshaped, num_new_t, eps_res, num_res_unconverged, & + i_iter, max_res_norm, unit_nr, Res_norms) + + davidson_converged = -1 + IF (num_res_unconverged == 0 .AND. bse_davidson_abort_cond /= 0) THEN + davidson_converged = 1 + success_abort_string = "RESIDUALS" + ELSE IF (num_en_unconverged == 0 .AND. (bse_davidson_abort_cond /= 1)) THEN + davidson_converged = 1 + success_abort_string = "ENERGIES" + ELSE IF (i_iter == num_davidson_iter) THEN + davidson_converged = -100 + success_abort_string = "-----" + ELSE + davidson_converged = -1 + END IF + + IF (bse_davidson_abort_cond == 0) THEN + bse_davidson_abort_cond_string = "ENERGY" + ELSE IF (bse_davidson_abort_cond == 1) THEN + bse_davidson_abort_cond_string = "RESIDUAL" + ELSE + bse_davidson_abort_cond_string = "EITHER" + END IF + + IF (davidson_converged == 1) THEN + CALL success_message(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, & + bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, & + num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space, max_res_norm, & + max_en_diff, num_en_unconverged, bse_davidson_abort_cond_string, & + eps_exc_en, success_abort_string, z_space_energy_cutoff) + + !Deallocate matrices, which are otherwise not cleared due to exiting the loop + DEALLOCATE (AZ, BZ, & + Z_vectors, M_ia_tmp, M_ji_tmp, RI_vector, Subspace_prev_eigenval, & + Subspace_new_eigenval, Subspace_new_eigenvec, Subspace_residuals_reshaped, & + Subspace_add_dir, AZ_reshaped, Z_vectors_reshaped, Subspace_ritzvec, Subspace_full_eigenval) + + EXIT + ELSE IF (davidson_converged < -1) THEN + CALL print_davidson_parameter(i_iter, num_davidson_iter, num_Z_vectors, num_res_unconverged, max_res_norm, & + eps_res, num_en_unconverged, max_en_diff, eps_exc_en, num_exc_en, & + num_Z_vectors_init, num_max_z_space, num_new_t, unit_nr, & + success_abort_string, bse_davidson_abort_cond_string, z_space_energy_cutoff) + + CALL cp_abort(__LOCATION__, "BSE/TDA-Davidson did not converge using "// & + bse_davidson_abort_cond_string//" threshold condition!") + END IF + + ! Calculate and add next orthonormal vector and update num_Z_vectors + CALL compute_new_directions(homo, virtual, Subspace_residuals_reshaped, Subspace_new_eigenval, Eigenval, & + num_new_t, Subspace_add_dir) + + !If exact-diag: compute difference to exact eigenvalues + IF (bse_full_diag_debug) THEN + ALLOCATE (En_diffs_exact(num_exc_en)) + num_exact_en_unconverged = 0 + DO j_print = 1, num_exc_en + En_diffs_exact(j_print) = ABS(Subspace_full_eigenval(j_print) - Full_exc_spectrum(j_print)) + IF (En_diffs_exact(j_print) > eps_exc_en) num_exact_en_unconverged = num_exact_en_unconverged + 1 + END DO + END IF + + !Check dimensions and orthonormalize vector system, depending on dimensionality + CALL check_Z_space_dimension(W_vectors, Z_vectors, Subspace_add_dir, Subspace_ritzvec, & + num_Z_vectors, num_new_t, num_max_z_space, homo, virtual, i_iter, unit_nr) + + !Copy eigenvalues for threshold + Subspace_prev_eigenval(:) = Subspace_full_eigenval(:num_exc_en) + + DEALLOCATE (AZ, & !BZ, + Z_vectors, M_ia_tmp, M_ji_tmp, RI_vector, & + Subspace_new_eigenval, Subspace_new_eigenvec, Subspace_residuals_reshaped, & + Subspace_add_dir, AZ_reshaped, Z_vectors_reshaped, Subspace_ritzvec, Subspace_full_eigenval, & + Res_norms, En_diffs) + + IF (bse_full_diag_debug) THEN + DEALLOCATE (En_diffs_exact) + END IF + + !Orthonorm: + CALL orthonormalize_W(W_vectors, num_Z_vectors, homo, virtual) + + END DO + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param W_vectors ... +!> \param Z_vectors ... +!> \param Subspace_add_dir ... +!> \param Subspace_ritzvec ... +!> \param num_Z_vectors ... +!> \param num_new_t ... +!> \param num_max_z_space ... +!> \param homo ... +!> \param virtual ... +!> \param i_iter ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE check_Z_space_dimension(W_vectors, Z_vectors, Subspace_add_dir, Subspace_ritzvec, & + num_Z_vectors, num_new_t, num_max_z_space, homo, virtual, i_iter, unit_nr) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: W_vectors, Z_vectors, Subspace_add_dir, & + Subspace_ritzvec + INTEGER :: num_Z_vectors, num_new_t, & + num_max_z_space, homo, virtual, & + i_iter, unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'check_Z_space_dimension' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (num_Z_vectors + num_new_t .LE. num_max_z_space) THEN + W_vectors(:, :, :num_Z_vectors) = Z_vectors(:, :, :) + W_vectors(:, :, num_Z_vectors + 1:) = Subspace_add_dir + num_Z_vectors = num_Z_vectors + num_new_t + ELSE + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Resetting dimension in i_iter=", i_iter + END IF + DEALLOCATE (W_vectors) + ALLOCATE (W_vectors(homo, virtual, 2*num_new_t)) + W_vectors(:, :, :num_new_t) = Subspace_ritzvec(:, :, :) + W_vectors(:, :, num_new_t + 1:) = Subspace_add_dir + num_Z_vectors = 2*num_new_t + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param AZ ... +!> \param Z_vectors_reshaped ... +!> \param AZ_reshaped ... +!> \param BZ ... +!> \param M_ia_tmp ... +!> \param M_ji_tmp ... +!> \param RI_vector ... +!> \param Subspace_new_eigenval ... +!> \param Subspace_full_eigenval ... +!> \param Subspace_new_eigenvec ... +!> \param Subspace_residuals_reshaped ... +!> \param Subspace_ritzvec ... +!> \param Subspace_add_dir ... +!> \param W_vectors ... +!> \param homo ... +!> \param virtual ... +!> \param num_Z_vectors ... +!> \param local_RI_size ... +!> \param num_new_t ... +! ************************************************************************************************** + SUBROUTINE create_bse_work_arrays(AZ, Z_vectors_reshaped, AZ_reshaped, BZ, M_ia_tmp, M_ji_tmp, & + RI_vector, Subspace_new_eigenval, Subspace_full_eigenval, Subspace_new_eigenvec, & + Subspace_residuals_reshaped, Subspace_ritzvec, Subspace_add_dir, W_vectors, & + homo, virtual, num_Z_vectors, local_RI_size, num_new_t) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Z_vectors_reshaped, AZ_reshaped + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BZ + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_ia_tmp, M_ji_tmp, RI_vector + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval, & + Subspace_full_eigenval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_new_eigenvec, & + Subspace_residuals_reshaped + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Subspace_ritzvec, Subspace_add_dir, & + W_vectors + INTEGER :: homo, virtual, num_Z_vectors, & + local_RI_size, num_new_t + + CHARACTER(LEN=*), PARAMETER :: routineN = 'create_bse_work_arrays' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + ALLOCATE (AZ(homo, virtual, num_Z_vectors)) + AZ = 0.0_dp + + ALLOCATE (Z_vectors_reshaped(homo*virtual, num_Z_vectors)) + Z_vectors_reshaped = 0.0_dp + + ALLOCATE (AZ_reshaped(homo*virtual, num_Z_vectors)) + AZ_reshaped = 0.0_dp + + ALLOCATE (BZ(homo, virtual, num_Z_vectors)) + BZ = 0.0_dp + + ALLOCATE (M_ia_tmp(homo, virtual)) + M_ia_tmp = 0.0_dp + + ALLOCATE (M_ji_tmp(homo, homo)) + M_ji_tmp = 0.0_dp + + ALLOCATE (RI_vector(local_RI_size, num_Z_vectors)) + RI_vector = 0.0_dp + + ALLOCATE (Subspace_new_eigenval(num_new_t)) + Subspace_new_eigenval = 0.0_dp + + ALLOCATE (Subspace_full_eigenval(num_Z_vectors)) + Subspace_full_eigenval = 0.0_dp + + ALLOCATE (Subspace_new_eigenvec(num_Z_vectors, num_new_t)) + Subspace_new_eigenvec = 0.0_dp + + ALLOCATE (subspace_residuals_reshaped(homo*virtual, num_new_t)) + subspace_residuals_reshaped = 0.0_dp + + ALLOCATE (Subspace_ritzvec(homo, virtual, num_new_t)) + Subspace_ritzvec = 0.0_dp + + ALLOCATE (Subspace_add_dir(homo, virtual, num_new_t)) + Subspace_add_dir = 0.0_dp + + ALLOCATE (W_vectors(homo, virtual, num_Z_vectors + num_new_t)) + W_vectors = 0.0_dp + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param Subspace_full_eigenval ... +!> \param num_new_t ... +!> \param eps_res ... +!> \param num_res_unconverged ... +!> \param bse_spin_config ... +!> \param unit_nr ... +!> \param num_exc_en ... +!> \param num_Z_vectors_init ... +!> \param num_davidson_iter ... +!> \param i_iter ... +!> \param num_Z_vectors ... +!> \param num_max_z_space ... +!> \param max_res_norm ... +!> \param max_en_diff ... +!> \param num_en_unconverged ... +!> \param bse_davidson_abort_cond_string ... +!> \param eps_exc_en ... +!> \param success_abort_string ... +!> \param z_space_energy_cutoff ... +! ************************************************************************************************** + SUBROUTINE success_message(Subspace_full_eigenval, num_new_t, eps_res, num_res_unconverged, & + bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, & + num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space, max_res_norm, & + max_en_diff, num_en_unconverged, bse_davidson_abort_cond_string, & + eps_exc_en, success_abort_string, z_space_energy_cutoff) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_full_eigenval + INTEGER :: num_new_t + REAL(KIND=dp) :: eps_res + INTEGER :: num_res_unconverged, bse_spin_config, unit_nr, num_exc_en, num_Z_vectors_init, & + num_davidson_iter, i_iter, num_Z_vectors, num_max_z_space + REAL(KIND=dp) :: max_res_norm, max_en_diff + INTEGER :: num_en_unconverged + CHARACTER(LEN=10) :: bse_davidson_abort_cond_string + REAL(KIND=dp) :: eps_exc_en + CHARACTER(LEN=10) :: success_abort_string + REAL(KIND=dp) :: z_space_energy_cutoff + + CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message' + + CHARACTER(LEN=10) :: multiplet + INTEGER :: handle, i + REAL(KIND=dp) :: alpha + + CALL timeset(routineN, handle) + + !Prepare variables for printing + SELECT CASE (bse_spin_config) + CASE (bse_singlet) + alpha = 2.0_dp + multiplet = "Singlet" + CASE (bse_triplet) + alpha = 0.0_dp + multiplet = "Triplet" + END SELECT + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) ' ' + WRITE (unit_nr, '(T3,A)') '******************************************************************************' + WRITE (unit_nr, '(T3,A)') '** **' + WRITE (unit_nr, '(T3,A)') '** BSE-TDA EXCITONIC ENERGIES **' + WRITE (unit_nr, '(T3,A)') '** **' + WRITE (unit_nr, '(T3,A)') '******************************************************************************' + WRITE (unit_nr, '(T3,A)') ' ' + WRITE (unit_nr, '(T3,A)') ' ' + WRITE (unit_nr, '(T3,A)') ' The excitation energies are calculated by iteratively diagonalizing: ' + WRITE (unit_nr, '(T3,A)') ' ' + WRITE (unit_nr, '(T3,A)') ' A_iajb = (E_a-E_i) delta_ij delta_ab + alpha * v_iajb - W_ijab ' + WRITE (unit_nr, '(T3,A)') ' ' + WRITE (unit_nr, '(T3,A48,A7,A12,F3.1)') & + ' The spin-dependent factor for the requested ', multiplet, " is alpha = ", alpha + WRITE (unit_nr, '(T3,A)') ' ' + WRITE (unit_nr, '(T3,A16,T50,A22)') & + ' Excitonic level', 'Excitation energy (eV)' + !prints actual energies values + DO i = 1, num_exc_en + WRITE (unit_nr, '(T3,I16,T50,F22.3)') i, Subspace_full_eigenval(i)*evolt + END DO + + WRITE (unit_nr, '(T3,A)') ' ' + + !prints parameters of Davidson algorithm + CALL print_davidson_parameter(i_iter, num_davidson_iter, num_Z_vectors, num_res_unconverged, max_res_norm, & + eps_res, num_en_unconverged, max_en_diff, eps_exc_en, num_exc_en, & + num_Z_vectors_init, num_max_z_space, num_new_t, unit_nr, & + success_abort_string, bse_davidson_abort_cond_string, z_space_energy_cutoff) + + !Insert warning if energies are not converged (could probably be the case if one uses residual threshold) + IF (num_en_unconverged > 0) THEN + WRITE (unit_nr, '(T3,A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + WRITE (unit_nr, '(T3,A2,T79,A2)') '!!', "!!" + WRITE (unit_nr, '(T3,A2,T8,A65,T79,A2)') '!!', "THERE ARE UNCONVERGED ENERGIES PRINTED OUT, SOMETHING WENT WRONG!", "!!" + WRITE (unit_nr, '(T3,A2,T79,A2)') '!!', "!!" + WRITE (unit_nr, '(T3,A)') '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + END IF + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param i_iter ... +!> \param num_davidson_iter ... +!> \param num_Z_vectors ... +!> \param num_res_unconverged ... +!> \param max_res_norm ... +!> \param eps_res ... +!> \param num_en_unconverged ... +!> \param max_en_diff ... +!> \param eps_exc_en ... +!> \param num_exc_en ... +!> \param num_Z_vectors_init ... +!> \param num_max_z_space ... +!> \param num_new_t ... +!> \param unit_nr ... +!> \param success_abort_string ... +!> \param bse_davidson_abort_cond_string ... +!> \param z_space_energy_cutoff ... +! ************************************************************************************************** + SUBROUTINE print_davidson_parameter(i_iter, num_davidson_iter, num_Z_vectors, num_res_unconverged, max_res_norm, & + eps_res, num_en_unconverged, max_en_diff, eps_exc_en, num_exc_en, & + num_Z_vectors_init, num_max_z_space, num_new_t, unit_nr, & + success_abort_string, bse_davidson_abort_cond_string, z_space_energy_cutoff) + + INTEGER :: i_iter, num_davidson_iter, & + num_Z_vectors, num_res_unconverged + REAL(KIND=dp) :: max_res_norm, eps_res + INTEGER :: num_en_unconverged + REAL(KIND=dp) :: max_en_diff, eps_exc_en + INTEGER :: num_exc_en, num_Z_vectors_init, & + num_max_z_space, num_new_t, unit_nr + CHARACTER(LEN=10) :: success_abort_string, & + bse_davidson_abort_cond_string + REAL(KIND=dp) :: z_space_energy_cutoff + + CHARACTER(LEN=*), PARAMETER :: routineN = 'print_davidson_parameter' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + WRITE (unit_nr, '(T3,A)') '******************************************************************************' + WRITE (unit_nr, '(T3,A2,T15,A49,T79,A2)') & + '**', "Parameters of the BSE-Davidson solver:", "**" + WRITE (unit_nr, '(T3,A2,T79,A2)') & + '**', "**" + WRITE (unit_nr, '(T3,A2,T79,A2)') & + '**', "**" + WRITE (unit_nr, '(T3,A2,T10,A16,I5,A12,I5,A8,T79,A2)') & + '**', "Converged after ", i_iter, " of maximal ", num_davidson_iter, " cycles,", "**" + WRITE (unit_nr, '(T3,A2,T20,A11,A9,A7,A8,A20,T79,A2)') & + '**', "because of ", success_abort_string, " using ", & + bse_davidson_abort_cond_string, " threshold condition", "**" + WRITE (unit_nr, '(T3,A2,T79,A2)') & + '**', "**" + WRITE (unit_nr, '(T3,A2,T10,A32,T65,I11,T79,A2)') & + '**', "The Z space has at the end dim. ", num_Z_vectors, "**" + WRITE (unit_nr, '(T3,A2,T10,A45,T65,I11,T79,A2)') & + '**', "Number of unconverged residuals in subspace: ", num_res_unconverged, "**" + WRITE (unit_nr, '(T3,A2,T10,A35,T65,E11.4,T79,A2)') & + '**', "largest unconverged residual (eV): ", max_res_norm*evolt, "**" + WRITE (unit_nr, '(T3,A2,T10,A45,T65,E11.4,T79,A2)') & + '**', "threshold for convergence of residuals (eV): ", eps_res*evolt, "**" + WRITE (unit_nr, '(T3,A2,T10,A45,T65,I11,T79,A2)') & + '**', "Number of desired, but unconverged energies: ", num_en_unconverged, "**" + WRITE (unit_nr, '(T3,A2,T10,A44,T65,E11.4,T79,A2)') & + '**', "largest unconverged energy difference (eV): ", max_en_diff*evolt, "**" + WRITE (unit_nr, '(T3,A2,T10,A44,T65,E11.4,T79,A2)') & + '**', "threshold for convergence of energies (eV): ", eps_exc_en*evolt, "**" + WRITE (unit_nr, '(T3,A2,T10,A40,T65,I11,T79,A2)') & + '**', "number of computed excitation energies: ", num_exc_en, "**" + + IF (z_space_energy_cutoff > 0) THEN + WRITE (unit_nr, '(T3,A2,T10,A37,T65,E11.4,T79,A2)') & + '**', "cutoff for excitation energies (eV): ", z_space_energy_cutoff*evolt, "**" + END IF + + WRITE (unit_nr, '(T3,A2,T10,A36,T65,I11,T79,A2)') & + '**', "number of Z space at the beginning: ", num_Z_vectors_init, "**" + WRITE (unit_nr, '(T3,A2,T10,A30,T65,I11,T79,A2)') & + '**', "maximal dimension of Z space: ", num_max_z_space, "**" + WRITE (unit_nr, '(T3,A2,T10,A31,T65,I11,T79,A2)') & + '**', "added directions per iteration: ", num_new_t, "**" + WRITE (unit_nr, '(T3,A2,T79,A2)') & + '**', "**" + WRITE (unit_nr, '(T3,A2,T79,A2)') & + '**', "**" + WRITE (unit_nr, '(T3,A)') '******************************************************************************' + WRITE (unit_nr, '(T3,A)') ' ' + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param Subspace_full_eigenval ... +!> \param Subspace_prev_eigenval ... +!> \param eps_exc_en ... +!> \param num_en_unconverged ... +!> \param num_exc_en ... +!> \param max_en_diff ... +!> \param En_diffs ... +! ************************************************************************************************** + SUBROUTINE check_en_convergence(Subspace_full_eigenval, Subspace_prev_eigenval, eps_exc_en, num_en_unconverged, & + num_exc_en, max_en_diff, En_diffs) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_full_eigenval, & + Subspace_prev_eigenval + REAL(KIND=dp) :: eps_exc_en + INTEGER :: num_en_unconverged, num_exc_en + REAL(KIND=dp) :: max_en_diff + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: En_diffs + + CHARACTER(LEN=*), PARAMETER :: routineN = 'check_en_convergence' + + INTEGER :: handle, mu_l + + CALL timeset(routineN, handle) + + num_en_unconverged = 0 + ALLOCATE (En_diffs(num_exc_en)) + DO mu_l = 1, num_exc_en + En_diffs(mu_l) = ABS(Subspace_full_eigenval(mu_l) - Subspace_prev_eigenval(mu_l)) + IF (En_diffs(mu_l) > eps_exc_en) num_en_unconverged = num_en_unconverged + 1 + END DO + max_en_diff = MAXVAL(En_diffs) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param Subspace_residuals_reshaped ... +!> \param num_new_t ... +!> \param eps_res ... +!> \param num_res_unconverged ... +!> \param i_iter ... +!> \param max_res_norm ... +!> \param unit_nr ... +!> \param Res_norms ... +! ************************************************************************************************** + SUBROUTINE check_res_convergence(Subspace_residuals_reshaped, num_new_t, eps_res, num_res_unconverged, & + i_iter, max_res_norm, unit_nr, Res_norms) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_residuals_reshaped + INTEGER :: num_new_t + REAL(KIND=dp) :: eps_res + INTEGER :: num_res_unconverged, i_iter + REAL(KIND=dp) :: max_res_norm + INTEGER :: unit_nr + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Res_norms + + CHARACTER(LEN=*), PARAMETER :: routineN = 'check_res_convergence' + + INTEGER :: handle, mu_l + + CALL timeset(routineN, handle) + + num_res_unconverged = 0 + ALLOCATE (Res_norms(num_new_t)) + DO mu_l = 1, num_new_t + Res_norms(mu_l) = NORM2(Subspace_residuals_reshaped(:, mu_l)) + IF (Res_norms(mu_l) > eps_res) THEN + num_res_unconverged = num_res_unconverged + 1 + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Unconverged res in i_iter=", i_iter, "is:", Res_norms(mu_l) + END IF + END IF + END DO + max_res_norm = MAXVAL(Res_norms) + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) "Maximal unconverged res (of ", num_res_unconverged, & + " unconverged res in this step) in i_iter=", i_iter, "is:", max_res_norm + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param W_vectors ... +!> \param num_Z_vectors ... +!> \param homo ... +!> \param virtual ... +! ************************************************************************************************** + SUBROUTINE orthonormalize_W(W_vectors, num_Z_vectors, homo, virtual) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: W_vectors + INTEGER :: num_Z_vectors, homo, virtual + + CHARACTER(LEN=*), PARAMETER :: routineN = 'orthonormalize_W' + + INTEGER :: handle, info_dor, info_orth, LWORK_dor, & + LWORK_W + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Tau_W, WORK_W, WORK_W_dor + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: W_vectors_reshaped + + CALL timeset(routineN, handle) + + ALLOCATE (W_vectors_reshaped(homo*virtual, num_Z_vectors)) + W_vectors_reshaped(:, :) = RESHAPE(W_vectors, (/homo*virtual, num_Z_vectors/)) + + ALLOCATE (Tau_W(MIN(homo*virtual, num_Z_vectors))) + Tau_W = 0.0_dp + + ALLOCATE (WORK_W(1)) + WORK_W = 0.0_dp + + ALLOCATE (WORK_W_dor(1)) + WORK_W_dor = 0.0_dp + + CALL DGEQRF(homo*virtual, num_Z_vectors, W_vectors_reshaped, homo*virtual, Tau_W, WORK_W, -1, info_orth) + LWORK_W = INT(WORK_W(1)) + DEALLOCATE (WORK_W) + ALLOCATE (WORK_W(LWORK_W)) + WORK_W = 0.0_dp + CALL DGEQRF(homo*virtual, num_Z_vectors, W_vectors_reshaped, homo*virtual, Tau_W, WORK_W, LWORK_W, info_orth) + IF (info_orth /= 0) THEN + CPABORT("QR Decomp Step 1 doesnt work") + END IF + CALL DORGQR(homo*virtual, num_Z_vectors, MIN(homo*virtual, num_Z_vectors), W_vectors_reshaped, homo*virtual, & + Tau_W, WORK_W_dor, -1, info_dor) + LWORK_dor = INT(WORK_W_dor(1)) + DEALLOCATE (WORK_W_dor) + ALLOCATE (WORK_W_dor(LWORK_dor)) + WORK_W_dor = 0.0_dp + CALL DORGQR(homo*virtual, num_Z_vectors, MIN(homo*virtual, num_Z_vectors), W_vectors_reshaped, homo*virtual, & + Tau_W, WORK_W_dor, LWORK_dor, info_dor) + IF (info_orth /= 0) THEN + CPABORT("QR Decomp Step 2 doesnt work") + END IF + + W_vectors(:, :, :) = RESHAPE(W_vectors_reshaped, (/homo, virtual, num_Z_vectors/)) + + DEALLOCATE (WORK_W, WORK_W_dor, Tau_W, W_vectors_reshaped) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param homo ... +!> \param virtual ... +!> \param Subspace_residuals_reshaped ... +!> \param Subspace_new_eigenval ... +!> \param Eigenval ... +!> \param num_new_t ... +!> \param Subspace_add_dir ... +! ************************************************************************************************** + SUBROUTINE compute_new_directions(homo, virtual, Subspace_residuals_reshaped, Subspace_new_eigenval, Eigenval, & + num_new_t, Subspace_add_dir) + + INTEGER :: homo, virtual + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_residuals_reshaped + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval + REAL(KIND=dp), DIMENSION(:) :: Eigenval + INTEGER :: num_new_t + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Subspace_add_dir + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_new_directions' + + INTEGER :: a_virt, handle, i_occ, mu_subspace, & + prec_neg + REAL(KIND=dp) :: prec_scalar + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_add_dir_reshaped + + CALL timeset(routineN, handle) + + ALLOCATE (Subspace_add_dir_reshaped(homo*virtual, num_new_t)) + + prec_neg = 0 + DO mu_subspace = 1, num_new_t + DO i_occ = 1, homo + DO a_virt = 1, virtual + !MG to check: Indexorder and range of indices + prec_scalar = -1/(Subspace_new_eigenval(mu_subspace) - (Eigenval(a_virt + homo) - Eigenval(i_occ))) + IF (prec_scalar < 0) THEN + prec_neg = prec_neg + 1 + !prec_scalar = - prec_scalar + END IF + Subspace_add_dir_reshaped((i_occ - 1)*virtual + a_virt, mu_subspace) = prec_scalar* & + Subspace_residuals_reshaped((i_occ - 1)*virtual + a_virt, mu_subspace) + END DO + END DO + END DO + + Subspace_add_dir(:, :, :) = RESHAPE(Subspace_add_dir_reshaped, (/homo, virtual, num_new_t/)) + + DEALLOCATE (Subspace_add_dir_reshaped) + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param AZ_reshaped ... +!> \param Z_vectors_reshaped ... +!> \param Subspace_new_eigenval ... +!> \param Subspace_new_eigenvec ... +!> \param Subspace_residuals_reshaped ... +!> \param homo ... +!> \param virtual ... +!> \param num_new_t ... +!> \param num_Z_vectors ... +!> \param Subspace_ritzvec ... +! ************************************************************************************************** + SUBROUTINE compute_residuals(AZ_reshaped, Z_vectors_reshaped, Subspace_new_eigenval, Subspace_new_eigenvec, & + Subspace_residuals_reshaped, homo, virtual, num_new_t, num_Z_vectors, Subspace_ritzvec) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: AZ_reshaped, Z_vectors_reshaped + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_new_eigenvec, & + Subspace_residuals_reshaped + INTEGER :: homo, virtual, num_new_t, num_Z_vectors + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Subspace_ritzvec + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_residuals' + + INTEGER :: handle, mu_subspace + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_res_A, Subspace_res_ev + + CALL timeset(routineN, handle) + + ALLOCATE (Subspace_res_ev(homo*virtual, num_new_t)) + Subspace_res_ev = 0.0_dp + + ALLOCATE (Subspace_res_A(homo*virtual, num_new_t)) + Subspace_res_A = 0.0_dp + + !Compute all residuals in one loop, iterating over number of new/added t per iteration + DO mu_subspace = 1, num_new_t + + CALL DGEMM("N", "N", homo*virtual, 1, num_Z_vectors, 1.0_dp, Z_vectors_reshaped, homo*virtual, & + Subspace_new_eigenvec(:, mu_subspace), num_Z_vectors, 0.0_dp, Subspace_res_ev(:, mu_subspace), homo*virtual) + + CALL DGEMM("N", "N", homo*virtual, 1, num_Z_vectors, 1.0_dp, AZ_reshaped, homo*virtual, & + Subspace_new_eigenvec(:, mu_subspace), num_Z_vectors, 0.0_dp, Subspace_res_A(:, mu_subspace), homo*virtual) + + Subspace_residuals_reshaped(:, mu_subspace) = Subspace_new_eigenval(mu_subspace)*Subspace_res_ev(:, mu_subspace) & + - Subspace_res_A(:, mu_subspace) + + END DO + Subspace_ritzvec(:, :, :) = RESHAPE(Subspace_res_ev, (/homo, virtual, num_new_t/)) + DEALLOCATE (Subspace_res_ev, Subspace_res_A) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param AZ_reshaped ... +!> \param Z_vectors_reshaped ... +!> \param num_Z_vectors ... +!> \param Subspace_new_eigenval ... +!> \param Subspace_new_eigenvec ... +!> \param num_new_t ... +!> \param Subspace_full_eigenval ... +!> \param para_env ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE compute_diagonalize_ZAZ(AZ_reshaped, Z_vectors_reshaped, num_Z_vectors, Subspace_new_eigenval, & + Subspace_new_eigenvec, num_new_t, Subspace_full_eigenval, para_env, unit_nr) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: AZ_reshaped, Z_vectors_reshaped + INTEGER, INTENT(in) :: num_Z_vectors + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_new_eigenval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Subspace_new_eigenvec + INTEGER, INTENT(in) :: num_new_t + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Subspace_full_eigenval + TYPE(mp_para_env_type), INTENT(IN) :: para_env + INTEGER, INTENT(in) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_diagonalize_ZAZ' + + INTEGER :: handle, i_Z_vector, j_Z_vector, LWORK, & + ZAZ_diag_info + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: WORK + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ZAZ + + CALL timeset(routineN, handle) + + ALLOCATE (ZAZ(num_Z_vectors, num_Z_vectors)) + ZAZ(:, :) = 0.0_dp + + !Flatten AZ and Z matrices of a certain j_Z_vector w.r.t. occ and virt indices + !Multiply for each j_Z_vec and write into matrix of dim (num_Z_vec, num_Z_vec) + DO i_Z_vector = 1, num_Z_vectors + DO j_Z_vector = 1, num_Z_vectors + ZAZ(j_Z_vector, i_Z_vector) = DOT_PRODUCT(Z_vectors_reshaped(:, j_Z_vector), AZ_reshaped(:, i_Z_vector)) + END DO + END DO + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "Before Diag" + END IF + + !MG to do: Check for symmetry of ZAZ! + ALLOCATE (WORK(1)) + WORK = 0.0_dp + CALL DSYEV("V", "U", num_Z_vectors, ZAZ, num_Z_vectors, Subspace_full_eigenval, WORK, -1, ZAZ_diag_info) + LWORK = INT(WORK(1)) + DEALLOCATE (WORK) + ALLOCATE (WORK(LWORK)) + WORK = 0.0_dp + !MG to check: Usage of symmetric routine okay? (Correct LWORK?) + CALL DSYEV("V", "U", num_Z_vectors, ZAZ, num_Z_vectors, Subspace_full_eigenval, WORK, LWORK, ZAZ_diag_info) + + IF (ZAZ_diag_info /= 0) THEN + CPABORT("ZAZ could not be diagonalized successfully.") + END IF + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) 'ProcNr', para_env%mepos, "After Diag" + END IF + + Subspace_new_eigenval(1:num_new_t) = Subspace_full_eigenval(1:num_new_t) + Subspace_new_eigenvec(:, 1:num_new_t) = ZAZ(:, 1:num_new_t) + DEALLOCATE (WORK) + DEALLOCATE (ZAZ) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param BZ ... +!> \param Z_vectors ... +!> \param B_iaQ_bse_local ... +!> \param B_bar_iaQ_bse_local ... +!> \param M_ji_tmp ... +!> \param homo ... +!> \param virtual ... +!> \param num_Z_vectors ... +!> \param local_RI_size ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, & + M_ji_tmp, homo, virtual, num_Z_vectors, local_RI_size, para_env) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BZ, Z_vectors, B_iaQ_bse_local, & + B_bar_iaQ_bse_local + REAL(KIND=dp), DIMENSION(:, :) :: M_ji_tmp + INTEGER :: homo, virtual, num_Z_vectors, & + local_RI_size + TYPE(mp_para_env_type), INTENT(IN) :: para_env + + INTEGER :: i_Z_vector, LLL + + BZ(:, :, :) = 0.0_dp + + !CALL compute_v_ia_jb_part(BZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, & + ! num_Z_vectors, homo, virtual) + + DO i_Z_vector = 1, num_Z_vectors + + DO LLL = 1, local_RI_size + + ! M_ji^P = sum_b Z_jb*B_bi^P + CALL DGEMM("N", "T", homo, homo, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, & + B_iaQ_bse_local(:, :, LLL), homo, 0.0_dp, M_ji_tmp, homo) + ! (BZ)_ia = sum_jP M_ij^P*B^bar_ja^P + CALL DGEMM("T", "N", homo, virtual, homo, 1.0_dp, M_ji_tmp, homo, & + B_bar_iaQ_bse_local, homo, 1.0_dp, BZ(:, :, i_Z_vector), homo) + + END DO + + END DO + + ! we make the sum to sum over all RI basis functions + CALL para_env%sum(BZ) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param AZ ... +!> \param Z_vectors ... +!> \param B_iaQ_bse_local ... +!> \param B_bar_ijQ_bse_local ... +!> \param B_abQ_bse_local ... +!> \param M_ia_tmp ... +!> \param RI_vector ... +!> \param Eigenval ... +!> \param homo ... +!> \param virtual ... +!> \param num_Z_vectors ... +!> \param local_RI_size ... +!> \param para_env ... +!> \param bse_spin_config ... +!> \param z_space_energy_cutoff ... +!> \param i_iter ... +!> \param bse_full_diag_debug ... +!> \param Full_exc_spectrum ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, M_ia_tmp, & + RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, & + para_env, bse_spin_config, z_space_energy_cutoff, i_iter, bse_full_diag_debug, & + Full_exc_spectrum, unit_nr) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, Z_vectors, B_iaQ_bse_local, & + B_bar_ijQ_bse_local, B_abQ_bse_local + REAL(KIND=dp), DIMENSION(:, :) :: M_ia_tmp, RI_vector + REAL(KIND=dp), DIMENSION(:) :: Eigenval + INTEGER :: homo, virtual, num_Z_vectors, & + local_RI_size + TYPE(mp_para_env_type), INTENT(IN) :: para_env + INTEGER :: bse_spin_config + REAL(KIND=dp) :: z_space_energy_cutoff + INTEGER :: i_iter + LOGICAL :: bse_full_diag_debug + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Full_exc_spectrum + INTEGER :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_AZ' + + INTEGER :: a, a_virt, b, diag_info, handle, i, & + i_occ, i_Z_vector, j, LLL, LWORK, m, n + REAL(KIND=dp) :: eigen_diff + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: WORK + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A_full_reshaped + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: A_full, v_iajb, W_ijab + + CALL timeset(routineN, handle) + AZ(:, :, :) = 0.0_dp + + IF (i_iter == 1 .AND. bse_full_diag_debug) THEN + ALLOCATE (W_ijab(homo, homo, virtual, virtual)) + ALLOCATE (A_full(homo, virtual, homo, virtual)) + ALLOCATE (A_full_reshaped(homo*virtual, homo*virtual)) + ALLOCATE (Full_exc_spectrum(homo*virtual)) + W_ijab = 0.0_dp + A_full = 0.0_dp + A_full_reshaped = 0.0_dp + Full_exc_spectrum = 0.0_dp + END IF + + CALL compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, & + num_Z_vectors, homo, virtual, bse_spin_config, v_iajb, bse_full_diag_debug, i_iter, & + para_env) + + DO i_Z_vector = 1, num_Z_vectors + + ! JW TO DO: OMP PARALLELIZATION + DO LLL = 1, local_RI_size + + ! M_ja^P = sum_b Z_jb*B_ba^P + CALL DGEMM("N", "N", homo, virtual, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, & + B_abQ_bse_local(:, :, LLL), virtual, 0.0_dp, M_ia_tmp, homo) + + ! (AZ)_ia = sum_jP B_bar_ij^P*M_ja^P + CALL DGEMM("N", "N", homo, virtual, homo, -1.0_dp, B_bar_ijQ_bse_local(:, :, LLL), homo, & + M_ia_tmp, homo, 1.0_dp, AZ(:, :, i_Z_vector), homo) + + END DO + END DO + + IF (i_iter == 1 .AND. bse_full_diag_debug) THEN + W_ijab = 0.0_dp + !Create screened 4c integrals for check + DO LLL = 1, local_RI_size + DO i = 1, homo + DO j = 1, homo + DO a = 1, virtual + DO b = 1, virtual + W_ijab(i, j, a, b) = W_ijab(i, j, a, b) + B_bar_ijQ_bse_local(i, j, LLL)*B_abQ_bse_local(a, b, LLL) + END DO + END DO + END DO + END DO + END DO + ! we make the mp_sum to sum over all RI basis functions + CALL para_env%sum(W_ijab) + END IF + + ! we make the mp_sum to sum over all RI basis functions + CALL para_env%sum(AZ) + + ! add (e_a-e_i)*Z_ia + DO i_occ = 1, homo + DO a_virt = 1, virtual + + eigen_diff = Eigenval(a_virt + homo) - Eigenval(i_occ) + IF (unit_nr > 0 .AND. i_iter == 1) THEN + WRITE (unit_nr, *) "Ediff at (i_occ,a_virt)=", i_occ, a_virt, " is: ", eigen_diff + END IF + + AZ(i_occ, a_virt, :) = AZ(i_occ, a_virt, :) + Z_vectors(i_occ, a_virt, :)*eigen_diff + + END DO + END DO + + !cut off contributions, which are too high in the excitation spectrum + IF (z_space_energy_cutoff > 0) THEN + DO i_occ = 1, homo + DO a_virt = 1, virtual + + IF (Eigenval(a_virt + homo) > z_space_energy_cutoff .OR. -Eigenval(i_occ) > z_space_energy_cutoff) THEN + AZ(i_occ, a_virt, :) = 0 + END IF + + END DO + END DO + END IF + + !Debugging purposes: full diagonalization of A + IF (i_iter == 1 .AND. bse_full_diag_debug) THEN + n = 0 + DO i = 1, homo + DO a = 1, virtual + n = n + 1 + m = 0 + DO j = 1, homo + DO b = 1, virtual + m = m + 1 + IF (a == b .AND. i == j) THEN + eigen_diff = Eigenval(a + homo) - Eigenval(i) + ELSE + eigen_diff = 0 + END IF + A_full_reshaped(n, m) = eigen_diff + 2*v_iajb(i, a, j, b) - W_ijab(i, j, a, b) + A_full(i, a, j, b) = eigen_diff + 2*v_iajb(i, a, j, b) - W_ijab(i, j, a, b) + END DO + END DO + END DO + END DO + + !MG to do: Check for symmetry of ZAZ! + ALLOCATE (WORK(1)) + WORK = 0.0_dp + CALL DSYEV("N", "U", homo*virtual, A_full_reshaped, homo*virtual, Full_exc_spectrum, WORK, -1, diag_info) + LWORK = INT(WORK(1)) + DEALLOCATE (WORK) + ALLOCATE (WORK(LWORK)) + WORK = 0.0_dp + !MG to check: Usage of symmetric routine okay? (Correct LWORK?) + CALL DSYEV("N", "U", homo*virtual, A_full_reshaped, homo*virtual, Full_exc_spectrum, WORK, LWORK, diag_info) + + DEALLOCATE (WORK) + + DEALLOCATE (W_ijab, v_iajb, A_full, A_full_reshaped) + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param AZ ... +!> \param Z_vectors ... +!> \param B_iaQ_bse_local ... +!> \param RI_vector ... +!> \param local_RI_size ... +!> \param num_Z_vectors ... +!> \param homo ... +!> \param virtual ... +!> \param bse_spin_config ... +!> \param v_iajb ... +!> \param bse_full_diag_debug ... +!> \param i_iter ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, & + num_Z_vectors, homo, virtual, bse_spin_config, v_iajb, bse_full_diag_debug, i_iter, & + para_env) + + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: AZ, Z_vectors, B_iaQ_bse_local + REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: RI_vector + INTEGER, INTENT(IN) :: local_RI_size, num_Z_vectors, homo, & + virtual, bse_spin_config + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: v_iajb + LOGICAL :: bse_full_diag_debug + INTEGER, INTENT(IN) :: i_iter + TYPE(mp_para_env_type), INTENT(IN) :: para_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_ia_jb_part' + + INTEGER :: a, a_virt, b, handle, i, i_occ, & + i_Z_vector, j, LLL + REAL(KIND=dp) :: alpha + +!debugging: + + CALL timeset(routineN, handle) + + !Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F) + SELECT CASE (bse_spin_config) + CASE (bse_singlet) + alpha = 2.0_dp + CASE (bse_triplet) + alpha = 0.0_dp + END SELECT + + RI_vector = 0.0_dp + + ! v_P = sum_jb B_jb^P Z_jb + DO LLL = 1, local_RI_size + DO i_Z_vector = 1, num_Z_vectors + DO i_occ = 1, homo + DO a_virt = 1, virtual + + RI_vector(LLL, i_Z_vector) = RI_vector(LLL, i_Z_vector) + & + Z_vectors(i_occ, a_virt, i_Z_vector)* & + B_iaQ_bse_local(i_occ, a_virt, LLL) + + END DO + END DO + END DO + END DO + + ! AZ = sum_P B_ia^P*v_P + ... + DO LLL = 1, local_RI_size + DO i_Z_vector = 1, num_Z_vectors + DO i_occ = 1, homo + DO a_virt = 1, virtual + !MG to check: Minus sign at v oder W? Factor for triplet/singlet + AZ(i_occ, a_virt, i_Z_vector) = AZ(i_occ, a_virt, i_Z_vector) + & + alpha*RI_vector(LLL, i_Z_vector)* & + B_iaQ_bse_local(i_occ, a_virt, LLL) + + END DO + END DO + END DO + END DO + IF (i_iter == 1 .AND. bse_full_diag_debug) THEN + ALLOCATE (v_iajb(homo, virtual, homo, virtual)) + v_iajb = 0.0_dp + !Create unscreened 4c integrals for check + DO LLL = 1, local_RI_size + DO i = 1, homo + DO j = 1, homo + DO a = 1, virtual + DO b = 1, virtual + v_iajb(i, a, j, b) = v_iajb(i, a, j, b) + B_iaQ_bse_local(i, a, LLL)*B_iaQ_bse_local(j, b, LLL) + END DO + END DO + END DO + END DO + END DO + ! we make the mp_sum to sum over all RI basis functions + CALL para_env%sum(v_iajb) + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ...Eigenval +!> \param Z_vectors ... +!> \param Eigenval ... +!> \param num_Z_vectors ... +!> \param homo ... +!> \param virtual ... +! ************************************************************************************************** + SUBROUTINE initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual) + + REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: Z_vectors + REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval + INTEGER, INTENT(IN) :: num_Z_vectors, homo, virtual + + CHARACTER(LEN=*), PARAMETER :: routineN = 'initial_guess_Z_vectors' + + INTEGER :: a_virt, handle, i_occ, i_Z_vector, & + min_loc(2) + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigen_diff_ia + + CALL timeset(routineN, handle) + + ALLOCATE (eigen_diff_ia(homo, virtual)) + + DO i_occ = 1, homo + DO a_virt = 1, virtual + eigen_diff_ia(i_occ, a_virt) = Eigenval(a_virt + homo) - Eigenval(i_occ) + END DO + END DO + + DO i_Z_vector = 1, num_Z_vectors + + min_loc = MINLOC(eigen_diff_ia) + + Z_vectors(min_loc(1), min_loc(2), i_Z_vector) = 1.0_dp + + eigen_diff_ia(min_loc(1), min_loc(2)) = 1.0E20_dp + + END DO + + DEALLOCATE (eigen_diff_ia) + + CALL timestop(handle) + + END SUBROUTINE + + ! ************************************************************************************************** +!> \brief ... +!> \param fm_mat_S_ab_bse ... +!> \param fm_mat_S ... +!> \param fm_mat_S_bar_ia_bse ... +!> \param fm_mat_S_bar_ij_bse ... +!> \param B_bar_ijQ_bse_local ... +!> \param B_abQ_bse_local ... +!> \param B_bar_iaQ_bse_local ... +!> \param B_iaQ_bse_local ... +!> \param dimen_RI ... +!> \param homo ... +!> \param virtual ... +!> \param gd_array ... +!> \param color_sub ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE fill_local_3c_arrays(fm_mat_S_ab_bse, fm_mat_S, & + fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, & + B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & + B_iaQ_bse_local, dimen_RI, homo, virtual, & + gd_array, color_sub, para_env) + + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ab_bse, fm_mat_S, & + fm_mat_S_bar_ia_bse, & + fm_mat_S_bar_ij_bse + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & + INTENT(OUT) :: B_bar_ijQ_bse_local, B_abQ_bse_local, & + B_bar_iaQ_bse_local, B_iaQ_bse_local + INTEGER, INTENT(IN) :: dimen_RI, homo, virtual + TYPE(group_dist_d1_type), INTENT(IN) :: gd_array + INTEGER, INTENT(IN) :: color_sub + TYPE(mp_para_env_type), INTENT(IN) :: para_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_local_3c_arrays' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + CALL allocate_and_fill_local_array(B_iaQ_bse_local, fm_mat_S, gd_array, color_sub, homo, virtual, dimen_RI, para_env) + + CALL allocate_and_fill_local_array(B_bar_iaQ_bse_local, fm_mat_S_bar_ia_bse, gd_array, color_sub, homo, virtual, & + dimen_RI, para_env) + + CALL allocate_and_fill_local_array(B_bar_ijQ_bse_local, fm_mat_S_bar_ij_bse, gd_array, color_sub, homo, homo, & + dimen_RI, para_env) + + CALL allocate_and_fill_local_array(B_abQ_bse_local, fm_mat_S_ab_bse, gd_array, color_sub, virtual, virtual, & + dimen_RI, para_env) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param B_local ... +!> \param fm_mat_S ... +!> \param gd_array ... +!> \param color_sub ... +!> \param small_size ... +!> \param big_size ... +!> \param dimen_RI ... +!> \param para_env ... +! ************************************************************************************************** + SUBROUTINE allocate_and_fill_local_array(B_local, fm_mat_S, gd_array, & + color_sub, small_size, big_size, dimen_RI, para_env) + + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & + INTENT(OUT) :: B_local + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S + TYPE(group_dist_d1_type), INTENT(IN) :: gd_array + INTEGER, INTENT(IN) :: color_sub, small_size, big_size, dimen_RI + TYPE(mp_para_env_type), INTENT(IN) :: para_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_and_fill_local_array' + + INTEGER :: combi_index, end_RI, handle, handle1, i_comm, i_entry, iiB, imepos, jjB, & + level_big_size, level_small_size, ncol_local, nrow_local, num_comm_cycles, RI_index, & + size_RI, start_RI + INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, mepos_from_RI_index, & + num_entries_rec, num_entries_send + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + REAL(KIND=dp) :: matrix_el + TYPE(integ_mat_buffer_type), ALLOCATABLE, & + DIMENSION(:) :: buffer_rec, buffer_send + TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array + + CALL timeset(routineN, handle) + + ALLOCATE (mepos_from_RI_index(dimen_RI)) + mepos_from_RI_index = 0 + + DO imepos = 0, para_env%num_pe - 1 + + CALL get_group_dist(gd_array, pos=imepos, starts=start_RI, ends=end_RI) + + mepos_from_RI_index(start_RI:end_RI) = imepos + + END DO + + ! color_sub is automatically the number of the process since every subgroup has only one MPI rank + CALL get_group_dist(gd_array, color_sub, start_RI, end_RI, size_RI) + + ALLOCATE (B_local(small_size, big_size, 1:size_RI)) + + ALLOCATE (num_entries_send(0:para_env%num_pe - 1)) + ALLOCATE (num_entries_rec(0:para_env%num_pe - 1)) + + ALLOCATE (req_array(1:para_env%num_pe, 4)) + + ALLOCATE (entry_counter(0:para_env%num_pe - 1)) + + CALL cp_fm_get_info(matrix=fm_mat_S, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + num_comm_cycles = 10 + + ! communicate not all due to huge memory overhead, since for every number in fm_mat_S, we store + ! three additional ones (RI index, first MO index, second MO index!!) + DO i_comm = 0, num_comm_cycles - 1 + + num_entries_send = 0 + num_entries_rec = 0 + + ! loop over RI index to get the number of sent entries + DO jjB = 1, nrow_local + + RI_index = row_indices(jjB) + + IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE + + imepos = mepos_from_RI_index(RI_index) + + num_entries_send(imepos) = num_entries_send(imepos) + ncol_local + + END DO + + CALL para_env%alltoall(num_entries_send, num_entries_rec, 1) + + ALLOCATE (buffer_rec(0:para_env%num_pe - 1)) + ALLOCATE (buffer_send(0:para_env%num_pe - 1)) + + ! allocate data message and corresponding indices + DO imepos = 0, para_env%num_pe - 1 + + ALLOCATE (buffer_rec(imepos)%msg(num_entries_rec(imepos))) + buffer_rec(imepos)%msg = 0.0_dp + + ALLOCATE (buffer_send(imepos)%msg(num_entries_send(imepos))) + buffer_send(imepos)%msg = 0.0_dp + + ALLOCATE (buffer_rec(imepos)%indx(num_entries_rec(imepos), 3)) + buffer_rec(imepos)%indx = 0 + + ALLOCATE (buffer_send(imepos)%indx(num_entries_send(imepos), 3)) + buffer_send(imepos)%indx = 0 + + END DO + + entry_counter(:) = 0 + + ! loop over RI index for filling the send-buffer + DO jjB = 1, nrow_local + + RI_index = row_indices(jjB) + + IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE + + imepos = mepos_from_RI_index(RI_index) + + DO iiB = 1, ncol_local + + combi_index = col_indices(iiB) + level_small_size = MAX(1, combi_index - 1)/MAX(big_size, 2) + 1 + level_big_size = combi_index - (level_small_size - 1)*big_size + + entry_counter(imepos) = entry_counter(imepos) + 1 + + buffer_send(imepos)%msg(entry_counter(imepos)) = fm_mat_S%local_data(jjB, iiB) + + buffer_send(imepos)%indx(entry_counter(imepos), 1) = RI_index + buffer_send(imepos)%indx(entry_counter(imepos), 2) = level_small_size + buffer_send(imepos)%indx(entry_counter(imepos), 3) = level_big_size + + END DO + + END DO + + CALL timeset("BSE_comm_data", handle1) + + CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array) + + CALL timestop(handle1) + + ! fill B_local + DO imepos = 0, para_env%num_pe - 1 + + DO i_entry = 1, num_entries_rec(imepos) + + RI_index = buffer_rec(imepos)%indx(i_entry, 1) - start_RI + 1 + level_small_size = buffer_rec(imepos)%indx(i_entry, 2) + level_big_size = buffer_rec(imepos)%indx(i_entry, 3) + + matrix_el = buffer_rec(imepos)%msg(i_entry) + + B_local(level_small_size, level_big_size, RI_index) = matrix_el + + END DO + + END DO + + DO imepos = 0, para_env%num_pe - 1 + DEALLOCATE (buffer_send(imepos)%msg) + DEALLOCATE (buffer_send(imepos)%indx) + DEALLOCATE (buffer_rec(imepos)%msg) + DEALLOCATE (buffer_rec(imepos)%indx) + END DO + + DEALLOCATE (buffer_rec, buffer_send) + + END DO + + DEALLOCATE (num_entries_send, num_entries_rec) + + DEALLOCATE (mepos_from_RI_index) + + DEALLOCATE (entry_counter, req_array) + + CALL timestop(handle) + + END SUBROUTINE + +END MODULE bse_iterative diff --git a/src/bse_main.F b/src/bse_main.F new file mode 100644 index 0000000000..a9688c6609 --- /dev/null +++ b/src/bse_main.F @@ -0,0 +1,221 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Main routines for GW + Bethe-Salpeter for computing electronic excitations +!> \par History +!> 04.2024 created [Maximilian Graml] +! ************************************************************************************************** + +MODULE bse_main + + USE bse_full_diag, ONLY: create_A,& + create_B,& + create_hermitian_form_of_ABBA,& + diagonalize_A,& + diagonalize_C + USE bse_iterative, ONLY: do_subspace_iterations,& + fill_local_3c_arrays + USE bse_util, ONLY: deallocate_matrices_bse,& + estimate_BSE_resources,& + mult_B_with_W,& + print_BSE_start_flag,& + truncate_BSE_matrices + USE cp_fm_types, ONLY: cp_fm_release,& + cp_fm_type + USE group_dist_types, ONLY: group_dist_d1_type + USE input_constants, ONLY: bse_abba,& + bse_both,& + bse_fulldiag,& + bse_iterdiag,& + bse_tda + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type + USE mp2_types, ONLY: mp2_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main' + + PUBLIC :: start_bse_calculation + +CONTAINS + +! ************************************************************************************************** +!> \brief Main subroutine managing BSE calculations +!> \param fm_mat_S_ia_bse ... +!> \param fm_mat_S_ij_bse ... +!> \param fm_mat_S_ab_bse ... +!> \param fm_mat_Q_static_bse ... +!> \param fm_mat_Q_static_bse_gemm ... +!> \param Eigenval ... +!> \param homo ... +!> \param virtual ... +!> \param dimen_RI ... +!> \param dimen_RI_red ... +!> \param gd_array ... +!> \param color_sub ... +!> \param mp2_env ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, & + fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, & + Eigenval, homo, virtual, dimen_RI, dimen_RI_red, & + gd_array, color_sub, mp2_env, unit_nr) + + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, & + fm_mat_S_ab_bse + TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q_static_bse, & + fm_mat_Q_static_bse_gemm + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), & + INTENT(IN) :: Eigenval + INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual + INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red + TYPE(group_dist_d1_type), INTENT(IN) :: gd_array + INTEGER, INTENT(IN) :: color_sub + TYPE(mp2_type) :: mp2_env + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation' + + INTEGER :: handle, homo_red, virtual_red + LOGICAL :: my_do_abba, my_do_fulldiag, & + my_do_iterat_diag, my_do_tda + REAL(KIND=dp) :: diag_runtime_est + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, & + B_bar_ijQ_bse_local, B_iaQ_bse_local + TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, fm_inv_sqrt_A_minus_B, fm_mat_S_ab_trunc, & + fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, & + fm_sqrt_A_minus_B + TYPE(mp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + para_env => fm_mat_S_ia_bse%matrix_struct%para_env + + my_do_fulldiag = .FALSE. + my_do_iterat_diag = .FALSE. + my_do_tda = .FALSE. + my_do_abba = .FALSE. + !Method: Iterative or full diagonalization + SELECT CASE (mp2_env%ri_g0w0%bse_diag_method) + CASE (bse_iterdiag) + my_do_iterat_diag = .TRUE. + !MG: Basics of the Davidson solver are implemented, but not rigorously checked. + CPABORT("Iterative BSE not yet implemented") + CASE (bse_fulldiag) + my_do_fulldiag = .TRUE. + END SELECT + !Approximation: TDA and/or full ABBA matrix + SELECT CASE (mp2_env%ri_g0w0%bse_approx) + CASE (bse_tda) + my_do_tda = .TRUE. + CASE (bse_abba) + my_do_abba = .TRUE. + CASE (bse_both) + my_do_tda = .TRUE. + my_do_abba = .TRUE. + END SELECT + + CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr) + + CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync() + + ! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices + CALL truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, & + fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, & + Eigenval(:, 1, 1), Eigenval_reduced, & + homo(1), virtual(1), dimen_RI, unit_nr, & + homo_red, virtual_red, & + mp2_env) + ! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs) + ! r,s: MO-index, P,R,T: RI-index + ! B: fm_mat_S_..., W: fm_mat_Q_... + CALL mult_B_with_W(fm_mat_S_ij_trunc, fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, & + fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, & + dimen_RI_red, homo_red, virtual_red) + + IF (my_do_iterat_diag) THEN + CALL fill_local_3c_arrays(fm_mat_S_ab_trunc, fm_mat_S_ia_trunc, & + fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, & + B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & + B_iaQ_bse_local, dimen_RI_red, homo_red, virtual_red, & + gd_array, color_sub, para_env) + END IF + + IF (my_do_fulldiag) THEN + ! Quick estimate of memory consumption and runtime of diagonalizations + CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, & + para_env, diag_runtime_est) + ! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W) + ! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab + ! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced + ! α is a spin-dependent factor + ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction) + ! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction) + CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, & + fm_A_BSE, Eigenval_reduced, unit_nr, & + homo_red, virtual_red, dimen_RI, mp2_env, & + para_env) + IF (my_do_abba) THEN + ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W) + ! B_ia,jb = α * v_ia,jb - W_ib,aj + ! α is a spin-dependent factor + ! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction) + ! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction) + CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, fm_B_BSE, & + homo_red, virtual_red, dimen_RI, unit_nr, mp2_env) + ! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem + ! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)). + ! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions + ! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001). + CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, & + fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, & + homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est) + END IF + CALL cp_fm_release(fm_B_BSE) + IF (my_do_tda) THEN + ! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n + CALL diagonalize_A(fm_A_BSE, homo_red, virtual_red, homo(1), & + unit_nr, diag_runtime_est, mp2_env) + END IF + ! Release to avoid faulty use of changed A matrix + CALL cp_fm_release(fm_A_BSE) + IF (my_do_abba) THEN + ! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n . + ! Here, the eigenvectors Z^n relate to X^n via + ! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001). + CALL diagonalize_C(fm_C_BSE, homo_red, virtual_red, homo(1), & + fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, & + unit_nr, diag_runtime_est, mp2_env) + END IF + ! Release to avoid faulty use of changed C matrix + CALL cp_fm_release(fm_C_BSE) + END IF + + CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, & + fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, & + fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm) + DEALLOCATE (Eigenval_reduced) + IF (my_do_iterat_diag) THEN + ! Contains untested Block-Davidson algorithm + CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & + B_iaQ_bse_local, homo(1), virtual(1), mp2_env%ri_g0w0%bse_spin_config, unit_nr, & + Eigenval(:, 1, 1), para_env, mp2_env) + ! Deallocate local 3c-B-matrices + DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local) + END IF + + CALL timestop(handle) + + END SUBROUTINE start_bse_calculation + +END MODULE bse_main diff --git a/src/bse_util.F b/src/bse_util.F new file mode 100644 index 0000000000..7b7546f93c --- /dev/null +++ b/src/bse_util.F @@ -0,0 +1,1351 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations +!> \par History +!> 11.2023 created [Maximilian Graml] +! ************************************************************************************************** +MODULE bse_util + USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_fm_basic_linalg, ONLY: cp_fm_upper_to_full + USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& + cp_fm_cholesky_invert + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type + USE cp_fm_types, ONLY: cp_fm_create,& + cp_fm_get_info,& + cp_fm_indxg2l,& + cp_fm_indxg2p,& + cp_fm_release,& + cp_fm_set_all,& + cp_fm_to_fm_submat_general,& + cp_fm_type + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type,& + mp_request_type + USE mp2_types, ONLY: integ_mat_buffer_type,& + mp2_type + USE parallel_gemm_api, ONLY: parallel_gemm + USE physcon, ONLY: evolt + USE rpa_communication, ONLY: communicate_buffer + USE util, ONLY: sort,& + sort_unique +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_util' + + PUBLIC :: mult_B_with_W, fm_general_add_bse, truncate_fm, fm_write_thresh, print_BSE_start_flag, & + deallocate_matrices_bse, comp_eigvec_coeff_BSE, sort_excitations, & + estimate_BSE_resources, filter_eigvec_contrib, truncate_BSE_matrices + +CONTAINS + +! ************************************************************************************************** +!> \brief Multiplies B-matrix (RI-3c-Integrals) with W (screening) to obtain \bar{B} +!> \param fm_mat_S_ij_bse ... +!> \param fm_mat_S_ia_bse ... +!> \param fm_mat_S_bar_ia_bse ... +!> \param fm_mat_S_bar_ij_bse ... +!> \param fm_mat_Q_static_bse_gemm ... +!> \param dimen_RI ... +!> \param homo ... +!> \param virtual ... +! ************************************************************************************************** + SUBROUTINE mult_B_with_W(fm_mat_S_ij_bse, fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, & + fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, & + dimen_RI, homo, virtual) + + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ij_bse, fm_mat_S_ia_bse + TYPE(cp_fm_type), INTENT(OUT) :: fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q_static_bse_gemm + INTEGER, INTENT(IN) :: dimen_RI, homo, virtual + + CHARACTER(LEN=*), PARAMETER :: routineN = 'mult_B_with_W' + + INTEGER :: handle, i_global, iiB, info_chol, & + j_global, jjB, ncol_local, nrow_local + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + TYPE(cp_fm_type) :: fm_work + + CALL timeset(routineN, handle) + + CALL cp_fm_create(fm_mat_S_bar_ia_bse, fm_mat_S_ia_bse%matrix_struct) + CALL cp_fm_set_all(fm_mat_S_bar_ia_bse, 0.0_dp) + + CALL cp_fm_create(fm_mat_S_bar_ij_bse, fm_mat_S_ij_bse%matrix_struct) + CALL cp_fm_set_all(fm_mat_S_bar_ij_bse, 0.0_dp) + + CALL cp_fm_create(fm_work, fm_mat_Q_static_bse_gemm%matrix_struct) + CALL cp_fm_set_all(fm_work, 0.0_dp) + + ! get info of fm_mat_Q_static_bse and compute ((1+Q(0))^-1-1) + CALL cp_fm_get_info(matrix=fm_mat_Q_static_bse_gemm, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + DO jjB = 1, ncol_local + j_global = col_indices(jjB) + DO iiB = 1, nrow_local + i_global = row_indices(iiB) + IF (j_global == i_global .AND. i_global <= dimen_RI) THEN + fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) = fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) + 1.0_dp + END IF + END DO + END DO + + ! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition + CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q_static_bse_gemm, n=dimen_RI, info_out=info_chol) + + CPASSERT(info_chol == 0) + + ! calculate [1+Q(i0)]^-1 + CALL cp_fm_cholesky_invert(fm_mat_Q_static_bse_gemm) + + ! symmetrize the result + CALL cp_fm_upper_to_full(fm_mat_Q_static_bse_gemm, fm_work) + + CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=homo**2, k=dimen_RI, alpha=1.0_dp, & + matrix_a=fm_mat_Q_static_bse_gemm, matrix_b=fm_mat_S_ij_bse, beta=0.0_dp, & + matrix_c=fm_mat_S_bar_ij_bse) + + ! fm_mat_S_bar_ia_bse has a different blacs_env as fm_mat_S_ij_bse since we take + ! fm_mat_S_ia_bse from RPA. Therefore, we also need a different fm_mat_Q_static_bse_gemm + CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=homo*virtual, k=dimen_RI, alpha=1.0_dp, & + matrix_a=fm_mat_Q_static_bse_gemm, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, & + matrix_c=fm_mat_S_bar_ia_bse) + + CALL cp_fm_release(fm_work) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,ab +!> to A_ia, which needs MPI communication. +!> \param fm_out ... +!> \param fm_in ... +!> \param beta ... +!> \param nrow_secidx_in ... +!> \param ncol_secidx_in ... +!> \param nrow_secidx_out ... +!> \param ncol_secidx_out ... +!> \param unit_nr ... +!> \param reordering ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, & + nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env) + + TYPE(cp_fm_type), INTENT(INOUT) :: fm_out, fm_in + REAL(kind=dp) :: beta + INTEGER, INTENT(IN) :: nrow_secidx_in, ncol_secidx_in, & + nrow_secidx_out, ncol_secidx_out + INTEGER :: unit_nr + INTEGER, DIMENSION(4) :: reordering + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_general_add_bse' + + INTEGER :: col_idx_loc, dummy, handle, handle2, i_entry_rec, idx_col_out, idx_row_out, ii, & + iproc, jj, ncol_block_in, ncol_block_out, ncol_local_in, ncol_local_out, npcol, nprocs, & + nprow, nrow_block_in, nrow_block_out, nrow_local_in, nrow_local_out, proc_send, & + row_idx_loc, send_pcol, send_prow + INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, & + num_entries_send + INTEGER, DIMENSION(4) :: indices_in + INTEGER, DIMENSION(:), POINTER :: col_indices_in, col_indices_out, & + row_indices_in, row_indices_out + TYPE(integ_mat_buffer_type), ALLOCATABLE, & + DIMENSION(:) :: buffer_rec, buffer_send + TYPE(mp_para_env_type), POINTER :: para_env_out + TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array + + CALL timeset(routineN, handle) + CALL timeset(routineN//"_1_setup", handle2) + + para_env_out => fm_out%matrix_struct%para_env + ! A_iajb + ! We start by moving data from local parts of W_ijab to the full matrix A_iajb using buffers + CALL cp_fm_get_info(matrix=fm_out, & + nrow_local=nrow_local_out, & + ncol_local=ncol_local_out, & + row_indices=row_indices_out, & + col_indices=col_indices_out, & + nrow_block=nrow_block_out, & + ncol_block=ncol_block_out) + + nprow = fm_out%matrix_struct%context%num_pe(1) + npcol = fm_out%matrix_struct%context%num_pe(2) + + ALLOCATE (num_entries_rec(0:para_env_out%num_pe - 1)) + ALLOCATE (num_entries_send(0:para_env_out%num_pe - 1)) + + num_entries_rec(:) = 0 + num_entries_send(:) = 0 + + dummy = 0 + + CALL cp_fm_get_info(matrix=fm_in, & + nrow_local=nrow_local_in, & + ncol_local=ncol_local_in, & + row_indices=row_indices_in, & + col_indices=col_indices_in, & + nrow_block=nrow_block_in, & + ncol_block=ncol_block_in) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, & + fm_out%matrix_struct%nrow_global + WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, & + fm_out%matrix_struct%ncol_global + + WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_out%name, nrow_block_out + WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_out%name, ncol_block_out + + WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_in%name, & + fm_in%matrix_struct%nrow_global + WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_in%name, & + fm_in%matrix_struct%ncol_global + + WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_in%name, nrow_block_in + WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_in%name, ncol_block_in + END IF + + ! Use scalapack wrapper to find process index in fm_out + ! To that end, we obtain the global index in fm_out from the level indices + indices_in(:) = 0 + DO row_idx_loc = 1, nrow_local_in + indices_in(1) = (row_indices_in(row_idx_loc) - 1)/nrow_secidx_in + 1 + indices_in(2) = MOD(row_indices_in(row_idx_loc) - 1, nrow_secidx_in) + 1 + DO col_idx_loc = 1, ncol_local_in + indices_in(3) = (col_indices_in(col_idx_loc) - 1)/ncol_secidx_in + 1 + indices_in(4) = MOD(col_indices_in(col_idx_loc) - 1, ncol_secidx_in) + 1 + + idx_row_out = indices_in(reordering(2)) + (indices_in(reordering(1)) - 1)*nrow_secidx_out + idx_col_out = indices_in(reordering(4)) + (indices_in(reordering(3)) - 1)*ncol_secidx_out + + send_prow = cp_fm_indxg2p(idx_row_out, nrow_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(1), nprow) + + send_pcol = cp_fm_indxg2p(idx_col_out, ncol_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(2), npcol) + + proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol) + + num_entries_send(proc_send) = num_entries_send(proc_send) + 1 + + END DO + END DO + + CALL timestop(handle2) + + CALL timeset(routineN//"_2_comm_entry_nums", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A27)') 'BSE|DEBUG|', 'Communicating entry numbers' + END IF + + CALL para_env_out%alltoall(num_entries_send, num_entries_rec, 1) + + CALL timestop(handle2) + + CALL timeset(routineN//"_3_alloc_buffer", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A18)') 'BSE|DEBUG|', 'Allocating buffers' + END IF + + ! Buffers for entries and their indices + ALLOCATE (buffer_rec(0:para_env_out%num_pe - 1)) + ALLOCATE (buffer_send(0:para_env_out%num_pe - 1)) + + ! allocate data message and corresponding indices + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc))) + buffer_rec(iproc)%msg = 0.0_dp + + END DO + + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc))) + buffer_send(iproc)%msg = 0.0_dp + + END DO + + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2)) + buffer_rec(iproc)%indx = 0 + + END DO + + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2)) + buffer_send(iproc)%indx = 0 + + END DO + + CALL timestop(handle2) + + CALL timeset(routineN//"_4_buf_from_fmin_"//fm_out%name, handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers' + END IF + + ALLOCATE (entry_counter(0:para_env_out%num_pe - 1)) + entry_counter(:) = 0 + + ! Now we can write the actual data and indices to the send-buffer + DO row_idx_loc = 1, nrow_local_in + indices_in(1) = (row_indices_in(row_idx_loc) - 1)/nrow_secidx_in + 1 + indices_in(2) = MOD(row_indices_in(row_idx_loc) - 1, nrow_secidx_in) + 1 + DO col_idx_loc = 1, ncol_local_in + indices_in(3) = (col_indices_in(col_idx_loc) - 1)/ncol_secidx_in + 1 + indices_in(4) = MOD(col_indices_in(col_idx_loc) - 1, ncol_secidx_in) + 1 + + idx_row_out = indices_in(reordering(2)) + (indices_in(reordering(1)) - 1)*nrow_secidx_out + idx_col_out = indices_in(reordering(4)) + (indices_in(reordering(3)) - 1)*ncol_secidx_out + + send_prow = cp_fm_indxg2p(idx_row_out, nrow_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(1), nprow) + + send_pcol = cp_fm_indxg2p(idx_col_out, ncol_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(2), npcol) + + proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol) + entry_counter(proc_send) = entry_counter(proc_send) + 1 + + buffer_send(proc_send)%msg(entry_counter(proc_send)) = & + fm_in%local_data(row_idx_loc, col_idx_loc) + + buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = idx_row_out + buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = idx_col_out + + END DO + END DO + + ALLOCATE (req_array(1:para_env_out%num_pe, 4)) + + CALL timestop(handle2) + + CALL timeset(routineN//"_5_comm_buffer", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A21)') 'BSE|DEBUG|', 'Communicating buffers' + END IF + + ! communicate the buffer + CALL communicate_buffer(para_env_out, num_entries_rec, num_entries_send, buffer_rec, & + buffer_send, req_array) + + CALL timestop(handle2) + + CALL timeset(routineN//"_6_buffer_to_fmout"//fm_out%name, handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name + END IF + + ! fill fm_out with the entries from buffer_rec, i.e. buffer_rec are parts of fm_in + nprocs = para_env_out%num_pe + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP SHARED(fm_out, nprocs, nrow_block_out, ncol_block_out, & +!$OMP num_entries_rec, buffer_rec, beta, dummy, nprow, npcol) & +!$OMP PRIVATE(iproc, i_entry_rec, ii, jj) + DO iproc = 0, nprocs - 1 + DO i_entry_rec = 1, num_entries_rec(iproc) + ii = cp_fm_indxg2l(buffer_rec(iproc)%indx(i_entry_rec, 1), nrow_block_out, & + dummy, dummy, nprow) + jj = cp_fm_indxg2l(buffer_rec(iproc)%indx(i_entry_rec, 2), ncol_block_out, & + dummy, dummy, npcol) + + fm_out%local_data(ii, jj) = fm_out%local_data(ii, jj) + beta*buffer_rec(iproc)%msg(i_entry_rec) + END DO + END DO +!$OMP END PARALLEL DO + + CALL timestop(handle2) + + CALL timeset(routineN//"_7_cleanup", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers' + END IF + + !Clean up all the arrays from the communication process + DO iproc = 0, para_env_out%num_pe - 1 + DEALLOCATE (buffer_rec(iproc)%msg) + DEALLOCATE (buffer_rec(iproc)%indx) + DEALLOCATE (buffer_send(iproc)%msg) + DEALLOCATE (buffer_send(iproc)%indx) + END DO + DEALLOCATE (buffer_rec, buffer_send) + DEALLOCATE (req_array) + DEALLOCATE (entry_counter) + DEALLOCATE (num_entries_rec, num_entries_send) + + CALL timestop(handle2) + CALL timestop(handle) + + END SUBROUTINE fm_general_add_bse + +! ************************************************************************************************** +!> \brief Routine for truncating a full matrix as given by the energy cutoffs in the input file. +!> Logic: Matrices have some dimension dimen_RI x nrow_in*ncol_in for the incoming (untruncated) matrix +!> and dimen_RI x nrow_out*ncol_out for the truncated matrix. The truncation is done by resorting the indices +!> via parallel communication. +!> \param fm_out ... +!> \param fm_in ... +!> \param ncol_in ... +!> \param nrow_out ... +!> \param ncol_out ... +!> \param unit_nr ... +!> \param mp2_env ... +!> \param nrow_offset ... +!> \param ncol_offset ... +! ************************************************************************************************** + SUBROUTINE truncate_fm(fm_out, fm_in, ncol_in, & + nrow_out, ncol_out, unit_nr, mp2_env, & + nrow_offset, ncol_offset) + + TYPE(cp_fm_type), INTENT(INOUT) :: fm_out + TYPE(cp_fm_type), INTENT(IN) :: fm_in + INTEGER :: ncol_in, nrow_out, ncol_out, unit_nr + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + INTEGER, INTENT(IN), OPTIONAL :: nrow_offset, ncol_offset + + CHARACTER(LEN=*), PARAMETER :: routineN = 'truncate_fm' + + INTEGER :: col_idx_loc, dummy, handle, handle2, i_entry_rec, idx_col_first, idx_col_in, & + idx_col_out, idx_col_sec, idx_row_in, ii, iproc, jj, ncol_block_in, ncol_block_out, & + ncol_local_in, ncol_local_out, npcol, nprocs, nprow, nrow_block_in, nrow_block_out, & + nrow_local_in, nrow_local_out, proc_send, row_idx_loc, send_pcol, send_prow + INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, & + num_entries_send + INTEGER, DIMENSION(:), POINTER :: col_indices_in, col_indices_out, & + row_indices_in, row_indices_out + LOGICAL :: correct_ncol, correct_nrow + TYPE(integ_mat_buffer_type), ALLOCATABLE, & + DIMENSION(:) :: buffer_rec, buffer_send + TYPE(mp_para_env_type), POINTER :: para_env_out + TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array + + CALL timeset(routineN, handle) + CALL timeset(routineN//"_1_setup", handle2) + + correct_nrow = .FALSE. + correct_ncol = .FALSE. + !In case of truncation in the occupied space, we need to correct the interval of indices + IF (PRESENT(nrow_offset)) THEN + correct_nrow = .TRUE. + END IF + IF (PRESENT(ncol_offset)) THEN + correct_ncol = .TRUE. + END IF + + para_env_out => fm_out%matrix_struct%para_env + + CALL cp_fm_get_info(matrix=fm_out, & + nrow_local=nrow_local_out, & + ncol_local=ncol_local_out, & + row_indices=row_indices_out, & + col_indices=col_indices_out, & + nrow_block=nrow_block_out, & + ncol_block=ncol_block_out) + + nprow = fm_out%matrix_struct%context%num_pe(1) + npcol = fm_out%matrix_struct%context%num_pe(2) + + ALLOCATE (num_entries_rec(0:para_env_out%num_pe - 1)) + ALLOCATE (num_entries_send(0:para_env_out%num_pe - 1)) + + num_entries_rec(:) = 0 + num_entries_send(:) = 0 + + dummy = 0 + + CALL cp_fm_get_info(matrix=fm_in, & + nrow_local=nrow_local_in, & + ncol_local=ncol_local_in, & + row_indices=row_indices_in, & + col_indices=col_indices_in, & + nrow_block=nrow_block_in, & + ncol_block=ncol_block_in) + + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_out%name, & + fm_out%matrix_struct%nrow_global + WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_out%name, & + fm_out%matrix_struct%ncol_global + + WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_out%name, nrow_block_out + WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_out%name, ncol_block_out + + WRITE (unit_nr, '(T2,A10,T13,A14,A10,T71,I10)') 'BSE|DEBUG|', 'Row number of ', fm_in%name, & + fm_in%matrix_struct%nrow_global + WRITE (unit_nr, '(T2,A10,T13,A17,A10,T71,I10)') 'BSE|DEBUG|', 'Column number of ', fm_in%name, & + fm_in%matrix_struct%ncol_global + + WRITE (unit_nr, '(T2,A10,T13,A18,A10,T71,I10)') 'BSE|DEBUG|', 'Row block size of ', fm_in%name, nrow_block_in + WRITE (unit_nr, '(T2,A10,T13,A21,A10,T71,I10)') 'BSE|DEBUG|', 'Column block size of ', fm_in%name, ncol_block_in + END IF + + ! We find global indices in S with nrow_in and ncol_in for truncation + DO col_idx_loc = 1, ncol_local_in + idx_col_in = col_indices_in(col_idx_loc) + + idx_col_first = (idx_col_in - 1)/ncol_in + 1 + idx_col_sec = MOD(idx_col_in - 1, ncol_in) + 1 + + ! If occupied orbitals are included, these have to be handled differently + ! due to their reversed indexing + IF (correct_nrow) THEN + idx_col_first = idx_col_first - nrow_offset + 1 + IF (idx_col_first .LE. 0) CYCLE + ELSE + IF (idx_col_first > nrow_out) EXIT + END IF + IF (correct_ncol) THEN + idx_col_sec = idx_col_sec - ncol_offset + 1 + IF (idx_col_sec .LE. 0) CYCLE + ELSE + IF (idx_col_sec > ncol_out) CYCLE + END IF + + idx_col_out = idx_col_sec + (idx_col_first - 1)*ncol_out + + DO row_idx_loc = 1, nrow_local_in + idx_row_in = row_indices_in(row_idx_loc) + + send_prow = cp_fm_indxg2p(idx_row_in, nrow_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(1), nprow) + + send_pcol = cp_fm_indxg2p(idx_col_out, ncol_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(2), npcol) + + proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol) + + num_entries_send(proc_send) = num_entries_send(proc_send) + 1 + + END DO + END DO + + CALL timestop(handle2) + + CALL timeset(routineN//"_2_comm_entry_nums", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A27)') 'BSE|DEBUG|', 'Communicating entry numbers' + END IF + + CALL para_env_out%alltoall(num_entries_send, num_entries_rec, 1) + + CALL timestop(handle2) + + CALL timeset(routineN//"_3_alloc_buffer", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A18)') 'BSE|DEBUG|', 'Allocating buffers' + END IF + + ! Buffers for entries and their indices + ALLOCATE (buffer_rec(0:para_env_out%num_pe - 1)) + ALLOCATE (buffer_send(0:para_env_out%num_pe - 1)) + + ! allocate data message and corresponding indices + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc))) + buffer_rec(iproc)%msg = 0.0_dp + + END DO + + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc))) + buffer_send(iproc)%msg = 0.0_dp + + END DO + + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2)) + buffer_rec(iproc)%indx = 0 + + END DO + + DO iproc = 0, para_env_out%num_pe - 1 + + ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2)) + buffer_send(iproc)%indx = 0 + + END DO + + CALL timestop(handle2) + + CALL timeset(routineN//"_4_buf_from_fmin_"//fm_out%name, handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A18,A10,A13)') 'BSE|DEBUG|', 'Writing data from ', fm_in%name, ' into buffers' + END IF + + ALLOCATE (entry_counter(0:para_env_out%num_pe - 1)) + entry_counter(:) = 0 + + ! Now we can write the actual data and indices to the send-buffer + DO col_idx_loc = 1, ncol_local_in + idx_col_in = col_indices_in(col_idx_loc) + + idx_col_first = (idx_col_in - 1)/ncol_in + 1 + idx_col_sec = MOD(idx_col_in - 1, ncol_in) + 1 + + ! If occupied orbitals are included, these have to be handled differently + ! due to their reversed indexing + IF (correct_nrow) THEN + idx_col_first = idx_col_first - nrow_offset + 1 + IF (idx_col_first .LE. 0) CYCLE + ELSE + IF (idx_col_first > nrow_out) EXIT + END IF + IF (correct_ncol) THEN + idx_col_sec = idx_col_sec - ncol_offset + 1 + IF (idx_col_sec .LE. 0) CYCLE + ELSE + IF (idx_col_sec > ncol_out) CYCLE + END IF + + idx_col_out = idx_col_sec + (idx_col_first - 1)*ncol_out + + DO row_idx_loc = 1, nrow_local_in + idx_row_in = row_indices_in(row_idx_loc) + + send_prow = cp_fm_indxg2p(idx_row_in, nrow_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(1), nprow) + + send_pcol = cp_fm_indxg2p(idx_col_out, ncol_block_out, dummy, & + fm_out%matrix_struct%first_p_pos(2), npcol) + + proc_send = fm_out%matrix_struct%context%blacs2mpi(send_prow, send_pcol) + entry_counter(proc_send) = entry_counter(proc_send) + 1 + + buffer_send(proc_send)%msg(entry_counter(proc_send)) = & + fm_in%local_data(row_idx_loc, col_idx_loc) + !No need to create row_out, since it is identical to incoming + !We dont change the RI index for any fm_mat_XX_BSE + buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = idx_row_in + buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = idx_col_out + + END DO + END DO + + ALLOCATE (req_array(1:para_env_out%num_pe, 4)) + + CALL timestop(handle2) + + CALL timeset(routineN//"_5_comm_buffer", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A21)') 'BSE|DEBUG|', 'Communicating buffers' + END IF + + ! communicate the buffer + CALL communicate_buffer(para_env_out, num_entries_rec, num_entries_send, buffer_rec, & + buffer_send, req_array) + + CALL timestop(handle2) + + CALL timeset(routineN//"_6_buffer_to_fmout"//fm_out%name, handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A24,A10)') 'BSE|DEBUG|', 'Writing from buffers to ', fm_out%name + END IF + + ! fill fm_out with the entries from buffer_rec, i.e. buffer_rec are parts of fm_in + nprocs = para_env_out%num_pe + +!$OMP PARALLEL DO DEFAULT(NONE) & +!$OMP SHARED(fm_out, nprocs, nrow_block_out, ncol_block_out, & +!$OMP num_entries_rec, buffer_rec, nprow, npcol, dummy) & +!$OMP PRIVATE(iproc, i_entry_rec, ii, jj) + DO iproc = 0, nprocs - 1 + DO i_entry_rec = 1, num_entries_rec(iproc) + ii = cp_fm_indxg2l(buffer_rec(iproc)%indx(i_entry_rec, 1), nrow_block_out, & + dummy, dummy, nprow) + jj = cp_fm_indxg2l(buffer_rec(iproc)%indx(i_entry_rec, 2), ncol_block_out, & + dummy, dummy, npcol) + + fm_out%local_data(ii, jj) = fm_out%local_data(ii, jj) + buffer_rec(iproc)%msg(i_entry_rec) + END DO + END DO +!$OMP END PARALLEL DO + + CALL timestop(handle2) + + CALL timeset(routineN//"_7_cleanup", handle2) + IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN + WRITE (unit_nr, '(T2,A10,T13,A41)') 'BSE|DEBUG|', 'Starting cleanup of communication buffers' + END IF + + !Clean up all the arrays from the communication process + DO iproc = 0, para_env_out%num_pe - 1 + DEALLOCATE (buffer_rec(iproc)%msg) + DEALLOCATE (buffer_rec(iproc)%indx) + DEALLOCATE (buffer_send(iproc)%msg) + DEALLOCATE (buffer_send(iproc)%indx) + END DO + DEALLOCATE (buffer_rec, buffer_send) + DEALLOCATE (req_array) + DEALLOCATE (entry_counter) + DEALLOCATE (num_entries_rec, num_entries_send) + + CALL timestop(handle2) + CALL timestop(handle) + + END SUBROUTINE truncate_fm + +! ************************************************************************************************** +!> \brief Debug function to write elements of a full matrix to file, if they are larger than a given threshold +!> \param fm ... +!> \param thresh ... +!> \param header ... +!> \param unit_nr ... +!> \param abs_vals ... +! ************************************************************************************************** + SUBROUTINE fm_write_thresh(fm, thresh, header, unit_nr, abs_vals) + + TYPE(cp_fm_type), INTENT(IN) :: fm + REAL(KIND=dp), INTENT(IN) :: thresh + CHARACTER(LEN=*), INTENT(IN) :: header + INTEGER, INTENT(IN) :: unit_nr + LOGICAL, OPTIONAL :: abs_vals + + CHARACTER(LEN=*), PARAMETER :: my_footer = " | ENDING WRITING OF MATRIX", & + routineN = 'fm_write_thresh' + + INTEGER :: handle, i, j, ncol_local, nrow_local + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + LOGICAL :: my_abs_vals + + CALL timeset(routineN, handle) + + IF (PRESENT(abs_vals)) THEN + my_abs_vals = abs_vals + ELSE + my_abs_vals = .FALSE. + END IF + + CALL cp_fm_get_info(matrix=fm, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) header + END IF + IF (my_abs_vals) THEN + DO i = 1, nrow_local + DO j = 1, ncol_local + IF (ABS(fm%local_data(i, j)) > thresh) THEN + WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), & + ABS(fm%local_data(i, j)) + END IF + END DO + END DO + ELSE + DO i = 1, nrow_local + DO j = 1, ncol_local + IF (ABS(fm%local_data(i, j)) > thresh) THEN + WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), & + fm%local_data(i, j) + END IF + END DO + END DO + END IF + CALL fm%matrix_struct%para_env%sync() + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) my_footer + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param bse_tda ... +!> \param bse_abba ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE print_BSE_start_flag(bse_tda, bse_abba, unit_nr) + + LOGICAL, INTENT(IN) :: bse_tda, bse_abba + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'print_BSE_start_flag' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) ' ' + WRITE (unit_nr, '(T2,A79)') '*******************************************************************************' + WRITE (unit_nr, '(T2,A79)') '** **' + WRITE (unit_nr, '(T2,A79)') '** Bethe Salpeter equation (BSE) for excitation energies **' + IF (bse_tda .AND. bse_abba) THEN + WRITE (unit_nr, '(T2,A79)') '** solved with and without Tamm-Dancoff approximation (TDA) **' + ELSE IF (bse_tda) THEN + WRITE (unit_nr, '(T2,A79)') '** solved with Tamm-Dancoff approximation (TDA) **' + ELSE + WRITE (unit_nr, '(T2,A79)') '** solved without Tamm-Dancoff approximation (TDA) **' + END IF + + WRITE (unit_nr, '(T2,A79)') '** **' + WRITE (unit_nr, '(T2,A79)') '*******************************************************************************' + WRITE (unit_nr, *) ' ' + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param fm_mat_S_bar_ia_bse ... +!> \param fm_mat_S_bar_ij_bse ... +!> \param fm_mat_S_trunc ... +!> \param fm_mat_S_ij_trunc ... +!> \param fm_mat_S_ab_trunc ... +!> \param fm_mat_Q_static_bse ... +!> \param fm_mat_Q_static_bse_gemm ... +! ************************************************************************************************** + SUBROUTINE deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, & + fm_mat_S_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, & + fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm) + + TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_trunc, & + fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm + + CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_matrices_bse' + + INTEGER :: handle + + CALL timeset(routineN, handle) + + CALL cp_fm_release(fm_mat_S_bar_ia_bse) + CALL cp_fm_release(fm_mat_S_bar_ij_bse) + CALL cp_fm_release(fm_mat_S_trunc) + CALL cp_fm_release(fm_mat_S_ij_trunc) + CALL cp_fm_release(fm_mat_S_ab_trunc) + CALL cp_fm_release(fm_mat_Q_static_bse) + CALL cp_fm_release(fm_mat_Q_static_bse_gemm) + + CALL timestop(handle) + + END SUBROUTINE deallocate_matrices_bse + +! ************************************************************************************************** +!> \brief Routine for computing the coefficients of the eigenvectors of the BSE matrix from a +!> multiplication with the eigenvalues +!> \param fm_work ... +!> \param eig_vals ... +!> \param beta ... +!> \param gamma ... +!> \param do_transpose ... +! ************************************************************************************************** + SUBROUTINE comp_eigvec_coeff_BSE(fm_work, eig_vals, beta, gamma, do_transpose) + + TYPE(cp_fm_type), INTENT(INOUT) :: fm_work + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & + INTENT(IN) :: eig_vals + REAL(KIND=dp), INTENT(IN) :: beta + REAL(KIND=dp), INTENT(IN), OPTIONAL :: gamma + LOGICAL, INTENT(IN), OPTIONAL :: do_transpose + + CHARACTER(LEN=*), PARAMETER :: routineN = 'comp_eigvec_coeff_BSE' + + INTEGER :: handle, i_row_global, ii, j_col_global, & + jj, ncol_local, nrow_local + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + LOGICAL :: my_do_transpose + REAL(KIND=dp) :: coeff, my_gamma + + CALL timeset(routineN, handle) + + IF (PRESENT(gamma)) THEN + my_gamma = gamma + ELSE + my_gamma = 2.0_dp + END IF + + IF (PRESENT(do_transpose)) THEN + my_do_transpose = do_transpose + ELSE + my_do_transpose = .FALSE. + END IF + + CALL cp_fm_get_info(matrix=fm_work, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + IF (my_do_transpose) THEN + DO jj = 1, ncol_local + j_col_global = col_indices(jj) + DO ii = 1, nrow_local + coeff = (eig_vals(j_col_global)**beta)/my_gamma + fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff + END DO + END DO + ELSE + DO jj = 1, ncol_local + DO ii = 1, nrow_local + i_row_global = row_indices(ii) + coeff = (eig_vals(i_row_global)**beta)/my_gamma + fm_work%local_data(ii, jj) = fm_work%local_data(ii, jj)*coeff + END DO + END DO + END IF + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param idx_prim ... +!> \param idx_sec ... +!> \param eigvec_entries ... +! ************************************************************************************************** + SUBROUTINE sort_excitations(idx_prim, idx_sec, eigvec_entries) + + INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim, idx_sec + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries + + CHARACTER(LEN=*), PARAMETER :: routineN = 'sort_excitations' + + INTEGER :: handle, ii, kk, num_entries, num_mults + INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_prim_work, idx_sec_work, tmp_index + LOGICAL :: unique_entries + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries_work + + CALL timeset(routineN, handle) + + num_entries = SIZE(idx_prim) + + ALLOCATE (tmp_index(num_entries)) + + CALL sort(idx_prim, num_entries, tmp_index) + + ALLOCATE (idx_sec_work(num_entries)) + ALLOCATE (eigvec_entries_work(num_entries)) + + DO ii = 1, num_entries + idx_sec_work(ii) = idx_sec(tmp_index(ii)) + eigvec_entries_work(ii) = eigvec_entries(tmp_index(ii)) + END DO + + DEALLOCATE (tmp_index) + DEALLOCATE (idx_sec) + DEALLOCATE (eigvec_entries) + + CALL MOVE_ALLOC(idx_sec_work, idx_sec) + CALL MOVE_ALLOC(eigvec_entries_work, eigvec_entries) + + !Now check for multiple entries in first idx to check necessity of sorting in second idx + CALL sort_unique(idx_prim, unique_entries) + IF (.NOT. unique_entries) THEN + ALLOCATE (idx_prim_work(num_entries)) + idx_prim_work(:) = idx_prim(:) + ! Find duplicate entries in idx_prim + DO ii = 1, num_entries + IF (idx_prim_work(ii) == 0) CYCLE + num_mults = COUNT(idx_prim_work == idx_prim_work(ii)) + IF (num_mults > 1) THEN + !Set all duplicate entries to 0 + idx_prim_work(ii:ii + num_mults - 1) = 0 + !Start sorting in secondary index + ALLOCATE (idx_sec_work(num_mults)) + ALLOCATE (eigvec_entries_work(num_mults)) + idx_sec_work(:) = idx_sec(ii:ii + num_mults - 1) + eigvec_entries_work(:) = eigvec_entries(ii:ii + num_mults - 1) + ALLOCATE (tmp_index(num_mults)) + CALL sort(idx_sec_work, num_mults, tmp_index) + + !Now write newly sorted indices to original arrays + DO kk = ii, ii + num_mults - 1 + idx_sec(kk) = idx_sec_work(kk - ii + 1) + eigvec_entries(kk) = eigvec_entries_work(tmp_index(kk - ii + 1)) + END DO + !Deallocate work arrays + DEALLOCATE (tmp_index) + DEALLOCATE (idx_sec_work) + DEALLOCATE (eigvec_entries_work) + END IF + idx_prim_work(ii) = idx_prim(ii) + END DO + DEALLOCATE (idx_prim_work) + END IF + + CALL timestop(handle) + + END SUBROUTINE sort_excitations + +! ************************************************************************************************** +!> \brief Roughly estimates the needed runtime and memory during the BSE run +!> \param homo_red ... +!> \param virtual_red ... +!> \param unit_nr ... +!> \param bse_abba ... +!> \param para_env ... +!> \param diag_runtime_est ... +! ************************************************************************************************** + SUBROUTINE estimate_BSE_resources(homo_red, virtual_red, unit_nr, bse_abba, & + para_env, diag_runtime_est) + + INTEGER :: homo_red, virtual_red, unit_nr + LOGICAL :: bse_abba + TYPE(mp_para_env_type), POINTER :: para_env + REAL(KIND=dp) :: diag_runtime_est + + CHARACTER(LEN=*), PARAMETER :: routineN = 'estimate_BSE_resources' + + INTEGER :: handle, num_BSE_matrices + REAL(KIND=dp) :: mem_est + + CALL timeset(routineN, handle) + + ! Number of matrices with size of A in TDA is 2 (A itself and W_ijab) + num_BSE_matrices = 2 + ! With the full diagonalization of ABBA, we need several auxiliary matrices in the process + ! The maximum number is 2 + 4 (additional B and C matrix as well as 4 matrices to create C) + IF (bse_abba) THEN + num_BSE_matrices = 2 + 4 + END IF + mem_est = (8*homo_red**2*virtual_red**2/para_env%num_pe*num_BSE_matrices)/(1024**3) + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A47,T73,F8.3)') 'BSE|', 'Peak memory estimate per MPI rank from BSE [GB]', & + mem_est + WRITE (unit_nr, '(T2,A4)') 'BSE|' + END IF + ! Rough estimation of diagonalization runtimes. Baseline was a full BSE Naphthalene + ! run with 11000x11000 entries in A/B/C, which took 10s on 32 ranks + diag_runtime_est = (homo_red*virtual_red/11000)**3*10*32/para_env%num_pe + + CALL timestop(handle) + + END SUBROUTINE estimate_BSE_resources + +! ************************************************************************************************** +!> \brief Filters eigenvector entries above a given threshold to describe excitations in the +!> singleparticle basis +!> \param fm_eigvec ... +!> \param idx_homo ... +!> \param idx_virt ... +!> \param eigvec_entries ... +!> \param i_exc ... +!> \param virtual ... +!> \param num_entries ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, & + i_exc, virtual, num_entries, mp2_env) + + TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec + INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries + INTEGER :: i_exc, virtual, num_entries + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'filter_eigvec_contrib' + + INTEGER :: eigvec_idx, handle, ii, iproc, jj, kk, & + ncol_local, nrow_local, & + num_entries_local + INTEGER, ALLOCATABLE, DIMENSION(:) :: num_entries_to_comm + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + REAL(KIND=dp) :: eigvec_entry + TYPE(integ_mat_buffer_type), ALLOCATABLE, & + DIMENSION(:) :: buffer_entries + TYPE(mp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + para_env => fm_eigvec%matrix_struct%para_env + + CALL cp_fm_get_info(matrix=fm_eigvec, & + nrow_local=nrow_local, & + ncol_local=ncol_local, & + row_indices=row_indices, & + col_indices=col_indices) + + ALLOCATE (num_entries_to_comm(0:para_env%num_pe - 1)) + num_entries_to_comm(:) = 0 + + DO jj = 1, ncol_local + !First check if i is localized on this proc + IF (col_indices(jj) /= i_exc) THEN + CYCLE + END IF + DO ii = 1, nrow_local + eigvec_idx = row_indices(ii) + eigvec_entry = fm_eigvec%local_data(ii, jj)/SQRT(2.0_dp) + IF (ABS(eigvec_entry) > mp2_env%ri_g0w0%eps_x) THEN + num_entries_to_comm(para_env%mepos) = num_entries_to_comm(para_env%mepos) + 1 + END IF + END DO + END DO + + !Gather number of entries of other processes + CALL para_env%sum(num_entries_to_comm) + + num_entries_local = num_entries_to_comm(para_env%mepos) + + ALLOCATE (buffer_entries(0:para_env%num_pe - 1)) + + DO iproc = 0, para_env%num_pe - 1 + ALLOCATE (buffer_entries(iproc)%msg(num_entries_to_comm(iproc))) + ALLOCATE (buffer_entries(iproc)%indx(num_entries_to_comm(iproc), 2)) + buffer_entries(iproc)%msg = 0.0_dp + buffer_entries(iproc)%indx = 0 + END DO + + kk = 1 + DO jj = 1, ncol_local + !First check if i is localized on this proc + IF (col_indices(jj) /= i_exc) THEN + CYCLE + END IF + DO ii = 1, nrow_local + eigvec_idx = row_indices(ii) + eigvec_entry = fm_eigvec%local_data(ii, jj)/SQRT(2.0_dp) + IF (ABS(eigvec_entry) > mp2_env%ri_g0w0%eps_x) THEN + buffer_entries(para_env%mepos)%indx(kk, 1) = (eigvec_idx - 1)/virtual + 1 + buffer_entries(para_env%mepos)%indx(kk, 2) = MOD(eigvec_idx - 1, virtual) + 1 + buffer_entries(para_env%mepos)%msg(kk) = eigvec_entry + kk = kk + 1 + END IF + END DO + END DO + + DO iproc = 0, para_env%num_pe - 1 + CALL para_env%sum(buffer_entries(iproc)%msg) + CALL para_env%sum(buffer_entries(iproc)%indx) + END DO + + !Now sum up gathered information + num_entries = SUM(num_entries_to_comm) + ALLOCATE (idx_homo(num_entries)) + ALLOCATE (idx_virt(num_entries)) + ALLOCATE (eigvec_entries(num_entries)) + + kk = 1 + DO iproc = 0, para_env%num_pe - 1 + IF (num_entries_to_comm(iproc) /= 0) THEN + DO ii = 1, num_entries_to_comm(iproc) + idx_homo(kk) = buffer_entries(iproc)%indx(ii, 1) + idx_virt(kk) = buffer_entries(iproc)%indx(ii, 2) + eigvec_entries(kk) = buffer_entries(iproc)%msg(ii) + kk = kk + 1 + END DO + END IF + END DO + + !Deallocate all the used arrays + DO iproc = 0, para_env%num_pe - 1 + DEALLOCATE (buffer_entries(iproc)%msg) + DEALLOCATE (buffer_entries(iproc)%indx) + END DO + DEALLOCATE (buffer_entries) + DEALLOCATE (num_entries_to_comm) + NULLIFY (row_indices) + NULLIFY (col_indices) + + !Now sort the results according to the involved singleparticle orbitals + ! (homo first, then virtual) + CALL sort_excitations(idx_homo, idx_virt, eigvec_entries) + + CALL timestop(handle) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief Determines indices within the given energy cutoffs and truncates Eigenvalues and matrices +!> \param fm_mat_S_ia_bse ... +!> \param fm_mat_S_ij_bse ... +!> \param fm_mat_S_ab_bse ... +!> \param fm_mat_S_trunc ... +!> \param fm_mat_S_ij_trunc ... +!> \param fm_mat_S_ab_trunc ... +!> \param Eigenval ... +!> \param Eigenval_reduced ... +!> \param homo ... +!> \param virtual ... +!> \param dimen_RI ... +!> \param unit_nr ... +!> \param homo_red ... +!> \param virt_red ... +!> \param mp2_env ... +! ************************************************************************************************** + SUBROUTINE truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, & + fm_mat_S_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, & + Eigenval, Eigenval_reduced, & + homo, virtual, dimen_RI, unit_nr, & + homo_red, virt_red, & + mp2_env) + + TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, & + fm_mat_S_ab_bse + TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_S_trunc, fm_mat_S_ij_trunc, & + fm_mat_S_ab_trunc + REAL(KIND=dp), DIMENSION(:) :: Eigenval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced + INTEGER, INTENT(IN) :: homo, virtual, dimen_RI, unit_nr + INTEGER, INTENT(OUT) :: homo_red, virt_red + TYPE(mp2_type), INTENT(INOUT) :: mp2_env + + CHARACTER(LEN=*), PARAMETER :: routineN = 'truncate_BSE_matrices' + + INTEGER :: handle, homo_incl, i_homo, j_virt, & + virt_incl + TYPE(cp_blacs_env_type), POINTER :: context + TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia, fm_struct_ij + TYPE(mp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + ! Determine index in homo and virtual for truncation + ! Uses indices of outermost orbitals within energy range (-mp2_env%ri_g0w0%bse_cutoff_occ,mp2_env%ri_g0w0%bse_cutoff_virt) + IF (mp2_env%ri_g0w0%bse_cutoff_occ > 0 .OR. mp2_env%ri_g0w0%bse_cutoff_virt > 0) THEN + IF (-mp2_env%ri_g0w0%bse_cutoff_occ .LT. Eigenval(1) - Eigenval(homo) & + .OR. mp2_env%ri_g0w0%bse_cutoff_occ < 0) THEN + homo_red = homo + homo_incl = 1 + ELSE + homo_incl = 1 + DO i_homo = 1, homo + IF (Eigenval(i_homo) - Eigenval(homo) .GT. -mp2_env%ri_g0w0%bse_cutoff_occ) THEN + homo_incl = i_homo + EXIT + END IF + END DO + homo_red = homo - homo_incl + 1 + END IF + + IF (mp2_env%ri_g0w0%bse_cutoff_virt .GT. Eigenval(homo + virtual) - Eigenval(homo + 1) & + .OR. mp2_env%ri_g0w0%bse_cutoff_virt < 0) THEN + virt_red = virtual + virt_incl = virtual + ELSE + virt_incl = homo + 1 + DO j_virt = 1, virtual + IF (Eigenval(homo + j_virt) - Eigenval(homo + 1) .GT. mp2_env%ri_g0w0%bse_cutoff_virt) THEN + virt_incl = j_virt - 1 + EXIT + END IF + END DO + virt_red = virt_incl + END IF + ELSE + homo_red = homo + virt_red = virtual + homo_incl = 1 + virt_incl = virtual + END IF + IF (unit_nr > 0) THEN + IF (mp2_env%ri_g0w0%bse_cutoff_occ > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A29,T71,F10.3)') 'BSE|', 'Cutoff occupied orbitals [eV]', & + mp2_env%ri_g0w0%bse_cutoff_occ*evolt + ELSE + WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No cutoff given for occupied orbitals' + END IF + IF (mp2_env%ri_g0w0%bse_cutoff_virt > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A28,T71,F10.3)') 'BSE|', 'Cutoff virtual orbitals [eV]', & + mp2_env%ri_g0w0%bse_cutoff_virt*evolt + ELSE + WRITE (unit_nr, '(T2,A4,T7,A36)') 'BSE|', 'No cutoff given for virtual orbitals' + END IF + WRITE (unit_nr, '(T2,A4,T7,A20,T71,I10)') 'BSE|', 'First occupied index', homo_incl + WRITE (unit_nr, '(T2,A4,T7,A18,T71,I10)') 'BSE|', 'Last virtual index', virt_incl + WRITE (unit_nr, '(T2,A4,T7,A35,T71,F10.3)') 'BSE|', 'Energy of first occupied index [eV]', Eigenval(homo_incl)*evolt + WRITE (unit_nr, '(T2,A4,T7,A33,T71,F10.3)') 'BSE|', 'Energy of last virtual index [eV]', Eigenval(homo + virt_incl)*evolt + WRITE (unit_nr, '(T2,A4,T7,A54,T71,F10.3)') 'BSE|', 'Energy difference of first occupied index to homo [eV]', & + -(Eigenval(homo_incl) - Eigenval(homo))*evolt + WRITE (unit_nr, '(T2,A4,T7,A52,T71,F10.3)') 'BSE|', 'Energy difference of last virtual index to LUMO [eV]', & + (Eigenval(homo + virt_incl) - Eigenval(homo + 1))*evolt + WRITE (unit_nr, '(T2,A4,T7,A35,T71,I10)') 'BSE|', 'Number of GW-corrected occupied MOs', mp2_env%ri_g0w0%corr_mos_occ + WRITE (unit_nr, '(T2,A4,T7,A34,T71,I10)') 'BSE|', 'Number of GW-corrected virtual MOs', mp2_env%ri_g0w0%corr_mos_virt + WRITE (unit_nr, '(T2,A4)') 'BSE|' + END IF + IF (unit_nr > 0) THEN + IF (homo - homo_incl + 1 > mp2_env%ri_g0w0%corr_mos_occ) THEN + CPABORT("Number of GW-corrected occupied MOs too small for chosen BSE cutoff") + END IF + IF (virt_incl > mp2_env%ri_g0w0%corr_mos_virt) THEN + CPABORT("Number of GW-corrected virtual MOs too small for chosen BSE cutoff") + END IF + END IF + !Truncate full fm_S matrices + !Allocate new truncated matrices of proper size + para_env => fm_mat_S_ia_bse%matrix_struct%para_env + context => fm_mat_S_ia_bse%matrix_struct%context + + CALL cp_fm_struct_create(fm_struct_ia, para_env, context, dimen_RI, homo_red*virt_red) + CALL cp_fm_struct_create(fm_struct_ij, para_env, context, dimen_RI, homo_red*homo_red) + CALL cp_fm_struct_create(fm_struct_ab, para_env, context, dimen_RI, virt_red*virt_red) + + CALL cp_fm_create(fm_mat_S_trunc, fm_struct_ia, "fm_S_trunc") + CALL cp_fm_create(fm_mat_S_ij_trunc, fm_struct_ij, "fm_S_ij_trunc") + CALL cp_fm_create(fm_mat_S_ab_trunc, fm_struct_ab, "fm_S_ab_trunc") + + !Copy parts of original matrices to truncated ones + IF (mp2_env%ri_g0w0%bse_cutoff_occ > 0 .OR. mp2_env%ri_g0w0%bse_cutoff_virt > 0) THEN + !Truncate eigenvals + ALLOCATE (Eigenval_reduced(homo_red + virt_red)) + Eigenval_reduced(:) = Eigenval(homo_incl:homo + virt_incl) + + CALL truncate_fm(fm_mat_S_trunc, fm_mat_S_ia_bse, virtual, & + homo_red, virt_red, unit_nr, mp2_env, & + nrow_offset=homo_incl) + CALL truncate_fm(fm_mat_S_ij_trunc, fm_mat_S_ij_bse, homo, & + homo_red, homo_red, unit_nr, mp2_env, & + homo_incl, homo_incl) + CALL truncate_fm(fm_mat_S_ab_trunc, fm_mat_S_ab_bse, virtual, & + virt_red, virt_red, unit_nr, mp2_env) + + ELSE + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T2,A4,T7,A37)') 'BSE|', 'No truncation of BSE matrices applied' + WRITE (unit_nr, '(T2,A4)') 'BSE|' + END IF + ALLOCATE (Eigenval_reduced(homo_red + virt_red)) + Eigenval_reduced(:) = Eigenval(:) + CALL cp_fm_to_fm_submat_general(fm_mat_S_ia_bse, fm_mat_S_trunc, dimen_RI, homo_red*virt_red, & + 1, 1, 1, 1, context) + CALL cp_fm_to_fm_submat_general(fm_mat_S_ij_bse, fm_mat_S_ij_trunc, dimen_RI, homo_red*homo_red, & + 1, 1, 1, 1, context) + CALL cp_fm_to_fm_submat_general(fm_mat_S_ab_bse, fm_mat_S_ab_trunc, dimen_RI, virt_red*virt_red, & + 1, 1, 1, 1, context) + END IF + + CALL cp_fm_struct_release(fm_struct_ia) + CALL cp_fm_struct_release(fm_struct_ij) + CALL cp_fm_struct_release(fm_struct_ab) + + NULLIFY (para_env) + NULLIFY (context) + + CALL timestop(handle) + + END SUBROUTINE truncate_BSE_matrices + +END MODULE diff --git a/src/input_constants.F b/src/input_constants.F index 6f1cd2bc5b..6d74a430e6 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1278,5 +1278,22 @@ MODULE input_constants INTEGER, PARAMETER, PUBLIC :: e_dens_total_hard_approx = 1, & e_dens_total_density = 2, & e_dens_soft_density = 3 + ! BSE spin configuration + INTEGER, PARAMETER, PUBLIC :: bse_singlet = 0, & + bse_triplet = 1 + + ! BSE method of diagonalization (full diagonalization or iteratively) + INTEGER, PARAMETER, PUBLIC :: bse_fulldiag = 0, & + bse_iterdiag = 1 + + ! BSE level of approximation to the ABBA matrix (i.e. with/without TDA) + INTEGER, PARAMETER, PUBLIC :: bse_tda = 0, & + bse_abba = 1, & + bse_both = 2 + + ! BSE iterative abortion condition + INTEGER, PARAMETER, PUBLIC :: bse_iter_en_cond = 0, & + bse_iter_res_cond = 1, & + bse_iter_both_cond = 2 END MODULE input_constants diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index 55150bb2ec..314bcb889e 100644 --- a/src/input_cp2k_mp2.F +++ b/src/input_cp2k_mp2.F @@ -12,45 +12,47 @@ !> \author MDB ! ************************************************************************************************** MODULE input_cp2k_mp2 - USE bibliography, ONLY: & - Bates2013, DelBen2012, DelBen2013, DelBen2015, DelBen2015b, Rybkin2016, Wilhelm2016a, & - Wilhelm2016b, Wilhelm2017, Wilhelm2018 - USE cp_eri_mme_interface, ONLY: create_eri_mme_section - USE cp_output_handling, ONLY: add_last_numeric,& - cp_print_key_section_create,& - debug_print_level,& - high_print_level,& - low_print_level,& - medium_print_level,& - silent_print_level - USE cp_units, ONLY: cp_unit_to_cp2k - USE input_constants, ONLY: & - do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, & - do_potential_long, do_potential_mix_cl, do_potential_short, do_potential_truncated, & - do_potential_tshpsc, eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, & - gw_read_exx, gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, & - kp_weights_W_tailored, kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, & - mp2_method_none, numerical, ot_precond_full_all, ot_precond_full_kinetic, & - ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, & - ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, & - ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, & - wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, & - z_solver_sd - USE input_cp2k_hfx, ONLY: create_hfx_section - USE input_cp2k_kpoints, ONLY: create_kpoint_set_section - USE input_keyword_types, ONLY: keyword_create,& - keyword_release,& - keyword_type - USE input_section_types, ONLY: section_add_keyword,& - section_add_subsection,& - section_create,& - section_release,& - section_type - USE input_val_types, ONLY: integer_t,& - real_t - USE kinds, ONLY: dp - USE string_utilities, ONLY: newline,& - s2a + USE bibliography, ONLY: & + Bates2013, DelBen2012, DelBen2013, DelBen2015, DelBen2015b, Rybkin2016, Wilhelm2016a, & + Wilhelm2016b, Wilhelm2017, Wilhelm2018 + USE cp_eri_mme_interface, ONLY: create_eri_mme_section + USE cp_output_handling, ONLY: add_last_numeric, & + cp_print_key_section_create, & + debug_print_level, & + high_print_level, & + low_print_level, & + medium_print_level, & + silent_print_level + USE cp_units, ONLY: cp_unit_to_cp2k + USE input_constants, ONLY: & + bse_fulldiag, bse_iterdiag, bse_tda, bse_abba, bse_both, & + bse_iter_both_cond, bse_iter_en_cond, bse_iter_res_cond, bse_singlet, & + bse_triplet, do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, & + do_potential_long, do_potential_mix_cl, do_potential_short, do_potential_truncated, & + do_potential_tshpsc, eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, & + gw_read_exx, gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, & + kp_weights_W_tailored, kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, & + mp2_method_none, numerical, ot_precond_full_all, ot_precond_full_kinetic, & + ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, & + ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, & + ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, & + wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, & + z_solver_sd + USE input_cp2k_hfx, ONLY: create_hfx_section + USE input_cp2k_kpoints, ONLY: create_kpoint_set_section + USE input_keyword_types, ONLY: keyword_create, & + keyword_release, & + keyword_type + USE input_section_types, ONLY: section_add_keyword, & + section_add_subsection, & + section_create, & + section_release, & + section_type + USE input_val_types, ONLY: integer_t, & + real_t + USE kinds, ONLY: dp + USE string_utilities, ONLY: newline, & + s2a #include "./base/base_uses.f90" IMPLICIT NONE @@ -852,16 +854,6 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="BSE", & - description="If true, electronic excitation energies are computed from the "// & - "Bethe-Salpeter equation on top of GW eigenvalues. Parameter of BSE can be adjusted in "// & - "the corresponding section.", & - usage="BSE", & - default_l_val=.FALSE., & - lone_keyword_l_val=.TRUE.) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="IMAGE_CHARGE_MODEL", & variants=(/"IC"/), & description="If true, an image charge model is applied to mimic the renormalization of "// & @@ -1215,39 +1207,205 @@ CONTAINS TYPE(section_type), POINTER :: section TYPE(keyword_type), POINTER :: keyword + TYPE(section_type), POINTER :: subsection CPASSERT(.NOT. ASSOCIATED(section)) CALL section_create(section, __LOCATION__, name="BSE", & - description="Parameters influencing the Bethe-Salpeter calculations on "// & - "top of GW eigenvalues.", & - n_keywords=12, n_subsections=1, repeats=.FALSE.) + description="Parameters for a calculation solving the Bethe-Salpeter equation "// & + "(BSE) for electronic excitations. The full BSE "// & + "$\left( \begin{array}{cc}A & B\\B & A\end{array} \right)\left( \begin{array}{cc}X^{(n)}\\Y^{(n)}\end{array} \right) = "// & + "\Omega^{(n)}\left(\begin{array}{cc}1&0\\0&-1\end{array}\right)\left(\begin{array}{cc}X^{(n)}\\Y^{(n)}\end{array}\right)$ "// & + "enables of electronic excitation energies $\Omega^{(n)}$. The BSE can be solved by diagonalizing "// & + "the full ABBA-matrix or by setting B=0, i.e. within the Tamm-Dancoff approximation (TDA). "// & + "Preliminary reference: Eq. (35) in PRB 92, 045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209", & + n_keywords=8, n_subsections=2, repeats=.FALSE.) NULLIFY (keyword) - CALL keyword_create(keyword, __LOCATION__, name="NUM_Z_VECTORS", & - description="Number Z vectors used in the subspace iterations. This is a convergence "// & - "parameter: Increasing NUM_Z_VECTORS should keep the result constant in case of convergence,.", & - usage="NUM_Z_VECTORS 50", & - default_i_val=20) + CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Activates BSE calculations.", & + usage="&BSE .TRUE.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="THRESHOLD_MIN_TRANS", & - variants=(/"EPS"/), & + CALL keyword_create(keyword, __LOCATION__, name="SPIN_CONFIG", & + description="Choose between calculation of singlet or triplet excitation (cf. given Reference above).", & + usage="SPIN_CONFIG TRIPLET", & + enum_c_vals=s2a("SINGLET", "TRIPLET"), & + enum_i_vals=(/bse_singlet, bse_triplet/), & + enum_desc=s2a("Computes singlet excitations.", & + "Computes triplet excitations."), & + default_i_val=bse_singlet) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="BSE_DIAG_METHOD", & + description="Method for BSE calculations. "// & + "Choose between full or iterative diagonalization.", & + usage="&BSE_DIAG_METHOD FULLDIAG", & + enum_c_vals=s2a("FULLDIAG", "ITERDIAG"), & + enum_i_vals=(/bse_fulldiag, bse_iterdiag/), & + enum_desc=s2a("Fully diagonalizes the BSE matrices within the chosen level of approximation.", & + "Iterative diagonalization has not been implemented yet."), & + default_i_val=bse_fulldiag) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="BSE_APPROX", & + description="Level of approximation applied to BSE calculations. "// & + "Choose between Tamm Dancoff approximation (TDA) and diagonalization of the full ABBA-matrix.", & + usage="&BSE_APPROX TDA", & + enum_c_vals=s2a("TDA", "ABBA", "BOTH"), & + enum_i_vals=(/bse_tda, bse_abba, bse_both/), & + enum_desc=s2a("The TDA is applied, i.e. B=0.", & + "The ABBA-matrix is diagonalized, i.e. the TDA is not applied.", & + "The BSE is solved within the TDA (B=0) as well as for the full ABBA-matrix."), & + default_i_val=bse_tda) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ENERGY_CUTOFF_OCC", & + description="Removing all orbitals with indices i,j from A_ia,jb and B_ia,jb with energy difference "// & + "to HOMO level larger than the given energy cutoff, i.e. "// & + "$\epsilon_{\mathrm{HOMO}}-\epsilon_a>\mathtt{ENERGY_CUTOFF_OCC}$. "// & + "Can be used to accelerate runtime and reduce memory consumption.", & + usage="ENERGY_CUTOFF_OCC 10.0", unit_str="eV", & + type_of_var=real_t, default_r_val=-1.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ENERGY_CUTOFF_VIRT", & + description="Removing all orbitals with indices a,b from A_ia,jb and B_ia,jb with energy difference "// & + "to LUMO level larger than the given energy cutoff, i.e. "// & + "$\epsilon_i-\epsilon_{\mathrm{LUMO}}>\mathtt{ENERGY_CUTOFF_VIRT}$. "// & + "Can be used to accelerate runtime and reduce memory consumption.", & + usage="ENERGY_CUTOFF_VIRT 10.0", unit_str="eV", & + type_of_var=real_t, default_r_val=-1.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="BSE_DEBUG_PRINT", & + description="Activates debug print statements in the BSE calculation.", & + usage="&BSE_DEBUG_PRINT .TRUE.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="NUM_PRINT_EXC", & + description="Number of printed excitation energies. Does not affect computation time.", & + usage="NUM_PRINT_EXC 10", & + default_i_val=10) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EPS_X", & + description="Threshold for printing contributions of singleparticle "// & + "transitions, i.e. elements of the eigenvector $X_{ia}^n$.", & + usage="EPS_X 0.1", & + type_of_var=real_t, default_r_val=0.1_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + NULLIFY (subsection) + CALL create_bse_iterat_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param section ... +! ************************************************************************************************** + SUBROUTINE create_bse_iterat_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + CALL section_create(section, __LOCATION__, name="BSE_ITERAT", & + description="Parameters influencing the iterative Bethe-Salpeter calculation. "// & + "The iterative solver has not been fully implemented yet.", & + n_keywords=9, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + + CALL keyword_create(keyword, __LOCATION__, name="DAVIDSON_ABORT_COND", & + description="Desired abortion condition for Davidson solver", & + usage="DAVIDSON_ABORT_COND OR", & + enum_c_vals=s2a("EN", "RES", "OR"), & + enum_i_vals=(/bse_iter_en_cond, bse_iter_res_cond, bse_iter_both_cond/), & + enum_desc=s2a("Uses energy threshold for successfully exiting solver.", & + "Uses residual threshold for successfully exiting solver.", & + "Uses either energy or residual threshold for successfully exiting solver."), & + default_i_val=bse_iter_en_cond) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="NUM_EXC_EN", & + description="Number of lowest excitation energies to be computed.", & + usage="NUM_EXC_EN 3", & + default_i_val=3) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="NUM_ADD_START_Z_SPACE", & + description="Determines the initial dimension of the subspace as "// & + "dim = (NUM_EXC_EN+NUM_ADD_START_Z_SPACE)", & + usage="NUM_ADD_START_Z_SPACE 1", & + default_i_val=0) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="FAC_MAX_Z_SPACE", & + description="Factor to determine maximum dimension of the Davidson subspace. "// & + "dimension = (NUM_EXC_EN+NUM_ADD_START_Z_SPACE)*FAC_MAX_Z_SPACE", & + usage="FAC_MAX_Z_SPACE 5", & + default_i_val=5) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="NUM_NEW_T", & + description="Number of new t vectors added. "// & + "Must be smaller/equals (NUM_EXC_EN+NUM_ADD_START_Z_SPACE)", & + usage="NUM_NEW_T 4", & + default_i_val=1) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="EPS_RES", & description="Threshold for stopping the iteration for computing the transition energies. "// & - "If the lowest excitation changes by less than THRESHOLD_MIN_TRANS (in eV), the ieration stops.", & - usage="THRESHOLD_MIN_TRANS 0.001", & - default_r_val=0.001_dp) + "If the residuals inside the Davidson space change by less than EPS_RES (in eV), the iteration "// & + "stops.", & + usage="EPS_RES 0.001", unit_str="eV", & + type_of_var=real_t, default_r_val=0.001_dp) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER", & + CALL keyword_create(keyword, __LOCATION__, name="EPS_EXC_EN", & + description="Threshold for stopping the iteration for computing the transition energies. "// & + "If the desired excitation energies change by less than EPS_EXC_EN (in eV), the iteration "// & + "stops.", & + usage="EPS_EXC_EN 0.001", unit_str="eV", & + type_of_var=real_t, default_r_val=0.001_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="NUM_DAVIDSON_ITER", & description="Maximum number of iterations for determining the transition energies.", & - usage="MAX_ITER 200", & - default_i_val=200) + usage="MAX_ITER 100", & + default_i_val=100) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="Z_SPACE_ENERGY_CUTOFF", & + description="Cutoff (in eV) for maximal energy difference entering the A matrix. "// & + "Per default and for negative values, there is no cutoff applied.", & + usage="Z_SPACE_ENERGY_CUTOFF 60", unit_str="eV", & + type_of_var=real_t, default_r_val=-1.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) END SUBROUTINE ! ************************************************************************************************** diff --git a/src/mp2_integrals.F b/src/mp2_integrals.F index 9d4b44a73a..237b71535b 100644 --- a/src/mp2_integrals.F +++ b/src/mp2_integrals.F @@ -1182,13 +1182,9 @@ CONTAINS CALL timeset(routineN//"_E_Ex_"//TRIM(intermed_mat%descr), handle) CALL copy_dbcsr_to_fm(intermed_mat%matrix_ia_jb, intermed_mat%fm_BIb_jb) - IF (.NOT. (TRIM(intermed_mat%descr) .EQ. "bse_ab")) THEN - - CALL grep_my_integrals(para_env, intermed_mat%fm_BIb_jb, BIb_jb, intermed_mat%max_row_col_local, & - intermed_mat%local_col_row_info, & - my_B_end, my_B_start) - - END IF + CALL grep_my_integrals(para_env, intermed_mat%fm_BIb_jb, BIb_jb, intermed_mat%max_row_col_local, & + intermed_mat%local_col_row_info, & + my_B_end, my_B_start) CALL timestop(handle) END SUBROUTINE ao_to_mo_and_store_B diff --git a/src/mp2_setup.F b/src/mp2_setup.F index 79fc7035d8..7d22124569 100644 --- a/src/mp2_setup.F +++ b/src/mp2_setup.F @@ -71,9 +71,9 @@ CONTAINS do_ri_sos_mp2, do_rpa REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals TYPE(cp_logger_type), POINTER :: logger - TYPE(section_vals_type), POINTER :: cphf_section, eri_mme_section, & - gw_section, low_scaling_section, & - mp2_section + TYPE(section_vals_type), POINTER :: bse_section, cphf_section, & + eri_mme_section, gw_section, & + low_scaling_section, mp2_section CALL timeset(routineN, handle) logger => cp_get_default_logger() @@ -168,14 +168,44 @@ CONTAINS CALL section_vals_val_get(mp2_section, "RI_RPA%GW%SOC_ENERGY_WINDOW", & r_val=mp2_env%ri_g0w0%soc_energy_window) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE", & + NULLIFY (bse_section) + bse_section => section_vals_get_subs_vals(mp2_section, "RI_RPA%GW%BSE") + CALL section_vals_val_get(bse_section, "_SECTION_PARAMETERS_", & l_val=mp2_env%ri_g0w0%do_bse) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%NUM_Z_VECTORS", & - i_val=mp2_env%ri_g0w0%num_z_vectors) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%THRESHOLD_MIN_TRANS", & - r_val=mp2_env%ri_g0w0%eps_min_trans) - CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%MAX_ITER", & - i_val=mp2_env%ri_g0w0%max_iter_bse) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%SPIN_CONFIG", & + i_val=mp2_env%ri_g0w0%bse_spin_config) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%ENERGY_CUTOFF_OCC", & + r_val=mp2_env%ri_g0w0%bse_cutoff_occ) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%ENERGY_CUTOFF_VIRT", & + r_val=mp2_env%ri_g0w0%bse_cutoff_virt) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%NUM_PRINT_EXC", & + i_val=mp2_env%ri_g0w0%num_print_exc) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%EPS_X", & + r_val=mp2_env%ri_g0w0%eps_x) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%DAVIDSON_ABORT_COND", & + i_val=mp2_env%ri_g0w0%davidson_abort_cond) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_EXC_EN", & + i_val=mp2_env%ri_g0w0%num_exc_en) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_ADD_START_Z_SPACE", & + i_val=mp2_env%ri_g0w0%num_add_start_z_space) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%FAC_MAX_Z_SPACE", & + i_val=mp2_env%ri_g0w0%fac_max_z_space) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_NEW_T", & + i_val=mp2_env%ri_g0w0%num_new_t) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%EPS_RES", & + r_val=mp2_env%ri_g0w0%eps_res) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%EPS_EXC_EN", & + r_val=mp2_env%ri_g0w0%eps_exc_en) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_DAVIDSON_ITER", & + i_val=mp2_env%ri_g0w0%num_davidson_iter) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%Z_SPACE_ENERGY_CUTOFF", & + r_val=mp2_env%ri_g0w0%z_space_energy_cutoff) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_DIAG_METHOD", & + i_val=mp2_env%ri_g0w0%bse_diag_method) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_APPROX", & + i_val=mp2_env%ri_g0w0%bse_approx) + CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_DEBUG_PRINT", & + l_val=mp2_env%ri_g0w0%bse_debug_print) CALL section_vals_val_get(mp2_section, "RI_RPA%GW%IMAGE_CHARGE_MODEL", & l_val=mp2_env%ri_g0w0%do_ic_model) diff --git a/src/mp2_types.F b/src/mp2_types.F index 32a1fb9e61..da1b794784 100644 --- a/src/mp2_types.F +++ b/src/mp2_types.F @@ -220,11 +220,25 @@ MODULE mp2_types do_aux_bas_gw = .FALSE. REAL(KIND=dp) :: frac_aux_mos = 0.0_dp INTEGER :: num_omega_points = 0 - LOGICAL :: do_bse = .FALSE. - INTEGER :: num_z_vectors = 0, & - max_iter_bse = 0 - REAL(KIND=dp) :: eps_min_trans = 0.0_dp - LOGICAL :: do_ic_model = .FALSE., & + INTEGER :: bse_spin_config = 0, & + bse_diag_method = 0, & + bse_approx = 0, & + num_exc_en = 0, & + num_print_exc = 0, & + num_add_start_z_space = 0, & + fac_max_z_space = 0, & + num_new_t = 0, & + num_davidson_iter = 0, & + davidson_abort_cond = 0 + REAL(KIND=dp) :: eps_res = 0.0_dp, & + eps_exc_en = 0.0_dp, & + eps_x = 0.0_dp, & + bse_cutoff_occ = 0.0_dp, & + bse_cutoff_virt = 0.0_dp, & + z_space_energy_cutoff = 0.0_dp + LOGICAL :: do_bse = .FALSE., & + bse_debug_print = .FALSE., & + do_ic_model = .FALSE., & print_ic_values = .FALSE. REAL(KIND=dp) :: eps_dist = 0.0_dp TYPE(one_dim_real_array), DIMENSION(2) :: ic_corr_list = one_dim_real_array(NULL()) diff --git a/src/rpa_gw.F b/src/rpa_gw.F index cb1c0fe3d2..cd3cab5b00 100644 --- a/src/rpa_gw.F +++ b/src/rpa_gw.F @@ -773,6 +773,7 @@ CONTAINS !> \param jquad ... !> \param nmo ... !> \param num_fit_points ... +!> \param num_integ_points ... !> \param do_bse ... !> \param do_im_time ... !> \param do_periodic ... @@ -800,8 +801,8 @@ CONTAINS ! ************************************************************************************************** SUBROUTINE compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, & gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, & - do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, & - fermi_level_offset, & + num_integ_points, do_bse, do_im_time, do_periodic, & + first_cycle_periodic_correction, fermi_level_offset, & omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, wj, & fm_mat_Q, fm_mat_Q_static_bse, fm_mat_R_gw, fm_mat_S_gw, & fm_mat_S_gw_work, mo_coeff, para_env, & @@ -812,7 +813,8 @@ CONTAINS DIMENSION(:, :, :, :), INTENT(INOUT) :: vec_Sigma_c_gw INTEGER, INTENT(IN) :: dimen_nm_gw, dimen_RI INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo - INTEGER, INTENT(IN) :: jquad, nmo, num_fit_points + INTEGER, INTENT(IN) :: jquad, nmo, num_fit_points, & + num_integ_points LOGICAL, INTENT(IN) :: do_bse, do_im_time, do_periodic LOGICAL, INTENT(INOUT) :: first_cycle_periodic_correction REAL(KIND=dp), INTENT(INOUT) :: fermi_level_offset, omega @@ -858,7 +860,8 @@ CONTAINS ! symmetrize the result, fm_mat_R_gw is only temporary work matrix CALL cp_fm_upper_to_full(fm_mat_Q, fm_mat_R_gw) - IF (do_bse .AND. jquad == 1) THEN + ! Omega=0 is at last index, not at jquad==1 + IF (do_bse .AND. jquad == num_integ_points) THEN CALL cp_fm_to_fm(fm_mat_Q, fm_mat_Q_static_bse) END IF diff --git a/src/rpa_main.F b/src/rpa_main.F index 41bcae16e8..eeea8b6582 100644 --- a/src/rpa_main.F +++ b/src/rpa_main.F @@ -18,8 +18,7 @@ MODULE rpa_main USE bibliography, ONLY: & Bates2013, DelBen2013, DelBen2015, Ren2011, Ren2013, Wilhelm2016a, Wilhelm2016b, & Wilhelm2017, Wilhelm2018, cite_reference - USE bse, ONLY: do_subspace_iterations,& - mult_B_with_W_and_fill_local_3c_arrays + USE bse_main, ONLY: start_bse_calculation USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_release,& cp_blacs_env_type @@ -585,7 +584,6 @@ CONTAINS ! for Bethe-Salpeter, we need other matrix fm_mat_S IF (do_bse) THEN - CALL create_integ_mat(BIb_C_2D_bse_ij, para_env, para_env_sub, color_sub, ngroup, integ_group_size, & dimen_RI_red, [dimen_homo_square], color_rpa_group, & mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, & @@ -860,7 +858,7 @@ CONTAINS ELSE grid_2D = 1 END IF - CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2d) + CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2D) IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & @@ -1102,9 +1100,9 @@ CONTAINS COMPLEX(KIND=dp), ALLOCATABLE, & DIMENSION(:, :, :, :) :: vec_Sigma_c_gw INTEGER :: count_ev_sc_GW, cut_memory, group_size_P, gw_corr_lev_tot, handle, handle3, i, & - ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, max_iter_bse, min_bsize, mm_style, & - nkp, nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, & - num_Z_vectors, Pspin, Qspin, size_P + ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, min_bsize, mm_style, nkp, & + nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, Pspin, Qspin, & + size_P INTEGER(int_8) :: dbcsr_nflop INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c @@ -1115,20 +1113,20 @@ CONTAINS first_cycle_periodic_correction, my_open_shell, print_ic_values LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_P_blocks REAL(KIND=dp) :: a_scaling, alpha, dbcsr_time, e_axk, e_axk_corr, eps_filter, & - eps_filter_im_time, eps_min_trans, ext_scaling, fermi_level_offset, & - fermi_level_offset_input, my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old + eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, & + my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, & trace_Qomega, vec_omega_fit_gw, wj, & wkp_W REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_W_gw, weights_cos_tf_t_to_w, & weights_cos_tf_w_to_t, & weights_sin_tf_t_to_w - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, & - B_bar_ijQ_bse_local, B_iaQ_bse_local, Eigenval_last, Eigenval_scf, vec_Sigma_x_gw + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_last, Eigenval_scf, & + vec_Sigma_x_gw TYPE(cp_cfm_type) :: cfm_mat_Q TYPE(cp_fm_type) :: fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, fm_mat_RI_global_work, & - fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & - fm_scaled_dm_virt_tau + fm_mat_S_ia_bse, fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, & + fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_gw_work, fm_mat_W, & fm_mo_coeff_occ, fm_mo_coeff_virt TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_L_kpoints, fm_mat_Minv_L_kpoints @@ -1311,16 +1309,12 @@ CONTAINS IF (do_bse) THEN - num_Z_vectors = mp2_env%ri_g0w0%num_z_vectors - eps_min_trans = mp2_env%ri_g0w0%eps_min_trans - max_iter_bse = mp2_env%ri_g0w0%max_iter_bse - CALL cp_fm_create(fm_mat_Q_static_bse_gemm, fm_mat_Q_gemm(1)%matrix_struct) CALL cp_fm_to_fm(fm_mat_Q_gemm(1), fm_mat_Q_static_bse_gemm) CALL cp_fm_set_all(fm_mat_Q_static_bse_gemm, 0.0_dp) CALL cp_fm_create(fm_mat_Q_static_bse, fm_mat_Q(1)%matrix_struct) - CALL cp_fm_to_fm(fm_mat_Q_gemm(1), fm_mat_Q_static_bse) + CALL cp_fm_to_fm(fm_mat_Q(1), fm_mat_Q_static_bse) CALL cp_fm_set_all(fm_mat_Q_static_bse, 0.0_dp) END IF @@ -1475,14 +1469,16 @@ CONTAINS CALL calc_mat_Q(fm_mat_S(1), do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual(1), & Eigenval(:, 1, 1), Eigenval_scf(:, 1, 1), & homo(1), omega, omega_old, jquad, mm_style, dimen_RI_red, dimen_ia(1), alpha, fm_mat_Q(1), & - fm_mat_Q_gemm(1), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter) + fm_mat_Q_gemm(1), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, & + num_integ_points) IF (my_open_shell) THEN CALL calc_mat_Q(fm_mat_S(2), do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual(2), & Eigenval(:, 1, 2), Eigenval_scf(:, 1, 2), & homo(2), omega, omega_old, jquad, mm_style, & dimen_RI_red, dimen_ia(2), alpha, fm_mat_Q(2), fm_mat_Q_gemm(2), do_bse, & - fm_mat_Q_static_bse_gemm, dgemm_counter) + fm_mat_Q_static_bse_gemm, dgemm_counter, & + num_integ_points) ! For SOS-MP2 we need both matrices separately IF (.NOT. do_ri_sos_laplace_mp2) THEN @@ -1551,7 +1547,7 @@ CONTAINS END IF ELSE CALL compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI_red, gw_corr_lev_occ, & - gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, & + gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, num_integ_points, & do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, & fermi_level_offset, & omega, Eigenval(:, 1, :), delta_corr, vec_omega_fit_gw, vec_W_gw, wj, & @@ -1712,16 +1708,24 @@ CONTAINS ! postprocessing after GW for Bethe-Salpeter IF (do_bse) THEN - CALL mult_B_with_W_and_fill_local_3c_arrays(fm_mat_S_ij_bse, fm_mat_S_ab_bse, fm_mat_S(1), fm_mat_Q_static_bse, & - fm_mat_Q_static_bse_gemm, & - B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & - B_iaQ_bse_local, dimen_RI_red, homo(1), virtual(1), dimen_ia(1), & - gd_array, color_sub, para_env) - - CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, & - B_iaQ_bse_local, homo(1), virtual(1), num_Z_vectors, & - max_iter_bse, eps_min_trans, Eigenval(:, 1, 1), para_env) - + ! Create a copy of fm_mat_S for usage in BSE + CALL cp_fm_create(fm_mat_S_ia_bse, fm_mat_S(1)%matrix_struct) + CALL cp_fm_to_fm(fm_mat_S(1), fm_mat_S_ia_bse) + ! Remove energy/frequency factor from 3c-Integral for BSE + IF (iter_sc_gw0 == 1) THEN + CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), & + Eigenval_last(:, 1, 1), homo(1), omega) + ELSE + CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), & + Eigenval_scf(:, 1, 1), homo(1), omega) + END IF + ! Main routine for all BSE postprocessing + CALL start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, & + fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, & + Eigenval, homo, virtual, dimen_RI, dimen_RI_red, & + gd_array, color_sub, mp2_env, unit_nr) + ! Release BSE-copy of fm_mat_S + CALL cp_fm_release(fm_mat_S_ia_bse) END IF IF (my_do_gw) THEN diff --git a/src/rpa_util.F b/src/rpa_util.F index 526066a953..551bc273d0 100644 --- a/src/rpa_util.F +++ b/src/rpa_util.F @@ -653,11 +653,12 @@ CONTAINS !> \param do_bse ... !> \param fm_mat_Q_static_bse_gemm ... !> \param dgemm_counter ... +!> \param num_integ_points ... ! ************************************************************************************************** SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual, & Eigenval, Eigenval_scf, & homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, & - do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter) + do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, num_integ_points) TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, first_cycle INTEGER, INTENT(IN) :: iter_sc_GW0, virtual @@ -670,6 +671,7 @@ CONTAINS LOGICAL, INTENT(IN) :: do_bse TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q_static_bse_gemm TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter + INTEGER, INTENT(IN) :: num_integ_points CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_Q' @@ -692,8 +694,10 @@ CONTAINS CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, & fm_mat_Q, dgemm_counter) - - IF (do_bse .AND. jquad == 1) THEN + ! fm_mat_Q_static_bse_gemm does not enter W_ijab (A matrix in TDA), but only full ABBA + ! (since only B_ij_bar enters W_ijab) + ! Changing jquad, since omega=0 is at last idx + IF (do_bse .AND. jquad == num_integ_points) THEN CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm) END IF CALL timestop(handle) diff --git a/tests/QS/regtest-bse/BSE_H2O_PBE.inp b/tests/QS/regtest-bse/BSE_H2O_PBE.inp new file mode 100644 index 0000000000..fef7ace9a4 --- /dev/null +++ b/tests/QS/regtest-bse/BSE_H2O_PBE.inp @@ -0,0 +1,80 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT BSE_H2O_PBE + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END TIMINGS +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-7 + MAX_SCF 100 + SCF_GUESS ATOMIC + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + &RI_RPA + QUADRATURE_POINTS 100 + &GW + CORR_MOS_OCC 1000 + CORR_MOS_VIRT 1000 + RI_SIGMA_X + &BSE + BSE_APPROX ABBA + BSE_DIAG_METHOD FULLDIAG + ENERGY_CUTOFF_OCC 60.0 + ENERGY_CUTOFF_VIRT 60.0 + &END BSE + &END GW + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.000 6.000 6.000 + PERIODIC NONE + &END CELL + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + COORD_FILE_FORMAT xyz + COORD_FILE_NAME H2O_gas.xyz + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-bse/H2O_gas.xyz b/tests/QS/regtest-bse/H2O_gas.xyz new file mode 100644 index 0000000000..d1b6342153 --- /dev/null +++ b/tests/QS/regtest-bse/H2O_gas.xyz @@ -0,0 +1,5 @@ + 3 + + O 0.000000 0.000000 -0.111000 + H 0.000000 -0.744000 0.495000 + H 0.000000 0.744000 0.495000 diff --git a/tests/QS/regtest-bse/TDA_H2O_PBE.inp b/tests/QS/regtest-bse/TDA_H2O_PBE.inp new file mode 100644 index 0000000000..15e5a5c06e --- /dev/null +++ b/tests/QS/regtest-bse/TDA_H2O_PBE.inp @@ -0,0 +1,80 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT BSE_H2O_PBE + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END TIMINGS +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 100 + REL_CUTOFF 20 + &END MGRID + &POISSON + PERIODIC NONE + POISSON_SOLVER WAVELET + &END POISSON + &QS + EPS_DEFAULT 1.0E-15 + EPS_PGF_ORB 1.0E-30 + METHOD GPW + &END QS + &SCF + EPS_SCF 1.0E-7 + MAX_SCF 100 + SCF_GUESS ATOMIC + &PRINT + &RESTART OFF + &END RESTART + &END PRINT + &END SCF + &XC + &WF_CORRELATION + &RI_RPA + QUADRATURE_POINTS 100 + &GW + CORR_MOS_OCC 1000 + CORR_MOS_VIRT 1000 + RI_SIGMA_X + &BSE + BSE_APPROX TDA + BSE_DIAG_METHOD FULLDIAG + ENERGY_CUTOFF_OCC 60.0 + ENERGY_CUTOFF_VIRT 60.0 + &END BSE + &END GW + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC [angstrom] 6.000 6.000 6.000 + PERIODIC NONE + &END CELL + &KIND H + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q1 + &END KIND + &KIND O + BASIS_SET DZVP-GTH + BASIS_SET RI_AUX RI_DZVP-GTH + POTENTIAL GTH-PBE-q6 + &END KIND + &TOPOLOGY + COORD_FILE_FORMAT xyz + COORD_FILE_NAME H2O_gas.xyz + &CENTER_COORDINATES + &END CENTER_COORDINATES + &END TOPOLOGY + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-bse/TEST_FILES b/tests/QS/regtest-bse/TEST_FILES new file mode 100644 index 0000000000..9301e4979f --- /dev/null +++ b/tests/QS/regtest-bse/TEST_FILES @@ -0,0 +1,3 @@ +TDA_H2O_PBE.inp 109 1e-05 7.7697 +BSE_H2O_PBE.inp 110 1e-05 7.7251 +#EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 74572880a2..178b79cc50 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -361,3 +361,4 @@ Fist/regtest-allegro libtorch Fist/regtest-deepmd deepmd Fist/regtest-16 DFTB/regtest-vdw +QS/regtest-bse libint diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index 7d5f195e46..1e08711e9e 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -1,4 +1,4 @@ -108 +110 Total energy:!3 MD| Potential energy!5 Total energy \[eV\]:!4 @@ -107,6 +107,8 @@ TDDFT+SOC !5 G0W0 direct band gap !6 SCF+SOC direct band gap !6 G0W0+SOC direct band gap !6 +BSE| 1 Singlet State -TDA- !6 +BSE| 1 Singlet State -full- !6 # # these are the tests the can be selected for regtesting. # do regtest will grep for test_grep (first column) and look if the numeric value