From fdff660ab3f8c59f8cf5ba5dac2a7e7ff49db258 Mon Sep 17 00:00:00 2001 From: Anna Hehn Date: Mon, 21 Jul 2025 14:57:46 +0200 Subject: [PATCH] Sternheimer BSE (1. Version with W matrix given) --- src/CMakeLists.txt | 1 + src/bse_full_diag.F | 32 +- src/bse_main.F | 4 +- src/cp_control_types.F | 2 + src/cp_control_utils.F | 1 + src/exstates_types.F | 12 + src/gw_utils.F | 4 +- src/input_cp2k_properties_dft.F | 6 + src/mp2.F | 20 +- src/qs_tddfpt2_bse_utils.F | 375 ++++++++++++++++++++ src/qs_tddfpt2_eigensolver.F | 66 ++-- src/qs_tddfpt2_fhxc.F | 82 ++--- src/qs_tddfpt2_methods.F | 5 +- src/qs_tddfpt2_utils.F | 61 +++- tests/QS/regtest-tddfpt-bse/H2O.inp | 94 +++++ tests/QS/regtest-tddfpt-bse/H2O_gas.xyz | 5 + tests/QS/regtest-tddfpt-bse/TEST_FILES.toml | 7 + tests/TEST_DIRS | 1 + 18 files changed, 688 insertions(+), 90 deletions(-) create mode 100644 src/qs_tddfpt2_bse_utils.F create mode 100644 tests/QS/regtest-tddfpt-bse/H2O.inp create mode 100644 tests/QS/regtest-tddfpt-bse/H2O_gas.xyz create mode 100644 tests/QS/regtest-tddfpt-bse/TEST_FILES.toml diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 19618a2661..d2baf2a32b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -723,6 +723,7 @@ list( qs_tddfpt2_soc_types.F qs_tddfpt2_stda_types.F qs_tddfpt2_stda_utils.F + qs_tddfpt2_bse_utils.F qs_tddfpt2_smearing_methods.F qs_tddfpt2_subgroups.F qs_tddfpt2_types.F diff --git a/src/bse_full_diag.F b/src/bse_full_diag.F index 9a486e5dbd..455e63f7b6 100644 --- a/src/bse_full_diag.F +++ b/src/bse_full_diag.F @@ -29,6 +29,8 @@ MODULE bse_full_diag USE cp_blacs_env, ONLY: cp_blacs_env_create,& cp_blacs_env_release,& cp_blacs_env_type + USE cp_control_types, ONLY: dft_control_type,& + tddfpt2_control_type USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add USE cp_fm_diag, ONLY: choose_eigv_solver,& cp_fm_power @@ -41,6 +43,7 @@ MODULE bse_full_diag cp_fm_set_all,& cp_fm_to_fm,& cp_fm_type + USE exstates_types, ONLY: excited_energy_type USE input_constants, ONLY: bse_screening_alpha,& bse_screening_rpa,& bse_singlet,& @@ -49,7 +52,8 @@ MODULE bse_full_diag USE message_passing, ONLY: mp_para_env_type USE mp2_types, ONLY: mp2_type USE parallel_gemm_api, ONLY: parallel_gemm - 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 @@ -81,11 +85,12 @@ CONTAINS !> \param dimen_RI ... !> \param mp2_env ... !> \param para_env ... +!> \param qs_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) + para_env, qs_env) TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, & fm_mat_S_ab_bse @@ -94,21 +99,29 @@ CONTAINS 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 + TYPE(qs_environment_type), POINTER :: qs_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 + ncol_local_A, nrow_local_A, sizeeigen INTEGER, DIMENSION(4) :: reordering INTEGER, DIMENSION(:), POINTER :: col_indices_A, row_indices_A 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(dft_control_type), POINTER :: dft_control + TYPE(excited_energy_type), POINTER :: ex_env + TYPE(tddfpt2_control_type), POINTER :: tddfpt_control CALL timeset(routineN, handle) + NULLIFY (dft_control, tddfpt_control) + CALL get_qs_env(qs_env, dft_control=dft_control) + tddfpt_control => dft_control%tddfpt2_control + IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A' END IF @@ -186,6 +199,13 @@ CONTAINS virtual, virtual, unit_nr, reordering, mp2_env) END IF !full matrix W is not needed anymore, release it to save memory + IF (tddfpt_control%do_bse) 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)) + END IF CALL cp_fm_release(fm_W) !Now add the energy differences (ε_a-ε_i) on the diagonal (i.e. δ_ij δ_ab) of A_ia,jb @@ -207,6 +227,12 @@ CONTAINS END DO END DO + IF (tddfpt_control%do_bse) THEN + sizeeigen = SIZE(Eigenval) + ALLOCATE (ex_env%gw_eigen(sizeeigen)) ! for now only closed-shell + ex_env%gw_eigen(:) = Eigenval(:) + END IF + CALL cp_fm_struct_release(fm_struct_A) CALL cp_fm_struct_release(fm_struct_W) diff --git a/src/bse_main.F b/src/bse_main.F index ff4cec873d..239bde441f 100644 --- a/src/bse_main.F +++ b/src/bse_main.F @@ -190,12 +190,12 @@ CONTAINS CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, & fm_A_BSE, Eigenval_reduced, unit_nr, & homo_red, virtual_red, dimen_RI, mp2_env, & - para_env) + para_env, qs_env) ELSE 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) + para_env, qs_env) END IF IF (my_do_abba) THEN ! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W) diff --git a/src/cp_control_types.F b/src/cp_control_types.F index 86a5f50e9e..6f2408fd1b 100644 --- a/src/cp_control_types.F +++ b/src/cp_control_types.F @@ -532,6 +532,8 @@ MODULE cp_control_types LOGICAL :: do_lrigpw = .FALSE. !> smeared occupation LOGICAL :: do_smearing = .FALSE. + !> dynamical correlation + LOGICAL :: do_bse = .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 91517280b0..d4775a3e3a 100644 --- a/src/cp_control_utils.F +++ b/src/cp_control_utils.F @@ -1687,6 +1687,7 @@ CONTAINS CALL section_vals_val_get(t_section, "RKS_TRIPLETS", l_val=t_control%rks_triplets) 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, "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 888afb687a..5f0fe3692b 100644 --- a/src/exstates_types.F +++ b/src/exstates_types.F @@ -75,6 +75,9 @@ MODULE exstates_types TYPE(local_rho_type), POINTER :: local_rho_set => NULL() 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() + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: gw_eigen + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks => NULL() END TYPE excited_energy_type CONTAINS @@ -91,6 +94,8 @@ CONTAINS CALL cp_fm_release(ex_env%evect) CALL cp_fm_release(ex_env%cpmos) ! + CALL cp_fm_release(ex_env%bse_w_matrix_MO) + ! CALL exstate_matrix_release(ex_env) ! CALL exstate_potential_release(ex_env) @@ -98,6 +103,9 @@ CONTAINS CALL cp_fm_release(ex_env%wfn_history%evect) CALL cp_fm_release(ex_env%wfn_history%cpmos) + IF (ALLOCATED(ex_env%gw_eigen)) & + DEALLOCATE (ex_env%gw_eigen) + DEALLOCATE (ex_env) END IF @@ -112,6 +120,8 @@ CONTAINS TYPE(excited_energy_type), POINTER :: ex_env IF (ASSOCIATED(ex_env)) THEN + IF (ASSOCIATED(ex_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_ks) + NULLIFY (ex_env%matrix_ks) IF (ASSOCIATED(ex_env%matrix_pe)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_pe) NULLIFY (ex_env%matrix_pe) IF (ASSOCIATED(ex_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_hz) @@ -196,6 +206,8 @@ CONTAINS ALLOCATE (ex_env) ex_env%evalue = 0.0_dp NULLIFY (ex_env%evect) + NULLIFY (ex_env%cpmos) + NULLIFY (ex_env%bse_w_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/gw_utils.F b/src/gw_utils.F index 781d5005b0..56327e72e1 100644 --- a/src/gw_utils.F +++ b/src/gw_utils.F @@ -1190,8 +1190,8 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'check_for_restart_files' CHARACTER(LEN=9) :: frmt - CHARACTER(LEN=default_path_length) :: project_name - CHARACTER(LEN=default_string_length) :: f_chi, f_S_n, f_S_p, f_S_x, f_W_t, prefix + CHARACTER(len=default_path_length) :: project_name + CHARACTER(len=default_string_length) :: f_chi, f_s_n, f_s_p, f_s_x, f_w_t, prefix INTEGER :: handle, i_spin, i_t_or_w, ind, n_spin, & num_time_freq_points LOGICAL :: chi_exists, Sigma_neg_time_exists, & diff --git a/src/input_cp2k_properties_dft.F b/src/input_cp2k_properties_dft.F index ced89c25a4..9a5601f381 100644 --- a/src/input_cp2k_properties_dft.F +++ b/src/input_cp2k_properties_dft.F @@ -1805,6 +1805,12 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL keyword_create(keyword, __LOCATION__, name="DO_BSE", & + description="Choosing BSE kernel.", & + usage="DO_BSE", 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 d768516cb3..66467e8667 100644 --- a/src/mp2.F +++ b/src/mp2.F @@ -26,9 +26,12 @@ MODULE mp2 cite_reference USE cp_blacs_env, ONLY: cp_blacs_env_type USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_api, ONLY: dbcsr_get_info,& + USE cp_dbcsr_api, ONLY: dbcsr_copy,& + dbcsr_create,& + dbcsr_get_info,& dbcsr_p_type - USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm + USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& + dbcsr_allocate_matrix_set USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& cp_fm_syrk,& cp_fm_triangular_invert,& @@ -48,6 +51,7 @@ MODULE mp2 cp_logger_type USE cp_output_handling, ONLY: cp_print_key_finished_output,& cp_print_key_unit_nr + USE exstates_types, ONLY: excited_energy_type USE hfx_exx, ONLY: calculate_exx USE hfx_types, ONLY: & alloc_containers, dealloc_containers, hfx_basis_info_type, hfx_basis_type, & @@ -152,6 +156,7 @@ CONTAINS TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_transl, matrix_s_kp TYPE(dft_control_type), POINTER :: dft_control + TYPE(excited_energy_type), POINTER :: ex_env TYPE(hfx_basis_info_type) :: basis_info TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers @@ -229,6 +234,17 @@ 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 + 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 new file mode 100644 index 0000000000..ee40a90b4e --- /dev/null +++ b/src/qs_tddfpt2_bse_utils.F @@ -0,0 +1,375 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +MODULE qs_tddfpt2_bse_utils + USE cp_blacs_env, ONLY: cp_blacs_env_type + USE cp_dbcsr_api, ONLY: dbcsr_create,& + dbcsr_p_type,& + dbcsr_release,& + dbcsr_set,& + dbcsr_type + USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl + USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& + copy_fm_to_dbcsr,& + cp_dbcsr_sm_fm_multiply + USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& + cp_fm_gemm,& + 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 exstates_types, ONLY: excited_energy_type + USE kinds, ONLY: dp + USE message_passing, ONLY: mp_para_env_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_bse_utils' + + LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. + ! number of first derivative components (3: d/dx, d/dy, d/dz) + INTEGER, PARAMETER, PRIVATE :: nderivs = 3 + INTEGER, PARAMETER, PRIVATE :: maxspins = 2 + + PUBLIC:: tddfpt_apply_bse + PUBLIC:: zeroth_order_gw + +CONTAINS +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param Aop_evects ... +!> \param evects ... +!> \param S_evects ... +!> \param gs_mos ... +!> \param matrix_s ... +!> \param matrix_ks ... +! ************************************************************************************************** + SUBROUTINE zeroth_order_gw(qs_env, Aop_evects, evects, S_evects, gs_mos, matrix_s, matrix_ks) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects + TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects, S_evects + TYPE(tddfpt_ground_state_mos), DIMENSION(:), & + INTENT(in) :: gs_mos + TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s + TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks + + CHARACTER(LEN=*), PARAMETER :: routineN = 'zeroth_order_gw' + + INTEGER :: handle, i, ispin, ivect, nactive, nao, & + nmo, nspins, nvects, occ, virt + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: gw_occ, gw_virt + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_fm_struct_type), POINTER :: fmstruct, matrix_struct + TYPE(cp_fm_type) :: fms, hevec, matrixtmp, matrixtmp2, & + matrixtmp3, Sweighted_vect, & + weighted_vect + TYPE(dbcsr_type) :: matrixf + TYPE(excited_energy_type), POINTER :: ex_env + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb + + CALL timeset(routineN, handle) + + NULLIFY (ex_env, sab_orb) + CALL get_qs_env(qs_env, exstate_env=ex_env, sab_orb=sab_orb) + + nspins = SIZE(matrix_ks, 1) + nspins = SIZE(evects, 1) + nvects = SIZE(evects, 2) + + DO ispin = 1, nspins + + CALL dbcsr_create(matrixf, template=matrix_s) + nmo = SIZE(ex_env%gw_eigen) + CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, & + nrow_global=nao, ncol_global=nactive) + NULLIFY (blacs_env, para_env) + CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env) + + occ = SIZE(gs_mos(ispin)%evals_occ) + nmo = SIZE(ex_env%gw_eigen) + virt = SIZE(gs_mos(ispin)%evals_virt) + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=virt, ncol_global=virt) + CALL cp_fm_create(matrixtmp, fmstruct) + CALL cp_fm_struct_release(fmstruct) + + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=virt, ncol_global=nao) + CALL cp_fm_create(matrixtmp2, fmstruct) + CALL cp_fm_struct_release(fmstruct) + + NULLIFY (fmstruct) + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=nao, ncol_global=nao) + CALL cp_fm_create(matrixtmp3, fmstruct) + CALL cp_fm_create(fms, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL cp_dbcsr_alloc_block_from_nbl(matrixf, sab_orb) + +!--add virt eigenvalues + CALL dbcsr_set(matrixf, 0.0_dp) + CALL cp_fm_create(weighted_vect, gs_mos(ispin)%mos_virt%matrix_struct) + CALL cp_fm_create(Sweighted_vect, gs_mos(ispin)%mos_virt%matrix_struct) + CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, weighted_vect) + CALL copy_dbcsr_to_fm(matrix_s, fms) + + ALLOCATE (gw_virt(virt)) + ALLOCATE (gw_occ(occ)) + gw_virt(1:virt) = ex_env%gw_eigen(occ + 1:nmo) + gw_occ(1:occ) = ex_env%gw_eigen(1:occ) + + CALL cp_fm_set_all(matrixtmp, 0.0_dp) + DO i = 1, virt + CALL cp_fm_set_element(matrixtmp, i, i, gw_virt(i)) + END DO + DEALLOCATE (gw_virt) + CALL cp_fm_gemm('N', 'N', nao, virt, nao, 1.0_dp, fms, weighted_vect, 0.0_dp, & + Sweighted_vect) + CALL cp_fm_gemm('N', 'T', virt, nao, virt, 1.0_dp, matrixtmp, Sweighted_vect, 0.0_dp, & + matrixtmp2) + CALL cp_fm_gemm('N', 'N', nao, nao, virt, 1.0_dp, Sweighted_vect, matrixtmp2, 0.0_dp, & + matrixtmp3) + CALL copy_fm_to_dbcsr(matrixtmp3, matrixf) + + CALL cp_fm_release(weighted_vect) + CALL cp_fm_release(Sweighted_vect) + CALL cp_fm_release(fmS) +!--add occ eigenvalues + CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, & + nrow_global=nao, ncol_global=nactive) + CALL cp_fm_create(hevec, matrix_struct) + + DO ivect = 1, nvects + CALL cp_dbcsr_sm_fm_multiply(matrixf, evects(ispin, ivect), & + Aop_evects(ispin, ivect), ncol=nactive, & + alpha=1.0_dp, beta=1.0_dp) + + CALL cp_fm_to_fm(S_evects(ispin, ivect), hevec) + CALL cp_fm_column_scale(hevec, gw_occ) + + CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, hevec) + END DO !ivect + DEALLOCATE (gw_occ) + + CALL cp_fm_release(matrixtmp) + CALL cp_fm_release(matrixtmp2) + CALL cp_fm_release(matrixtmp3) + + CALL dbcsr_release(matrixf) + CALL cp_fm_release(hevec) + END DO !ispin + + virt = SIZE(Aop_evects, 2) + CALL timestop(handle) + + END SUBROUTINE zeroth_order_gw + +! ************************************************************************************************** +!> \brief Update action of TDDFPT operator on trial vectors by adding BSE W term. +!> \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) + + 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 + + CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse' + + 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 + INTEGER, DIMENSION(2) :: nactive + INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices + REAL(KIND=dp) :: alpha + 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(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_w_matrix_MO) + bse_w_matrix_MO => ex_env%bse_w_matrix_MO(1, 1) + + DO ivect = 1, nvects + DO ispin = 1, nspins + NULLIFY (matrix_struct, fmstruct) + CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, & + 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_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=nactive(ispin)*nactive(ispin), & + ncol_global=nvirt*nao) + CALL cp_fm_create(WXvirtao, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL cp_fm_set_all(WXvirtao, 0.0_dp) + + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=nactive(ispin)*nactive(ispin), & + ncol_global=nao*nao) + CALL cp_fm_create(WXaoao, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL cp_fm_set_all(WXaoao, 0.0_dp) + + CALL cp_fm_gemm('T', 'N', nvirt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, fms, 0.0_dp, & + CSvirt) + 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) + + 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 + j_col_global = col_indices(jj) + + 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 + + DO c_virt_col = 1, nvirt + mu_col_global = (a_virt_col - 1)*nvirt + c_virt_col + + WXvirtao%local_data(i_row_global, j_col_global) = WXvirtao%local_data(i_row_global, j_col_global) + & + bse_w_matrix_MO%local_data(i_row_global, mu_col_global)*CSvirt%local_data(c_virt_col, b_nao_col) + + END DO + END DO + END DO + + NULLIFY (row_indices, col_indices) ! redefine indices + CALL cp_fm_get_info(matrix=WXaoao, nrow_local=nrow_local, ncol_local=ncol_local, & + row_indices=row_indices, col_indices=col_indices, & + nrow_block=nrow_block, ncol_block=ncol_block) + + CALL cp_fm_set_all(WXaoao, 0.0_dp) + DO ii = 1, nrow_local + i_row_global = row_indices(ii) + DO jj = 1, ncol_local + j_col_global = col_indices(jj) + + 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 + + DO k_occ_col = 1, nvirt + mu_col_global = (k_occ_col - 1)*nao + a_nao_col + + WXaoao%local_data(i_row_global, j_col_global) = WXaoao%local_data(i_row_global, j_col_global) + & + WXvirtao%local_data(i_row_global, mu_col_global)*CSvirt%local_data(k_occ_col, b_nao_col) + + END DO + END DO + END DO + + CALL cp_fm_release(WXvirtao) + CALL cp_fm_release(CSvirt) + + CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, & + context=blacs_env, nrow_global=nao, ncol_global=nactive(ispin)) + CALL cp_fm_create(WXmat2, fmstruct) + CALL cp_fm_struct_release(fmstruct) + CALL cp_fm_set_all(WXmat2, 0.0_dp) + + DO ii = 1, nrow_local + i_row_global = row_indices(ii) + DO jj = 1, ncol_local + j_col_global = col_indices(jj) + + 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 + IF (a_nao_col .EQ. 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 .NE. b_nao_col) THEN + 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) + + CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, WXmat2) + + CALL cp_fm_release(WXmat2) + + END DO! ispin + END DO !ivect + + CALL cp_fm_release(fms) + + CALL timestop(handle) + + END SUBROUTINE tddfpt_apply_bse +END MODULE qs_tddfpt2_bse_utils diff --git a/src/qs_tddfpt2_eigensolver.F b/src/qs_tddfpt2_eigensolver.F index 19ab144676..7bd318964c 100644 --- a/src/qs_tddfpt2_eigensolver.F +++ b/src/qs_tddfpt2_eigensolver.F @@ -45,6 +45,8 @@ MODULE qs_tddfpt2_eigensolver USE qs_kernel_types, ONLY: full_kernel_env_type,& kernel_env_type USE qs_scf_methods, ONLY: eigensolver + USE qs_tddfpt2_bse_utils, ONLY: tddfpt_apply_bse,& + zeroth_order_gw USE qs_tddfpt2_fhxc, ONLY: fhxc_kernel,& stda_kernel USE qs_tddfpt2_operators, ONLY: tddfpt_apply_energy_diff,& @@ -373,8 +375,9 @@ CONTAINS INTEGER :: handle, ispin, ivect, nspins, nvects INTEGER, DIMENSION(maxspins) :: nmo_occ - LOGICAL :: do_admm, do_hfx, do_lri_response, & - is_rks_triplets, re_int + LOGICAL :: do_admm, do_bse, do_hfx, & + do_lri_response, is_rks_triplets, & + re_int REAL(KIND=dp) :: rcut, scale TYPE(cp_fm_type) :: fm_dummy TYPE(full_kernel_env_type), POINTER :: kernel_env_admm_aux @@ -393,6 +396,8 @@ CONTAINS END IF is_rks_triplets = tddfpt_control%rks_triplets do_lri_response = tddfpt_control%do_lrigpw + do_bse = tddfpt_control%do_bse + IF (do_bse) do_hfx = .FALSE. IF (debug_this_module) THEN CPASSERT(nspins > 0) @@ -413,17 +418,17 @@ CONTAINS DO ivect = 1, nvects DO ispin = 1, SIZE(evects, 1) ASSOCIATE (evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect)) - IF (ASSOCIATED(evect%matrix_struct)) THEN - IF (ASSOCIATED(work_matrix%matrix_struct)) THEN - CALL cp_fm_copy_general(evect, work_matrix, para_env) - ELSE - CALL cp_fm_copy_general(evect, fm_dummy, para_env) - END IF - ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN - CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env) - ELSE - CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env) - END IF + IF (ASSOCIATED(evect%matrix_struct)) THEN + IF (ASSOCIATED(work_matrix%matrix_struct)) THEN + CALL cp_fm_copy_general(evect, work_matrix, para_env) + ELSE + CALL cp_fm_copy_general(evect, fm_dummy, para_env) + END IF + ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN + CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env) + ELSE + CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env) + END IF END ASSOCIATE END DO END DO @@ -455,25 +460,30 @@ CONTAINS DO ispin = 1, SIZE(evects, 1) ASSOCIATE (Aop_evect => Aop_evects(ispin, ivect), & work_matrix => work_matrices%Aop_evects_sub(ispin, ivect)) - IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN - IF (ASSOCIATED(work_matrix%matrix_struct)) THEN - CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env) - ELSE - CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env) - END IF - ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN - CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env) - ELSE - CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env) - END IF + IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN + IF (ASSOCIATED(work_matrix%matrix_struct)) THEN + CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env) + ELSE + CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env) + END IF + ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN + CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env) + ELSE + CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env) + END IF END ASSOCIATE END DO END DO END IF ! orbital energy difference term - 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) + IF (.NOT. do_bse) 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 + CALL zeroth_order_gw(qs_env=qs_env, Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, & + gs_mos=gs_mos, matrix_s=matrix_s, matrix_ks=matrix_ks) + END IF ! if smeared occupation, then add aCCSX here IF (tddfpt_control%do_smearing) THEN @@ -547,6 +557,10 @@ CONTAINS END DO END IF END IF + IF (do_bse) THEN + ! add dynamical screening + CALL tddfpt_apply_bse(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env) + END IF END IF diff --git a/src/qs_tddfpt2_fhxc.F b/src/qs_tddfpt2_fhxc.F index 5252d78f94..7158b0347f 100644 --- a/src/qs_tddfpt2_fhxc.F +++ b/src/qs_tddfpt2_fhxc.F @@ -242,50 +242,52 @@ CONTAINS ! Skip kernel if collinear xc-kernel for spin-flip is requested IF (spinflip /= tddfpt_sf_col) 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 - IF (kernel_env%do_exck) THEN - CPABORT("NYA") + IF (.NOT. dft_control%tddfpt2_control%do_bse) 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 + IF (kernel_env%do_exck) THEN + CPABORT("NYA") + ELSE + CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, & + rho_ia_struct=work_matrices%rho_xc_struct_sub, & + is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, & + work_v_xc=work_matrices%wpw_rspace_sub, & + work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, & + spinflip=spinflip) + END IF + DO ispin = 1, nspins + CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), & + work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol) + CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), & + hmat=work_matrices%A_ia_munu_sub(ispin), & + qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw_xc, & + pw_env_external=sub_env%pw_env, & + task_list_external=sub_env%task_list_orb_soft) + CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin)) + END DO ELSE - CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, & - rho_ia_struct=work_matrices%rho_xc_struct_sub, & - is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, & - work_v_xc=work_matrices%wpw_rspace_sub, & - work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, & - spinflip=spinflip) + IF (kernel_env%do_exck) THEN + CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, & + work_matrices%rho_orb_struct_sub, is_rks_triplets) + ELSE + CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, & + rho_ia_struct=work_matrices%rho_orb_struct_sub, & + is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, & + work_v_xc=work_matrices%wpw_rspace_sub, & + work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, & + spinflip=spinflip) + END IF END IF - DO ispin = 1, nspins - CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), & - work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol) - CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), & - hmat=work_matrices%A_ia_munu_sub(ispin), & - qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw_xc, & - pw_env_external=sub_env%pw_env, & - task_list_external=sub_env%task_list_orb_soft) - CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin)) - END DO - ELSE - IF (kernel_env%do_exck) THEN - CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, & - work_matrices%rho_orb_struct_sub, is_rks_triplets) - ELSE - CALL tddfpt_apply_xc(A_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, & - rho_ia_struct=work_matrices%rho_orb_struct_sub, & - is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, & - work_v_xc=work_matrices%wpw_rspace_sub, & - work_v_xc_tau=work_matrices%wpw_tau_rspace_sub, & - spinflip=spinflip) + IF (gapw .OR. gapw_xc) THEN + rho_atom_set => sub_env%local_rho_set%rho_atom_set + rho1_atom_set => work_matrices%local_rho_set%rho_atom_set + CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, kernel_env%xc_section, & + sub_env%para_env, do_tddfpt2=.TRUE., do_triplet=is_rks_triplets, & + do_sf=do_noncol) END IF - END IF - IF (gapw .OR. gapw_xc) THEN - rho_atom_set => sub_env%local_rho_set%rho_atom_set - rho1_atom_set => work_matrices%local_rho_set%rho_atom_set - CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, kernel_env%xc_section, & - sub_env%para_env, do_tddfpt2=.TRUE., do_triplet=is_rks_triplets, & - do_sf=do_noncol) - END IF + END IF ! do_bse END IF ! spin-flip ! ADMM correction diff --git a/src/qs_tddfpt2_methods.F b/src/qs_tddfpt2_methods.F index 3cdd0d68be..c59a468c0c 100644 --- a/src/qs_tddfpt2_methods.F +++ b/src/qs_tddfpt2_methods.F @@ -1087,8 +1087,9 @@ CONTAINS ! according to their excitation energies log_unit = cp_print_key_unit_nr(logger, tddfpt_print_section, & "GUESS_VECTORS", extension=".tddfptLog") - CALL tddfpt_guess_vectors(tddfpt_control, nspins, evects=evects, evals=evals, & - gs_mos=gs_mos, log_unit=log_unit) + CALL tddfpt_guess_vectors(evects=evects, evals=evals, & + gs_mos=gs_mos, log_unit=log_unit, tddfpt_control=tddfpt_control, & + qs_env=qs_env, nspins=nspins) CALL cp_print_key_finished_output(log_unit, logger, & tddfpt_print_section, "GUESS_VECTORS") diff --git a/src/qs_tddfpt2_utils.F b/src/qs_tddfpt2_utils.F index 5dc6c53c30..b7ff19b5fa 100644 --- a/src/qs_tddfpt2_utils.F +++ b/src/qs_tddfpt2_utils.F @@ -37,6 +37,7 @@ MODULE qs_tddfpt2_utils 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 input_constants, ONLY: & cholesky_dbcsr, cholesky_inverse, cholesky_off, cholesky_restore, no_sf_tddfpt, oe_gllb, & oe_lb, oe_none, oe_saop, oe_shift @@ -56,6 +57,8 @@ MODULE qs_tddfpt2_utils USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix + USE qs_ks_types, ONLY: qs_ks_env_type,& + set_ks_env USE qs_mo_types, ONLY: allocate_mo_set,& deallocate_mo_set,& get_mo_set,& @@ -117,7 +120,9 @@ CONTAINS TYPE(cp_fm_type), POINTER :: mos_virt_spin TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(dft_control_type), POINTER :: dft_control + TYPE(excited_energy_type), POINTER :: ex_env TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(qs_ks_env_type), POINTER :: ks_env TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(section_vals_type), POINTER :: print_section TYPE(tddfpt2_control_type), POINTER :: tddfpt_control @@ -127,6 +132,12 @@ 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 + 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) + CALL set_ks_env(ks_env, matrix_ks=matrix_ks) + END IF CPASSERT(.NOT. ASSOCIATED(gs_mos)) ! obtain occupied and virtual (unoccupied) ground-state Kohn-Sham orbitals @@ -918,12 +929,13 @@ CONTAINS ! ************************************************************************************************** !> \brief Generate missed guess vectors. -!> \param tddfpt_control ... -!> \param nspins ... !> \param evects guess vectors distributed across all processors (initialised on exit) !> \param evals guessed transition energies (initialised on exit) !> \param gs_mos occupied and virtual molecular orbitals optimised for the ground state !> \param log_unit output unit +!> \param tddfpt_control ... +!> \param qs_env ... +!> \param nspins ... !> \par History !> * 05.2016 created as tddfpt_guess() [Sergey Chulkov] !> * 06.2016 renamed, altered prototype, supports spin-polarised density [Sergey Chulkov] @@ -938,28 +950,28 @@ CONTAINS !> a restart file). !> \endparblock ! ************************************************************************************************** - SUBROUTINE tddfpt_guess_vectors(tddfpt_control, nspins, evects, evals, gs_mos, log_unit) - TYPE(tddfpt2_control_type), POINTER :: tddfpt_control - INTEGER :: nspins + SUBROUTINE tddfpt_guess_vectors(evects, evals, gs_mos, log_unit, tddfpt_control, qs_env, nspins) TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals + INTEGER, INTENT(in) :: nspins + TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env + TYPE(tddfpt2_control_type), INTENT(in), POINTER :: tddfpt_control + INTEGER, INTENT(in) :: log_unit TYPE(tddfpt_ground_state_mos), DIMENSION(nspins), & INTENT(in) :: gs_mos - INTEGER, INTENT(in) :: log_unit CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_guess_vectors' CHARACTER(len=5) :: spin_label1, spin_label2 - INTEGER :: handle, imo_occ, imo_virt, ind, ispin, & - istate, jspin, nstates, & - nstates_occ_virt_alpha, & - nstates_selected, spin1, spin2 + INTEGER :: handle, imo_occ, imo_virt, ind, ispin, istate, jspin, nstates, & + nstates_occ_virt_alpha, nstates_selected, spin1, spin2, total INTEGER, ALLOCATABLE, DIMENSION(:) :: inds INTEGER, DIMENSION(maxspins) :: nmo_occ_avail, nmo_occ_selected, & nmo_virt_selected REAL(kind=dp) :: e_occ - REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ + REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_virt_minus_occ, gw_occ, gw_virt TYPE(cp_fm_struct_p_type), DIMENSION(maxspins) :: fm_struct_evects + TYPE(excited_energy_type), POINTER :: ex_env CALL timeset(routineN, handle) @@ -970,6 +982,9 @@ CONTAINS CPASSERT(nspins == 1 .OR. nspins == 2) END IF + NULLIFY (ex_env) + CALL get_qs_env(qs_env, exstate_env=ex_env) + DO ispin = 1, nspins ! number of occupied orbitals for each spin component nmo_occ_avail(ispin) = SIZE(gs_mos(ispin)%evals_occ) @@ -996,15 +1011,35 @@ CONTAINS IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN ! Guess for spin-conserving TDDFT DO ispin = 1, nspins + ! if do_bse and do_gw, take gw zeroth order + IF (tddfpt_control%do_bse) THEN + ALLOCATE (gw_virt(nmo_virt_selected(ispin)), gw_occ(nmo_occ_selected(ispin))) + total = nmo_virt_selected(ispin) + nmo_occ_selected(ispin) + gw_virt(1:nmo_virt_selected(ispin)) = ex_env%gw_eigen((nmo_occ_selected(ispin) + 1):total) + gw_occ(1:nmo_occ_selected(ispin)) = ex_env%gw_eigen(1:nmo_occ_selected(ispin)) + END IF + DO imo_occ = 1, nmo_occ_selected(ispin) ! Here imo_occ enumerate Occupied orbitals in inverse order (from the last to the first element) - e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1) + IF (.NOT. tddfpt_control%do_bse) THEN + e_occ = gs_mos(ispin)%evals_occ(nmo_occ_avail(ispin) - imo_occ + 1) + ELSE + e_occ = gw_occ(nmo_occ_avail(ispin) - imo_occ + 1) + END IF DO imo_virt = 1, nmo_virt_selected(ispin) istate = istate + 1 - e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ + IF (.NOT. tddfpt_control%do_bse) THEN + e_virt_minus_occ(istate) = gs_mos(ispin)%evals_virt(imo_virt) - e_occ + ELSE + e_virt_minus_occ(istate) = gw_virt(imo_virt) - e_occ + END IF END DO END DO + + IF (tddfpt_control%do_bse) THEN + DEALLOCATE (gw_virt, gw_occ) + END IF END DO ELSE ! Guess for spin-flip TDDFT diff --git a/tests/QS/regtest-tddfpt-bse/H2O.inp b/tests/QS/regtest-tddfpt-bse/H2O.inp new file mode 100644 index 0000000000..daaad581fc --- /dev/null +++ b/tests/QS/regtest-tddfpt-bse/H2O.inp @@ -0,0 +1,94 @@ +&GLOBAL + PRINT_LEVEL MEDIUM + PROJECT BSE_H2O + RUN_TYPE ENERGY + &TIMINGS + THRESHOLD 0.01 + &END TIMINGS +&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 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 + 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-5 + DO_BSE T + 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/H2O_gas.xyz b/tests/QS/regtest-tddfpt-bse/H2O_gas.xyz new file mode 100644 index 0000000000..d1b6342153 --- /dev/null +++ b/tests/QS/regtest-tddfpt-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-tddfpt-bse/TEST_FILES.toml b/tests/QS/regtest-tddfpt-bse/TEST_FILES.toml new file mode 100644 index 0000000000..d686bb0363 --- /dev/null +++ b/tests/QS/regtest-tddfpt-bse/TEST_FILES.toml @@ -0,0 +1,7 @@ +# runs are executed in the same order as in this file +# the second field tells which test should be run in order to compare with the last available output +# 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=1.0E-05, ref=2.1463}] +#EOF diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index b270ab92d9..8764f08f5a 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -23,6 +23,7 @@ QS/regtest-dlaf-2 dlaf libint !ifx QS/regtest-pao-2 QS/regtest-tddfpt-force-gapw libint libxc !ifx QS/regtest-tddfpt-4 libint libxc !ifx +QS/regtest-tddfpt-bse libint libxc QS/regtest-tddfpt-force-3 libint !ifx QS/regtest-debug-5 libint !ifx QS/regtest-gpw-4