mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 22:25:32 -04:00
Bethe Salpeter equation for molecules (#3308)
This commit is contained in:
Max Graml
GitHub
parent
cdee4f21b0
commit
6a30819246
20 changed files with 4398 additions and 773 deletions
|
|
@ -69,7 +69,10 @@ list(
|
|||
basis_set_output.F
|
||||
beta_gamma_psi.F
|
||||
block_p_types.F
|
||||
bse.F
|
||||
bse_main.F
|
||||
bse_full_diag.F
|
||||
bse_iterative.F
|
||||
bse_util.F
|
||||
bsse.F
|
||||
cell_methods.F
|
||||
colvar_methods.F
|
||||
|
|
|
|||
646
src/bse.F
646
src/bse.F
|
|
@ -1,646 +0,0 @@
|
|||
!--------------------------------------------------------------------------------------------------!
|
||||
! CP2K: A general program to perform molecular dynamics simulations !
|
||||
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
|
||||
! !
|
||||
! SPDX-License-Identifier: GPL-2.0-or-later !
|
||||
!--------------------------------------------------------------------------------------------------!
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Routines for GW + Bethe-Salpeter for computing electronic excitations
|
||||
!> \par History
|
||||
!> 04.2017 created [Jan Wilhelm]
|
||||
! **************************************************************************************************
|
||||
MODULE bse
|
||||
USE cp_fm_basic_linalg, ONLY: cp_fm_upper_to_full
|
||||
USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
|
||||
cp_fm_cholesky_invert
|
||||
USE cp_fm_types, ONLY: cp_fm_create,&
|
||||
cp_fm_get_info,&
|
||||
cp_fm_release,&
|
||||
cp_fm_set_all,&
|
||||
cp_fm_to_fm,&
|
||||
cp_fm_type
|
||||
USE group_dist_types, ONLY: get_group_dist,&
|
||||
group_dist_d1_type
|
||||
USE kinds, ONLY: dp
|
||||
USE message_passing, ONLY: mp_para_env_type,&
|
||||
mp_request_type
|
||||
USE mp2_types, ONLY: integ_mat_buffer_type
|
||||
USE parallel_gemm_api, ONLY: parallel_gemm
|
||||
USE rpa_communication, ONLY: communicate_buffer
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
PRIVATE
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse'
|
||||
|
||||
PUBLIC :: mult_B_with_W_and_fill_local_3c_arrays, do_subspace_iterations
|
||||
|
||||
CONTAINS
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param B_bar_ijQ_bse_local ...
|
||||
!> \param B_abQ_bse_local ...
|
||||
!> \param B_bar_iaQ_bse_local ...
|
||||
!> \param B_iaQ_bse_local ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param num_Z_vectors ...
|
||||
!> \param max_iter ...
|
||||
!> \param threshold_min_trans ...
|
||||
!> \param Eigenval ...
|
||||
!> \param para_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_iaQ_bse_local, homo, virtual, num_Z_vectors, &
|
||||
max_iter, threshold_min_trans, Eigenval, para_env)
|
||||
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_bar_ijQ_bse_local, B_abQ_bse_local, &
|
||||
B_bar_iaQ_bse_local, B_iaQ_bse_local
|
||||
INTEGER :: homo, virtual, num_Z_vectors, max_iter
|
||||
REAL(KIND=dp) :: threshold_min_trans
|
||||
REAL(KIND=dp), DIMENSION(:) :: Eigenval
|
||||
TYPE(mp_para_env_type), INTENT(IN) :: para_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'do_subspace_iterations'
|
||||
|
||||
INTEGER :: handle, i_iter, local_RI_size
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: M_ia_tmp, M_ji_tmp, RI_vector
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, BZ, Z_vectors
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
! JW hack 2del
|
||||
threshold_min_trans = 0.01_dp
|
||||
|
||||
ALLOCATE (Z_vectors(homo, virtual, num_Z_vectors))
|
||||
Z_vectors = 0.0_dp
|
||||
|
||||
ALLOCATE (AZ(homo, virtual, num_Z_vectors))
|
||||
AZ = 0.0_dp
|
||||
|
||||
ALLOCATE (BZ(homo, virtual, num_Z_vectors))
|
||||
BZ = 0.0_dp
|
||||
|
||||
local_RI_size = SIZE(B_iaQ_bse_local, 3)
|
||||
|
||||
ALLOCATE (M_ia_tmp(homo, virtual))
|
||||
M_ia_tmp = 0.0_dp
|
||||
|
||||
ALLOCATE (M_ji_tmp(homo, homo))
|
||||
M_ji_tmp = 0.0_dp
|
||||
|
||||
ALLOCATE (RI_vector(local_RI_size, num_Z_vectors))
|
||||
RI_vector = 0.0_dp
|
||||
|
||||
CALL initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual)
|
||||
|
||||
DO i_iter = 1, max_iter
|
||||
|
||||
CALL compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, &
|
||||
M_ia_tmp, RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, &
|
||||
para_env)
|
||||
|
||||
CALL compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
M_ji_tmp, RI_vector, homo, virtual, num_Z_vectors, local_RI_size, &
|
||||
para_env)
|
||||
|
||||
END DO
|
||||
|
||||
DEALLOCATE (AZ, BZ, Z_vectors, M_ia_tmp, M_ji_tmp, RI_vector)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param BZ ...
|
||||
!> \param Z_vectors ...
|
||||
!> \param B_iaQ_bse_local ...
|
||||
!> \param B_bar_iaQ_bse_local ...
|
||||
!> \param M_ji_tmp ...
|
||||
!> \param RI_vector ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param num_Z_vectors ...
|
||||
!> \param local_RI_size ...
|
||||
!> \param para_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE compute_BZ(BZ, Z_vectors, B_iaQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
M_ji_tmp, RI_vector, homo, virtual, num_Z_vectors, local_RI_size, para_env)
|
||||
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BZ, Z_vectors, B_iaQ_bse_local, &
|
||||
B_bar_iaQ_bse_local
|
||||
REAL(KIND=dp), DIMENSION(:, :) :: M_ji_tmp, RI_vector
|
||||
INTEGER :: homo, virtual, num_Z_vectors, &
|
||||
local_RI_size
|
||||
TYPE(mp_para_env_type), INTENT(IN) :: para_env
|
||||
|
||||
INTEGER :: i_Z_vector, LLL
|
||||
|
||||
BZ(:, :, :) = 0.0_dp
|
||||
|
||||
CALL compute_v_ia_jb_part(BZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
|
||||
num_Z_vectors, homo, virtual)
|
||||
|
||||
DO i_Z_vector = 1, num_Z_vectors
|
||||
|
||||
DO LLL = 1, local_RI_size
|
||||
|
||||
! M_ji^P = sum_b Z_jb*B_bi^P
|
||||
CALL DGEMM("N", "T", homo, homo, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, &
|
||||
B_iaQ_bse_local(:, :, LLL), homo, 0.0_dp, M_ji_tmp, homo)
|
||||
! (BZ)_ia = sum_jP M_ij^P*B^bar_ja^P
|
||||
CALL DGEMM("T", "N", homo, virtual, homo, 1.0_dp, M_ji_tmp, homo, &
|
||||
B_bar_iaQ_bse_local, homo, 1.0_dp, BZ(:, :, i_Z_vector), homo)
|
||||
|
||||
END DO
|
||||
|
||||
END DO
|
||||
|
||||
! we make the sum over all RI basis functions
|
||||
CALL para_env%sum(BZ)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param AZ ...
|
||||
!> \param Z_vectors ...
|
||||
!> \param B_iaQ_bse_local ...
|
||||
!> \param B_bar_ijQ_bse_local ...
|
||||
!> \param B_abQ_bse_local ...
|
||||
!> \param M_ia_tmp ...
|
||||
!> \param RI_vector ...
|
||||
!> \param Eigenval ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param num_Z_vectors ...
|
||||
!> \param local_RI_size ...
|
||||
!> \param para_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE compute_AZ(AZ, Z_vectors, B_iaQ_bse_local, B_bar_ijQ_bse_local, B_abQ_bse_local, M_ia_tmp, &
|
||||
RI_vector, Eigenval, homo, virtual, num_Z_vectors, local_RI_size, para_env)
|
||||
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: AZ, Z_vectors, B_iaQ_bse_local, &
|
||||
B_bar_ijQ_bse_local, B_abQ_bse_local
|
||||
REAL(KIND=dp), DIMENSION(:, :) :: M_ia_tmp, RI_vector
|
||||
REAL(KIND=dp), DIMENSION(:) :: Eigenval
|
||||
INTEGER :: homo, virtual, num_Z_vectors, &
|
||||
local_RI_size
|
||||
TYPE(mp_para_env_type), INTENT(IN) :: para_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_AZ'
|
||||
|
||||
INTEGER :: a_virt, handle, i_occ, i_Z_vector, LLL
|
||||
REAL(KIND=dp) :: eigen_diff
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
AZ(:, :, :) = 0.0_dp
|
||||
|
||||
CALL compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
|
||||
num_Z_vectors, homo, virtual)
|
||||
|
||||
DO i_Z_vector = 1, num_Z_vectors
|
||||
|
||||
! JW TO DO: OMP PARALLELIZATION
|
||||
DO LLL = 1, local_RI_size
|
||||
|
||||
! M_ja^P = sum_j Z_jb*B_ba^P
|
||||
CALL DGEMM("N", "N", homo, virtual, virtual, 1.0_dp, Z_vectors(:, :, i_Z_vector), homo, &
|
||||
B_abQ_bse_local(:, :, LLL), virtual, 0.0_dp, M_ia_tmp, homo)
|
||||
|
||||
! (AZ)_ia = sum_jP B_bar_ij^P*M_ja^P
|
||||
CALL DGEMM("N", "N", homo, virtual, homo, 1.0_dp, B_bar_ijQ_bse_local(:, :, LLL), homo, &
|
||||
M_ia_tmp, homo, 1.0_dp, AZ(:, :, i_Z_vector), homo)
|
||||
|
||||
END DO
|
||||
|
||||
END DO
|
||||
|
||||
! we make the sum over all RI basis functions
|
||||
CALL para_env%sum(AZ)
|
||||
|
||||
! add (e_a-e_i)*Z_ia
|
||||
DO i_occ = 1, homo
|
||||
DO a_virt = 1, virtual
|
||||
|
||||
eigen_diff = Eigenval(a_virt + homo) - Eigenval(i_occ)
|
||||
|
||||
AZ(i_occ, a_virt, :) = AZ(i_occ, a_virt, :) + Z_vectors(i_occ, a_virt, :)*eigen_diff
|
||||
|
||||
END DO
|
||||
END DO
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param AZ ...
|
||||
!> \param Z_vectors ...
|
||||
!> \param B_iaQ_bse_local ...
|
||||
!> \param RI_vector ...
|
||||
!> \param local_RI_size ...
|
||||
!> \param num_Z_vectors ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE compute_v_ia_jb_part(AZ, Z_vectors, B_iaQ_bse_local, RI_vector, local_RI_size, &
|
||||
num_Z_vectors, homo, virtual)
|
||||
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: AZ, Z_vectors, B_iaQ_bse_local
|
||||
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: RI_vector
|
||||
INTEGER, INTENT(IN) :: local_RI_size, num_Z_vectors, homo, &
|
||||
virtual
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_v_ia_jb_part'
|
||||
|
||||
INTEGER :: a_virt, handle, i_occ, i_Z_vector, LLL
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
RI_vector = 0.0_dp
|
||||
|
||||
! v_P = sum_jb B_jb^P Z_jb
|
||||
DO LLL = 1, local_RI_size
|
||||
DO i_Z_vector = 1, num_Z_vectors
|
||||
DO i_occ = 1, homo
|
||||
DO a_virt = 1, virtual
|
||||
|
||||
RI_vector(LLL, i_Z_vector) = RI_vector(LLL, i_Z_vector) + &
|
||||
Z_vectors(i_occ, a_virt, i_Z_vector)* &
|
||||
B_iaQ_bse_local(i_occ, a_virt, LLL)
|
||||
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
|
||||
! AZ = sum_P B_ia^P*v_P + ...
|
||||
DO LLL = 1, local_RI_size
|
||||
DO i_Z_vector = 1, num_Z_vectors
|
||||
DO i_occ = 1, homo
|
||||
DO a_virt = 1, virtual
|
||||
|
||||
AZ(i_occ, a_virt, i_Z_vector) = AZ(i_occ, a_virt, i_Z_vector) + &
|
||||
RI_vector(LLL, i_Z_vector)* &
|
||||
B_iaQ_bse_local(i_occ, a_virt, LLL)
|
||||
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
END DO
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param Z_vectors ...
|
||||
!> \param Eigenval ...
|
||||
!> \param num_Z_vectors ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE initial_guess_Z_vectors(Z_vectors, Eigenval, num_Z_vectors, homo, virtual)
|
||||
|
||||
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: Z_vectors
|
||||
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eigenval
|
||||
INTEGER, INTENT(IN) :: num_Z_vectors, homo, virtual
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'initial_guess_Z_vectors'
|
||||
|
||||
INTEGER :: a_virt, handle, i_occ, i_Z_vector, &
|
||||
min_loc(2)
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigen_diff_ia
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
ALLOCATE (eigen_diff_ia(homo, virtual))
|
||||
|
||||
DO i_occ = 1, homo
|
||||
DO a_virt = 1, virtual
|
||||
eigen_diff_ia(i_occ, a_virt) = Eigenval(a_virt + homo) - Eigenval(i_occ)
|
||||
END DO
|
||||
END DO
|
||||
|
||||
DO i_Z_vector = 1, num_Z_vectors
|
||||
|
||||
min_loc = MINLOC(eigen_diff_ia)
|
||||
|
||||
Z_vectors(min_loc(1), min_loc(2), i_Z_vector) = 1.0_dp
|
||||
|
||||
eigen_diff_ia(min_loc(1), min_loc(2)) = 1.0E20_dp
|
||||
|
||||
END DO
|
||||
|
||||
DEALLOCATE (eigen_diff_ia)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param fm_mat_S_ij_bse ...
|
||||
!> \param fm_mat_S_ab_bse ...
|
||||
!> \param fm_mat_S ...
|
||||
!> \param fm_mat_Q_static_bse ...
|
||||
!> \param fm_mat_Q_static_bse_gemm ...
|
||||
!> \param B_bar_ijQ_bse_local ...
|
||||
!> \param B_abQ_bse_local ...
|
||||
!> \param B_bar_iaQ_bse_local ...
|
||||
!> \param B_iaQ_bse_local ...
|
||||
!> \param dimen_RI ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param dimen_ia ...
|
||||
!> \param gd_array ...
|
||||
!> \param color_sub ...
|
||||
!> \param para_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE mult_B_with_W_and_fill_local_3c_arrays(fm_mat_S_ij_bse, fm_mat_S_ab_bse, fm_mat_S, fm_mat_Q_static_bse, &
|
||||
fm_mat_Q_static_bse_gemm, &
|
||||
B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_iaQ_bse_local, dimen_RI, homo, virtual, dimen_ia, &
|
||||
gd_array, color_sub, para_env)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
|
||||
fm_mat_S, fm_mat_Q_static_bse, &
|
||||
fm_mat_Q_static_bse_gemm
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
|
||||
INTENT(OUT) :: B_bar_ijQ_bse_local, B_abQ_bse_local, &
|
||||
B_bar_iaQ_bse_local, B_iaQ_bse_local
|
||||
INTEGER, INTENT(IN) :: dimen_RI, homo, virtual, dimen_ia
|
||||
TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
|
||||
INTEGER, INTENT(IN) :: color_sub
|
||||
TYPE(mp_para_env_type), INTENT(IN) :: para_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'mult_B_with_W_and_fill_local_3c_arrays'
|
||||
|
||||
INTEGER :: handle, i_global, iiB, info_chol, &
|
||||
j_global, jjB, ncol_local, nrow_local
|
||||
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
|
||||
TYPE(cp_fm_type) :: fm_mat_S_bar_ia_bse, &
|
||||
fm_mat_S_bar_ij_bse, fm_mat_work
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
CALL cp_fm_create(fm_mat_S_bar_ia_bse, fm_mat_S%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_mat_S, fm_mat_S_bar_ia_bse)
|
||||
CALL cp_fm_set_all(fm_mat_S_bar_ia_bse, 0.0_dp)
|
||||
|
||||
CALL cp_fm_create(fm_mat_S_bar_ij_bse, fm_mat_S_ij_bse%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_mat_S_ij_bse, fm_mat_S_bar_ij_bse)
|
||||
CALL cp_fm_set_all(fm_mat_S_bar_ij_bse, 0.0_dp)
|
||||
|
||||
CALL cp_fm_create(fm_mat_work, fm_mat_Q_static_bse_gemm%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_mat_Q_static_bse_gemm, fm_mat_work)
|
||||
CALL cp_fm_set_all(fm_mat_work, 0.0_dp)
|
||||
|
||||
! get info of fm_mat_Q_static_bse and compute ((1+Q(0))^-1-1)
|
||||
CALL cp_fm_get_info(matrix=fm_mat_Q_static_bse_gemm, &
|
||||
nrow_local=nrow_local, &
|
||||
ncol_local=ncol_local, &
|
||||
row_indices=row_indices, &
|
||||
col_indices=col_indices)
|
||||
|
||||
DO jjB = 1, ncol_local
|
||||
j_global = col_indices(jjB)
|
||||
DO iiB = 1, nrow_local
|
||||
i_global = row_indices(iiB)
|
||||
IF (j_global == i_global .AND. i_global <= dimen_RI) THEN
|
||||
fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) = fm_mat_Q_static_bse_gemm%local_data(iiB, jjB) + 1.0_dp
|
||||
END IF
|
||||
END DO
|
||||
END DO
|
||||
|
||||
! calculate Trace(Log(Matrix)) as Log(DET(Matrix)) via cholesky decomposition
|
||||
CALL cp_fm_cholesky_decompose(matrix=fm_mat_Q_static_bse_gemm, n=dimen_RI, info_out=info_chol)
|
||||
CPASSERT(info_chol == 0)
|
||||
|
||||
! calculate [1+Q(i0)]^-1
|
||||
CALL cp_fm_cholesky_invert(fm_mat_Q_static_bse_gemm)
|
||||
! symmetrize the result
|
||||
CALL cp_fm_upper_to_full(fm_mat_Q_static_bse_gemm, fm_mat_work)
|
||||
|
||||
CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=homo**2, k=dimen_RI, alpha=1.0_dp, &
|
||||
matrix_a=fm_mat_Q_static_bse, matrix_b=fm_mat_S_ij_bse, beta=0.0_dp, &
|
||||
matrix_c=fm_mat_S_bar_ij_bse)
|
||||
|
||||
! fm_mat_S_bar_ia_bse has a different blacs_env as fm_mat_S_ij_bse since we take
|
||||
! fm_mat_S from RPA. Therefore, we also need a different fm_mat_Q_static_bse_gemm
|
||||
CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=dimen_ia, k=dimen_RI, alpha=1.0_dp, &
|
||||
matrix_a=fm_mat_Q_static_bse_gemm, matrix_b=fm_mat_S, beta=0.0_dp, &
|
||||
matrix_c=fm_mat_S_bar_ia_bse)
|
||||
|
||||
CALL allocate_and_fill_local_array(B_iaQ_bse_local, fm_mat_S, gd_array, color_sub, homo, virtual, dimen_RI, para_env)
|
||||
|
||||
CALL allocate_and_fill_local_array(B_bar_iaQ_bse_local, fm_mat_S_bar_ia_bse, gd_array, color_sub, homo, virtual, &
|
||||
dimen_RI, para_env)
|
||||
|
||||
CALL allocate_and_fill_local_array(B_bar_ijQ_bse_local, fm_mat_S_bar_ij_bse, gd_array, color_sub, homo, homo, &
|
||||
dimen_RI, para_env)
|
||||
|
||||
CALL allocate_and_fill_local_array(B_abQ_bse_local, fm_mat_S_ab_bse, gd_array, color_sub, virtual, virtual, &
|
||||
dimen_RI, para_env)
|
||||
|
||||
CALL cp_fm_release(fm_mat_S_bar_ia_bse)
|
||||
CALL cp_fm_release(fm_mat_S_bar_ij_bse)
|
||||
CALL cp_fm_release(fm_mat_work)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param B_local ...
|
||||
!> \param fm_mat_S ...
|
||||
!> \param gd_array ...
|
||||
!> \param color_sub ...
|
||||
!> \param small_size ...
|
||||
!> \param big_size ...
|
||||
!> \param dimen_RI ...
|
||||
!> \param para_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE allocate_and_fill_local_array(B_local, fm_mat_S, gd_array, &
|
||||
color_sub, small_size, big_size, dimen_RI, para_env)
|
||||
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
|
||||
INTENT(OUT) :: B_local
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
|
||||
TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
|
||||
INTEGER, INTENT(IN) :: color_sub, small_size, big_size, dimen_RI
|
||||
TYPE(mp_para_env_type), INTENT(IN) :: para_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_and_fill_local_array'
|
||||
|
||||
INTEGER :: combi_index, end_RI, handle, handle1, i_comm, i_entry, iiB, imepos, jjB, &
|
||||
level_big_size, level_small_size, ncol_local, nrow_local, num_comm_cycles, RI_index, &
|
||||
size_RI, start_RI
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, mepos_from_RI_index, &
|
||||
num_entries_rec, num_entries_send
|
||||
INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
|
||||
REAL(KIND=dp) :: matrix_el
|
||||
TYPE(integ_mat_buffer_type), ALLOCATABLE, &
|
||||
DIMENSION(:) :: buffer_rec, buffer_send
|
||||
TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
ALLOCATE (mepos_from_RI_index(dimen_RI))
|
||||
mepos_from_RI_index = 0
|
||||
|
||||
DO imepos = 0, para_env%num_pe - 1
|
||||
|
||||
CALL get_group_dist(gd_array, pos=imepos, starts=start_RI, ends=end_RI)
|
||||
|
||||
mepos_from_RI_index(start_RI:end_RI) = imepos
|
||||
|
||||
END DO
|
||||
|
||||
! color_sub is automatically the number of the process since every subgroup has only one MPI rank
|
||||
CALL get_group_dist(gd_array, color_sub, start_RI, end_RI, size_RI)
|
||||
|
||||
ALLOCATE (B_local(small_size, big_size, 1:size_RI))
|
||||
|
||||
ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
|
||||
ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
|
||||
|
||||
ALLOCATE (req_array(1:para_env%num_pe, 4))
|
||||
|
||||
ALLOCATE (entry_counter(0:para_env%num_pe - 1))
|
||||
|
||||
CALL cp_fm_get_info(matrix=fm_mat_S, &
|
||||
nrow_local=nrow_local, &
|
||||
ncol_local=ncol_local, &
|
||||
row_indices=row_indices, &
|
||||
col_indices=col_indices)
|
||||
|
||||
num_comm_cycles = 10
|
||||
|
||||
! communicate not all due to huge memory overhead, since for every number in fm_mat_S, we store
|
||||
! three additional ones (RI index, first MO index, second MO index!!)
|
||||
DO i_comm = 0, num_comm_cycles - 1
|
||||
|
||||
num_entries_send = 0
|
||||
num_entries_rec = 0
|
||||
|
||||
! loop over RI index to get the number of sent entries
|
||||
DO jjB = 1, nrow_local
|
||||
|
||||
RI_index = row_indices(jjB)
|
||||
|
||||
IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE
|
||||
|
||||
imepos = mepos_from_RI_index(RI_index)
|
||||
|
||||
num_entries_send(imepos) = num_entries_send(imepos) + ncol_local
|
||||
|
||||
END DO
|
||||
|
||||
CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
|
||||
|
||||
ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
|
||||
ALLOCATE (buffer_send(0:para_env%num_pe - 1))
|
||||
|
||||
! allocate data message and corresponding indices
|
||||
DO imepos = 0, para_env%num_pe - 1
|
||||
|
||||
ALLOCATE (buffer_rec(imepos)%msg(num_entries_rec(imepos)))
|
||||
buffer_rec(imepos)%msg = 0.0_dp
|
||||
|
||||
ALLOCATE (buffer_send(imepos)%msg(num_entries_send(imepos)))
|
||||
buffer_send(imepos)%msg = 0.0_dp
|
||||
|
||||
ALLOCATE (buffer_rec(imepos)%indx(num_entries_rec(imepos), 3))
|
||||
buffer_rec(imepos)%indx = 0
|
||||
|
||||
ALLOCATE (buffer_send(imepos)%indx(num_entries_send(imepos), 3))
|
||||
buffer_send(imepos)%indx = 0
|
||||
|
||||
END DO
|
||||
|
||||
entry_counter(:) = 0
|
||||
|
||||
! loop over RI index for filling the send-buffer
|
||||
DO jjB = 1, nrow_local
|
||||
|
||||
RI_index = row_indices(jjB)
|
||||
|
||||
IF (MODULO(RI_index, num_comm_cycles) /= i_comm) CYCLE
|
||||
|
||||
imepos = mepos_from_RI_index(RI_index)
|
||||
|
||||
DO iiB = 1, ncol_local
|
||||
|
||||
combi_index = col_indices(iiB)
|
||||
level_small_size = MAX(1, combi_index - 1)/big_size + 1
|
||||
level_big_size = combi_index - (level_small_size - 1)*big_size
|
||||
|
||||
entry_counter(imepos) = entry_counter(imepos) + 1
|
||||
|
||||
buffer_send(imepos)%msg(entry_counter(imepos)) = fm_mat_S%local_data(jjB, iiB)
|
||||
|
||||
buffer_send(imepos)%indx(entry_counter(imepos), 1) = RI_index
|
||||
buffer_send(imepos)%indx(entry_counter(imepos), 2) = level_small_size
|
||||
buffer_send(imepos)%indx(entry_counter(imepos), 3) = level_big_size
|
||||
|
||||
END DO
|
||||
|
||||
END DO
|
||||
|
||||
CALL timeset("BSE_comm_data", handle1)
|
||||
|
||||
CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, buffer_send, req_array)
|
||||
|
||||
CALL timestop(handle1)
|
||||
|
||||
! fill B_local
|
||||
DO imepos = 0, para_env%num_pe - 1
|
||||
|
||||
DO i_entry = 1, num_entries_rec(imepos)
|
||||
|
||||
RI_index = buffer_rec(imepos)%indx(i_entry, 1) - start_RI + 1
|
||||
level_small_size = buffer_rec(imepos)%indx(i_entry, 2)
|
||||
level_big_size = buffer_rec(imepos)%indx(i_entry, 3)
|
||||
|
||||
matrix_el = buffer_rec(imepos)%msg(i_entry)
|
||||
|
||||
B_local(level_small_size, level_big_size, RI_index) = matrix_el
|
||||
|
||||
END DO
|
||||
|
||||
END DO
|
||||
|
||||
DO imepos = 0, para_env%num_pe - 1
|
||||
DEALLOCATE (buffer_send(imepos)%msg)
|
||||
DEALLOCATE (buffer_send(imepos)%indx)
|
||||
DEALLOCATE (buffer_rec(imepos)%msg)
|
||||
DEALLOCATE (buffer_rec(imepos)%indx)
|
||||
END DO
|
||||
|
||||
DEALLOCATE (buffer_rec, buffer_send)
|
||||
|
||||
END DO
|
||||
|
||||
DEALLOCATE (num_entries_send, num_entries_rec)
|
||||
|
||||
DEALLOCATE (mepos_from_RI_index)
|
||||
|
||||
DEALLOCATE (entry_counter, req_array)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
END MODULE bse
|
||||
712
src/bse_full_diag.F
Normal file
712
src/bse_full_diag.F
Normal file
|
|
@ -0,0 +1,712 @@
|
|||
!--------------------------------------------------------------------------------------------------!
|
||||
! CP2K: A general program to perform molecular dynamics simulations !
|
||||
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
|
||||
! !
|
||||
! SPDX-License-Identifier: GPL-2.0-or-later !
|
||||
!--------------------------------------------------------------------------------------------------!
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Routines for the full diagonalization of GW + Bethe-Salpeter for computing
|
||||
!> electronic excitations
|
||||
!> \par History
|
||||
!> 10.2023 created [Maximilian Graml]
|
||||
! **************************************************************************************************
|
||||
MODULE bse_full_diag
|
||||
|
||||
USE bse_util, ONLY: comp_eigvec_coeff_BSE,&
|
||||
filter_eigvec_contrib,&
|
||||
fm_general_add_bse
|
||||
USE cp_blacs_env, ONLY: cp_blacs_env_create,&
|
||||
cp_blacs_env_release,&
|
||||
cp_blacs_env_type
|
||||
USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
|
||||
USE cp_fm_diag, ONLY: choose_eigv_solver,&
|
||||
cp_fm_power
|
||||
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
|
||||
cp_fm_struct_release,&
|
||||
cp_fm_struct_type
|
||||
USE cp_fm_types, ONLY: cp_fm_create,&
|
||||
cp_fm_get_info,&
|
||||
cp_fm_release,&
|
||||
cp_fm_set_all,&
|
||||
cp_fm_to_fm,&
|
||||
cp_fm_type
|
||||
USE input_constants, ONLY: bse_singlet,&
|
||||
bse_triplet
|
||||
USE kinds, ONLY: dp
|
||||
USE message_passing, ONLY: mp_para_env_type
|
||||
USE mp2_types, ONLY: mp2_type
|
||||
USE parallel_gemm_api, ONLY: parallel_gemm
|
||||
USE physcon, ONLY: evolt
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
PRIVATE
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_full_diag'
|
||||
|
||||
PUBLIC :: create_A, diagonalize_A, create_B, create_hermitian_form_of_ABBA, &
|
||||
diagonalize_C
|
||||
|
||||
CONTAINS
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
|
||||
!> A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
|
||||
!> ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
|
||||
!> α is a spin-dependent factor
|
||||
!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
|
||||
!> W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
|
||||
!> \param fm_mat_S_ia_bse ...
|
||||
!> \param fm_mat_S_bar_ij_bse ...
|
||||
!> \param fm_mat_S_ab_bse ...
|
||||
!> \param fm_A ...
|
||||
!> \param Eigenval ...
|
||||
!> \param unit_nr ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param dimen_RI ...
|
||||
!> \param mp2_env ...
|
||||
!> \param para_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
|
||||
fm_A, Eigenval, unit_nr, &
|
||||
homo, virtual, dimen_RI, mp2_env, &
|
||||
para_env)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, &
|
||||
fm_mat_S_ab_bse
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
|
||||
REAL(KIND=dp), DIMENSION(:) :: Eigenval
|
||||
INTEGER, INTENT(IN) :: unit_nr, homo, virtual, dimen_RI
|
||||
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
|
||||
TYPE(mp_para_env_type), INTENT(INOUT) :: para_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'create_A'
|
||||
|
||||
INTEGER :: a_virt_row, handle, i_occ_row, &
|
||||
i_row_global, ii, j_col_global, jj, &
|
||||
ncol_local_A, nrow_local_A
|
||||
INTEGER, DIMENSION(4) :: reordering
|
||||
INTEGER, DIMENSION(:), POINTER :: col_indices_A, row_indices_A
|
||||
REAL(KIND=dp) :: alpha, eigen_diff
|
||||
TYPE(cp_blacs_env_type), POINTER :: blacs_env
|
||||
TYPE(cp_fm_struct_type), POINTER :: fm_struct_A, fm_struct_W
|
||||
TYPE(cp_fm_type) :: fm_W
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating A'
|
||||
END IF
|
||||
|
||||
!Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
|
||||
SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
|
||||
CASE (bse_singlet)
|
||||
alpha = 2.0_dp
|
||||
CASE (bse_triplet)
|
||||
alpha = 0.0_dp
|
||||
END SELECT
|
||||
|
||||
! create the blacs env for ij matrices (NOT fm_mat_S_ia_bse%matrix_struct related parallel_gemms!)
|
||||
NULLIFY (blacs_env)
|
||||
CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env)
|
||||
|
||||
! We have to use the same blacs_env for A as for the matrices fm_mat_S_ia_bse from RPA
|
||||
! Logic: A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
|
||||
! We create v_ia,jb and W_ij,ab, then we communicate entries from local W_ij,ab
|
||||
! to the full matrix v_ia,jb. By adding these and the energy diffenences: v_ia,jb -> A_ia,jb
|
||||
! We use the A matrix already from the start instead of v
|
||||
CALL cp_fm_struct_create(fm_struct_A, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, &
|
||||
ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env)
|
||||
CALL cp_fm_create(fm_A, fm_struct_A, name="fm_A_iajb")
|
||||
CALL cp_fm_set_all(fm_A, 0.0_dp)
|
||||
|
||||
CALL cp_fm_struct_create(fm_struct_W, context=fm_mat_S_ab_bse%matrix_struct%context, nrow_global=homo**2, &
|
||||
ncol_global=virtual**2, para_env=fm_mat_S_ab_bse%matrix_struct%para_env)
|
||||
CALL cp_fm_create(fm_W, fm_struct_W, name="fm_W_ijab")
|
||||
CALL cp_fm_set_all(fm_W, 0.0_dp)
|
||||
|
||||
! Create A matrix from GW Energies, v_ia,jb and W_ij,ab (different blacs_env!)
|
||||
! v_ia,jb, which is directly initialized in A (with a factor of alpha)
|
||||
! v_ia,jb = \sum_P B^P_ia B^P_jb
|
||||
CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, &
|
||||
matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
|
||||
matrix_c=fm_A)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated A_iajb'
|
||||
END IF
|
||||
|
||||
!W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab
|
||||
CALL parallel_gemm(transa="T", transb="N", m=homo**2, n=virtual**2, k=dimen_RI, alpha=1.0_dp, &
|
||||
matrix_a=fm_mat_S_bar_ij_bse, matrix_b=fm_mat_S_ab_bse, beta=0.0_dp, &
|
||||
matrix_c=fm_W)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
|
||||
END IF
|
||||
|
||||
! We start by moving data from local parts of W_ij,ab to the full matrix A_ia,jb using buffers
|
||||
CALL cp_fm_get_info(matrix=fm_A, &
|
||||
nrow_local=nrow_local_A, &
|
||||
ncol_local=ncol_local_A, &
|
||||
row_indices=row_indices_A, &
|
||||
col_indices=col_indices_A)
|
||||
! Writing -1.0_dp * W_ij,ab to A_ia,jb, i.e. beta = -1.0_dp,
|
||||
! W_ij,ab: nrow_secidx_in = homo, ncol_secidx_in = virtual
|
||||
! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual
|
||||
reordering = (/1, 3, 2, 4/)
|
||||
CALL fm_general_add_bse(fm_A, fm_W, -1.0_dp, homo, virtual, &
|
||||
virtual, virtual, unit_nr, reordering, mp2_env)
|
||||
!full matrix W is not needed anymore, release it to save memory
|
||||
CALL cp_fm_struct_release(fm_struct_W)
|
||||
CALL cp_fm_release(fm_W)
|
||||
CALL cp_fm_struct_release(fm_struct_A)
|
||||
|
||||
!Now add the energy differences (ε_a-ε_i) on the diagonal (i.e. δ_ij δ_ab) of A_ia,jb
|
||||
DO ii = 1, nrow_local_A
|
||||
|
||||
i_row_global = row_indices_A(ii)
|
||||
|
||||
DO jj = 1, ncol_local_A
|
||||
|
||||
j_col_global = col_indices_A(jj)
|
||||
|
||||
IF (i_row_global == j_col_global) THEN
|
||||
i_occ_row = (i_row_global - 1)/virtual + 1
|
||||
a_virt_row = MOD(i_row_global - 1, virtual) + 1
|
||||
eigen_diff = Eigenval(a_virt_row + homo) - Eigenval(i_occ_row)
|
||||
fm_A%local_data(ii, jj) = fm_A%local_data(ii, jj) + eigen_diff
|
||||
|
||||
END IF
|
||||
END DO
|
||||
END DO
|
||||
|
||||
CALL cp_fm_struct_release(fm_struct_A)
|
||||
CALL cp_fm_struct_release(fm_struct_W)
|
||||
|
||||
CALL cp_blacs_env_release(blacs_env)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE create_A
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
|
||||
!> B_ia,jb = α * v_ia,jb - W_ib,aj
|
||||
!> α is a spin-dependent factor
|
||||
!> v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
|
||||
!> W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
|
||||
!> \param fm_mat_S_ia_bse ...
|
||||
!> \param fm_mat_S_bar_ia_bse ...
|
||||
!> \param fm_B ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param dimen_RI ...
|
||||
!> \param unit_nr ...
|
||||
!> \param mp2_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
|
||||
homo, virtual, dimen_RI, unit_nr, mp2_env)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_B
|
||||
INTEGER, INTENT(IN) :: homo, virtual, dimen_RI, unit_nr
|
||||
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'create_B'
|
||||
|
||||
INTEGER :: handle
|
||||
INTEGER, DIMENSION(4) :: reordering
|
||||
REAL(KIND=dp) :: alpha
|
||||
TYPE(cp_fm_struct_type), POINTER :: fm_struct_v
|
||||
TYPE(cp_fm_type) :: fm_W
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A10)') 'BSE|DEBUG|', 'Creating B'
|
||||
END IF
|
||||
|
||||
! Determines factor of exchange term, depending on requested spin configuration (cf. input_constants.F)
|
||||
SELECT CASE (mp2_env%ri_g0w0%bse_spin_config)
|
||||
CASE (bse_singlet)
|
||||
alpha = 2.0_dp
|
||||
CASE (bse_triplet)
|
||||
alpha = 0.0_dp
|
||||
END SELECT
|
||||
|
||||
CALL cp_fm_struct_create(fm_struct_v, context=fm_mat_S_ia_bse%matrix_struct%context, nrow_global=homo*virtual, &
|
||||
ncol_global=homo*virtual, para_env=fm_mat_S_ia_bse%matrix_struct%para_env)
|
||||
CALL cp_fm_create(fm_B, fm_struct_v, name="fm_B_iajb")
|
||||
CALL cp_fm_set_all(fm_B, 0.0_dp)
|
||||
|
||||
CALL cp_fm_create(fm_W, fm_struct_v, name="fm_W_ibaj")
|
||||
CALL cp_fm_set_all(fm_W, 0.0_dp)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated B_iajb'
|
||||
END IF
|
||||
! v_ia,jb = \sum_P B^P_ia B^P_jb
|
||||
CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=alpha, &
|
||||
matrix_a=fm_mat_S_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
|
||||
matrix_c=fm_B)
|
||||
|
||||
! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj
|
||||
CALL parallel_gemm(transa="T", transb="N", m=homo*virtual, n=homo*virtual, k=dimen_RI, alpha=1.0_dp, &
|
||||
matrix_a=fm_mat_S_bar_ia_bse, matrix_b=fm_mat_S_ia_bse, beta=0.0_dp, &
|
||||
matrix_c=fm_W)
|
||||
! from W_ib,ja to A_ia,jb (formally: W_ib,aj, but our internal indexorder is different)
|
||||
! Writing -1.0_dp * W_ib,ja to A_ia,jb, i.e. beta = -1.0_dp,
|
||||
! W_ib,ja: nrow_secidx_in = virtual, ncol_secidx_in = virtual
|
||||
! A_ia,jb: nrow_secidx_out = virtual, ncol_secidx_out = virtual
|
||||
reordering = (/1, 4, 3, 2/)
|
||||
CALL fm_general_add_bse(fm_B, fm_W, -1.0_dp, virtual, virtual, &
|
||||
virtual, virtual, unit_nr, reordering, mp2_env)
|
||||
|
||||
CALL cp_fm_release(fm_W)
|
||||
CALL cp_fm_struct_release(fm_struct_v)
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE create_B
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
|
||||
!> (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
|
||||
!> We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
|
||||
!> of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
|
||||
!> \param fm_A ...
|
||||
!> \param fm_B ...
|
||||
!> \param fm_C ...
|
||||
!> \param fm_sqrt_A_minus_B ...
|
||||
!> \param fm_inv_sqrt_A_minus_B ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param unit_nr ...
|
||||
!> \param mp2_env ...
|
||||
!> \param diag_est ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE create_hermitian_form_of_ABBA(fm_A, fm_B, fm_C, &
|
||||
fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
|
||||
homo, virtual, unit_nr, mp2_env, diag_est)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_A, fm_B
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_C, fm_sqrt_A_minus_B, &
|
||||
fm_inv_sqrt_A_minus_B
|
||||
INTEGER, INTENT(IN) :: homo, virtual, unit_nr
|
||||
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
|
||||
REAL(KIND=dp), INTENT(IN) :: diag_est
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'create_hermitian_form_of_ABBA'
|
||||
|
||||
INTEGER :: dim_mat, handle, n_dependent
|
||||
REAL(KIND=dp), DIMENSION(2) :: eigvals_AB_diff
|
||||
TYPE(cp_fm_type) :: fm_A_minus_B, fm_A_plus_B, fm_dummy, &
|
||||
fm_work_product
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4,T7,A25,A39,ES6.0,A3)') 'BSE|', 'Diagonalizing aux. matrix', &
|
||||
' with size of A. This will take around ', diag_est, " s."
|
||||
END IF
|
||||
|
||||
! Create work matrices, which will hold A+B and A-B and their powers
|
||||
! C is created afterwards to save memory
|
||||
! Final result: C = (A-B)^0.5 (A+B) (A-B)^0.5 EQ.I
|
||||
! \_______/ \___/ \______/
|
||||
! fm_sqrt_A_minus_B fm_A_plus_B fm_sqrt_A_minus_B
|
||||
! (EQ.Ia) (EQ.Ib) (EQ.Ia)
|
||||
! Intermediate work matrices:
|
||||
! fm_inv_sqrt_A_minus_B: (A-B)^-0.5 EQ.II
|
||||
! fm_A_minus_B: (A-B) EQ.III
|
||||
! fm_work_product: (A-B)^0.5 (A+B) from (EQ.Ia) and (EQ.Ib) EQ.IV
|
||||
CALL cp_fm_create(fm_A_plus_B, fm_A%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_A, fm_A_plus_B)
|
||||
CALL cp_fm_create(fm_A_minus_B, fm_A%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_A, fm_A_minus_B)
|
||||
CALL cp_fm_create(fm_sqrt_A_minus_B, fm_A%matrix_struct)
|
||||
CALL cp_fm_set_all(fm_sqrt_A_minus_B, 0.0_dp)
|
||||
CALL cp_fm_create(fm_inv_sqrt_A_minus_B, fm_A%matrix_struct)
|
||||
CALL cp_fm_set_all(fm_inv_sqrt_A_minus_B, 0.0_dp)
|
||||
|
||||
CALL cp_fm_create(fm_work_product, fm_A%matrix_struct)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A19)') 'BSE|DEBUG|', 'Created work arrays'
|
||||
END IF
|
||||
|
||||
! Add/Substract B (cf. EQs. Ib and III)
|
||||
CALL cp_fm_scale_and_add(1.0_dp, fm_A_plus_B, 1.0_dp, fm_B)
|
||||
CALL cp_fm_scale_and_add(1.0_dp, fm_A_minus_B, -1.0_dp, fm_B)
|
||||
|
||||
! cp_fm_power will overwrite matrix, therefore we create copies
|
||||
CALL cp_fm_to_fm(fm_A_minus_B, fm_inv_sqrt_A_minus_B)
|
||||
|
||||
! In order to avoid a second diagonalization (cp_fm_power), we create (A-B)^0.5 (EQ.Ia)
|
||||
! from (A-B)^-0.5 (EQ.II) by multiplication with (A-B) (EQ.III) afterwards.
|
||||
|
||||
! Raise A-B to -0.5_dp, no quenching of eigenvectors, hence threshold=0.0_dp
|
||||
CALL cp_fm_create(fm_dummy, fm_A%matrix_struct)
|
||||
! Create (A-B)^-0.5 (cf. EQ.II)
|
||||
CALL cp_fm_power(fm_inv_sqrt_A_minus_B, fm_dummy, -0.5_dp, 0.0_dp, n_dependent, eigvals=eigvals_AB_diff)
|
||||
CALL cp_fm_release(fm_dummy)
|
||||
! Raise an error in case the the matrix A-B is not positive definite (i.e. negative eigenvalues)
|
||||
! In this case, the procedure for hermitian form of ABBA is not applicable
|
||||
IF (eigvals_AB_diff(1) < 0) THEN
|
||||
CALL cp_abort(__LOCATION__, &
|
||||
"Matrix (A-B) is not positive definite. "// &
|
||||
"Hermitian diagonalization of full BSE matrix is ill-defined.")
|
||||
END IF
|
||||
|
||||
! We keep fm_inv_sqrt_A_minus_B for print of singleparticle transitions of ABBA
|
||||
! We further create (A-B)^0.5 for the singleparticle transitions of ABBA
|
||||
! Create (A-B)^0.5= (A-B)^-0.5 * (A-B) (EQ.Ia)
|
||||
dim_mat = homo*virtual
|
||||
CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_inv_sqrt_A_minus_B, fm_A_minus_B, 0.0_dp, &
|
||||
fm_sqrt_A_minus_B)
|
||||
|
||||
! Compute and store LHS of C, i.e. (A-B)^0.5 (A+B) (EQ.IV)
|
||||
CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_sqrt_A_minus_B, fm_A_plus_B, 0.0_dp, &
|
||||
fm_work_product)
|
||||
|
||||
! Release to save memory
|
||||
CALL cp_fm_release(fm_A_plus_B)
|
||||
CALL cp_fm_release(fm_A_minus_B)
|
||||
|
||||
! Now create full
|
||||
CALL cp_fm_create(fm_C, fm_A%matrix_struct)
|
||||
CALL cp_fm_set_all(fm_C, 0.0_dp)
|
||||
! Compute C=(A-B)^0.5 (A+B) (A-B)^0.5 (EQ.I)
|
||||
CALL parallel_gemm("N", "N", dim_mat, dim_mat, dim_mat, 1.0_dp, fm_work_product, fm_sqrt_A_minus_B, 0.0_dp, &
|
||||
fm_C)
|
||||
CALL cp_fm_release(fm_work_product)
|
||||
|
||||
IF (unit_nr > 0 .AND. mp2_env%ri_g0w0%bse_debug_print) THEN
|
||||
WRITE (unit_nr, '(T2,A10,T13,A36)') 'BSE|DEBUG|', 'Filled C=(A-B)^0.5 (A+B) (A-B)^0.5'
|
||||
END IF
|
||||
|
||||
CALL timestop(handle)
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
|
||||
!> Here, the eigenvectors Z^n relate to X^n via
|
||||
!> Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
|
||||
!> \param fm_C ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param homo_irred ...
|
||||
!> \param fm_sqrt_A_minus_B ...
|
||||
!> \param fm_inv_sqrt_A_minus_B ...
|
||||
!> \param unit_nr ...
|
||||
!> \param diag_est ...
|
||||
!> \param mp2_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, &
|
||||
fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
|
||||
unit_nr, diag_est, mp2_env)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_C
|
||||
INTEGER, INTENT(IN) :: homo, virtual, homo_irred
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B
|
||||
INTEGER, INTENT(IN) :: unit_nr
|
||||
REAL(KIND=dp), INTENT(IN) :: diag_est
|
||||
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_C'
|
||||
|
||||
INTEGER :: diag_info, handle
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
|
||||
TYPE(cp_fm_type) :: fm_eigvec, fm_mat_eigvec_transform, &
|
||||
fm_mat_eigvec_transform_neg
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing C. ', &
|
||||
'This will take around ', diag_est, ' s.'
|
||||
END IF
|
||||
|
||||
!We have now the full matrix C=(A-B)^0.5 (A+B) (A-B)^0.5
|
||||
!Now: Diagonalize it
|
||||
CALL cp_fm_create(fm_eigvec, fm_C%matrix_struct)
|
||||
|
||||
ALLOCATE (Exc_ens(homo*virtual))
|
||||
|
||||
CALL choose_eigv_solver(fm_C, fm_eigvec, Exc_ens, diag_info)
|
||||
CPASSERT(diag_info == 0)
|
||||
Exc_ens = SQRT(Exc_ens)
|
||||
|
||||
! Prepare eigenvector for interpretation of singleparticle transitions
|
||||
! Compare: F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)
|
||||
! We aim for the upper part of the vector (X,Y) for a direct comparison with the TDA result
|
||||
|
||||
! Following Furche, we basically use Eqs. (A10): First, we multiply
|
||||
! the (A-B)^+-0.5 with eigenvectors and then the eigenvalues
|
||||
! One has to be careful about the index structure, since the eigenvector matrix is not symmetric anymore!
|
||||
|
||||
! First, Eq. I from (A10) from Furche: (X+Y)_n = (Ω_n)^0.5 (A-B)^0.5 T_n
|
||||
CALL cp_fm_create(fm_mat_eigvec_transform, fm_C%matrix_struct)
|
||||
CALL cp_fm_set_all(fm_mat_eigvec_transform, 0.0_dp)
|
||||
CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
|
||||
matrix_a=fm_sqrt_A_minus_B, matrix_b=fm_eigvec, beta=0.0_dp, &
|
||||
matrix_c=fm_mat_eigvec_transform)
|
||||
CALL cp_fm_release(fm_sqrt_A_minus_B)
|
||||
CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform, Exc_ens, -0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
|
||||
|
||||
! Second, Eq. II from (A10) from Furche: (X-Y)_n = (Ω_n)^0.5 (A-B)^-0.5 T_n
|
||||
CALL cp_fm_create(fm_mat_eigvec_transform_neg, fm_C%matrix_struct)
|
||||
CALL cp_fm_set_all(fm_mat_eigvec_transform_neg, 0.0_dp)
|
||||
CALL parallel_gemm(transa="N", transb="N", m=homo*virtual, n=homo*virtual, k=homo*virtual, alpha=1.0_dp, &
|
||||
matrix_a=fm_inv_sqrt_A_minus_B, matrix_b=fm_eigvec, beta=0.0_dp, &
|
||||
matrix_c=fm_mat_eigvec_transform_neg)
|
||||
CALL cp_fm_release(fm_inv_sqrt_A_minus_B)
|
||||
CALL cp_fm_release(fm_eigvec)
|
||||
CALL comp_eigvec_coeff_BSE(fm_mat_eigvec_transform_neg, Exc_ens, 0.5_dp, gamma=2.0_dp, do_transpose=.TRUE.)
|
||||
|
||||
! Now, we add the two equations to obtain X_n
|
||||
CALL cp_fm_scale_and_add(1.0_dp, fm_mat_eigvec_transform, 1.0_dp, fm_mat_eigvec_transform_neg)
|
||||
|
||||
!Cleanup
|
||||
CALL cp_fm_release(fm_mat_eigvec_transform_neg)
|
||||
|
||||
CALL success_message(Exc_ens, fm_mat_eigvec_transform, mp2_env, &
|
||||
homo, virtual, homo_irred, unit_nr, .FALSE.)
|
||||
|
||||
DEALLOCATE (Exc_ens)
|
||||
CALL cp_fm_release(fm_mat_eigvec_transform)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE diagonalize_C
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Solving hermitian eigenvalue equation A X^n = Ω^n X^n
|
||||
!> \param fm_A ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param homo_irred ...
|
||||
!> \param unit_nr ...
|
||||
!> \param diag_est ...
|
||||
!> \param mp2_env ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE diagonalize_A(fm_A, homo, virtual, homo_irred, &
|
||||
unit_nr, diag_est, mp2_env)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_A
|
||||
INTEGER, INTENT(IN) :: homo, virtual, homo_irred, unit_nr
|
||||
REAL(KIND=dp), INTENT(IN) :: diag_est
|
||||
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_A'
|
||||
|
||||
INTEGER :: diag_info, handle
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
|
||||
TYPE(cp_fm_type) :: fm_eigvec
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
!Continue with formatting of subroutine create_A
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4,T7,A17,A22,ES6.0,A3)') 'BSE|', 'Diagonalizing A. ', &
|
||||
'This will take around ', diag_est, ' s.'
|
||||
END IF
|
||||
|
||||
!We have now the full matrix A_iajb, distributed over all ranks
|
||||
!Now: Diagonalize it
|
||||
CALL cp_fm_create(fm_eigvec, fm_A%matrix_struct)
|
||||
|
||||
ALLOCATE (Exc_ens(homo*virtual))
|
||||
|
||||
CALL choose_eigv_solver(fm_A, fm_eigvec, Exc_ens, diag_info)
|
||||
CPASSERT(diag_info == 0)
|
||||
CALL success_message(Exc_ens, fm_eigvec, mp2_env, &
|
||||
homo, virtual, homo_irred, unit_nr, .TRUE.)
|
||||
|
||||
CALL cp_fm_release(fm_eigvec)
|
||||
DEALLOCATE (Exc_ens)
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE diagonalize_A
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Prints the success message (incl. energies) for full diag of BSE (TDA/full ABBA via flag)
|
||||
!> \param Exc_ens ...
|
||||
!> \param fm_eigvec ...
|
||||
!> \param mp2_env ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param homo_irred ...
|
||||
!> \param unit_nr ...
|
||||
!> \param flag_TDA ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE success_message(Exc_ens, fm_eigvec, mp2_env, &
|
||||
homo, virtual, homo_irred, unit_nr, flag_TDA)
|
||||
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Exc_ens
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec
|
||||
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
|
||||
INTEGER :: homo, virtual, homo_irred, unit_nr
|
||||
LOGICAL, OPTIONAL :: flag_TDA
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'success_message'
|
||||
|
||||
CHARACTER(LEN=10) :: info_approximation, multiplet
|
||||
INTEGER :: handle, i_exc, k, num_entries
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:) :: idx_homo, idx_virt
|
||||
REAL(KIND=dp) :: alpha
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
!Prepare variables for printing
|
||||
IF (mp2_env%ri_g0w0%bse_spin_config == 0) THEN
|
||||
multiplet = "Singlet"
|
||||
alpha = 2.0_dp
|
||||
ELSE
|
||||
multiplet = "Triplet"
|
||||
alpha = 0.0_dp
|
||||
END IF
|
||||
IF (.NOT. PRESENT(flag_TDA)) THEN
|
||||
flag_TDA = .FALSE.
|
||||
END IF
|
||||
IF (flag_TDA) THEN
|
||||
info_approximation = " -TDA- "
|
||||
ELSE
|
||||
info_approximation = "-full-"
|
||||
END IF
|
||||
|
||||
!Get information about eigenvector
|
||||
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
IF (flag_TDA) THEN
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA) *'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
|
||||
'the BSE within the TDA:'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb X_jb^n ) = Ω^n X_ia^n'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
|
||||
ELSE
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '* Full Bethe Salpeter equation (BSE) (i.e. without TDA) *'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
|
||||
'the BSE without the TDA:'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n| |1 0| |X^n|'
|
||||
WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', ' sum_jb ( A_ia,jb X_jb^n + B_ia,jb Y_jb^n ) = Ω^n X_ia^n'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A62)') 'BSE|', '- sum_jb ( B_ia,jb X_jb^n + A_ia,jb Y_jb^n ) = Ω^n Y_ia^n'
|
||||
END IF
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
|
||||
'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
|
||||
'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab'
|
||||
IF (.NOT. flag_TDA) THEN
|
||||
WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
|
||||
END IF
|
||||
IF (.NOT. flag_TDA) THEN
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', 'The BSE is solved for Ω^n and X_ia^n as a hermitian problem, e.g. Eq.(42)'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
|
||||
! WRITE (unit_nr, '(T2,A4,T7,A69)') 'BSE|', 'C_ia,jb = sum_kc,ld ((A-B)^0.5)_ia,kc (A+B)_kc,ld ((A-B)^0.5)_ld,jb'
|
||||
! WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X+Y)_ia,n = sum_jb,m ((A-B)^0.5)ia,jb Z_jb,m (Ω_m)^-0.5'
|
||||
! WRITE (unit_nr, '(T2,A4,T7,A58)') 'BSE|', '(X-Y)_ia,n = sum_jb,m ((A-B)^-0.5)ia,jb Z_jb,m (Ω_m)^0.5'
|
||||
END IF
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A7,T19,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A7,T19,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A3,T19,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A6,T18,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A6,T18,A27)') 'BSE|', 'W_...:', 'Screened direct interaction'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1))') 'BSE|', &
|
||||
'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
IF (flag_TDA) THEN
|
||||
WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
|
||||
ELSE
|
||||
WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
|
||||
END IF
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T13,A10,T27,A13,T42,A12,T59,A22)') 'BSE|', &
|
||||
'Excitation', "Multiplet", 'TDA/full BSE', 'Excitation energy (eV)'
|
||||
END IF
|
||||
!prints actual energies values
|
||||
IF (unit_nr > 0) THEN
|
||||
DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
|
||||
WRITE (unit_nr, '(T2,A4,T7,I16,T27,A7,A6,T48,A6,T59,F22.4)') &
|
||||
'BSE|', i_exc, multiplet, " State", info_approximation, Exc_ens(i_exc)*evolt
|
||||
END DO
|
||||
END IF
|
||||
!prints single particle transitions
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
|
||||
WRITE (unit_nr, '(T2,A4,T7,A70)') &
|
||||
'BSE|', "Excitations are built up by the following single-particle transitions,"
|
||||
WRITE (unit_nr, '(T2,A4,T7,A42,F5.2,A2)') &
|
||||
'BSE|', "neglecting contributions where |X_ia^n| < ", mp2_env%ri_g0w0%eps_x, " :"
|
||||
|
||||
WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
|
||||
homo_irred, ' and LUMO a =', homo_irred + 1, " --"
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A9,T20,A1,T22,A2,T29,A1,T42,A12,T71,A10)') &
|
||||
"BSE|", "n-th exc.", "i", "=>", "a", 'TDA/full BSE', "|X_ia^n|"
|
||||
END IF
|
||||
DO i_exc = 1, MIN(homo*virtual, mp2_env%ri_g0w0%num_print_exc)
|
||||
!Iterate through eigenvector and print values above threshold
|
||||
CALL filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
|
||||
i_exc, virtual, num_entries, mp2_env)
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
DO k = 1, num_entries
|
||||
WRITE (unit_nr, '(T2,A4,T11,I5,T16,I5,T22,A2,T25,I5,T48,A6,T59,F22.4)') &
|
||||
"BSE|", i_exc, homo_irred - homo + idx_homo(k), "=>", &
|
||||
homo_irred + idx_virt(k), info_approximation, ABS(eigvec_entries(k))
|
||||
END DO
|
||||
END IF
|
||||
|
||||
DEALLOCATE (idx_homo)
|
||||
DEALLOCATE (idx_virt)
|
||||
DEALLOCATE (eigvec_entries)
|
||||
END DO
|
||||
IF (unit_nr > 0) THEN
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4)') 'BSE|'
|
||||
WRITE (unit_nr, '(T2,A4,T7,A53)') 'BSE|', 'The BSE was successfully calculated. Have a nice day!'
|
||||
END IF
|
||||
|
||||
CALL timestop(handle)
|
||||
END SUBROUTINE success_message
|
||||
|
||||
END MODULE bse_full_diag
|
||||
1587
src/bse_iterative.F
Normal file
1587
src/bse_iterative.F
Normal file
File diff suppressed because it is too large
Load diff
221
src/bse_main.F
Normal file
221
src/bse_main.F
Normal file
|
|
@ -0,0 +1,221 @@
|
|||
!--------------------------------------------------------------------------------------------------!
|
||||
! CP2K: A general program to perform molecular dynamics simulations !
|
||||
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
|
||||
! !
|
||||
! SPDX-License-Identifier: GPL-2.0-or-later !
|
||||
!--------------------------------------------------------------------------------------------------!
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Main routines for GW + Bethe-Salpeter for computing electronic excitations
|
||||
!> \par History
|
||||
!> 04.2024 created [Maximilian Graml]
|
||||
! **************************************************************************************************
|
||||
|
||||
MODULE bse_main
|
||||
|
||||
USE bse_full_diag, ONLY: create_A,&
|
||||
create_B,&
|
||||
create_hermitian_form_of_ABBA,&
|
||||
diagonalize_A,&
|
||||
diagonalize_C
|
||||
USE bse_iterative, ONLY: do_subspace_iterations,&
|
||||
fill_local_3c_arrays
|
||||
USE bse_util, ONLY: deallocate_matrices_bse,&
|
||||
estimate_BSE_resources,&
|
||||
mult_B_with_W,&
|
||||
print_BSE_start_flag,&
|
||||
truncate_BSE_matrices
|
||||
USE cp_fm_types, ONLY: cp_fm_release,&
|
||||
cp_fm_type
|
||||
USE group_dist_types, ONLY: group_dist_d1_type
|
||||
USE input_constants, ONLY: bse_abba,&
|
||||
bse_both,&
|
||||
bse_fulldiag,&
|
||||
bse_iterdiag,&
|
||||
bse_tda
|
||||
USE kinds, ONLY: dp
|
||||
USE message_passing, ONLY: mp_para_env_type
|
||||
USE mp2_types, ONLY: mp2_type
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
||||
PRIVATE
|
||||
|
||||
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_main'
|
||||
|
||||
PUBLIC :: start_bse_calculation
|
||||
|
||||
CONTAINS
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief Main subroutine managing BSE calculations
|
||||
!> \param fm_mat_S_ia_bse ...
|
||||
!> \param fm_mat_S_ij_bse ...
|
||||
!> \param fm_mat_S_ab_bse ...
|
||||
!> \param fm_mat_Q_static_bse ...
|
||||
!> \param fm_mat_Q_static_bse_gemm ...
|
||||
!> \param Eigenval ...
|
||||
!> \param homo ...
|
||||
!> \param virtual ...
|
||||
!> \param dimen_RI ...
|
||||
!> \param dimen_RI_red ...
|
||||
!> \param gd_array ...
|
||||
!> \param color_sub ...
|
||||
!> \param mp2_env ...
|
||||
!> \param unit_nr ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
|
||||
fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, &
|
||||
Eigenval, homo, virtual, dimen_RI, dimen_RI_red, &
|
||||
gd_array, color_sub, mp2_env, unit_nr)
|
||||
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S_ia_bse, fm_mat_S_ij_bse, &
|
||||
fm_mat_S_ab_bse
|
||||
TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q_static_bse, &
|
||||
fm_mat_Q_static_bse_gemm
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
|
||||
INTENT(IN) :: Eigenval
|
||||
INTEGER, DIMENSION(:), INTENT(IN) :: homo, virtual
|
||||
INTEGER, INTENT(IN) :: dimen_RI, dimen_RI_red
|
||||
TYPE(group_dist_d1_type), INTENT(IN) :: gd_array
|
||||
INTEGER, INTENT(IN) :: color_sub
|
||||
TYPE(mp2_type) :: mp2_env
|
||||
INTEGER, INTENT(IN) :: unit_nr
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'start_bse_calculation'
|
||||
|
||||
INTEGER :: handle, homo_red, virtual_red
|
||||
LOGICAL :: my_do_abba, my_do_fulldiag, &
|
||||
my_do_iterat_diag, my_do_tda
|
||||
REAL(KIND=dp) :: diag_runtime_est
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: Eigenval_reduced
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_bar_ijQ_bse_local, B_iaQ_bse_local
|
||||
TYPE(cp_fm_type) :: fm_A_BSE, fm_B_BSE, fm_C_BSE, fm_inv_sqrt_A_minus_B, fm_mat_S_ab_trunc, &
|
||||
fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, &
|
||||
fm_sqrt_A_minus_B
|
||||
TYPE(mp_para_env_type), POINTER :: para_env
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
|
||||
para_env => fm_mat_S_ia_bse%matrix_struct%para_env
|
||||
|
||||
my_do_fulldiag = .FALSE.
|
||||
my_do_iterat_diag = .FALSE.
|
||||
my_do_tda = .FALSE.
|
||||
my_do_abba = .FALSE.
|
||||
!Method: Iterative or full diagonalization
|
||||
SELECT CASE (mp2_env%ri_g0w0%bse_diag_method)
|
||||
CASE (bse_iterdiag)
|
||||
my_do_iterat_diag = .TRUE.
|
||||
!MG: Basics of the Davidson solver are implemented, but not rigorously checked.
|
||||
CPABORT("Iterative BSE not yet implemented")
|
||||
CASE (bse_fulldiag)
|
||||
my_do_fulldiag = .TRUE.
|
||||
END SELECT
|
||||
!Approximation: TDA and/or full ABBA matrix
|
||||
SELECT CASE (mp2_env%ri_g0w0%bse_approx)
|
||||
CASE (bse_tda)
|
||||
my_do_tda = .TRUE.
|
||||
CASE (bse_abba)
|
||||
my_do_abba = .TRUE.
|
||||
CASE (bse_both)
|
||||
my_do_tda = .TRUE.
|
||||
my_do_abba = .TRUE.
|
||||
END SELECT
|
||||
|
||||
CALL print_BSE_start_flag(my_do_tda, my_do_abba, unit_nr)
|
||||
|
||||
CALL fm_mat_S_ia_bse%matrix_struct%para_env%sync()
|
||||
|
||||
! Reduce matrices in case of energy cutoff for occupied and unoccupied in A/B-BSE-matrices
|
||||
CALL truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
|
||||
fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
|
||||
Eigenval(:, 1, 1), Eigenval_reduced, &
|
||||
homo(1), virtual(1), dimen_RI, unit_nr, &
|
||||
homo_red, virtual_red, &
|
||||
mp2_env)
|
||||
! \bar{B}^P_rs = \sum_R W_PR B^R_rs where B^R_rs = \sum_T [1/sqrt(v)]_RT (T|rs)
|
||||
! r,s: MO-index, P,R,T: RI-index
|
||||
! B: fm_mat_S_..., W: fm_mat_Q_...
|
||||
CALL mult_B_with_W(fm_mat_S_ij_trunc, fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, &
|
||||
fm_mat_S_bar_ij_bse, fm_mat_Q_static_bse_gemm, &
|
||||
dimen_RI_red, homo_red, virtual_red)
|
||||
|
||||
IF (my_do_iterat_diag) THEN
|
||||
CALL fill_local_3c_arrays(fm_mat_S_ab_trunc, fm_mat_S_ia_trunc, &
|
||||
fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
|
||||
B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_iaQ_bse_local, dimen_RI_red, homo_red, virtual_red, &
|
||||
gd_array, color_sub, para_env)
|
||||
END IF
|
||||
|
||||
IF (my_do_fulldiag) THEN
|
||||
! Quick estimate of memory consumption and runtime of diagonalizations
|
||||
CALL estimate_BSE_resources(homo_red, virtual_red, unit_nr, my_do_abba, &
|
||||
para_env, diag_runtime_est)
|
||||
! Matrix A constructed from GW energies and 3c-B-matrices (cf. subroutine mult_B_with_W)
|
||||
! A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab
|
||||
! ε_a, ε_i are GW singleparticle energies from Eigenval_reduced
|
||||
! α is a spin-dependent factor
|
||||
! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
|
||||
! W_ij,ab = \sum_P \bar{B}^P_ij B^P_ab (screened Coulomb interaction)
|
||||
CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, &
|
||||
fm_A_BSE, Eigenval_reduced, unit_nr, &
|
||||
homo_red, virtual_red, dimen_RI, mp2_env, &
|
||||
para_env)
|
||||
IF (my_do_abba) THEN
|
||||
! Matrix B constructed from 3c-B-matrices (cf. subroutine mult_B_with_W)
|
||||
! B_ia,jb = α * v_ia,jb - W_ib,aj
|
||||
! α is a spin-dependent factor
|
||||
! v_ia,jb = \sum_P B^P_ia B^P_jb (unscreened Coulomb interaction)
|
||||
! W_ib,aj = \sum_P \bar{B}^P_ib B^P_aj (screened Coulomb interaction)
|
||||
CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_bar_ia_bse, fm_B_BSE, &
|
||||
homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
|
||||
! Construct Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 to solve full BSE matrix as a hermitian problem
|
||||
! (cf. Eq. (A7) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001)).
|
||||
! We keep fm_sqrt_A_minus_B and fm_inv_sqrt_A_minus_B for print of singleparticle transitions
|
||||
! of ABBA as described in Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
|
||||
CALL create_hermitian_form_of_ABBA(fm_A_BSE, fm_B_BSE, fm_C_BSE, &
|
||||
fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
|
||||
homo_red, virtual_red, unit_nr, mp2_env, diag_runtime_est)
|
||||
END IF
|
||||
CALL cp_fm_release(fm_B_BSE)
|
||||
IF (my_do_tda) THEN
|
||||
! Solving the hermitian eigenvalue equation A X^n = Ω^n X^n
|
||||
CALL diagonalize_A(fm_A_BSE, homo_red, virtual_red, homo(1), &
|
||||
unit_nr, diag_runtime_est, mp2_env)
|
||||
END IF
|
||||
! Release to avoid faulty use of changed A matrix
|
||||
CALL cp_fm_release(fm_A_BSE)
|
||||
IF (my_do_abba) THEN
|
||||
! Solving eigenvalue equation C Z^n = (Ω^n)^2 Z^n .
|
||||
! Here, the eigenvectors Z^n relate to X^n via
|
||||
! Eq. (A10) in F. Furche J. Chem. Phys., Vol. 114, No. 14, (2001).
|
||||
CALL diagonalize_C(fm_C_BSE, homo_red, virtual_red, homo(1), &
|
||||
fm_sqrt_A_minus_B, fm_inv_sqrt_A_minus_B, &
|
||||
unit_nr, diag_runtime_est, mp2_env)
|
||||
END IF
|
||||
! Release to avoid faulty use of changed C matrix
|
||||
CALL cp_fm_release(fm_C_BSE)
|
||||
END IF
|
||||
|
||||
CALL deallocate_matrices_bse(fm_mat_S_bar_ia_bse, fm_mat_S_bar_ij_bse, &
|
||||
fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
|
||||
fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm)
|
||||
DEALLOCATE (Eigenval_reduced)
|
||||
IF (my_do_iterat_diag) THEN
|
||||
! Contains untested Block-Davidson algorithm
|
||||
CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_iaQ_bse_local, homo(1), virtual(1), mp2_env%ri_g0w0%bse_spin_config, unit_nr, &
|
||||
Eigenval(:, 1, 1), para_env, mp2_env)
|
||||
! Deallocate local 3c-B-matrices
|
||||
DEALLOCATE (B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, B_iaQ_bse_local)
|
||||
END IF
|
||||
|
||||
CALL timestop(handle)
|
||||
|
||||
END SUBROUTINE start_bse_calculation
|
||||
|
||||
END MODULE bse_main
|
||||
1351
src/bse_util.F
Normal file
1351
src/bse_util.F
Normal file
File diff suppressed because it is too large
Load diff
|
|
@ -1278,5 +1278,22 @@ MODULE input_constants
|
|||
INTEGER, PARAMETER, PUBLIC :: e_dens_total_hard_approx = 1, &
|
||||
e_dens_total_density = 2, &
|
||||
e_dens_soft_density = 3
|
||||
! BSE spin configuration
|
||||
INTEGER, PARAMETER, PUBLIC :: bse_singlet = 0, &
|
||||
bse_triplet = 1
|
||||
|
||||
! BSE method of diagonalization (full diagonalization or iteratively)
|
||||
INTEGER, PARAMETER, PUBLIC :: bse_fulldiag = 0, &
|
||||
bse_iterdiag = 1
|
||||
|
||||
! BSE level of approximation to the ABBA matrix (i.e. with/without TDA)
|
||||
INTEGER, PARAMETER, PUBLIC :: bse_tda = 0, &
|
||||
bse_abba = 1, &
|
||||
bse_both = 2
|
||||
|
||||
! BSE iterative abortion condition
|
||||
INTEGER, PARAMETER, PUBLIC :: bse_iter_en_cond = 0, &
|
||||
bse_iter_res_cond = 1, &
|
||||
bse_iter_both_cond = 2
|
||||
|
||||
END MODULE input_constants
|
||||
|
|
|
|||
|
|
@ -12,45 +12,47 @@
|
|||
!> \author MDB
|
||||
! **************************************************************************************************
|
||||
MODULE input_cp2k_mp2
|
||||
USE bibliography, ONLY: &
|
||||
Bates2013, DelBen2012, DelBen2013, DelBen2015, DelBen2015b, Rybkin2016, Wilhelm2016a, &
|
||||
Wilhelm2016b, Wilhelm2017, Wilhelm2018
|
||||
USE cp_eri_mme_interface, ONLY: create_eri_mme_section
|
||||
USE cp_output_handling, ONLY: add_last_numeric,&
|
||||
cp_print_key_section_create,&
|
||||
debug_print_level,&
|
||||
high_print_level,&
|
||||
low_print_level,&
|
||||
medium_print_level,&
|
||||
silent_print_level
|
||||
USE cp_units, ONLY: cp_unit_to_cp2k
|
||||
USE input_constants, ONLY: &
|
||||
do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, &
|
||||
do_potential_long, do_potential_mix_cl, do_potential_short, do_potential_truncated, &
|
||||
do_potential_tshpsc, eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, &
|
||||
gw_read_exx, gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, &
|
||||
kp_weights_W_tailored, kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, &
|
||||
mp2_method_none, numerical, ot_precond_full_all, ot_precond_full_kinetic, &
|
||||
ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
|
||||
ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, &
|
||||
ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, &
|
||||
wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, &
|
||||
z_solver_sd
|
||||
USE input_cp2k_hfx, ONLY: create_hfx_section
|
||||
USE input_cp2k_kpoints, ONLY: create_kpoint_set_section
|
||||
USE input_keyword_types, ONLY: keyword_create,&
|
||||
keyword_release,&
|
||||
keyword_type
|
||||
USE input_section_types, ONLY: section_add_keyword,&
|
||||
section_add_subsection,&
|
||||
section_create,&
|
||||
section_release,&
|
||||
section_type
|
||||
USE input_val_types, ONLY: integer_t,&
|
||||
real_t
|
||||
USE kinds, ONLY: dp
|
||||
USE string_utilities, ONLY: newline,&
|
||||
s2a
|
||||
USE bibliography, ONLY: &
|
||||
Bates2013, DelBen2012, DelBen2013, DelBen2015, DelBen2015b, Rybkin2016, Wilhelm2016a, &
|
||||
Wilhelm2016b, Wilhelm2017, Wilhelm2018
|
||||
USE cp_eri_mme_interface, ONLY: create_eri_mme_section
|
||||
USE cp_output_handling, ONLY: add_last_numeric, &
|
||||
cp_print_key_section_create, &
|
||||
debug_print_level, &
|
||||
high_print_level, &
|
||||
low_print_level, &
|
||||
medium_print_level, &
|
||||
silent_print_level
|
||||
USE cp_units, ONLY: cp_unit_to_cp2k
|
||||
USE input_constants, ONLY: &
|
||||
bse_fulldiag, bse_iterdiag, bse_tda, bse_abba, bse_both, &
|
||||
bse_iter_both_cond, bse_iter_en_cond, bse_iter_res_cond, bse_singlet, &
|
||||
bse_triplet, do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, &
|
||||
do_potential_long, do_potential_mix_cl, do_potential_short, do_potential_truncated, &
|
||||
do_potential_tshpsc, eri_default, gaussian, gw_no_print_exx, gw_pade_approx, gw_print_exx, &
|
||||
gw_read_exx, gw_skip_for_regtest, gw_two_pole_model, kp_weights_W_auto, &
|
||||
kp_weights_W_tailored, kp_weights_W_uniform, mp2_method_direct, mp2_method_gpw, &
|
||||
mp2_method_none, numerical, ot_precond_full_all, ot_precond_full_kinetic, &
|
||||
ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, &
|
||||
ot_precond_s_inverse, ri_default, ri_rpa_g0w0_crossing_bisection, &
|
||||
ri_rpa_g0w0_crossing_newton, ri_rpa_g0w0_crossing_z_shot, soc_lda, soc_none, soc_pbe, &
|
||||
wfc_mm_style_gemm, wfc_mm_style_syrk, z_solver_cg, z_solver_pople, z_solver_richardson, &
|
||||
z_solver_sd
|
||||
USE input_cp2k_hfx, ONLY: create_hfx_section
|
||||
USE input_cp2k_kpoints, ONLY: create_kpoint_set_section
|
||||
USE input_keyword_types, ONLY: keyword_create, &
|
||||
keyword_release, &
|
||||
keyword_type
|
||||
USE input_section_types, ONLY: section_add_keyword, &
|
||||
section_add_subsection, &
|
||||
section_create, &
|
||||
section_release, &
|
||||
section_type
|
||||
USE input_val_types, ONLY: integer_t, &
|
||||
real_t
|
||||
USE kinds, ONLY: dp
|
||||
USE string_utilities, ONLY: newline, &
|
||||
s2a
|
||||
#include "./base/base_uses.f90"
|
||||
|
||||
IMPLICIT NONE
|
||||
|
|
@ -852,16 +854,6 @@ CONTAINS
|
|||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="BSE", &
|
||||
description="If true, electronic excitation energies are computed from the "// &
|
||||
"Bethe-Salpeter equation on top of GW eigenvalues. Parameter of BSE can be adjusted in "// &
|
||||
"the corresponding section.", &
|
||||
usage="BSE", &
|
||||
default_l_val=.FALSE., &
|
||||
lone_keyword_l_val=.TRUE.)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="IMAGE_CHARGE_MODEL", &
|
||||
variants=(/"IC"/), &
|
||||
description="If true, an image charge model is applied to mimic the renormalization of "// &
|
||||
|
|
@ -1215,39 +1207,205 @@ CONTAINS
|
|||
TYPE(section_type), POINTER :: section
|
||||
|
||||
TYPE(keyword_type), POINTER :: keyword
|
||||
TYPE(section_type), POINTER :: subsection
|
||||
|
||||
CPASSERT(.NOT. ASSOCIATED(section))
|
||||
CALL section_create(section, __LOCATION__, name="BSE", &
|
||||
description="Parameters influencing the Bethe-Salpeter calculations on "// &
|
||||
"top of GW eigenvalues.", &
|
||||
n_keywords=12, n_subsections=1, repeats=.FALSE.)
|
||||
description="Parameters for a calculation solving the Bethe-Salpeter equation "// &
|
||||
"(BSE) for electronic excitations. The full BSE "// &
|
||||
"$\left( \begin{array}{cc}A & B\\B & A\end{array} \right)\left( \begin{array}{cc}X^{(n)}\\Y^{(n)}\end{array} \right) = "// &
|
||||
"\Omega^{(n)}\left(\begin{array}{cc}1&0\\0&-1\end{array}\right)\left(\begin{array}{cc}X^{(n)}\\Y^{(n)}\end{array}\right)$ "// &
|
||||
"enables of electronic excitation energies $\Omega^{(n)}$. The BSE can be solved by diagonalizing "// &
|
||||
"the full ABBA-matrix or by setting B=0, i.e. within the Tamm-Dancoff approximation (TDA). "// &
|
||||
"Preliminary reference: Eq. (35) in PRB 92, 045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209", &
|
||||
n_keywords=8, n_subsections=2, repeats=.FALSE.)
|
||||
|
||||
NULLIFY (keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="NUM_Z_VECTORS", &
|
||||
description="Number Z vectors used in the subspace iterations. This is a convergence "// &
|
||||
"parameter: Increasing NUM_Z_VECTORS should keep the result constant in case of convergence,.", &
|
||||
usage="NUM_Z_VECTORS 50", &
|
||||
default_i_val=20)
|
||||
CALL keyword_create(keyword, __LOCATION__, name="_SECTION_PARAMETERS_", &
|
||||
description="Activates BSE calculations.", &
|
||||
usage="&BSE .TRUE.", &
|
||||
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="THRESHOLD_MIN_TRANS", &
|
||||
variants=(/"EPS"/), &
|
||||
CALL keyword_create(keyword, __LOCATION__, name="SPIN_CONFIG", &
|
||||
description="Choose between calculation of singlet or triplet excitation (cf. given Reference above).", &
|
||||
usage="SPIN_CONFIG TRIPLET", &
|
||||
enum_c_vals=s2a("SINGLET", "TRIPLET"), &
|
||||
enum_i_vals=(/bse_singlet, bse_triplet/), &
|
||||
enum_desc=s2a("Computes singlet excitations.", &
|
||||
"Computes triplet excitations."), &
|
||||
default_i_val=bse_singlet)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="BSE_DIAG_METHOD", &
|
||||
description="Method for BSE calculations. "// &
|
||||
"Choose between full or iterative diagonalization.", &
|
||||
usage="&BSE_DIAG_METHOD FULLDIAG", &
|
||||
enum_c_vals=s2a("FULLDIAG", "ITERDIAG"), &
|
||||
enum_i_vals=(/bse_fulldiag, bse_iterdiag/), &
|
||||
enum_desc=s2a("Fully diagonalizes the BSE matrices within the chosen level of approximation.", &
|
||||
"Iterative diagonalization has not been implemented yet."), &
|
||||
default_i_val=bse_fulldiag)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="BSE_APPROX", &
|
||||
description="Level of approximation applied to BSE calculations. "// &
|
||||
"Choose between Tamm Dancoff approximation (TDA) and diagonalization of the full ABBA-matrix.", &
|
||||
usage="&BSE_APPROX TDA", &
|
||||
enum_c_vals=s2a("TDA", "ABBA", "BOTH"), &
|
||||
enum_i_vals=(/bse_tda, bse_abba, bse_both/), &
|
||||
enum_desc=s2a("The TDA is applied, i.e. B=0.", &
|
||||
"The ABBA-matrix is diagonalized, i.e. the TDA is not applied.", &
|
||||
"The BSE is solved within the TDA (B=0) as well as for the full ABBA-matrix."), &
|
||||
default_i_val=bse_tda)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_CUTOFF_OCC", &
|
||||
description="Removing all orbitals with indices i,j from A_ia,jb and B_ia,jb with energy difference "// &
|
||||
"to HOMO level larger than the given energy cutoff, i.e. "// &
|
||||
"$\epsilon_{\mathrm{HOMO}}-\epsilon_a>\mathtt{ENERGY_CUTOFF_OCC}$. "// &
|
||||
"Can be used to accelerate runtime and reduce memory consumption.", &
|
||||
usage="ENERGY_CUTOFF_OCC 10.0", unit_str="eV", &
|
||||
type_of_var=real_t, default_r_val=-1.0_dp)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="ENERGY_CUTOFF_VIRT", &
|
||||
description="Removing all orbitals with indices a,b from A_ia,jb and B_ia,jb with energy difference "// &
|
||||
"to LUMO level larger than the given energy cutoff, i.e. "// &
|
||||
"$\epsilon_i-\epsilon_{\mathrm{LUMO}}>\mathtt{ENERGY_CUTOFF_VIRT}$. "// &
|
||||
"Can be used to accelerate runtime and reduce memory consumption.", &
|
||||
usage="ENERGY_CUTOFF_VIRT 10.0", unit_str="eV", &
|
||||
type_of_var=real_t, default_r_val=-1.0_dp)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="BSE_DEBUG_PRINT", &
|
||||
description="Activates debug print statements in the BSE calculation.", &
|
||||
usage="&BSE_DEBUG_PRINT .TRUE.", &
|
||||
default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="NUM_PRINT_EXC", &
|
||||
description="Number of printed excitation energies. Does not affect computation time.", &
|
||||
usage="NUM_PRINT_EXC 10", &
|
||||
default_i_val=10)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="EPS_X", &
|
||||
description="Threshold for printing contributions of singleparticle "// &
|
||||
"transitions, i.e. elements of the eigenvector $X_{ia}^n$.", &
|
||||
usage="EPS_X 0.1", &
|
||||
type_of_var=real_t, default_r_val=0.1_dp)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
NULLIFY (subsection)
|
||||
CALL create_bse_iterat_section(subsection)
|
||||
CALL section_add_subsection(section, subsection)
|
||||
CALL section_release(subsection)
|
||||
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
!> \brief ...
|
||||
!> \param section ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE create_bse_iterat_section(section)
|
||||
TYPE(section_type), POINTER :: section
|
||||
|
||||
TYPE(keyword_type), POINTER :: keyword
|
||||
|
||||
CPASSERT(.NOT. ASSOCIATED(section))
|
||||
CALL section_create(section, __LOCATION__, name="BSE_ITERAT", &
|
||||
description="Parameters influencing the iterative Bethe-Salpeter calculation. "// &
|
||||
"The iterative solver has not been fully implemented yet.", &
|
||||
n_keywords=9, n_subsections=0, repeats=.FALSE.)
|
||||
|
||||
NULLIFY (keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="DAVIDSON_ABORT_COND", &
|
||||
description="Desired abortion condition for Davidson solver", &
|
||||
usage="DAVIDSON_ABORT_COND OR", &
|
||||
enum_c_vals=s2a("EN", "RES", "OR"), &
|
||||
enum_i_vals=(/bse_iter_en_cond, bse_iter_res_cond, bse_iter_both_cond/), &
|
||||
enum_desc=s2a("Uses energy threshold for successfully exiting solver.", &
|
||||
"Uses residual threshold for successfully exiting solver.", &
|
||||
"Uses either energy or residual threshold for successfully exiting solver."), &
|
||||
default_i_val=bse_iter_en_cond)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="NUM_EXC_EN", &
|
||||
description="Number of lowest excitation energies to be computed.", &
|
||||
usage="NUM_EXC_EN 3", &
|
||||
default_i_val=3)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="NUM_ADD_START_Z_SPACE", &
|
||||
description="Determines the initial dimension of the subspace as "// &
|
||||
"dim = (NUM_EXC_EN+NUM_ADD_START_Z_SPACE)", &
|
||||
usage="NUM_ADD_START_Z_SPACE 1", &
|
||||
default_i_val=0)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="FAC_MAX_Z_SPACE", &
|
||||
description="Factor to determine maximum dimension of the Davidson subspace. "// &
|
||||
"dimension = (NUM_EXC_EN+NUM_ADD_START_Z_SPACE)*FAC_MAX_Z_SPACE", &
|
||||
usage="FAC_MAX_Z_SPACE 5", &
|
||||
default_i_val=5)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="NUM_NEW_T", &
|
||||
description="Number of new t vectors added. "// &
|
||||
"Must be smaller/equals (NUM_EXC_EN+NUM_ADD_START_Z_SPACE)", &
|
||||
usage="NUM_NEW_T 4", &
|
||||
default_i_val=1)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="EPS_RES", &
|
||||
description="Threshold for stopping the iteration for computing the transition energies. "// &
|
||||
"If the lowest excitation changes by less than THRESHOLD_MIN_TRANS (in eV), the ieration stops.", &
|
||||
usage="THRESHOLD_MIN_TRANS 0.001", &
|
||||
default_r_val=0.001_dp)
|
||||
"If the residuals inside the Davidson space change by less than EPS_RES (in eV), the iteration "// &
|
||||
"stops.", &
|
||||
usage="EPS_RES 0.001", unit_str="eV", &
|
||||
type_of_var=real_t, default_r_val=0.001_dp)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="MAX_ITER", &
|
||||
CALL keyword_create(keyword, __LOCATION__, name="EPS_EXC_EN", &
|
||||
description="Threshold for stopping the iteration for computing the transition energies. "// &
|
||||
"If the desired excitation energies change by less than EPS_EXC_EN (in eV), the iteration "// &
|
||||
"stops.", &
|
||||
usage="EPS_EXC_EN 0.001", unit_str="eV", &
|
||||
type_of_var=real_t, default_r_val=0.001_dp)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="NUM_DAVIDSON_ITER", &
|
||||
description="Maximum number of iterations for determining the transition energies.", &
|
||||
usage="MAX_ITER 200", &
|
||||
default_i_val=200)
|
||||
usage="MAX_ITER 100", &
|
||||
default_i_val=100)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
|
||||
CALL keyword_create(keyword, __LOCATION__, name="Z_SPACE_ENERGY_CUTOFF", &
|
||||
description="Cutoff (in eV) for maximal energy difference entering the A matrix. "// &
|
||||
"Per default and for negative values, there is no cutoff applied.", &
|
||||
usage="Z_SPACE_ENERGY_CUTOFF 60", unit_str="eV", &
|
||||
type_of_var=real_t, default_r_val=-1.0_dp)
|
||||
CALL section_add_keyword(section, keyword)
|
||||
CALL keyword_release(keyword)
|
||||
END SUBROUTINE
|
||||
|
||||
! **************************************************************************************************
|
||||
|
|
|
|||
|
|
@ -1182,13 +1182,9 @@ CONTAINS
|
|||
CALL timeset(routineN//"_E_Ex_"//TRIM(intermed_mat%descr), handle)
|
||||
CALL copy_dbcsr_to_fm(intermed_mat%matrix_ia_jb, intermed_mat%fm_BIb_jb)
|
||||
|
||||
IF (.NOT. (TRIM(intermed_mat%descr) .EQ. "bse_ab")) THEN
|
||||
|
||||
CALL grep_my_integrals(para_env, intermed_mat%fm_BIb_jb, BIb_jb, intermed_mat%max_row_col_local, &
|
||||
intermed_mat%local_col_row_info, &
|
||||
my_B_end, my_B_start)
|
||||
|
||||
END IF
|
||||
CALL grep_my_integrals(para_env, intermed_mat%fm_BIb_jb, BIb_jb, intermed_mat%max_row_col_local, &
|
||||
intermed_mat%local_col_row_info, &
|
||||
my_B_end, my_B_start)
|
||||
|
||||
CALL timestop(handle)
|
||||
END SUBROUTINE ao_to_mo_and_store_B
|
||||
|
|
|
|||
|
|
@ -71,9 +71,9 @@ CONTAINS
|
|||
do_ri_sos_mp2, do_rpa
|
||||
REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals
|
||||
TYPE(cp_logger_type), POINTER :: logger
|
||||
TYPE(section_vals_type), POINTER :: cphf_section, eri_mme_section, &
|
||||
gw_section, low_scaling_section, &
|
||||
mp2_section
|
||||
TYPE(section_vals_type), POINTER :: bse_section, cphf_section, &
|
||||
eri_mme_section, gw_section, &
|
||||
low_scaling_section, mp2_section
|
||||
|
||||
CALL timeset(routineN, handle)
|
||||
logger => cp_get_default_logger()
|
||||
|
|
@ -168,14 +168,44 @@ CONTAINS
|
|||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%SOC_ENERGY_WINDOW", &
|
||||
r_val=mp2_env%ri_g0w0%soc_energy_window)
|
||||
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE", &
|
||||
NULLIFY (bse_section)
|
||||
bse_section => section_vals_get_subs_vals(mp2_section, "RI_RPA%GW%BSE")
|
||||
CALL section_vals_val_get(bse_section, "_SECTION_PARAMETERS_", &
|
||||
l_val=mp2_env%ri_g0w0%do_bse)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%NUM_Z_VECTORS", &
|
||||
i_val=mp2_env%ri_g0w0%num_z_vectors)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%THRESHOLD_MIN_TRANS", &
|
||||
r_val=mp2_env%ri_g0w0%eps_min_trans)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%MAX_ITER", &
|
||||
i_val=mp2_env%ri_g0w0%max_iter_bse)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%SPIN_CONFIG", &
|
||||
i_val=mp2_env%ri_g0w0%bse_spin_config)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%ENERGY_CUTOFF_OCC", &
|
||||
r_val=mp2_env%ri_g0w0%bse_cutoff_occ)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%ENERGY_CUTOFF_VIRT", &
|
||||
r_val=mp2_env%ri_g0w0%bse_cutoff_virt)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%NUM_PRINT_EXC", &
|
||||
i_val=mp2_env%ri_g0w0%num_print_exc)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%EPS_X", &
|
||||
r_val=mp2_env%ri_g0w0%eps_x)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%DAVIDSON_ABORT_COND", &
|
||||
i_val=mp2_env%ri_g0w0%davidson_abort_cond)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_EXC_EN", &
|
||||
i_val=mp2_env%ri_g0w0%num_exc_en)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_ADD_START_Z_SPACE", &
|
||||
i_val=mp2_env%ri_g0w0%num_add_start_z_space)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%FAC_MAX_Z_SPACE", &
|
||||
i_val=mp2_env%ri_g0w0%fac_max_z_space)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_NEW_T", &
|
||||
i_val=mp2_env%ri_g0w0%num_new_t)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%EPS_RES", &
|
||||
r_val=mp2_env%ri_g0w0%eps_res)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%EPS_EXC_EN", &
|
||||
r_val=mp2_env%ri_g0w0%eps_exc_en)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%NUM_DAVIDSON_ITER", &
|
||||
i_val=mp2_env%ri_g0w0%num_davidson_iter)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_ITERAT%Z_SPACE_ENERGY_CUTOFF", &
|
||||
r_val=mp2_env%ri_g0w0%z_space_energy_cutoff)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_DIAG_METHOD", &
|
||||
i_val=mp2_env%ri_g0w0%bse_diag_method)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_APPROX", &
|
||||
i_val=mp2_env%ri_g0w0%bse_approx)
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%BSE%BSE_DEBUG_PRINT", &
|
||||
l_val=mp2_env%ri_g0w0%bse_debug_print)
|
||||
|
||||
CALL section_vals_val_get(mp2_section, "RI_RPA%GW%IMAGE_CHARGE_MODEL", &
|
||||
l_val=mp2_env%ri_g0w0%do_ic_model)
|
||||
|
|
|
|||
|
|
@ -220,11 +220,25 @@ MODULE mp2_types
|
|||
do_aux_bas_gw = .FALSE.
|
||||
REAL(KIND=dp) :: frac_aux_mos = 0.0_dp
|
||||
INTEGER :: num_omega_points = 0
|
||||
LOGICAL :: do_bse = .FALSE.
|
||||
INTEGER :: num_z_vectors = 0, &
|
||||
max_iter_bse = 0
|
||||
REAL(KIND=dp) :: eps_min_trans = 0.0_dp
|
||||
LOGICAL :: do_ic_model = .FALSE., &
|
||||
INTEGER :: bse_spin_config = 0, &
|
||||
bse_diag_method = 0, &
|
||||
bse_approx = 0, &
|
||||
num_exc_en = 0, &
|
||||
num_print_exc = 0, &
|
||||
num_add_start_z_space = 0, &
|
||||
fac_max_z_space = 0, &
|
||||
num_new_t = 0, &
|
||||
num_davidson_iter = 0, &
|
||||
davidson_abort_cond = 0
|
||||
REAL(KIND=dp) :: eps_res = 0.0_dp, &
|
||||
eps_exc_en = 0.0_dp, &
|
||||
eps_x = 0.0_dp, &
|
||||
bse_cutoff_occ = 0.0_dp, &
|
||||
bse_cutoff_virt = 0.0_dp, &
|
||||
z_space_energy_cutoff = 0.0_dp
|
||||
LOGICAL :: do_bse = .FALSE., &
|
||||
bse_debug_print = .FALSE., &
|
||||
do_ic_model = .FALSE., &
|
||||
print_ic_values = .FALSE.
|
||||
REAL(KIND=dp) :: eps_dist = 0.0_dp
|
||||
TYPE(one_dim_real_array), DIMENSION(2) :: ic_corr_list = one_dim_real_array(NULL())
|
||||
|
|
|
|||
11
src/rpa_gw.F
11
src/rpa_gw.F
|
|
@ -773,6 +773,7 @@ CONTAINS
|
|||
!> \param jquad ...
|
||||
!> \param nmo ...
|
||||
!> \param num_fit_points ...
|
||||
!> \param num_integ_points ...
|
||||
!> \param do_bse ...
|
||||
!> \param do_im_time ...
|
||||
!> \param do_periodic ...
|
||||
|
|
@ -800,8 +801,8 @@ CONTAINS
|
|||
! **************************************************************************************************
|
||||
SUBROUTINE compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI, gw_corr_lev_occ, &
|
||||
gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
|
||||
do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, &
|
||||
fermi_level_offset, &
|
||||
num_integ_points, do_bse, do_im_time, do_periodic, &
|
||||
first_cycle_periodic_correction, fermi_level_offset, &
|
||||
omega, Eigenval, delta_corr, vec_omega_fit_gw, vec_W_gw, wj, &
|
||||
fm_mat_Q, fm_mat_Q_static_bse, fm_mat_R_gw, fm_mat_S_gw, &
|
||||
fm_mat_S_gw_work, mo_coeff, para_env, &
|
||||
|
|
@ -812,7 +813,8 @@ CONTAINS
|
|||
DIMENSION(:, :, :, :), INTENT(INOUT) :: vec_Sigma_c_gw
|
||||
INTEGER, INTENT(IN) :: dimen_nm_gw, dimen_RI
|
||||
INTEGER, DIMENSION(:), INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt, homo
|
||||
INTEGER, INTENT(IN) :: jquad, nmo, num_fit_points
|
||||
INTEGER, INTENT(IN) :: jquad, nmo, num_fit_points, &
|
||||
num_integ_points
|
||||
LOGICAL, INTENT(IN) :: do_bse, do_im_time, do_periodic
|
||||
LOGICAL, INTENT(INOUT) :: first_cycle_periodic_correction
|
||||
REAL(KIND=dp), INTENT(INOUT) :: fermi_level_offset, omega
|
||||
|
|
@ -858,7 +860,8 @@ CONTAINS
|
|||
! symmetrize the result, fm_mat_R_gw is only temporary work matrix
|
||||
CALL cp_fm_upper_to_full(fm_mat_Q, fm_mat_R_gw)
|
||||
|
||||
IF (do_bse .AND. jquad == 1) THEN
|
||||
! Omega=0 is at last index, not at jquad==1
|
||||
IF (do_bse .AND. jquad == num_integ_points) THEN
|
||||
CALL cp_fm_to_fm(fm_mat_Q, fm_mat_Q_static_bse)
|
||||
END IF
|
||||
|
||||
|
|
|
|||
|
|
@ -18,8 +18,7 @@ MODULE rpa_main
|
|||
USE bibliography, ONLY: &
|
||||
Bates2013, DelBen2013, DelBen2015, Ren2011, Ren2013, Wilhelm2016a, Wilhelm2016b, &
|
||||
Wilhelm2017, Wilhelm2018, cite_reference
|
||||
USE bse, ONLY: do_subspace_iterations,&
|
||||
mult_B_with_W_and_fill_local_3c_arrays
|
||||
USE bse_main, ONLY: start_bse_calculation
|
||||
USE cp_blacs_env, ONLY: cp_blacs_env_create,&
|
||||
cp_blacs_env_release,&
|
||||
cp_blacs_env_type
|
||||
|
|
@ -585,7 +584,6 @@ CONTAINS
|
|||
|
||||
! for Bethe-Salpeter, we need other matrix fm_mat_S
|
||||
IF (do_bse) THEN
|
||||
|
||||
CALL create_integ_mat(BIb_C_2D_bse_ij, para_env, para_env_sub, color_sub, ngroup, integ_group_size, &
|
||||
dimen_RI_red, [dimen_homo_square], color_rpa_group, &
|
||||
mp2_env%block_size_row, mp2_env%block_size_col, unit_nr, &
|
||||
|
|
@ -860,7 +858,7 @@ CONTAINS
|
|||
ELSE
|
||||
grid_2D = 1
|
||||
END IF
|
||||
CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2d)
|
||||
CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env_RPA, grid_2d=grid_2D)
|
||||
|
||||
IF (unit_nr > 0 .AND. .NOT. my_do_im_time) THEN
|
||||
WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
|
||||
|
|
@ -1102,9 +1100,9 @@ CONTAINS
|
|||
COMPLEX(KIND=dp), ALLOCATABLE, &
|
||||
DIMENSION(:, :, :, :) :: vec_Sigma_c_gw
|
||||
INTEGER :: count_ev_sc_GW, cut_memory, group_size_P, gw_corr_lev_tot, handle, handle3, i, &
|
||||
ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, max_iter_bse, min_bsize, mm_style, &
|
||||
nkp, nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, &
|
||||
num_Z_vectors, Pspin, Qspin, size_P
|
||||
ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, min_bsize, mm_style, nkp, &
|
||||
nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, Pspin, Qspin, &
|
||||
size_P
|
||||
INTEGER(int_8) :: dbcsr_nflop
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c
|
||||
INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c
|
||||
|
|
@ -1115,20 +1113,20 @@ CONTAINS
|
|||
first_cycle_periodic_correction, my_open_shell, print_ic_values
|
||||
LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_P_blocks
|
||||
REAL(KIND=dp) :: a_scaling, alpha, dbcsr_time, e_axk, e_axk_corr, eps_filter, &
|
||||
eps_filter_im_time, eps_min_trans, ext_scaling, fermi_level_offset, &
|
||||
fermi_level_offset_input, my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old
|
||||
eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, &
|
||||
my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, &
|
||||
trace_Qomega, vec_omega_fit_gw, wj, &
|
||||
wkp_W
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vec_W_gw, weights_cos_tf_t_to_w, &
|
||||
weights_cos_tf_w_to_t, &
|
||||
weights_sin_tf_t_to_w
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_bar_ijQ_bse_local, B_iaQ_bse_local, Eigenval_last, Eigenval_scf, vec_Sigma_x_gw
|
||||
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Eigenval_last, Eigenval_scf, &
|
||||
vec_Sigma_x_gw
|
||||
TYPE(cp_cfm_type) :: cfm_mat_Q
|
||||
TYPE(cp_fm_type) :: fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, fm_mat_RI_global_work, &
|
||||
fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, &
|
||||
fm_scaled_dm_virt_tau
|
||||
fm_mat_S_ia_bse, fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
|
||||
fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau
|
||||
TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_gw_work, fm_mat_W, &
|
||||
fm_mo_coeff_occ, fm_mo_coeff_virt
|
||||
TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_L_kpoints, fm_mat_Minv_L_kpoints
|
||||
|
|
@ -1311,16 +1309,12 @@ CONTAINS
|
|||
|
||||
IF (do_bse) THEN
|
||||
|
||||
num_Z_vectors = mp2_env%ri_g0w0%num_z_vectors
|
||||
eps_min_trans = mp2_env%ri_g0w0%eps_min_trans
|
||||
max_iter_bse = mp2_env%ri_g0w0%max_iter_bse
|
||||
|
||||
CALL cp_fm_create(fm_mat_Q_static_bse_gemm, fm_mat_Q_gemm(1)%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_mat_Q_gemm(1), fm_mat_Q_static_bse_gemm)
|
||||
CALL cp_fm_set_all(fm_mat_Q_static_bse_gemm, 0.0_dp)
|
||||
|
||||
CALL cp_fm_create(fm_mat_Q_static_bse, fm_mat_Q(1)%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_mat_Q_gemm(1), fm_mat_Q_static_bse)
|
||||
CALL cp_fm_to_fm(fm_mat_Q(1), fm_mat_Q_static_bse)
|
||||
CALL cp_fm_set_all(fm_mat_Q_static_bse, 0.0_dp)
|
||||
|
||||
END IF
|
||||
|
|
@ -1475,14 +1469,16 @@ CONTAINS
|
|||
CALL calc_mat_Q(fm_mat_S(1), do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual(1), &
|
||||
Eigenval(:, 1, 1), Eigenval_scf(:, 1, 1), &
|
||||
homo(1), omega, omega_old, jquad, mm_style, dimen_RI_red, dimen_ia(1), alpha, fm_mat_Q(1), &
|
||||
fm_mat_Q_gemm(1), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter)
|
||||
fm_mat_Q_gemm(1), do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, &
|
||||
num_integ_points)
|
||||
|
||||
IF (my_open_shell) THEN
|
||||
CALL calc_mat_Q(fm_mat_S(2), do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual(2), &
|
||||
Eigenval(:, 1, 2), Eigenval_scf(:, 1, 2), &
|
||||
homo(2), omega, omega_old, jquad, mm_style, &
|
||||
dimen_RI_red, dimen_ia(2), alpha, fm_mat_Q(2), fm_mat_Q_gemm(2), do_bse, &
|
||||
fm_mat_Q_static_bse_gemm, dgemm_counter)
|
||||
fm_mat_Q_static_bse_gemm, dgemm_counter, &
|
||||
num_integ_points)
|
||||
|
||||
! For SOS-MP2 we need both matrices separately
|
||||
IF (.NOT. do_ri_sos_laplace_mp2) THEN
|
||||
|
|
@ -1551,7 +1547,7 @@ CONTAINS
|
|||
END IF
|
||||
ELSE
|
||||
CALL compute_GW_self_energy(vec_Sigma_c_gw, dimen_nm_gw, dimen_RI_red, gw_corr_lev_occ, &
|
||||
gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, &
|
||||
gw_corr_lev_virt, homo, jquad, nmo, num_fit_points, num_integ_points, &
|
||||
do_bse, do_im_time, do_periodic, first_cycle_periodic_correction, &
|
||||
fermi_level_offset, &
|
||||
omega, Eigenval(:, 1, :), delta_corr, vec_omega_fit_gw, vec_W_gw, wj, &
|
||||
|
|
@ -1712,16 +1708,24 @@ CONTAINS
|
|||
|
||||
! postprocessing after GW for Bethe-Salpeter
|
||||
IF (do_bse) THEN
|
||||
CALL mult_B_with_W_and_fill_local_3c_arrays(fm_mat_S_ij_bse, fm_mat_S_ab_bse, fm_mat_S(1), fm_mat_Q_static_bse, &
|
||||
fm_mat_Q_static_bse_gemm, &
|
||||
B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_iaQ_bse_local, dimen_RI_red, homo(1), virtual(1), dimen_ia(1), &
|
||||
gd_array, color_sub, para_env)
|
||||
|
||||
CALL do_subspace_iterations(B_bar_ijQ_bse_local, B_abQ_bse_local, B_bar_iaQ_bse_local, &
|
||||
B_iaQ_bse_local, homo(1), virtual(1), num_Z_vectors, &
|
||||
max_iter_bse, eps_min_trans, Eigenval(:, 1, 1), para_env)
|
||||
|
||||
! Create a copy of fm_mat_S for usage in BSE
|
||||
CALL cp_fm_create(fm_mat_S_ia_bse, fm_mat_S(1)%matrix_struct)
|
||||
CALL cp_fm_to_fm(fm_mat_S(1), fm_mat_S_ia_bse)
|
||||
! Remove energy/frequency factor from 3c-Integral for BSE
|
||||
IF (iter_sc_gw0 == 1) THEN
|
||||
CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), &
|
||||
Eigenval_last(:, 1, 1), homo(1), omega)
|
||||
ELSE
|
||||
CALL remove_scaling_factor_rpa(fm_mat_S_ia_bse, virtual(1), &
|
||||
Eigenval_scf(:, 1, 1), homo(1), omega)
|
||||
END IF
|
||||
! Main routine for all BSE postprocessing
|
||||
CALL start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_bse, &
|
||||
fm_mat_Q_static_bse, fm_mat_Q_static_bse_gemm, &
|
||||
Eigenval, homo, virtual, dimen_RI, dimen_RI_red, &
|
||||
gd_array, color_sub, mp2_env, unit_nr)
|
||||
! Release BSE-copy of fm_mat_S
|
||||
CALL cp_fm_release(fm_mat_S_ia_bse)
|
||||
END IF
|
||||
|
||||
IF (my_do_gw) THEN
|
||||
|
|
|
|||
|
|
@ -653,11 +653,12 @@ CONTAINS
|
|||
!> \param do_bse ...
|
||||
!> \param fm_mat_Q_static_bse_gemm ...
|
||||
!> \param dgemm_counter ...
|
||||
!> \param num_integ_points ...
|
||||
! **************************************************************************************************
|
||||
SUBROUTINE calc_mat_Q(fm_mat_S, do_ri_sos_laplace_mp2, first_cycle, iter_sc_GW0, virtual, &
|
||||
Eigenval, Eigenval_scf, &
|
||||
homo, omega, omega_old, jquad, mm_style, dimen_RI, dimen_ia, alpha, fm_mat_Q, fm_mat_Q_gemm, &
|
||||
do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter)
|
||||
do_bse, fm_mat_Q_static_bse_gemm, dgemm_counter, num_integ_points)
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
|
||||
LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, first_cycle
|
||||
INTEGER, INTENT(IN) :: iter_sc_GW0, virtual
|
||||
|
|
@ -670,6 +671,7 @@ CONTAINS
|
|||
LOGICAL, INTENT(IN) :: do_bse
|
||||
TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q_static_bse_gemm
|
||||
TYPE(dgemm_counter_type), INTENT(INOUT) :: dgemm_counter
|
||||
INTEGER, INTENT(IN) :: num_integ_points
|
||||
|
||||
CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_mat_Q'
|
||||
|
||||
|
|
@ -692,8 +694,10 @@ CONTAINS
|
|||
|
||||
CALL contract_S_to_Q(mm_style, dimen_RI, dimen_ia, alpha, fm_mat_S, fm_mat_Q_gemm, &
|
||||
fm_mat_Q, dgemm_counter)
|
||||
|
||||
IF (do_bse .AND. jquad == 1) THEN
|
||||
! fm_mat_Q_static_bse_gemm does not enter W_ijab (A matrix in TDA), but only full ABBA
|
||||
! (since only B_ij_bar enters W_ijab)
|
||||
! Changing jquad, since omega=0 is at last idx
|
||||
IF (do_bse .AND. jquad == num_integ_points) THEN
|
||||
CALL cp_fm_to_fm(fm_mat_Q_gemm, fm_mat_Q_static_bse_gemm)
|
||||
END IF
|
||||
CALL timestop(handle)
|
||||
|
|
|
|||
80
tests/QS/regtest-bse/BSE_H2O_PBE.inp
Normal file
80
tests/QS/regtest-bse/BSE_H2O_PBE.inp
Normal file
|
|
@ -0,0 +1,80 @@
|
|||
&GLOBAL
|
||||
PRINT_LEVEL MEDIUM
|
||||
PROJECT BSE_H2O_PBE
|
||||
RUN_TYPE ENERGY
|
||||
&TIMINGS
|
||||
THRESHOLD 0.01
|
||||
&END TIMINGS
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
METHOD Quickstep
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME HFX_BASIS
|
||||
POTENTIAL_FILE_NAME GTH_POTENTIALS
|
||||
&MGRID
|
||||
CUTOFF 100
|
||||
REL_CUTOFF 20
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
POISSON_SOLVER WAVELET
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-15
|
||||
EPS_PGF_ORB 1.0E-30
|
||||
METHOD GPW
|
||||
&END QS
|
||||
&SCF
|
||||
EPS_SCF 1.0E-7
|
||||
MAX_SCF 100
|
||||
SCF_GUESS ATOMIC
|
||||
&PRINT
|
||||
&RESTART OFF
|
||||
&END RESTART
|
||||
&END PRINT
|
||||
&END SCF
|
||||
&XC
|
||||
&WF_CORRELATION
|
||||
&RI_RPA
|
||||
QUADRATURE_POINTS 100
|
||||
&GW
|
||||
CORR_MOS_OCC 1000
|
||||
CORR_MOS_VIRT 1000
|
||||
RI_SIGMA_X
|
||||
&BSE
|
||||
BSE_APPROX ABBA
|
||||
BSE_DIAG_METHOD FULLDIAG
|
||||
ENERGY_CUTOFF_OCC 60.0
|
||||
ENERGY_CUTOFF_VIRT 60.0
|
||||
&END BSE
|
||||
&END GW
|
||||
&END RI_RPA
|
||||
&END WF_CORRELATION
|
||||
&XC_FUNCTIONAL PBE
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC [angstrom] 6.000 6.000 6.000
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&KIND H
|
||||
BASIS_SET DZVP-GTH
|
||||
BASIS_SET RI_AUX RI_DZVP-GTH
|
||||
POTENTIAL GTH-PBE-q1
|
||||
&END KIND
|
||||
&KIND O
|
||||
BASIS_SET DZVP-GTH
|
||||
BASIS_SET RI_AUX RI_DZVP-GTH
|
||||
POTENTIAL GTH-PBE-q6
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
COORD_FILE_FORMAT xyz
|
||||
COORD_FILE_NAME H2O_gas.xyz
|
||||
&CENTER_COORDINATES
|
||||
&END CENTER_COORDINATES
|
||||
&END TOPOLOGY
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
5
tests/QS/regtest-bse/H2O_gas.xyz
Normal file
5
tests/QS/regtest-bse/H2O_gas.xyz
Normal 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
|
||||
80
tests/QS/regtest-bse/TDA_H2O_PBE.inp
Normal file
80
tests/QS/regtest-bse/TDA_H2O_PBE.inp
Normal file
|
|
@ -0,0 +1,80 @@
|
|||
&GLOBAL
|
||||
PRINT_LEVEL MEDIUM
|
||||
PROJECT BSE_H2O_PBE
|
||||
RUN_TYPE ENERGY
|
||||
&TIMINGS
|
||||
THRESHOLD 0.01
|
||||
&END TIMINGS
|
||||
&END GLOBAL
|
||||
|
||||
&FORCE_EVAL
|
||||
METHOD Quickstep
|
||||
&DFT
|
||||
BASIS_SET_FILE_NAME HFX_BASIS
|
||||
POTENTIAL_FILE_NAME GTH_POTENTIALS
|
||||
&MGRID
|
||||
CUTOFF 100
|
||||
REL_CUTOFF 20
|
||||
&END MGRID
|
||||
&POISSON
|
||||
PERIODIC NONE
|
||||
POISSON_SOLVER WAVELET
|
||||
&END POISSON
|
||||
&QS
|
||||
EPS_DEFAULT 1.0E-15
|
||||
EPS_PGF_ORB 1.0E-30
|
||||
METHOD GPW
|
||||
&END QS
|
||||
&SCF
|
||||
EPS_SCF 1.0E-7
|
||||
MAX_SCF 100
|
||||
SCF_GUESS ATOMIC
|
||||
&PRINT
|
||||
&RESTART OFF
|
||||
&END RESTART
|
||||
&END PRINT
|
||||
&END SCF
|
||||
&XC
|
||||
&WF_CORRELATION
|
||||
&RI_RPA
|
||||
QUADRATURE_POINTS 100
|
||||
&GW
|
||||
CORR_MOS_OCC 1000
|
||||
CORR_MOS_VIRT 1000
|
||||
RI_SIGMA_X
|
||||
&BSE
|
||||
BSE_APPROX TDA
|
||||
BSE_DIAG_METHOD FULLDIAG
|
||||
ENERGY_CUTOFF_OCC 60.0
|
||||
ENERGY_CUTOFF_VIRT 60.0
|
||||
&END BSE
|
||||
&END GW
|
||||
&END RI_RPA
|
||||
&END WF_CORRELATION
|
||||
&XC_FUNCTIONAL PBE
|
||||
&END XC_FUNCTIONAL
|
||||
&END XC
|
||||
&END DFT
|
||||
&SUBSYS
|
||||
&CELL
|
||||
ABC [angstrom] 6.000 6.000 6.000
|
||||
PERIODIC NONE
|
||||
&END CELL
|
||||
&KIND H
|
||||
BASIS_SET DZVP-GTH
|
||||
BASIS_SET RI_AUX RI_DZVP-GTH
|
||||
POTENTIAL GTH-PBE-q1
|
||||
&END KIND
|
||||
&KIND O
|
||||
BASIS_SET DZVP-GTH
|
||||
BASIS_SET RI_AUX RI_DZVP-GTH
|
||||
POTENTIAL GTH-PBE-q6
|
||||
&END KIND
|
||||
&TOPOLOGY
|
||||
COORD_FILE_FORMAT xyz
|
||||
COORD_FILE_NAME H2O_gas.xyz
|
||||
&CENTER_COORDINATES
|
||||
&END CENTER_COORDINATES
|
||||
&END TOPOLOGY
|
||||
&END SUBSYS
|
||||
&END FORCE_EVAL
|
||||
3
tests/QS/regtest-bse/TEST_FILES
Normal file
3
tests/QS/regtest-bse/TEST_FILES
Normal file
|
|
@ -0,0 +1,3 @@
|
|||
TDA_H2O_PBE.inp 109 1e-05 7.7697
|
||||
BSE_H2O_PBE.inp 110 1e-05 7.7251
|
||||
#EOF
|
||||
|
|
@ -361,3 +361,4 @@ Fist/regtest-allegro libtorch
|
|||
Fist/regtest-deepmd deepmd
|
||||
Fist/regtest-16
|
||||
DFTB/regtest-vdw
|
||||
QS/regtest-bse libint
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
108
|
||||
110
|
||||
Total energy:!3
|
||||
MD| Potential energy!5
|
||||
Total energy \[eV\]:!4
|
||||
|
|
@ -107,6 +107,8 @@ TDDFT+SOC !5
|
|||
G0W0 direct band gap !6
|
||||
SCF+SOC direct band gap !6
|
||||
G0W0+SOC direct band gap !6
|
||||
BSE| 1 Singlet State -TDA- !6
|
||||
BSE| 1 Singlet State -full- !6
|
||||
#
|
||||
# these are the tests the can be selected for regtesting.
|
||||
# do regtest will grep for test_grep (first column) and look if the numeric value
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue