From ed08826bade181817ca71e324003a8c2a892c363 Mon Sep 17 00:00:00 2001 From: annahehn <36073704+annahehn@users.noreply.github.com> Date: Mon, 30 Mar 2026 10:15:54 +0200 Subject: [PATCH] Sternheimer BSE 1.Version parallel (#4832) --- src/bse_full_diag.F | 41 +++- src/bse_main.F | 13 +- src/cp_control_types.F | 2 + src/cp_control_utils.F | 2 + src/exstates_types.F | 3 + src/input_cp2k_properties_dft.F | 12 + src/mp2.F | 21 +- src/qs_tddfpt2_bse_utils.F | 255 ++++++++++++++++++-- src/qs_tddfpt2_eigensolver.F | 15 +- src/qs_tddfpt2_fhxc.F | 12 +- src/qs_tddfpt2_utils.F | 6 +- tests/QS/regtest-tddfpt-bse/H2O.inp | 3 - tests/QS/regtest-tddfpt-bse/H2O_W_only.inp | 91 +++++++ tests/QS/regtest-tddfpt-bse/TEST_FILES.toml | 3 +- tests/TEST_DIRS | 4 +- 15 files changed, 423 insertions(+), 60 deletions(-) create mode 100644 tests/QS/regtest-tddfpt-bse/H2O_W_only.inp diff --git a/src/bse_full_diag.F b/src/bse_full_diag.F index 9b9d5952e4..b2d0ad2626 100644 --- a/src/bse_full_diag.F +++ b/src/bse_full_diag.F @@ -111,7 +111,7 @@ CONTAINS REAL(KIND=dp) :: alpha, alpha_screening, 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 + TYPE(cp_fm_type) :: fm_A_copy, fm_W TYPE(dft_control_type), POINTER :: dft_control TYPE(excited_energy_type), POINTER :: ex_env TYPE(tddfpt2_control_type), POINTER :: tddfpt_control @@ -153,6 +153,10 @@ CONTAINS 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) + IF (tddfpt_control%do_bse_w_only) THEN + CALL cp_fm_create(fm_A_copy, fm_struct_A, name="fm_A_iajb") + CALL cp_fm_set_all(fm_A_copy, 0.0_dp) + END IF 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) @@ -162,9 +166,11 @@ CONTAINS ! 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 ((.NOT. tddfpt_control%do_bse) .AND. (.NOT. tddfpt_control%do_bse_w_only)) THEN + 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) + END IF IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb' @@ -197,18 +203,34 @@ CONTAINS 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) + IF (tddfpt_control%do_bse_w_only) THEN + CALL fm_general_add_bse(fm_A_copy, fm_W, -1.0_dp, homo, virtual, & + virtual, virtual, unit_nr, reordering, mp2_env) + END IF END IF !full matrix W is not needed anymore, release it to save memory - IF (tddfpt_control%do_bse) THEN + IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. & + tddfpt_control%do_bse_gw_only) THEN NULLIFY (ex_env) CALL get_qs_env(qs_env, exstate_env=ex_env) - ALLOCATE (ex_env%bse_w_matrix_MO(1, 1)) ! for now only closed-shell - CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_W) - CALL cp_fm_to_fm(fm_W, ex_env%bse_w_matrix_MO(1, 1)) + IF (.NOT. tddfpt_control%do_bse_gw_only) THEN + ALLOCATE (ex_env%bse_w_matrix_MO(1, 1)) ! for now only closed-shell + ALLOCATE (ex_env%bse_a_matrix_MO(1, 1)) ! for now only closed-shell + CALL cp_fm_create(ex_env%bse_w_matrix_MO(1, 1), fm_struct_W) + CALL cp_fm_create(ex_env%bse_a_matrix_MO(1, 1), fm_struct_A) + CALL cp_fm_to_fm(fm_W, ex_env%bse_w_matrix_MO(1, 1)) + IF (tddfpt_control%do_bse_w_only) THEN + CALL cp_fm_to_fm(fm_A_copy, ex_env%bse_a_matrix_MO(1, 1)) + ELSE + CALL cp_fm_to_fm(fm_A, ex_env%bse_a_matrix_MO(1, 1)) + END IF + END IF END IF CALL cp_fm_release(fm_W) + IF (tddfpt_control%do_bse_w_only) CALL cp_fm_release(fm_A_copy) !Now add the energy differences (ε_a-ε_i) on the diagonal (i.e. δ_ij δ_ab) of A_ia,jb + IF (.NOT. tddfpt_control%do_bse) THEN DO ii = 1, nrow_local_A i_row_global = row_indices_A(ii) @@ -226,8 +248,9 @@ CONTAINS END IF END DO END DO + END IF - IF (tddfpt_control%do_bse) THEN + IF (tddfpt_control%do_bse .OR. tddfpt_control%do_bse_w_only .OR. tddfpt_control%do_bse_gw_only) THEN sizeeigen = SIZE(Eigenval) ALLOCATE (ex_env%gw_eigen(sizeeigen)) ! for now only closed-shell ex_env%gw_eigen(:) = Eigenval(:) diff --git a/src/bse_main.F b/src/bse_main.F index bcac5df4e3..0aea47722c 100644 --- a/src/bse_main.F +++ b/src/bse_main.F @@ -26,6 +26,8 @@ MODULE bse_main estimate_BSE_resources,& mult_B_with_W,& truncate_BSE_matrices + USE cp_control_types, ONLY: dft_control_type,& + tddfpt2_control_type USE cp_fm_types, ONLY: cp_fm_release,& cp_fm_type USE cp_log_handling, ONLY: cp_get_default_logger,& @@ -42,7 +44,8 @@ MODULE bse_main USE kinds, ONLY: dp USE message_passing, ONLY: mp_para_env_type USE mp2_types, ONLY: mp2_type - USE qs_environment_types, ONLY: qs_environment_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type #include "./base/base_uses.f90" IMPLICIT NONE @@ -108,7 +111,9 @@ CONTAINS 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(cp_logger_type), POINTER :: logger + TYPE(dft_control_type), POINTER :: dft_control TYPE(mp_para_env_type), POINTER :: para_env + TYPE(tddfpt2_control_type), POINTER :: tddfpt_control CALL timeset(routineN, handle) @@ -222,7 +227,11 @@ CONTAINS 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 + + NULLIFY (dft_control, tddfpt_control) + CALL get_qs_env(qs_env, dft_control=dft_control) + tddfpt_control => dft_control%tddfpt2_control + IF ((my_do_tda) .AND. (.NOT. tddfpt_control%do_bse)) 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, qs_env, mo_coeff) diff --git a/src/cp_control_types.F b/src/cp_control_types.F index bf1090fc6d..ef0ab9dbc1 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -559,6 +559,8 @@ MODULE cp_control_types LOGICAL :: do_smearing = .FALSE. !> dynamical correlation LOGICAL :: do_bse = .FALSE. + LOGICAL :: do_bse_w_only = .FALSE. + LOGICAL :: do_bse_gw_only = .FALSE. ! automatic generation of auxiliary basis for LRI-TDDFT INTEGER :: auto_basis_p_lri_aux = 1 !> use symmetric definition of ADMM Kernel correction diff --git a/src/cp_control_utils.F b/src/cp_control_utils.F index f2db225bed..d21b6d7008 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -1707,6 +1707,8 @@ CONTAINS CALL section_vals_val_get(t_section, "DO_LRIGPW", l_val=t_control%do_lrigpw) CALL section_vals_val_get(t_section, "DO_SMEARING", l_val=t_control%do_smearing) CALL section_vals_val_get(t_section, "DO_BSE", l_val=t_control%do_bse) + CALL section_vals_val_get(t_section, "DO_BSE_W_ONLY", l_val=t_control%do_bse_w_only) + CALL section_vals_val_get(t_section, "DO_BSE_GW_ONLY", l_val=t_control%do_bse_gw_only) CALL section_vals_val_get(t_section, "ADMM_KERNEL_CORRECTION_SYMMETRIC", l_val=t_control%admm_symm) CALL section_vals_val_get(t_section, "ADMM_KERNEL_XC_CORRECTION", l_val=t_control%admm_xc_correction) CALL section_vals_val_get(t_section, "EXCITON_DESCRIPTORS", l_val=t_control%do_exciton_descriptors) diff --git a/src/exstates_types.F b/src/exstates_types.F index 776e66dd8b..482c7abaff 100644 --- a/src/exstates_types.F +++ b/src/exstates_types.F @@ -76,6 +76,7 @@ MODULE exstates_types TYPE(local_rho_type), POINTER :: local_rho_set_admm => NULL() TYPE(wfn_history_type) :: wfn_history = wfn_history_type() TYPE(cp_fm_type), POINTER, DIMENSION(:, :) :: bse_w_matrix_MO => NULL() + TYPE(cp_fm_type), POINTER, DIMENSION(:, :) :: bse_a_matrix_MO => NULL() REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: gw_eigen TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks => NULL() END TYPE excited_energy_type @@ -95,6 +96,7 @@ CONTAINS CALL cp_fm_release(ex_env%cpmos) ! CALL cp_fm_release(ex_env%bse_w_matrix_MO) + CALL cp_fm_release(ex_env%bse_a_matrix_MO) ! CALL exstate_matrix_release(ex_env) ! @@ -208,6 +210,7 @@ CONTAINS NULLIFY (ex_env%evect) NULLIFY (ex_env%cpmos) NULLIFY (ex_env%bse_w_matrix_MO) + NULLIFY (ex_env%bse_a_matrix_MO) IF (excited_state) THEN CALL section_vals_val_get(dft_section, "EXCITED_STATES%STATE", i_val=ex_env%state) CALL section_vals_val_get(dft_section, "EXCITED_STATES%XC_KERNEL_METHOD", & diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index 471bd7915c..a107c7c91c 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -1848,6 +1848,18 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="DO_BSE_W_ONLY", & + description="Debug option for BSE kernel.", & + usage="DO_BSE_W_ONLY", 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="DO_BSE_GW_ONLY", & + description="Debug option for BSE kernel.", & + usage="DO_BSE_GW_ONLY", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + ! LRI subsection CALL create_lrigpw_section(subsection) CALL section_add_subsection(section, subsection) diff --git a/src/mp2.F b/src/mp2.F index 47b5b961a7..779267f31b 100644 --- a/src/mp2.F +++ b/src/mp2.F @@ -235,15 +235,18 @@ CONTAINS END IF ! IF DO_BSE In TDDFT, SAVE ks_matrix to ex_env - NULLIFY (ex_env) - CALL get_qs_env(qs_env, exstate_env=ex_env) - nspins = 1 ! for now only open-shell - CALL dbcsr_allocate_matrix_set(ex_env%matrix_ks, nspins) - DO ispin = 1, nspins - ALLOCATE (ex_env%matrix_ks(ispin)%matrix) - CALL dbcsr_create(ex_env%matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix) - CALL dbcsr_copy(ex_env%matrix_ks(ispin)%matrix, matrix_ks(ispin)%matrix) - END DO + IF (dft_control%tddfpt2_control%do_bse .OR. dft_control%tddfpt2_control%do_bse_w_only .OR. & + dft_control%tddfpt2_control%do_bse_gw_only) THEN + NULLIFY (ex_env) + CALL get_qs_env(qs_env, exstate_env=ex_env) + nspins = 1 ! for now only open-shell + CALL dbcsr_allocate_matrix_set(ex_env%matrix_ks, nspins) + DO ispin = 1, nspins + ALLOCATE (ex_env%matrix_ks(ispin)%matrix) + CALL dbcsr_create(ex_env%matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix) + CALL dbcsr_copy(ex_env%matrix_ks(ispin)%matrix, matrix_ks(ispin)%matrix) + END DO + END IF unit_nr = cp_print_key_unit_nr(logger, input, "DFT%XC%WF_CORRELATION%PRINT", & extension=".mp2Log") diff --git a/src/qs_tddfpt2_bse_utils.F b/src/qs_tddfpt2_bse_utils.F index b5270352b7..e5ee023cbc 100644 --- a/src/qs_tddfpt2_bse_utils.F +++ b/src/qs_tddfpt2_bse_utils.F @@ -17,17 +17,18 @@ MODULE qs_tddfpt2_bse_utils copy_fm_to_dbcsr,& cp_dbcsr_sm_fm_multiply USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& + cp_fm_gemm,& + cp_fm_matvec,& cp_fm_scale_and_add USE cp_fm_struct, ONLY: cp_fm_struct_create,& cp_fm_struct_release,& cp_fm_struct_type - USE cp_fm_types, ONLY: cp_fm_create,& - cp_fm_get_info,& - cp_fm_release,& - cp_fm_set_all,& - cp_fm_set_element,& - cp_fm_to_fm,& - cp_fm_type + USE cp_fm_types, ONLY: & + cp_fm_create, cp_fm_get_info, cp_fm_release, cp_fm_set_all, cp_fm_set_element, & + cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, cp_fm_vectorssum, cp_fm_write_formatted + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_io_unit,& + cp_logger_type USE exstates_types, ONLY: excited_energy_type USE kinds, ONLY: dp USE message_passing, ONLY: mp_para_env_type @@ -49,7 +50,7 @@ MODULE qs_tddfpt2_bse_utils INTEGER, PARAMETER, PRIVATE :: nderivs = 3 INTEGER, PARAMETER, PRIVATE :: maxspins = 2 - PUBLIC:: tddfpt_apply_bse + PUBLIC:: tddfpt_apply_bse_debug, tddfpt_apply_bse PUBLIC:: zeroth_order_gw CONTAINS @@ -192,14 +193,13 @@ CONTAINS ! ************************************************************************************************** !> \brief Update action of TDDFPT operator on trial vectors by adding BSE W term. +!> \brief debug version !> \param Aop_evects ... !> \param evects ... !> \param gs_mos ... !> \param qs_env ... -!> \note Based on the subroutine tddfpt_apply_hfx() which was originally created by -!> Mohamed Fawzi on 10.2002. ! ************************************************************************************************** - SUBROUTINE tddfpt_apply_bse(Aop_evects, evects, gs_mos, qs_env) + SUBROUTINE tddfpt_apply_bse_debug(Aop_evects, evects, gs_mos, qs_env) TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects @@ -207,22 +207,33 @@ CONTAINS INTENT(in) :: gs_mos TYPE(qs_environment_type), POINTER :: qs_env - CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse' + CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse_debug' INTEGER :: a_nao_col, a_virt_col, b_nao_col, c_virt_col, handle, i_occ_row, i_row_global, & - ii, ispin, ivect, j_col_global, j_occ_row, jj, k_occ_col, mu_col_global, nao, ncol_block, & - ncol_local, nrow_block, nrow_local, nspins, nvects, nvirt + ii, iounit, ispin, ivect, j_col_global, j_occ_row, jj, k_occ_col, mu_col_global, nao, & + ncol_block, ncol_block_bse, ncol_block_cs, ncol_local, ncol_local_bse, ncol_local_cs, & + nrow_block, nrow_block_bse, nrow_block_cs, nrow_local, nrow_local_bse, nrow_local_cs, & + nspins, nvects, nvirt INTEGER, DIMENSION(2) :: nactive - INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_bse, & + col_indices_cs, row_indices, & + row_indices_bse, row_indices_cs REAL(KIND=dp) :: alpha + REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), & + POINTER :: my_block, my_bse_w_matrix_MO, my_CSvirt TYPE(cp_blacs_env_type), POINTER :: blacs_env TYPE(cp_fm_struct_type), POINTER :: fmstruct, matrix_struct TYPE(cp_fm_type) :: CSvirt, fms, WXaoao, WXmat2, WXvirtao TYPE(cp_fm_type), POINTER :: bse_w_matrix_MO + TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s TYPE(excited_energy_type), POINTER :: ex_env TYPE(mp_para_env_type), POINTER :: para_env + NULLIFY (logger) !get output_unit + logger => cp_get_default_logger() + iounit = cp_logger_get_default_io_unit(logger) + CALL timeset(routineN, handle) nspins = SIZE(evects, 1) @@ -280,9 +291,16 @@ CONTAINS NULLIFY (row_indices, col_indices) CALL cp_fm_get_info(matrix=WXvirtao, nrow_local=nrow_local, ncol_local=ncol_local, & row_indices=row_indices, col_indices=col_indices, & - nrow_block=nrow_block, ncol_block=ncol_block) + nrow_block=nrow_block, ncol_block=ncol_block, local_data=my_block) + CALL cp_fm_get_info(matrix=CSvirt, nrow_local=nrow_local_cs, ncol_local=ncol_local_cs, & + row_indices=row_indices_cs, col_indices=col_indices_cs, & + nrow_block=nrow_block_cs, ncol_block=ncol_block_cs, local_data=my_CSvirt) + CALL cp_fm_get_info(matrix=bse_w_matrix_MO, nrow_local=nrow_local_bse, ncol_local=ncol_local_bse, & + row_indices=row_indices_bse, col_indices=col_indices_bse, & + nrow_block=nrow_block_bse, ncol_block=ncol_block_bse, local_data=my_bse_w_matrix_MO) CALL cp_fm_set_all(WXvirtao, 0.0_dp) + DO ii = 1, nrow_local i_row_global = row_indices(ii) DO jj = 1, ncol_local @@ -290,7 +308,6 @@ CONTAINS i_occ_row = (i_row_global - 1)/nactive(ispin) + 1 j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1 - a_virt_col = (j_col_global - 1)/nao + 1 b_nao_col = MOD(j_col_global - 1, nao) + 1 @@ -317,7 +334,6 @@ CONTAINS i_occ_row = (i_row_global - 1)/nactive(ispin) + 1 j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1 - a_nao_col = (j_col_global - 1)/nao + 1 b_nao_col = MOD(j_col_global - 1, nao) + 1 @@ -349,19 +365,18 @@ CONTAINS j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1 a_nao_col = (j_col_global - 1)/nao + 1 b_nao_col = MOD(j_col_global - 1, nao) + 1 - IF (a_nao_col == b_nao_col) THEN - WXmat2%local_data(b_nao_col, i_occ_row) = WXmat2%local_data(b_nao_col, i_occ_row) + & - WXaoao%local_data(i_row_global, j_col_global)*evects(ispin, ivect)%local_data(a_nao_col, j_occ_row) - END IF - IF (a_nao_col /= b_nao_col) THEN - WXmat2%local_data(a_nao_col, i_occ_row) = WXmat2%local_data(a_nao_col, i_occ_row) + & + + WXmat2%local_data(a_nao_col, i_occ_row) = WXmat2%local_data(a_nao_col, i_occ_row) + & WXaoao%local_data(i_row_global, j_col_global)*evects(ispin, ivect)%local_data(b_nao_col, j_occ_row) - END IF END DO END DO CALL cp_fm_release(WXaoao) + IF (iounit > 0) THEN + CALL cp_fm_write_formatted(WXmat2, iounit, "WXmat2") + END IF + CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, WXmat2) CALL cp_fm_release(WXmat2) @@ -373,5 +388,195 @@ CONTAINS CALL timestop(handle) + END SUBROUTINE tddfpt_apply_bse_debug + +! ************************************************************************************************** +!> \brief ... +!> \param Aop_evects ... +!> \param evects ... +!> \param gs_mos ... +!> \param qs_env ... +!> \param S_evects ... +! ************************************************************************************************** + SUBROUTINE tddfpt_apply_bse(Aop_evects, evects, gs_mos, qs_env, S_evects) + + TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects + TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects + TYPE(tddfpt_ground_state_mos), DIMENSION(:), & + INTENT(in) :: gs_mos + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: S_evects + + CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse' + + INTEGER :: handle, ispin, ivect, nao, nspins, & + nvects, nvirt + INTEGER, DIMENSION(2) :: nactive + REAL(KIND=dp) :: alpha + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries, rho + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_fm_struct_type), POINTER :: fmstruct + TYPE(cp_fm_type) :: CSvirt, evects_mo, evects_unsplit, fms, & + rhomu, rhosplit + TYPE(cp_fm_type), POINTER :: bse_a_matrix_MO + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s + TYPE(excited_energy_type), POINTER :: ex_env + TYPE(mp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + nspins = SIZE(evects, 1) + nvects = SIZE(evects, 2) + IF (nspins > 1) THEN + alpha = 1.0_dp + ELSE + alpha = 2.0_dp + END IF + CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao) + DO ispin = 1, nspins + CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin)) + END DO + + NULLIFY (ex_env, para_env, blacs_env, matrix_s) + CALL get_qs_env(qs_env, exstate_env=ex_env, para_env=para_env, blacs_env=blacs_env, & + matrix_s=matrix_s) + + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=nao, ncol_global=nao) + CALL cp_fm_create(fms, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fms) + + NULLIFY (bse_a_matrix_MO) + bse_a_matrix_MO => ex_env%bse_a_matrix_MO(1, 1) + para_env => bse_a_matrix_MO%matrix_struct%para_env + + DO ivect = 1, nvects + DO ispin = 1, nspins + NULLIFY (fmstruct) + CALL cp_fm_get_info(matrix=evects(ispin, 1), & + nrow_global=nao, ncol_global=nactive(ispin)) + nvirt = SIZE(gs_mos(ispin)%evals_virt) + + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=nvirt, ncol_global=nao) + CALL cp_fm_create(CSvirt, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL cp_fm_gemm('T', 'N', nvirt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, fms, 0.0_dp, & + CSvirt) + + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, & + nrow_global=nvirt, ncol_global=nactive(ispin)) + CALL cp_fm_create(evects_mo, fmstruct) + CALL cp_fm_struct_release(fmstruct) + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, & + nrow_global=nvirt, ncol_global=nactive(ispin)) + CALL cp_fm_create(rhosplit, fmstruct) + CALL cp_fm_struct_release(fmstruct) + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, & + nrow_global=nao, ncol_global=nactive(ispin)) + CALL cp_fm_create(rhomu, fmstruct) + CALL cp_fm_struct_release(fmstruct) + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, & + nrow_global=nactive(ispin)*nvirt, ncol_global=1) + CALL cp_fm_create(evects_unsplit, fmstruct) + CALL cp_fm_struct_release(fmstruct) + +! get X_jb + CALL parallel_gemm("T", "N", nvirt, nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, & + S_evects(ispin, ivect), 0.0_dp, evects_mo) +! rearrange X_jb + CALL contract_bse(evects_mo, evects_unsplit, nactive(ispin), nvirt, eigvec_entries) +! contract A_iajb X_jb + ALLOCATE (rho(nvirt*nactive(ispin))) + CALL cp_fm_matvec(bse_a_matrix_MO, eigvec_entries, rho, 1.0_dp, 0.0_dp) +! rearrange rho_ia + CALL split_bse(rho, rhosplit, nvirt) +! get rho_imu + CALL parallel_gemm("T", "N", nao, nactive(ispin), nvirt, 1.0_dp, CSvirt, & + rhosplit, 0.0_dp, rhomu) + + CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), 1.0_dp, rhomu) + + CALL cp_fm_release(rhomu) + CALL cp_fm_release(CSvirt) + CALL cp_fm_release(evects_mo) + CALL cp_fm_release(rhosplit) + CALL cp_fm_release(evects_unsplit) + DEALLOCATE (rho) + DEALLOCATE (eigvec_entries) + + END DO! ispin + END DO !ivect + + CALL cp_fm_release(fms) + + CALL timestop(handle) + END SUBROUTINE tddfpt_apply_bse + +! ************************************************************************************************** +!> \brief ... +!> \param evects_mo ... +!> \param evects_unsplit ... +!> \param nactive ... +!> \param nvirt ... +!> \param eigvec_entries ... +! ************************************************************************************************** + SUBROUTINE contract_bse(evects_mo, evects_unsplit, nactive, nvirt, eigvec_entries) + TYPE(cp_fm_type), INTENT(IN) :: evects_mo + TYPE(cp_fm_type), INTENT(INOUT) :: evects_unsplit + INTEGER :: nactive, nvirt + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries + + CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_bse' + + INTEGER :: handle, jj + + CALL timeset(routineN, handle) + + ALLOCATE (eigvec_entries(nvirt*nactive)) + DO jj = 1, nactive + CALL cp_fm_to_fm_submat(msource=evects_mo, mtarget=evects_unsplit, & + nrow=nvirt, ncol=1, s_firstrow=1, s_firstcol=jj, & + t_firstrow=(jj - 1)*nvirt + 1, t_firstcol=1) + END DO + CALL cp_fm_vectorssum(evects_unsplit, eigvec_entries, "R") + + CALL timestop(handle) + END SUBROUTINE contract_bse + +! ************************************************************************************************** +!> \brief ... +!> \param rho ... +!> \param rhosplit ... +!> \param nvirt ... +! ************************************************************************************************** + SUBROUTINE split_bse(rho, rhosplit, nvirt) + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & + INTENT(IN) :: rho + TYPE(cp_fm_type), INTENT(INOUT) :: rhosplit + INTEGER, INTENT(IN) :: nvirt + + CHARACTER(LEN=*), PARAMETER :: routineN = 'split_bse' + + INTEGER :: handle, i_occ_row, ii, j_occ_row + + CALL timeset(routineN, handle) + + CALL cp_fm_set_all(rhosplit, 0.0_dp) + + DO ii = 1, SIZE(rho) + i_occ_row = (ii - 1)/nvirt + 1 + j_occ_row = MOD(ii - 1, nvirt) + 1 + CALL cp_fm_set_element(rhosplit, j_occ_row, i_occ_row, rho(ii)) + END DO + + CALL timestop(handle) + END SUBROUTINE split_bse + END MODULE qs_tddfpt2_bse_utils diff --git a/src/qs_tddfpt2_eigensolver.F b/src/qs_tddfpt2_eigensolver.F index ca7225cd61..dcfa59166d 100644 --- a/src/qs_tddfpt2_eigensolver.F +++ b/src/qs_tddfpt2_eigensolver.F @@ -375,7 +375,8 @@ CONTAINS INTEGER :: handle, ispin, ivect, nspins, nvects INTEGER, DIMENSION(maxspins) :: nmo_occ - LOGICAL :: do_admm, do_bse, do_hfx, & + LOGICAL :: do_admm, do_bse, do_bse_gw_only, & + do_bse_w_only, do_hfx, & do_lri_response, is_rks_triplets, & re_int REAL(KIND=dp) :: rcut, scale @@ -397,7 +398,11 @@ CONTAINS is_rks_triplets = tddfpt_control%rks_triplets do_lri_response = tddfpt_control%do_lrigpw do_bse = tddfpt_control%do_bse + do_bse_w_only = tddfpt_control%do_bse_w_only + do_bse_gw_only = tddfpt_control%do_bse_gw_only IF (do_bse) do_hfx = .FALSE. + IF (do_bse_w_only) do_hfx = .FALSE. + IF (do_bse_gw_only) do_hfx = .FALSE. IF (debug_this_module) THEN CPASSERT(nspins > 0) @@ -477,7 +482,7 @@ CONTAINS END IF ! orbital energy difference term - IF (.NOT. do_bse) THEN + IF ((.NOT. do_bse) .AND. (.NOT. do_bse_gw_only)) THEN CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, & gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control) ELSE @@ -557,9 +562,11 @@ CONTAINS END DO END IF END IF - IF (do_bse) THEN + IF ((do_bse) .OR. (do_bse_w_only)) THEN ! add dynamical screening - CALL tddfpt_apply_bse(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env) + ! CALL tddfpt_apply_bse_debug(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env) + CALL tddfpt_apply_bse(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env, & + S_evects=S_evects) END IF END IF diff --git a/src/qs_tddfpt2_fhxc.F b/src/qs_tddfpt2_fhxc.F index 0e2ed36b5a..6c4ba207e8 100644 --- a/src/qs_tddfpt2_fhxc.F +++ b/src/qs_tddfpt2_fhxc.F @@ -245,8 +245,8 @@ CONTAINS ! Skip kernel if collinear xc-kernel for spin-flip is requested IF (spinflip /= tddfpt_sf_col) THEN - - IF (.NOT. dft_control%tddfpt2_control%do_bse) THEN + IF ((.NOT. dft_control%tddfpt2_control%do_bse) .AND. (.NOT. dft_control%tddfpt2_control%do_bse_w_only)) THEN + IF ((.NOT. dft_control%tddfpt2_control%do_bse_gw_only)) THEN ! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2 ! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation IF (gapw_xc) THEN @@ -293,10 +293,13 @@ CONTAINS do_sf=do_noncol) END IF + END IF ! do_bse END IF ! do_bse END IF ! spin-flip ! ADMM correction + IF ((.NOT. dft_control%tddfpt2_control%do_bse) .AND. (.NOT. dft_control%tddfpt2_control%do_bse_w_only) & + .AND. (.NOT. dft_control%tddfpt2_control%do_bse_gw_only)) THEN IF (do_admm .AND. admm_xc_correction) THEN IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, & @@ -410,8 +413,11 @@ CONTAINS END IF END IF END IF + END IF ! electron-hole Coulomb interaction + IF (.NOT. dft_control%tddfpt2_control%do_bse_w_only) THEN + IF (.NOT. dft_control%tddfpt2_control%do_bse_gw_only) THEN IF ((.NOT. is_rks_triplets) .AND. (spinflip == no_sf_tddfpt)) THEN ! a sum J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu can be computed by solving ! the Poisson equation for combined density (rho_{ia,alpha} + rho_{ia,beta}) . @@ -499,6 +505,8 @@ CONTAINS ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp) END IF END DO + END IF + END IF END DO CALL timestop(handle) diff --git a/src/qs_tddfpt2_utils.F b/src/qs_tddfpt2_utils.F index 6236e8beaf..9af2361c7e 100644 --- a/src/qs_tddfpt2_utils.F +++ b/src/qs_tddfpt2_utils.F @@ -133,7 +133,8 @@ CONTAINS CALL get_qs_env(qs_env, blacs_env=blacs_env, cell=cell, dft_control=dft_control, & matrix_ks=matrix_ks, matrix_s=matrix_s, mos=mos, scf_env=scf_env) tddfpt_control => dft_control%tddfpt2_control - IF (tddfpt_control%do_bse) THEN + IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. & + (tddfpt_control%do_bse_gw_only)) THEN NULLIFY (ks_env, ex_env) CALL get_qs_env(qs_env, exstate_env=ex_env, ks_env=ks_env) CALL dbcsr_copy(matrix_ks(1)%matrix, ex_env%matrix_ks(1)%matrix) @@ -1031,7 +1032,8 @@ CONTAINS nv = nmo_virt_selected(ispin) ALLOCATE (ev_virt(nv), ev_occ(no)) ! if do_bse and do_gw, take gw zeroth order - IF (tddfpt_control%do_bse) THEN + IF ((tddfpt_control%do_bse) .OR. (tddfpt_control%do_bse_w_only) .OR. & + (tddfpt_control%do_bse_gw_only)) THEN ev_virt(1:nv) = ex_env%gw_eigen(nmo(ispin) + 1:nmo(ispin) + nv) DO i = 1, no j = nmo_occ_avail(ispin) - i + 1 diff --git a/tests/QS/regtest-tddfpt-bse/H2O.inp b/tests/QS/regtest-tddfpt-bse/H2O.inp index daaad581fc..ba15c110f1 100644 --- a/tests/QS/regtest-tddfpt-bse/H2O.inp +++ b/tests/QS/regtest-tddfpt-bse/H2O.inp @@ -2,9 +2,6 @@ PRINT_LEVEL MEDIUM PROJECT BSE_H2O RUN_TYPE ENERGY - &TIMINGS - THRESHOLD 0.01 - &END TIMINGS &END GLOBAL &FORCE_EVAL diff --git a/tests/QS/regtest-tddfpt-bse/H2O_W_only.inp b/tests/QS/regtest-tddfpt-bse/H2O_W_only.inp new file mode 100644 index 0000000000..b2a2d109b5 --- /dev/null +++ b/tests/QS/regtest-tddfpt-bse/H2O_W_only.inp @@ -0,0 +1,91 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT BSE_H2O + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME HFX_BASIS + BASIS_SET_FILE_NAME BASIS_RI_cc-TZ + POTENTIAL_FILE_NAME GTH_POTENTIALS + &MGRID + CUTOFF 500 + REL_CUTOFF 60 + &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 + ## SELF_CONSISTENCY G0W0 + &BSE + BSE_DIAG_METHOD FULLDIAG + TDA ON + &END BSE + &END GW + &END RI_RPA + &END WF_CORRELATION + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &PROPERTIES + &TDDFPT + CONVERGENCE [eV] 1.0e-8 + DO_BSE_W_ONLY T ! THIS IS JUST FOR DEBUGGING ! + KERNEL FULL + MAX_ITER 20 + MAX_KV 900 + NLUMO 19 + NSTATES 10 + &DIPOLE_MOMENTS + DIPOLE_FORM LENGTH + &END DIPOLE_MOMENTS + &END TDDFPT + &END PROPERTIES + &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-tddfpt-bse/TEST_FILES.toml b/tests/QS/regtest-tddfpt-bse/TEST_FILES.toml index bb89667919..eb4773a25b 100644 --- a/tests/QS/regtest-tddfpt-bse/TEST_FILES.toml +++ b/tests/QS/regtest-tddfpt-bse/TEST_FILES.toml @@ -3,5 +3,6 @@ # e.g. 0 means do not compare anything, running is enough # 1 compares the last total energy in the file # for details see cp2k/tools/do_regtest -"H2O.inp" = [{matcher="M037", tol=2.0E-05, ref=2.14632}] +"H2O.inp" = [{matcher="M037", tol=2.0E-04, ref=2.14632}] +"H2O_W_only.inp" = [{matcher="M037", tol=2.0E-05, ref=0.952977}] #EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 83ad361b90..8b1928b076 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -410,6 +410,4 @@ QS/regtest-trexio trexio QS/regtest-trexio-2 trexio QS/regtest-rtbse-gxac libint greenx !ifx QS/regtest-cneo - -# TODO: Re-enable once problems with MPI are resolved (https://github.com/cp2k/cp2k/pull/4445) -# QS/regtest-tddfpt-bse libint libxc +QS/regtest-tddfpt-bse libint libxc