Skip to content

Commit

Permalink
BSE: Add more options for setup of Casida equation
Browse files Browse the repository at this point in the history
  • Loading branch information
@MGraml
MGraml authored Dec 6, 2024
1 parent 1018eda commit 15f7139
Show file tree
Hide file tree
Showing 14 changed files with 507 additions and 56 deletions.
14 changes: 7 additions & 7 deletions docs/methods/properties/optical/bethe-salpeter.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ Diagonalizing $A$ in TDA, or the full block-matrix $ABBA$, takes in the order of
$(N_\mathrm{occ} N_\mathrm{empty})^3$ floating point operations. This translates to a computational
scaling of $O(N^6)$ in the system size $N$.

### 1.2 Optical properties
### 1.2 Optical absorption spectrum

The BSE further allows the investigation of optical properties. For example, the optical absorption
spectrum can be computed as the imaginary part of the dynamical dipole polarizability tensor
Expand Down Expand Up @@ -132,7 +132,7 @@ $$

### 1.3 Visualizing excited states using Natural Transition Orbitals (NTOs)

In order to represent the exciton wave function independent of a specific choice of the molecular
In order to analyse the excitation wave function independent of a specific choice of the molecular
orbitals $\varphi_p(\mathbf{r})$, we can rewrite it as

$$
Expand Down Expand Up @@ -224,7 +224,7 @@ $$
where we drop the excitation index $n$ from now on for better readability.

For each excitation level $n$, we have then several quantities, which allow us to quantify the
spatial extend of the electron, the hole and their combined two-particle character as combined
spatial extent of the electron, the hole and their combined two-particle character as combined
electron-hole pair, i.e. as an exciton. By that, we can often determine the type of the excited
state, i.e. distinguish between, e.g., valence, Rydberg or charge-transfer states
\[[](#Mewes2018)\].
Expand Down Expand Up @@ -401,8 +401,8 @@ is well parallelized, i.e. you can use several nodes that can provide the memory
We have benchmarked the numerical precision of our BSE implementation and compared its results to
the BSE implementation in FHI aims \[[](#Liu2020)\]. For our recommended settings, i.e.
BSE@ev*GW*<sub>0</sub>@PBE with the aug-cc-pVDZ basis set, we have found excellent agreement with
less than 5 meV mean absolute deviation over the first 10 excitation levels and the 28 molecules in
*Thiel's set*.
less than 5 meV mean absolute deviation averaged over the first 10 excitation levels and the 28
molecules in *Thiel's set* for Singlet excitations.

The current BSE implementation in CP2K works for molecules. The inclusion of periodic boundary
conditions in a Γ-only approach and with full *k*-point sampling is work in progress.
Expand Down Expand Up @@ -552,8 +552,8 @@ BSE| 2 -ABBA- 1 1.00210
In the input file, we have specified a list of three excitation levels, no constraint on the
oscillator strengths $f^{(n)}$ and used the default threshold on the weights, given by
$\lambda_I^2 \leq 0.010$, which is repeated at the beginning of the section. Therefore, we obtain
two NTO pairs for $n=1$ with significant weight and one NTO pair for $n=2$ and $n=8$, respectively.
For $n=1$, we can verify that the NTO weights satisfy
two NTO pairs for $n=1$ with significant weight and one NTO pair for $n=2$. For $n=1$, we can verify
that the NTO weights satisfy
$1.01320 + 0.01320 \approx \sum_I {\left( \lambda_I^{(n)}\right)}^2 = c_n \approx 1.026$ when
comparing to the output of the exciton descriptor section.

Expand Down
77 changes: 52 additions & 25 deletions src/bse_full_diag.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,9 @@ MODULE bse_full_diag
cp_fm_set_all,&
cp_fm_to_fm,&
cp_fm_type
USE input_constants, ONLY: bse_singlet,&
USE input_constants, ONLY: bse_screening_alpha,&
bse_screening_rpa,&
bse_singlet,&
bse_triplet
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_para_env_type
Expand Down Expand Up @@ -99,7 +101,7 @@ SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
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
REAL(KIND=dp) :: alpha, alpha_screening, eigen_diff
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(cp_fm_struct_type), POINTER :: fm_struct_A, fm_struct_W
TYPE(cp_fm_type) :: fm_W
Expand All @@ -118,6 +120,12 @@ SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
alpha = 0.0_dp
END SELECT

IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
alpha_screening = mp2_env%bse%screening_factor
ELSE
alpha_screening = 1.0_dp
END IF

! 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)
Expand Down Expand Up @@ -148,10 +156,13 @@ SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
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 infinite screening is applied, fm_W is simply 0 - Otherwise it needs to be computed from 3c integrals
IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
!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=alpha_screening, &
matrix_a=fm_mat_S_bar_ij_bse, matrix_b=fm_mat_S_ab_bse, beta=0.0_dp, &
matrix_c=fm_W)
END IF

IF (unit_nr > 0 .AND. mp2_env%bse%bse_debug_print) THEN
WRITE (unit_nr, '(T2,A10,T13,A16)') 'BSE|DEBUG|', 'Allocated W_ijab'
Expand All @@ -166,9 +177,13 @@ SUBROUTINE create_A(fm_mat_S_ia_bse, fm_mat_S_bar_ij_bse, fm_mat_S_ab_bse, &
! 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)

! If infinite screening is applied, fm_W is simply 0 - Otherwise it needs to be computed from 3c integrals
IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
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)
END IF
!full matrix W is not needed anymore, release it to save memory
CALL cp_fm_release(fm_W)

Expand Down Expand Up @@ -227,7 +242,7 @@ SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &

INTEGER :: handle
INTEGER, DIMENSION(4) :: reordering
REAL(KIND=dp) :: alpha
REAL(KIND=dp) :: alpha, alpha_screening
TYPE(cp_fm_struct_type), POINTER :: fm_struct_v
TYPE(cp_fm_type) :: fm_W

Expand All @@ -245,6 +260,12 @@ SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
alpha = 0.0_dp
END SELECT

IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
alpha_screening = mp2_env%bse%screening_factor
ELSE
alpha_screening = 1.0_dp
END IF

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")
Expand All @@ -261,17 +282,21 @@ SUBROUTINE create_B(fm_mat_S_ia_bse, fm_mat_S_bar_ia_bse, fm_B, &
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)
! If infinite screening is applied, fm_W is simply 0 - Otherwise it needs to be computed from 3c integrals
IF (mp2_env%bse%screening_method /= bse_screening_rpa) THEN
! 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=alpha_screening, &
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)
END IF

CALL cp_fm_release(fm_W)
CALL cp_fm_struct_release(fm_struct_v)
Expand Down Expand Up @@ -454,9 +479,11 @@ SUBROUTINE diagonalize_C(fm_C, homo, virtual, homo_irred, &
! for positive definiteness (would make another O(N^6) Diagon. necessary)
! Instead, we include a check here
IF (Exc_ens(1) < 0) THEN
CALL cp_abort(__LOCATION__, &
"Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
"(A+B) is not positive definite.")
IF (unit_nr > 0) THEN
CALL cp_abort(__LOCATION__, &
"Matrix C=(A-B)^0.5 (A+B) (A-B)^0.5 has negative eigenvalues, i.e. "// &
"(A+B) is not positive definite.")
END IF
END IF
Exc_ens = SQRT(Exc_ens)

Expand Down Expand Up @@ -669,7 +696,7 @@ SUBROUTINE postprocess_bse(Exc_ens, fm_eigvec_X, mp2_env, qs_env, mo_coeff, &

! Prints basic definitions used in BSE calculation
CALL print_output_header(homo, virtual, homo_irred, flag_TDA, &
multiplet, alpha, unit_nr)
multiplet, alpha, mp2_env, unit_nr)

! Prints excitation energies up to user-specified number
CALL print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
Expand Down
32 changes: 26 additions & 6 deletions src/bse_main.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,8 @@ MODULE bse_main
bse_both,&
bse_fulldiag,&
bse_iterdiag,&
bse_screening_alpha,&
bse_screening_tdhf,&
bse_tda
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_para_env_type
Expand Down Expand Up @@ -181,18 +183,36 @@ SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_b
! α 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)

! For unscreened W matrix, we need fm_mat_S_ij_trunc
IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
mp2_env%bse%screening_method == bse_screening_alpha) THEN
CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_ij_trunc, fm_mat_S_ab_trunc, &
fm_A_BSE, Eigenval_reduced, unit_nr, &
homo_red, virtual_red, dimen_RI, mp2_env, &
para_env)
ELSE
CALL create_A(fm_mat_S_ia_trunc, fm_mat_S_bar_ij_bse, fm_mat_S_ab_trunc, &
fm_A_BSE, Eigenval_reduced, unit_nr, &
homo_red, virtual_red, dimen_RI, mp2_env, &
para_env)
END IF
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)

! For unscreened W matrix, we need fm_mat_S_ia_trunc
IF (mp2_env%bse%screening_method == bse_screening_tdhf .OR. &
mp2_env%bse%screening_method == bse_screening_alpha) THEN
CALL create_B(fm_mat_S_ia_trunc, fm_mat_S_ia_trunc, fm_B_BSE, &
homo_red, virtual_red, dimen_RI, unit_nr, mp2_env)
ELSE
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)
END IF
! 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
Expand Down
40 changes: 32 additions & 8 deletions src/bse_print.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,10 @@ MODULE bse_print
USE bse_util, ONLY: filter_eigvec_contrib
USE cp_fm_types, ONLY: cp_fm_get_info,&
cp_fm_type
USE input_constants, ONLY: bse_screening_alpha,&
bse_screening_rpa,&
bse_screening_tdhf,&
bse_screening_w0
USE kinds, ONLY: dp
USE mp2_types, ONLY: mp2_type
USE particle_types, ONLY: particle_type
Expand Down Expand Up @@ -86,15 +90,17 @@ SUBROUTINE print_BSE_start_flag(bse_tda, bse_abba, unit_nr)
!> \param flag_TDA ...
!> \param multiplet ...
!> \param alpha ...
!> \param mp2_env ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE print_output_header(homo, virtual, homo_irred, flag_TDA, &
multiplet, alpha, unit_nr)
multiplet, alpha, mp2_env, unit_nr)

INTEGER, INTENT(IN) :: homo, virtual, homo_irred
LOGICAL, INTENT(IN) :: flag_TDA
CHARACTER(LEN=10), INTENT(IN) :: multiplet
REAL(KIND=dp), INTENT(IN) :: alpha
TYPE(mp2_type), INTENT(INOUT) :: mp2_env
INTEGER, INTENT(IN) :: unit_nr

CHARACTER(LEN=*), PARAMETER :: routineN = 'print_output_header'
Expand Down Expand Up @@ -146,9 +152,17 @@ SUBROUTINE print_output_header(homo, virtual, homo_irred, flag_TDA, &
'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 (mp2_env%bse%screening_method == bse_screening_w0) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab'
ELSE IF (mp2_env%bse%screening_method == bse_screening_rpa) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb'
END IF
IF (.NOT. flag_TDA) THEN
WRITE (unit_nr, '(T2,A4,T7,A32)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
ELSE IF (mp2_env%bse%screening_method == bse_screening_rpa) THEN
WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'B_ia,jb = α * v_ia,jb'
END IF
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 .'
Expand All @@ -159,11 +173,21 @@ SUBROUTINE print_output_header(homo, virtual, homo_irred, flag_TDA, &
WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
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,T7,A7,T31,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
WRITE (unit_nr, '(T2,A4,T7,A7,T31,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
WRITE (unit_nr, '(T2,A4,T7,A3,T31,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
WRITE (unit_nr, '(T2,A4,T7,A6,T30,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
WRITE (unit_nr, '(T2,A4,T7,A,T31,A)') 'BSE|', 'W_... = 1/ϵ v_...:', &
'Direct interaction screened by '
WRITE (unit_nr, '(T2,A4,T30,A)') 'BSE|', &
'dielectric function ϵ(ω=0)'
ELSE IF (mp2_env%bse%screening_method == bse_screening_tdhf) THEN
WRITE (unit_nr, '(T2,A4,T7,A,T30,A)') 'BSE|', 'W_... = v_...:', 'Direct interaction without screening'
ELSE IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
WRITE (unit_nr, '(T2,A4,T7,A,T31,A,F5.2)') 'BSE|', 'W_... = γ v_...:', &
'Direct interaction with artificial screening γ=', mp2_env%bse%screening_factor
END IF
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4)') 'BSE|'
WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
Expand Down
Loading

0 comments on commit 15f7139

Please sign in to comment.