comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-27 21:55:16 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Sternheimer BSE (1. Version with W matrix given)

Browse source
This commit is contained in:
Anna Hehn 2025-07-21 14:57:46 +02:00 • committed by Ole Schütt
parent 67ae283c5d
commit fdff660ab3
18 changed files with 688 additions and 90 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

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

32
src/bse_full_diag.F
View file

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

4
src/bse_main.F
View file

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

2
src/cp_control_types.F
View file

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

1
src/cp_control_utils.F
View file

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

12
src/exstates_types.F
View file

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

4
src/gw_utils.F
View file

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

6
src/input_cp2k_properties_dft.F
View file

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

20
src/mp2.F
View file

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

375
src/qs_tddfpt2_bse_utils.F Normal file
View file

@ -0,0 +1,375 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! 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

66
src/qs_tddfpt2_eigensolver.F
View file

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

82
src/qs_tddfpt2_fhxc.F
View file

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

5
src/qs_tddfpt2_methods.F
View file

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

61
src/qs_tddfpt2_utils.F
View file

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

94
tests/QS/regtest-tddfpt-bse/H2O.inp Normal file
View file

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

5
tests/QS/regtest-tddfpt-bse/H2O_gas.xyz Normal file
View file

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

7
tests/QS/regtest-tddfpt-bse/TEST_FILES.toml Normal file
View file

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

1
tests/TEST_DIRS
View file

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