Sternheimer BSE 1.Version parallel (#4832)

This commit is contained in:
annahehn 2026-03-30 10:15:54 +02:00 committed by GitHub
parent b4df9bbab3
commit ed08826bad
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
15 changed files with 423 additions and 60 deletions

View file

@ -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(:)

View file

@ -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)

View file

@ -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

View file

@ -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)

View file

@ -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", &

View file

@ -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)

View file

@ -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")

View file

@ -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

View file

@ -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

View file

@ -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)

View file

@ -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

View file

@ -2,9 +2,6 @@
PRINT_LEVEL MEDIUM
PROJECT BSE_H2O
RUN_TYPE ENERGY
&TIMINGS
THRESHOLD 0.01
&END TIMINGS
&END GLOBAL
&FORCE_EVAL

View file

@ -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

View file

@ -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

View file

@ -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