From 15f7139ffc6a4f83bf4d3b38c4ac5b15c6d6b709 Mon Sep 17 00:00:00 2001 From: Max Graml <78024843+MGraml@users.noreply.github.com> Date: Fri, 6 Dec 2024 11:52:36 +0100 Subject: [PATCH] BSE: Add more options for setup of Casida equation --- .../properties/optical/bethe-salpeter.md | 14 ++-- src/bse_full_diag.F | 77 +++++++++++------ src/bse_main.F | 32 +++++-- src/bse_print.F | 40 +++++++-- src/bse_properties.F | 3 + src/bse_util.F | 55 +++++++++++- src/input_constants.F | 6 ++ src/input_cp2k_mp2.F | 67 ++++++++++++++- src/mp2_setup.F | 6 ++ src/mp2_types.F | 5 +- .../BSE_H2O_PBE_evGW0_RPA_screening.inp | 83 ++++++++++++++++++ .../BSE_H2O_PBE_evGW0_TDHF_screening.inp | 83 ++++++++++++++++++ .../BSE_H2O_PBE_evGW0_custom_screening.inp | 84 +++++++++++++++++++ tests/QS/regtest-bse/TEST_FILES | 8 +- 14 files changed, 507 insertions(+), 56 deletions(-) create mode 100644 tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_RPA_screening.inp create mode 100644 tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_TDHF_screening.inp create mode 100644 tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_custom_screening.inp diff --git a/docs/methods/properties/optical/bethe-salpeter.md b/docs/methods/properties/optical/bethe-salpeter.md index ec600ca573..4e4f6c3040 100644 --- a/docs/methods/properties/optical/bethe-salpeter.md +++ b/docs/methods/properties/optical/bethe-salpeter.md @@ -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 @@ -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 $$ @@ -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)\]. @@ -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*0@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. @@ -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. diff --git a/src/bse_full_diag.F b/src/bse_full_diag.F index 5facff40ae..4030c5729d 100644 --- a/src/bse_full_diag.F +++ b/src/bse_full_diag.F @@ -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 @@ -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 @@ -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) @@ -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' @@ -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) @@ -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 @@ -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") @@ -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) @@ -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) @@ -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, & diff --git a/src/bse_main.F b/src/bse_main.F index 7c42cc70b9..ce8550b998 100644 --- a/src/bse_main.F +++ b/src/bse_main.F @@ -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 @@ -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 diff --git a/src/bse_print.F b/src/bse_print.F index 0aa569b690..54a3ff7f68 100644 --- a/src/bse_print.F +++ b/src/bse_print.F @@ -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 @@ -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' @@ -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 .' @@ -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|', & diff --git a/src/bse_properties.F b/src/bse_properties.F index 2cacdf5d13..9fe7b8432a 100644 --- a/src/bse_properties.F +++ b/src/bse_properties.F @@ -275,6 +275,7 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_ ! where f_n are the oscillator strengths and E_exc the excitation energies ! For the full polarizability tensor, we have ! α_{µ µ'}(ω) = \sum_n [2 Ω^n d^n_µ d^n_µ'] / [(ω+iη)^2- (Ω^n)^2] + ! = \sum_n "polarizability_residues" / [(ω+iη)^2- (Ω^n)^2] DO i = 1, num_steps omega = freq_start + (i - 1)*freq_step abs_spectrum(i, 1) = omega @@ -283,6 +284,8 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_ AIMAG(1/((omega + CMPLX(0.0, eta, kind=dp))**2 - Exc_ens(j)**2)) DO idir = 1, 3 DO jdir = 1, 3 + ! Factor 2 from formula for tensor is already in the polarizability_residues + ! to follow the same convention as the oscillator strengths abs_spectrum(i, 2 + (idir - 1)*3 + jdir) = abs_spectrum(i, 2 + (idir - 1)*3 + jdir) & - polarizability_residues(idir, jdir, j)* & AIMAG(1/((omega + CMPLX(0.0, eta, kind=dp))**2 - Exc_ens(j)**2)) diff --git a/src/bse_util.F b/src/bse_util.F index 0dc59265e9..39fea374e3 100644 --- a/src/bse_util.F +++ b/src/bse_util.F @@ -43,7 +43,10 @@ MODULE bse_util USE cp_output_handling, ONLY: cp_print_key_finished_output,& cp_print_key_unit_nr USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube - USE input_constants, ONLY: use_mom_ref_coac + USE input_constants, ONLY: bse_screening_alpha,& + bse_screening_rpa,& + bse_screening_tdhf,& + use_mom_ref_coac USE input_section_types, ONLY: section_vals_type USE kinds, ONLY: default_path_length,& dp,& @@ -1282,7 +1285,12 @@ SUBROUTINE truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_b IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN !Truncate eigenvals ALLOCATE (Eigenval_reduced(homo_red + virt_red)) - Eigenval_reduced(:) = Eigenval(homo_incl:homo + virt_incl) + ! Include USE_KS_ENERGIES input + IF (mp2_env%bse%use_ks_energies) THEN + Eigenval_reduced(:) = Eigenval_scf(homo_incl:homo + virt_incl) + ELSE + Eigenval_reduced(:) = Eigenval(homo_incl:homo + virt_incl) + END IF CALL truncate_fm(fm_mat_S_trunc, fm_mat_S_ia_bse, virtual, & homo_red, virt_red, unit_nr, mp2_env, & @@ -1299,7 +1307,12 @@ SUBROUTINE truncate_BSE_matrices(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_b WRITE (unit_nr, '(T2,A4)') 'BSE|' END IF ALLOCATE (Eigenval_reduced(homo_red + virt_red)) - Eigenval_reduced(:) = Eigenval(:) + ! Include USE_KS_ENERGIES input + IF (mp2_env%bse%use_ks_energies) THEN + Eigenval_reduced(:) = Eigenval_scf(:) + ELSE + Eigenval_reduced(:) = Eigenval(:) + END IF CALL cp_fm_to_fm_submat_general(fm_mat_S_ia_bse, fm_mat_S_trunc, dimen_RI, homo_red*virt_red, & 1, 1, 1, 1, context) CALL cp_fm_to_fm_submat_general(fm_mat_S_ij_bse, fm_mat_S_ij_trunc, dimen_RI, homo_red*homo_red, & @@ -1600,6 +1613,42 @@ SUBROUTINE adapt_BSE_input_params(homo, virtual, unit_nr, mp2_env) mp2_env%bse%num_print_exc_descr = mp2_env%bse%num_print_exc END IF + ! Handle screening factor options + IF (mp2_env%BSE%screening_factor > 0.0_dp) THEN + IF (mp2_env%BSE%screening_method /= bse_screening_alpha) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Screening factor is only supported for &SCREENING_IN_W ALPHA. "// & + "Resetting SCREENING_IN_W to ALPHA.") + END IF + mp2_env%BSE%screening_method = bse_screening_alpha + END IF + IF (mp2_env%BSE%screening_factor > 1.0_dp) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Screening factor is larger than 1.0. ") + END IF + END IF + END IF + + IF (mp2_env%BSE%screening_factor < 0.0_dp .AND. & + mp2_env%BSE%screening_method == bse_screening_alpha) THEN + IF (unit_nr > 0) THEN + CALL cp_warn(__LOCATION__, & + "Screening factor is negative. Defaulting to 0.25") + END IF + mp2_env%BSE%screening_factor = 0.25_dp + END IF + + IF (mp2_env%BSE%screening_factor == 0.0_dp) THEN + ! Use RPA internally in this case + mp2_env%BSE%screening_method = bse_screening_rpa + END IF + IF (mp2_env%BSE%screening_factor == 1.0_dp) THEN + ! Use TDHF internally in this case + mp2_env%BSE%screening_method = bse_screening_tdhf + END IF + CALL timestop(handle) END SUBROUTINE adapt_BSE_input_params diff --git a/src/input_constants.F b/src/input_constants.F index 9198b3cd72..fc063b7db9 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1321,6 +1321,12 @@ MODULE input_constants bse_abba = 1, & bse_both = 2 + ! BSE level of approximation to the screening + INTEGER, PARAMETER, PUBLIC :: bse_screening_w0 = 0, & + bse_screening_tdhf = 1, & + bse_screening_rpa = 2, & + bse_screening_alpha = 3 + ! BSE iterative abortion condition INTEGER, PARAMETER, PUBLIC :: bse_iter_en_cond = 0, & bse_iter_res_cond = 1, & diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F index fc7b0e9c53..91839b1ee9 100644 --- a/src/input_cp2k_mp2.F +++ b/src/input_cp2k_mp2.F @@ -26,6 +26,7 @@ MODULE input_cp2k_mp2 USE cp_units, ONLY: cp_unit_to_cp2k USE input_constants, ONLY: & bse_fulldiag, bse_iterdiag, bse_tda, bse_abba, bse_both, & + bse_screening_w0, bse_screening_tdhf, bse_screening_rpa, bse_screening_alpha, & 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, & @@ -1330,7 +1331,7 @@ SUBROUTINE create_bse_section(section) 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", & + 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.", & @@ -1342,7 +1343,7 @@ SUBROUTINE create_bse_section(section) CALL keyword_create(keyword, __LOCATION__, name="TDA", & description="Level of approximation applied to BSE calculations. "// & "Choose between Tamm Dancoff approximation (TDA) and/or diagonalization of the full ABBA-matrix.", & - usage="&TDA ON", & + usage="TDA ON", & enum_c_vals=s2a("ON", "OFF", "TDA+ABBA"), & enum_i_vals=(/bse_tda, bse_abba, bse_both/), & enum_desc=s2a("The TDA is applied, i.e. B=0.", & @@ -1374,7 +1375,7 @@ SUBROUTINE create_bse_section(section) CALL keyword_create(keyword, __LOCATION__, name="BSE_DEBUG_PRINT", & description="Activates debug print statements in the BSE calculation.", & - usage="&BSE_DEBUG_PRINT .TRUE.", & + usage="BSE_DEBUG_PRINT .TRUE.", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -1398,7 +1399,7 @@ SUBROUTINE create_bse_section(section) CALL keyword_create(keyword, __LOCATION__, name="PRINT_DIRECTIONAL_EXC_DESCR", & description="Activates printing of exciton descriptors per direction.", & - usage="&PRINT_DIRECTIONAL_EXC_DESCR .TRUE.", & + usage="PRINT_DIRECTIONAL_EXC_DESCR .TRUE.", & default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -1411,6 +1412,18 @@ SUBROUTINE create_bse_section(section) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="USE_KS_ENERGIES", & + description="Uses KS energies instead of GW quasiparticle energies.", & + usage="USE_KS_ENERGIES .TRUE.", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + NULLIFY (subsection) + CALL create_bse_screening_section(subsection) + CALL section_add_subsection(section, subsection) + CALL section_release(subsection) + NULLIFY (subsection) CALL create_bse_iterat_section(subsection) CALL section_add_subsection(section, subsection) @@ -1428,6 +1441,52 @@ SUBROUTINE create_bse_section(section) END SUBROUTINE + SUBROUTINE create_bse_screening_section(section) + TYPE(section_type), POINTER :: section + + TYPE(keyword_type), POINTER :: keyword + + CPASSERT(.NOT. ASSOCIATED(section)) + + CALL section_create(section, __LOCATION__, name="SCREENING_IN_W", & + description="Screening $\epsilon$ applied to $W(\omega=0)=\epsilon^{-1}(\omega=0) v $ "// & + "in the BSE calculation. Besides default BSE, i.e. $W_0$ (screening with DFT energies), "// & + "a fixed $\alpha = \epsilon^{-1}(\omega)$ can be applied, which is similar to the mixing "// & + "parameter for hybrid functionals in LR-TDDFT. In addition, the keywords TDHF "// & + "(no screening - $\alpha = 1$) and RPA (infinite screening - $\alpha = 0$) can be applied.", & + n_keywords=2, n_subsections=0, repeats=.FALSE.) + + NULLIFY (keyword) + + CALL keyword_create( & + keyword, __LOCATION__, name="_SECTION_PARAMETERS_", & + description="Shortcut for the most common functional combinations.", & + usage="&xc_functional BLYP", & + enum_c_vals=s2a("W_0", "TDHF", "RPA", "ALPHA"), & + enum_i_vals=(/bse_screening_w0, bse_screening_tdhf, bse_screening_rpa, bse_screening_alpha/), & + enum_desc=s2a("The Coulomb interaction is screened by applying DFT energies "// & + "$\varepsilon_p^{DFT}$, which is typically used for GW-BSE and "// & + "often labeled as $W_0$.", & + "The Coulomb interaction is not screened, i.e. $W_{pq,rs}(\omega=0) "// & + "\rightarrow v_{pq,rs}$ enters.", & + "Infinite screening is applied, i.e. $W_{pq,rs}(\omega=0) \rightarrow 0$.", & + "Arbitrary screening parameter. Specify within section."), & + default_i_val=bse_screening_w0, & + lone_keyword_i_val=bse_screening_w0) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, __LOCATION__, name="ALPHA", & + description="Screening parameter similar to the mixing in hybrid functionals used in TDDFT. "// & + "$\alpha$ mimicks the screening $\epsilon^{-1}(\omega)$ and enforces $W = \alpha v$ "// & + "in the BSE calculation.", & + usage="ALPHA 0.25", & + type_of_var=real_t, default_r_val=-1.00_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + END SUBROUTINE create_bse_screening_section + SUBROUTINE create_bse_nto_section(print_key) TYPE(section_type), POINTER :: print_key diff --git a/src/mp2_setup.F b/src/mp2_setup.F index b3a708453f..dadffb74d5 100644 --- a/src/mp2_setup.F +++ b/src/mp2_setup.F @@ -228,6 +228,12 @@ SUBROUTINE read_mp2_section(input, mp2_env) i_val=mp2_env%bse%bse_diag_method) CALL section_vals_val_get(bse_section, "TDA", & i_val=mp2_env%bse%flag_tda) + CALL section_vals_val_get(bse_section, "USE_KS_ENERGIES", & + l_val=mp2_env%bse%use_ks_energies) + CALL section_vals_val_get(bse_section, "SCREENING_IN_W%_SECTION_PARAMETERS_", & + i_val=mp2_env%bse%screening_method) + CALL section_vals_val_get(bse_section, "SCREENING_IN_W%ALPHA", & + r_val=mp2_env%bse%screening_factor) CALL section_vals_val_get(bse_section, "BSE_DEBUG_PRINT", & l_val=mp2_env%bse%bse_debug_print) CALL section_vals_val_get(bse_section, "BSE_SPECTRUM%_SECTION_PARAMETERS_", & diff --git a/src/mp2_types.F b/src/mp2_types.F index f355eede2e..ff4cd19253 100644 --- a/src/mp2_types.F +++ b/src/mp2_types.F @@ -308,6 +308,7 @@ MODULE mp2_types num_exc_en = 0, & num_print_exc = 0, & num_print_exc_descr = 0, & + screening_method = 0, & num_add_start_z_space = 0, & fac_max_z_space = 0, & num_new_t = 0, & @@ -316,12 +317,14 @@ MODULE mp2_types REAL(KIND=dp) :: eps_res = 0.0_dp, & eps_exc_en = 0.0_dp, & eps_x = 0.0_dp, & + screening_factor = 0.0_dp, & bse_cutoff_occ = 0.0_dp, & bse_cutoff_empty = 0.0_dp, & z_space_energy_cutoff = 0.0_dp LOGICAL :: do_bse = .FALSE., & bse_debug_print = .FALSE., & - print_directional_exc_descr = .FALSE. + print_directional_exc_descr = .FALSE., & + use_ks_energies = .FALSE. !BSE optical spectrum REAL(KIND=dp) :: bse_spectrum_freq_step_size = 0.0_dp, & bse_spectrum_freq_start = 0.0_dp, & diff --git a/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_RPA_screening.inp b/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_RPA_screening.inp new file mode 100644 index 0000000000..89c302b4c4 --- /dev/null +++ b/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_RPA_screening.inp @@ -0,0 +1,83 @@ +&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 + SELF_CONSISTENCY EVGW0 + &BSE + BSE_DIAG_METHOD FULLDIAG + ENERGY_CUTOFF_EMPTY 50.0 + ENERGY_CUTOFF_OCC 60.0 + TDA OFF + &SCREENING_IN_W RPA + &END SCREENING_IN_W + &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 diff --git a/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_TDHF_screening.inp b/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_TDHF_screening.inp new file mode 100644 index 0000000000..bce874faca --- /dev/null +++ b/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_TDHF_screening.inp @@ -0,0 +1,83 @@ +&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 + SELF_CONSISTENCY EVGW0 + &BSE + BSE_DIAG_METHOD FULLDIAG + ENERGY_CUTOFF_EMPTY 50.0 + ENERGY_CUTOFF_OCC 60.0 + TDA OFF + &SCREENING_IN_W TDHF + &END SCREENING_IN_W + &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 diff --git a/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_custom_screening.inp b/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_custom_screening.inp new file mode 100644 index 0000000000..36dee74748 --- /dev/null +++ b/tests/QS/regtest-bse/BSE_H2O_PBE_evGW0_custom_screening.inp @@ -0,0 +1,84 @@ +&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 + SELF_CONSISTENCY EVGW0 + &BSE + BSE_DIAG_METHOD FULLDIAG + ENERGY_CUTOFF_EMPTY 50.0 + ENERGY_CUTOFF_OCC 60.0 + TDA OFF + &SCREENING_IN_W ALPHA + ALPHA 0.25 + &END SCREENING_IN_W + &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 diff --git a/tests/QS/regtest-bse/TEST_FILES b/tests/QS/regtest-bse/TEST_FILES index 2ccc430922..5dca13db46 100644 --- a/tests/QS/regtest-bse/TEST_FILES +++ b/tests/QS/regtest-bse/TEST_FILES @@ -14,11 +14,15 @@ TDA_H2O_PBE_G0W0_NTOs.inp 121 2e-03 BSE_H2O_PBE_G0W0_NTOs.inp 122 2e-03 0.01164 TDA_H2O_PBE_G0W0_EXC_DESCR.inp 123 5e-04 1.8290 BSE_H2O_PBE_G0W0_EXC_DESCR.inp 124 5e-04 1.8258 +BSE_H2O_PBE_evGW0_RPA_screening.inp 110 2e-05 17.5569 +BSE_H2O_PBE_evGW0_TDHF_screening.inp 110 2e-05 6.4895 +BSE_H2O_PBE_evGW0_custom_screening.inp 110 2e-05 14.9538 #EOF #Low convergence criteria for evGW make a higher threshold necessary for these tests #Logic: Absolute errors (for G0W0/evGW0) should be within 1e-4, i.e. relative errors are # in the range of 2e-5 for energies. # For oscillator strengths, they should be within 1e-3, so relative errors # are in the range of 7e-2 -# First 6 check excitation energies, then 6 checking oscillator strengths, two check NTOs -# and the last one checks exciton descriptors in TDA +# First 6 check excitation energies, then 6 checking oscillator strengths, two check NTOs, +# two checks for exciton descriptors in TDA and ABBA and custom screening for +# RPA, TDHF and ALPHA=0.25 is checked in the last three tests