BSE: Minor changes for spectra

- add computation and printing of photoabsorption cross section
- improve handling of BSE calculations where PERIODIC XYZ is invoked

docs/methods/properties/optical/bethe-salpeter.md

@@ -84,7 +84,7 @@ $\alpha_{\mu,\mu'}(\omega) $ with $(\mu,\mu'\in\{x,y,z\})$:
 $$
 \begin{align}
 \alpha_{\mu,\mu'}(\omega) 
-= \sum_n \frac{2 \Omega^{(n)} d^{(n)}_{\mu} d^{(n)}_{\mu'}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2}
+= - \sum_n \frac{2 \Omega^{(n)} d^{(n)}_{\mu} d^{(n)}_{\mu'}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2}
 \quad ,
 \end{align}
 $$
@@ -114,7 +114,7 @@ We can rewrite the last equation as
 $$
 \begin{align}
 \mathrm{Im}\left[\bar{\alpha}(\omega)\right] 
-= \mathrm{Im}\left[
+= - \mathrm{Im}\left[
   \sum_n \frac{f^{(n)}}{(\omega+i\eta)^2-\left(\Omega^{(n)}\right)^2}
   \right]
 \quad .
@@ -130,6 +130,17 @@ f^{(n)} = \frac{2}{3} \Omega^{(n)} \sum_{\mu\in\{x,y,z\}} | d^{(n)}_{\mu} |^2
 \end{align}
 $$
 
+Additionally, the photoabsorption cross section tensor
+
+$$
+\begin{align}
+\sigma_{\mu,\mu'}(\omega)  = \frac{4 \pi \omega}{c} \mathrm{Im}\left[\alpha_{\mu,\mu'}(\omega) \right]
+\quad .
+\end{align}
+$$
+
+is printed, where $c$ denotes the speed of light.
+
 ### 1.3 Visualizing excited states using Natural Transition Orbitals (NTOs)
 
 In order to analyse the excitation wave function independent of a specific choice of the molecular
@@ -353,11 +364,12 @@ In the upper GW/BSE section, the following keywords have been used:
 
 - [BSE_SPECTRUM](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM): Activates
   computation and printing of the optical absorption spectrum. For each chosen option in
-  [TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA), a file is printed with
+  [TDA](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.TDA), one file is printed with
   columns frequency $\omega$, $\mathrm{Im}\left[\bar{\alpha}(\omega)\right]$ and the imaginary part
   of the elements of the dynamical dipole polarizability tensor
-  $\mathrm{Im}\left[{\alpha_{\mu,\mu'}}(\omega)\right]$. The frequency range, step size and the
-  broadening $\eta$ can be specified by the user (cf. keywords in
+  $\mathrm{Im}\left[{\alpha_{\mu,\mu'}}(\omega)\right]$ as well as another file with the respective
+  entries for the photoabsorption cross section tensor ${\sigma_{\mu,\mu'}}(\omega)$. The frequency
+  range, step size and the broadening $\eta$ can be specified by the user (cf. keywords in
   [BSE_SPECTRUM](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM)). Further,
   multiple broadenings $\eta$ can be given for one cp2k run
   ([ETA_LIST](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.BSE_SPECTRUM.ETA_LIST)),

src/bse_main.F

@@ -171,7 +171,7 @@ SUBROUTINE start_bse_calculation(fm_mat_S_ia_bse, fm_mat_S_ij_bse, fm_mat_S_ab_b
                                    gd_array, color_sub, para_env)
       END IF
 
-      CALL adapt_BSE_input_params(homo_red, virtual_red, unit_nr, mp2_env)
+      CALL adapt_BSE_input_params(homo_red, virtual_red, unit_nr, mp2_env, qs_env)
 
       IF (my_do_fulldiag) THEN
          ! Quick estimate of memory consumption and runtime of diagonalizations

src/bse_properties.F

@@ -38,11 +38,13 @@ MODULE bse_properties
                                               section_vals_type,&
                                               section_vals_val_get
    USE kinds,                           ONLY: dp
+   USE mathconstants,                   ONLY: pi
    USE mp2_types,                       ONLY: mp2_type
    USE parallel_gemm_api,               ONLY: parallel_gemm
    USE particle_methods,                ONLY: write_qs_particle_coordinates
    USE particle_types,                  ONLY: particle_type
-   USE physcon,                         ONLY: evolt
+   USE physcon,                         ONLY: c_light_au,&
+                                              evolt
    USE qs_environment_types,            ONLY: get_qs_env,&
                                               qs_environment_type
    USE qs_kind_types,                   ONLY: qs_kind_type
@@ -233,12 +235,13 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_
       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_absorption_spectrum'
 
       CHARACTER(LEN=10)                                  :: eta_str, width_eta_format_str
-      CHARACTER(LEN=30)                                  :: file_name_spectrum
+      CHARACTER(LEN=40)                                  :: file_name_crosssection, &
+                                                            file_name_spectrum
       INTEGER                                            :: handle, i, idir, j, jdir, k, num_steps, &
                                                             unit_nr_file, width_eta
       REAL(KIND=dp)                                      :: eta, freq_end, freq_start, freq_step, &
                                                             omega
-      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: abs_spectrum
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: abs_cross_section, abs_spectrum
       REAL(KIND=dp), DIMENSION(:), POINTER               :: eta_list
       TYPE(cp_logger_type), POINTER                      :: logger
 
@@ -261,21 +264,29 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_
          WRITE (eta_str, width_eta_format_str) eta*evolt
          ! Filename itself
          file_name_spectrum = 'BSE'//TRIM(ADJUSTL(info_approximation))//'eta='//TRIM(eta_str)//'.spectrum'
+         file_name_crosssection = 'BSE'//TRIM(ADJUSTL(info_approximation))//'eta='//TRIM(eta_str)//'.crosssection'
 
          ! First column is frequency in eV, second column is imaginary part of the trace of the polarizability
          ! The following 9 columns are the entries of the polarizability tensor
          ALLOCATE (abs_spectrum(num_steps, 11))
          abs_spectrum(:, :) = 0.0_dp
+         ! Also calculate and print the photoabsorption cross section tensor
+         ! σ_{µ µ'}(ω) = 4πω Im{α_{µ µ'}(ω)} / c
+         ALLOCATE (abs_cross_section(num_steps, 11))
+         abs_cross_section(:, :) = 0.0_dp
 
          ! Calculate the imaginary part of the mean dipole polarizability α_{avg}(ω)
          ! which is given by (cf. C. Ullrichs Book on TDDFT, Eq. 7.51)
-         ! α_{avg}(ω) = \sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}
+         ! We introduce an additional - due to his convention for charge vs particle density, see also:
+         ! Computer Physics Communications, 208:149–161, November 2016
+         ! https://doi.org/10.1016/j.cpc.2016.06.019
+         ! α_{avg}(ω) = - \sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}
          ! and then the imaginary part is (in the limit η -> 0)
-         ! Im[α_{avg}(ω)] = \sum_{n=1}^{N_exc} f_n * η / ((ω - Ω^n)² + η²)
+         ! Im[α_{avg}(ω)] = - \sum_{n=1}^{N_exc} f_n * η / ((ω - Ω^n)² + η²)
          ! 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]
+         ! α_{µ µ'}(ω) = - \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
@@ -294,6 +305,25 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_
             END DO
          END DO
 
+         ! Extract cross section σ from polarizability tensor
+         DO i = 1, num_steps
+            omega = abs_spectrum(i, 1)
+            abs_cross_section(i, 1) = omega
+            abs_cross_section(i, 2:) = 4.0_dp*pi*abs_spectrum(i, 2:)*omega/c_light_au
+         END DO
+
+         !For debug runs: Export an entry of the two tensors to allow regtests on spectra
+         IF (mp2_env%bse%bse_debug_print) THEN
+            IF (unit_nr > 0) THEN
+               WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
+                  'Averaged dynamical dipole polarizability at 8.2 eV:', &
+                  abs_spectrum(83, 2)
+               WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
+                  'Averaged photoabsorption cross section at 8.2 eV:', &
+                  abs_cross_section(83, 2)
+            END IF
+         END IF
+
          ! Print it to file
          logger => cp_get_default_logger()
          IF (logger%para_env%is_source()) THEN
@@ -302,10 +332,26 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_
             unit_nr_file = -1
          END IF
 
+         IF (unit_nr_file > 0) THEN
+            CALL open_file(file_name_crosssection, unit_number=unit_nr_file, &
+                           file_status="UNKNOWN", file_action="WRITE")
+            WRITE (unit_nr_file, '(A,A6)') "# Photoabsorption cross section  σ_{µ µ'}(ω) =  -4πω/c * Im[ \sum_n "// &
+               "[2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] ] from Bethe Salpeter equation for method ", &
+               TRIM(ADJUSTL(info_approximation))
+            WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "#     Frequency (eV)", "σ_{avg}(ω)", "σ_xx(ω)", &
+               "σ_xy(ω)", "σ_xz(ω)", "σ_yx(ω)", "σ_yy(ω)", "σ_yz(ω)", "σ_zx(ω)", &
+               "σ_zy(ω)", "σ_zz(ω)"
+            DO i = 1, num_steps
+               WRITE (unit_nr_file, '(11(F20.8,1X))') abs_cross_section(i, 1)*evolt, abs_cross_section(i, 2:11)
+            END DO
+            CALL close_file(unit_nr_file)
+         END IF
+         DEALLOCATE (abs_cross_section)
+
          IF (unit_nr_file > 0) THEN
             CALL open_file(file_name_spectrum, unit_number=unit_nr_file, &
                            file_status="UNKNOWN", file_action="WRITE")
-            WRITE (unit_nr_file, '(A,A6)') "# Imaginary part of polarizability α_{µ µ'}(ω) =  \sum_n "// &
+            WRITE (unit_nr_file, '(A,A6)') "# Imaginary part of polarizability α_{µ µ'}(ω) = -\sum_n "// &
                "[2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] from Bethe Salpeter equation for method ", &
                TRIM(ADJUSTL(info_approximation))
             WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "#     Frequency (eV)", "Im{α_{avg}(ω)}", "Im{α_xx(ω)}", &
@@ -323,21 +369,35 @@ SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_
          WRITE (unit_nr, '(T2,A4)') 'BSE|'
          WRITE (unit_nr, '(T2,A4,T7,A50,A)') &
             'BSE|', "Printed optical absorption spectrum to local file ", file_name_spectrum
-         WRITE (unit_nr, '(T2,A4,T7,A50)') &
-            'BSE|', "using the Eq. 7.51 from C. Ullrichs Book on TDDFT:"
+         WRITE (unit_nr, '(T2,A4,T7,A44,A)') &
+            'BSE|', "as well as photoabsorption cross section to ", file_name_crosssection
+         WRITE (unit_nr, '(T2,A4,T7,A52)') &
+            'BSE|', "using the Eq. (7.51) from C. Ullrichs Book on TDDFT:"
          WRITE (unit_nr, '(T2,A4)') 'BSE|'
          WRITE (unit_nr, '(T2,A4,T10,A73)') &
-            'BSE|', "Im{α_{avg}(ω)} = Im{\sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}}"
+            'BSE|', "Im{α_{avg}(ω)} = -Im{\sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}}"
          WRITE (unit_nr, '(T2,A4)') 'BSE|'
          WRITE (unit_nr, '(T2,A4,T7,A)') &
             'BSE|', "or for the full polarizability tensor:"
          WRITE (unit_nr, '(T2,A4,T10,A)') &
-            'BSE|', "α_{µ µ'}(ω) =  \sum_n [2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²]"
+            'BSE|', "α_{µ µ'}(ω) =  -\sum_n [2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²]"
+         WRITE (unit_nr, '(T2,A4)') 'BSE|'
+         WRITE (unit_nr, '(T2,A4,T7,A)') &
+            'BSE|', "as well as Eq. (7.48):"
+         WRITE (unit_nr, '(T2,A4,T10,A)') &
+            'BSE|', "σ_{µ µ'}(ω) = 4πω Im{α_{µ µ'}(ω)} / c"
          WRITE (unit_nr, '(T2,A4)') 'BSE|'
          WRITE (unit_nr, '(T2,A4,T7,A)') &
             'BSE|', "with transition moments d_µ^n, oscillator strengths f_n,"
          WRITE (unit_nr, '(T2,A4,T7,A)') &
-            'BSE|', "and excitation energies Ω^n."
+            'BSE|', "excitation energies Ω^n and the speed of light c."
+         WRITE (unit_nr, '(T2,A4)') 'BSE|'
+         WRITE (unit_nr, '(T2,A4,T7,A)') &
+            'BSE|', "Please note that we adopt an additional minus sign for both quantities,"
+         WRITE (unit_nr, '(T2,A4,T7,A)') &
+            'BSE|', "due to the convention for charge vs particle density as done in MolGW:"
+         WRITE (unit_nr, '(T2,A4,T7,A)') &
+            'BSE|', "https://doi.org/10.1016/j.cpc.2016.06.019."
          WRITE (unit_nr, '(T2,A4)') 'BSE|'
       END IF
 

src/bse_util.F

@@ -62,6 +62,7 @@ MODULE bse_util
    USE physcon,                         ONLY: evolt
    USE pw_env_types,                    ONLY: pw_env_get,&
                                               pw_env_type
+   USE pw_poisson_types,                ONLY: pw_poisson_type
    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
                                               pw_pool_type
    USE pw_types,                        ONLY: pw_c1d_gs_type,&
@@ -1506,23 +1507,35 @@ SUBROUTINE print_bse_nto_cubes(qs_env, mos, istate, info_approximation, &
    END SUBROUTINE print_bse_nto_cubes
 
 ! **************************************************************************************************
-!> \brief Borrowed from the tddfpt module with slight adaptions
+!> \brief Checks BSE input section and adapts them if necessary
 !> \param homo ...
 !> \param virtual ...
 !> \param unit_nr ...
 !> \param mp2_env ...
+!> \param qs_env ...
 ! **************************************************************************************************
-   SUBROUTINE adapt_BSE_input_params(homo, virtual, unit_nr, mp2_env)
+   SUBROUTINE adapt_BSE_input_params(homo, virtual, unit_nr, mp2_env, qs_env)
 
       INTEGER, INTENT(IN)                                :: homo, virtual, unit_nr
       TYPE(mp2_type)                                     :: mp2_env
+      TYPE(qs_environment_type), POINTER                 :: qs_env
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'adapt_BSE_input_params'
 
-      INTEGER                                            :: handle, i, j, n, &
+      INTEGER                                            :: handle, i, j, n, ndim_periodic_cell, &
+                                                            ndim_periodic_poisson, &
                                                             num_state_list_exceptions
+      TYPE(cell_type), POINTER                           :: cell_ref
+      TYPE(pw_env_type), POINTER                         :: pw_env
+      TYPE(pw_poisson_type), POINTER                     :: poisson_env
 
       CALL timeset(routineN, handle)
+      ! Get environment infos for later usage
+      NULLIFY (pw_env, cell_ref, poisson_env)
+      CALL get_qs_env(qs_env, pw_env=pw_env, cell_ref=cell_ref)
+      CALL pw_env_get(pw_env, poisson_env=poisson_env)
+      ndim_periodic_poisson = COUNT(poisson_env%parameters%periodic == 1)
+      ndim_periodic_cell = SUM(cell_ref%perd(1:3)) ! Borrowed from cell_methods.F/write_cell_low
 
       ! Handle negative NUM_PRINT_EXC
       IF (mp2_env%bse%num_print_exc < 0 .OR. &
@@ -1649,6 +1662,33 @@ SUBROUTINE adapt_BSE_input_params(homo, virtual, unit_nr, mp2_env)
          mp2_env%BSE%screening_method = bse_screening_tdhf
       END IF
 
+      ! Add warning for usage of KS energies
+      IF (mp2_env%bse%use_ks_energies) THEN
+         IF (unit_nr > 0) THEN
+            CALL cp_warn(__LOCATION__, &
+                         "Using KS energies for BSE calculations. Therefore, no quantities "// &
+                         "of the preceeding GW calculation enter the BSE.")
+         END IF
+      END IF
+
+      ! Add warning if periodic calculation is invoked
+      IF (ndim_periodic_poisson /= 0) THEN
+         IF (unit_nr > 0) THEN
+            CALL cp_warn(__LOCATION__, &
+                         "Poisson solver should be invoked by PERIODIC NONE. "// &
+                         "The applied length gauge might give misleading results for "// &
+                         "oscillator strengths.")
+         END IF
+      END IF
+      IF (ndim_periodic_cell /= 0) THEN
+         IF (unit_nr > 0) THEN
+            CALL cp_warn(__LOCATION__, &
+                         "CELL in SUBSYS should be invoked with PERIODIC NONE. "// &
+                         "The applied length gauge might give misleading results for "// &
+                         "oscillator strengths.")
+         END IF
+      END IF
+
       CALL timestop(handle)
    END SUBROUTINE adapt_BSE_input_params
 

src/mp2_integrals.F

@@ -335,7 +335,12 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
          CPASSERT(my_do_gw)
          CPASSERT(.NOT. do_im_time)
          ! GPW integrals have to be implemented later
-         CPASSERT(.NOT. (eri_method == do_eri_gpw))
+         IF (eri_method == do_eri_gpw) THEN
+            CALL cp_abort(__LOCATION__, &
+                          "BSE calculations are not implemented for GPW integrals. "// &
+                          "This is probably caused by invoking a periodic calculation. "// &
+                          "Use PERIODIC NONE for BSE calculations.")
+         END IF
       END IF
 
       ngroup = para_env%num_pe/para_env_sub%num_pe

src/rpa_main.F

@@ -1749,8 +1749,7 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
          IF (mp2_env%ri_g0w0%iter_evGW > 1) THEN
             IF (unit_nr > 0) THEN
                CALL cp_warn(__LOCATION__, &
-                            "BSE@evGW applies W0, i.e. screening with DFT energies! "// &
-                            "Use BSE@evGW0 to avoid this warning!")
+                            "BSE@evGW applies W0, i.e. screening with DFT energies to the BSE!")
             END IF
          END IF
          ! Create a copy of fm_mat_S for usage in BSE