LCOV - code coverage report
Current view: top level - src - bse_print.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:2fce0f8) Lines: 332 407 81.6 %
Date: 2024-12-21 06:28:57 Functions: 8 9 88.9 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines for printing information in context of the BSE calculation
      10             : !> \par History
      11             : !>      10.2024 created [Maximilian Graml]
      12             : ! **************************************************************************************************
      13             : MODULE bse_print
      14             : 
      15             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      16             :    USE bse_properties,                  ONLY: compute_absorption_spectrum,&
      17             :                                               exciton_descr_type
      18             :    USE bse_util,                        ONLY: filter_eigvec_contrib
      19             :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      20             :                                               cp_fm_type
      21             :    USE input_constants,                 ONLY: bse_screening_alpha,&
      22             :                                               bse_screening_rpa,&
      23             :                                               bse_screening_tdhf,&
      24             :                                               bse_screening_w0
      25             :    USE kinds,                           ONLY: dp
      26             :    USE mp2_types,                       ONLY: mp2_type
      27             :    USE particle_types,                  ONLY: particle_type
      28             :    USE physcon,                         ONLY: angstrom,&
      29             :                                               evolt
      30             :    USE qs_environment_types,            ONLY: get_qs_env,&
      31             :                                               qs_environment_type
      32             : #include "./base/base_uses.f90"
      33             : 
      34             :    IMPLICIT NONE
      35             : 
      36             :    PRIVATE
      37             : 
      38             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_print'
      39             : 
      40             :    PUBLIC :: print_BSE_start_flag, fm_write_thresh, print_excitation_energies, &
      41             :              print_output_header, print_transition_amplitudes, print_optical_properties, &
      42             :              print_exciton_descriptors
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! **************************************************************************************************
      47             : !> \brief ...
      48             : !> \param bse_tda ...
      49             : !> \param bse_abba ...
      50             : !> \param unit_nr ...
      51             : ! **************************************************************************************************
      52          42 :    SUBROUTINE print_BSE_start_flag(bse_tda, bse_abba, unit_nr)
      53             : 
      54             :       LOGICAL, INTENT(IN)                                :: bse_tda, bse_abba
      55             :       INTEGER, INTENT(IN)                                :: unit_nr
      56             : 
      57             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_BSE_start_flag'
      58             : 
      59             :       INTEGER                                            :: handle
      60             : 
      61          42 :       CALL timeset(routineN, handle)
      62             : 
      63          42 :       IF (unit_nr > 0) THEN
      64          21 :          WRITE (unit_nr, *) ' '
      65          21 :          WRITE (unit_nr, '(T2,A79)') '*******************************************************************************'
      66          21 :          WRITE (unit_nr, '(T2,A79)') '**                                                                           **'
      67          21 :          WRITE (unit_nr, '(T2,A79)') '**           Bethe Salpeter equation (BSE) for excitation energies           **'
      68          21 :          IF (bse_tda .AND. bse_abba) THEN
      69           0 :             WRITE (unit_nr, '(T2,A79)') '**          solved with and without Tamm-Dancoff approximation (TDA)         **'
      70          21 :          ELSE IF (bse_tda) THEN
      71           8 :             WRITE (unit_nr, '(T2,A79)') '**                solved with Tamm-Dancoff approximation (TDA)               **'
      72             :          ELSE
      73          13 :             WRITE (unit_nr, '(T2,A79)') '**               solved without Tamm-Dancoff approximation (TDA)             **'
      74             :          END IF
      75             : 
      76          21 :          WRITE (unit_nr, '(T2,A79)') '**                                                                           **'
      77          21 :          WRITE (unit_nr, '(T2,A79)') '*******************************************************************************'
      78          21 :          WRITE (unit_nr, *) ' '
      79             :       END IF
      80             : 
      81          42 :       CALL timestop(handle)
      82             : 
      83          42 :    END SUBROUTINE
      84             : 
      85             : ! **************************************************************************************************
      86             : !> \brief ...
      87             : !> \param homo ...
      88             : !> \param virtual ...
      89             : !> \param homo_irred ...
      90             : !> \param flag_TDA ...
      91             : !> \param multiplet ...
      92             : !> \param alpha ...
      93             : !> \param mp2_env ...
      94             : !> \param unit_nr ...
      95             : ! **************************************************************************************************
      96          42 :    SUBROUTINE print_output_header(homo, virtual, homo_irred, flag_TDA, &
      97             :                                   multiplet, alpha, mp2_env, unit_nr)
      98             : 
      99             :       INTEGER, INTENT(IN)                                :: homo, virtual, homo_irred
     100             :       LOGICAL, INTENT(IN)                                :: flag_TDA
     101             :       CHARACTER(LEN=10), INTENT(IN)                      :: multiplet
     102             :       REAL(KIND=dp), INTENT(IN)                          :: alpha
     103             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     104             :       INTEGER, INTENT(IN)                                :: unit_nr
     105             : 
     106             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_output_header'
     107             : 
     108             :       INTEGER                                            :: handle
     109             : 
     110          42 :       CALL timeset(routineN, handle)
     111             : 
     112          42 :       IF (unit_nr > 0) THEN
     113          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     114          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     115          21 :          IF (flag_TDA) THEN
     116           8 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     117           8 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '*   Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA)  *'
     118           8 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     119           8 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     120           8 :             WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
     121          16 :                'the BSE within the TDA:'
     122           8 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     123           8 :             WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n'
     124           8 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     125           8 :             WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
     126           8 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     127           8 :             WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb   X_jb^n ) = Ω^n X_ia^n'
     128           8 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     129           8 :             WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
     130           8 :             WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
     131             :          ELSE
     132          13 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     133          13 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '*      Full ("ABBA") Bethe Salpeter equation (BSE) (i.e. without TDA)    *'
     134          13 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     135          13 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     136          13 :             WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
     137          26 :                'the BSE without the TDA:'
     138          13 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     139          13 :             WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n|       |1  0| |X^n|'
     140          13 :             WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|'
     141          13 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     142          13 :             WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
     143          13 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     144          13 :             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'
     145          13 :             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'
     146             :          END IF
     147          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     148          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     149          21 :          WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
     150          42 :             'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
     151          21 :          WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
     152          42 :             'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
     153          21 :          WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
     154          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     155          21 :          IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
     156          18 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb - W_ij,ab'
     157           3 :          ELSE IF (mp2_env%bse%screening_method == bse_screening_rpa) THEN
     158           1 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'A_ia,jb = (ε_a-ε_i) δ_ij δ_ab + α * v_ia,jb'
     159             :          END IF
     160          21 :          IF (.NOT. flag_TDA) THEN
     161          13 :             IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
     162          10 :                WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'B_ia,jb = α * v_ia,jb - W_ib,aj'
     163           3 :             ELSE IF (mp2_env%bse%screening_method == bse_screening_rpa) THEN
     164           1 :                WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', 'B_ia,jb = α * v_ia,jb'
     165             :             END IF
     166          13 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     167          13 :             WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
     168          13 :             WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
     169             :          END IF
     170          21 :          IF (.NOT. flag_TDA) THEN
     171          13 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     172          13 :             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)'
     173          13 :             WRITE (unit_nr, '(T2,A4,T7,A71)') 'BSE|', 'in PRB 92,045209 (2015); http://dx.doi.org/10.1103/PhysRevB.92.045209 .'
     174             :          END IF
     175          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     176          21 :          WRITE (unit_nr, '(T2,A4,T7,A7,T31,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
     177          21 :          WRITE (unit_nr, '(T2,A4,T7,A7,T31,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
     178          21 :          WRITE (unit_nr, '(T2,A4,T7,A3,T31,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
     179          21 :          WRITE (unit_nr, '(T2,A4,T7,A6,T30,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
     180          21 :          IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
     181          18 :             WRITE (unit_nr, '(T2,A4,T7,A,T31,A)') 'BSE|', 'W_... = 1/ϵ v_...:', &
     182          36 :                'Direct interaction screened by '
     183          18 :             WRITE (unit_nr, '(T2,A4,T30,A)') 'BSE|', &
     184          36 :                'dielectric function ϵ(ω=0)'
     185           3 :          ELSE IF (mp2_env%bse%screening_method == bse_screening_tdhf) THEN
     186           1 :             WRITE (unit_nr, '(T2,A4,T7,A,T30,A)') 'BSE|', 'W_... = v_...:', 'Direct interaction without screening'
     187           2 :          ELSE IF (mp2_env%bse%screening_method == bse_screening_alpha) THEN
     188           1 :             WRITE (unit_nr, '(T2,A4,T7,A,T31,A,F5.2)') 'BSE|', 'W_... = γ v_...:', &
     189           2 :                'Direct interaction with artificial screening γ=', mp2_env%bse%screening_factor
     190             :          END IF
     191          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     192          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     193          21 :          WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
     194          42 :             'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha
     195          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     196             :       END IF
     197             : 
     198          42 :       CALL timestop(handle)
     199             : 
     200          42 :    END SUBROUTINE
     201             : 
     202             : ! **************************************************************************************************
     203             : !> \brief ...
     204             : !> \param Exc_ens ...
     205             : !> \param homo ...
     206             : !> \param virtual ...
     207             : !> \param flag_TDA ...
     208             : !> \param multiplet ...
     209             : !> \param info_approximation ...
     210             : !> \param mp2_env ...
     211             : !> \param unit_nr ...
     212             : ! **************************************************************************************************
     213          42 :    SUBROUTINE print_excitation_energies(Exc_ens, homo, virtual, flag_TDA, multiplet, &
     214             :                                         info_approximation, mp2_env, unit_nr)
     215             : 
     216             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Exc_ens
     217             :       INTEGER, INTENT(IN)                                :: homo, virtual
     218             :       LOGICAL, INTENT(IN)                                :: flag_TDA
     219             :       CHARACTER(LEN=10), INTENT(IN)                      :: multiplet, info_approximation
     220             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     221             :       INTEGER, INTENT(IN)                                :: unit_nr
     222             : 
     223             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_excitation_energies'
     224             : 
     225             :       INTEGER                                            :: handle, i_exc
     226             : 
     227          42 :       CALL timeset(routineN, handle)
     228             : 
     229          42 :       IF (unit_nr > 0) THEN
     230          21 :          IF (flag_TDA) THEN
     231           8 :             WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
     232             :          ELSE
     233          13 :             WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
     234             :          END IF
     235          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     236          21 :          WRITE (unit_nr, '(T2,A4,T11,A12,T26,A11,T44,A8,T55,A27)') 'BSE|', &
     237          42 :             'Excitation n', "Spin Config", 'TDA/ABBA', 'Excitation energy Ω^n (eV)'
     238             :       END IF
     239             :       !prints actual energies values
     240          42 :       IF (unit_nr > 0) THEN
     241         546 :          DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
     242             :             WRITE (unit_nr, '(T2,A4,T7,I16,T30,A7,T46,A6,T59,F22.4)') &
     243         546 :                'BSE|', i_exc, multiplet, info_approximation, Exc_ens(i_exc)*evolt
     244             :          END DO
     245             :       END IF
     246             : 
     247          42 :       CALL timestop(handle)
     248             : 
     249          42 :    END SUBROUTINE print_excitation_energies
     250             : 
     251             : ! **************************************************************************************************
     252             : !> \brief ...
     253             : !> \param fm_eigvec_X ...
     254             : !> \param homo ...
     255             : !> \param virtual ...
     256             : !> \param homo_irred ...
     257             : !> \param info_approximation ...
     258             : !> \param mp2_env ...
     259             : !> \param unit_nr ...
     260             : !> \param fm_eigvec_Y ...
     261             : ! **************************************************************************************************
     262          42 :    SUBROUTINE print_transition_amplitudes(fm_eigvec_X, homo, virtual, homo_irred, &
     263             :                                           info_approximation, mp2_env, unit_nr, fm_eigvec_Y)
     264             : 
     265             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_eigvec_X
     266             :       INTEGER, INTENT(IN)                                :: homo, virtual, homo_irred
     267             :       CHARACTER(LEN=10), INTENT(IN)                      :: info_approximation
     268             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     269             :       INTEGER, INTENT(IN)                                :: unit_nr
     270             :       TYPE(cp_fm_type), INTENT(IN), OPTIONAL             :: fm_eigvec_Y
     271             : 
     272             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_transition_amplitudes'
     273             : 
     274             :       INTEGER                                            :: handle, i_exc
     275             : 
     276          42 :       CALL timeset(routineN, handle)
     277             : 
     278          42 :       IF (unit_nr > 0) THEN
     279          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     280             :          WRITE (unit_nr, '(T2,A4,T7,A61)') &
     281          21 :             'BSE|', "Single-particle transitions are built up by (de-)excitations,"
     282             :          WRITE (unit_nr, '(T2,A4,T7,A18)') &
     283          21 :             'BSE|', "which we denote by"
     284          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     285             :          WRITE (unit_nr, '(T2,A4,T20,A2,T30,A40)') &
     286          21 :             'BSE|', "=>", "for excitations, i.e. entries of X_ia^n,"
     287             :          WRITE (unit_nr, '(T2,A4,T20,A2,T30,A42)') &
     288          21 :             'BSE|', "<=", "for deexcitations, i.e. entries of Y_ia^n."
     289             :          WRITE (unit_nr, '(T2,A4)') &
     290          21 :             'BSE|'
     291             :          WRITE (unit_nr, '(T2,A4,T7,A73)') &
     292          21 :             'BSE|', "The following single-particle transitions have significant contributions,"
     293             :          WRITE (unit_nr, '(T2,A4,T7,A16,F5.3,A15,F5.3,A16)') &
     294          21 :             'BSE|', "i.e. |X_ia^n| > ", mp2_env%bse%eps_x, " or |Y_ia^n| > ", &
     295          42 :             mp2_env%bse%eps_x, ", respectively :"
     296             : 
     297          21 :          WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
     298          42 :             homo_irred, ' and LUMO a =', homo_irred + 1, " --"
     299          21 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     300             :          WRITE (unit_nr, '(T2,A4,T7,A12,T30,A1,T32,A5,T42,A1,T49,A8,T64,A17)') &
     301          21 :             "BSE|", "Excitation n", "i", "=>/<=", "a", 'TDA/ABBA', "|X_ia^n|/|Y_ia^n|"
     302             :       END IF
     303        1092 :       DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
     304        1050 :          IF (unit_nr > 0) THEN
     305         525 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     306             :          END IF
     307             :          !Iterate through eigenvector and print values above threshold
     308             :          CALL print_transition_amplitudes_core(fm_eigvec_X, "=>", info_approximation, &
     309             :                                                i_exc, virtual, homo, homo_irred, &
     310        1050 :                                                unit_nr, mp2_env)
     311        1092 :          IF (PRESENT(fm_eigvec_Y)) THEN
     312             :             CALL print_transition_amplitudes_core(fm_eigvec_Y, "<=", info_approximation, &
     313             :                                                   i_exc, virtual, homo, homo_irred, &
     314         650 :                                                   unit_nr, mp2_env)
     315             :          END IF
     316             :       END DO
     317             : 
     318          42 :       CALL timestop(handle)
     319             : 
     320          42 :    END SUBROUTINE print_transition_amplitudes
     321             : 
     322             : ! **************************************************************************************************
     323             : !> \brief ...
     324             : !> \param Exc_ens ...
     325             : !> \param oscill_str ...
     326             : !> \param trans_mom_bse ...
     327             : !> \param polarizability_residues ...
     328             : !> \param homo ...
     329             : !> \param virtual ...
     330             : !> \param homo_irred ...
     331             : !> \param flag_TDA ...
     332             : !> \param info_approximation ...
     333             : !> \param mp2_env ...
     334             : !> \param unit_nr ...
     335             : ! **************************************************************************************************
     336          42 :    SUBROUTINE print_optical_properties(Exc_ens, oscill_str, trans_mom_bse, polarizability_residues, &
     337             :                                        homo, virtual, homo_irred, flag_TDA, &
     338             :                                        info_approximation, mp2_env, unit_nr)
     339             : 
     340             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: Exc_ens, oscill_str
     341             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: trans_mom_bse, polarizability_residues
     342             :       INTEGER, INTENT(IN)                                :: homo, virtual, homo_irred
     343             :       LOGICAL, INTENT(IN)                                :: flag_TDA
     344             :       CHARACTER(LEN=10), INTENT(IN)                      :: info_approximation
     345             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     346             :       INTEGER, INTENT(IN)                                :: unit_nr
     347             : 
     348             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_optical_properties'
     349             : 
     350             :       INTEGER                                            :: handle, i_exc
     351             : 
     352          42 :       CALL timeset(routineN, handle)
     353             : 
     354             :       ! Discriminate between singlet and triplet, since triplet state can't couple to light
     355             :       ! and therefore calculations of dipoles etc are not necessary
     356          42 :       IF (mp2_env%bse%bse_spin_config == 0) THEN
     357          42 :          IF (unit_nr > 0) THEN
     358          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     359             :             WRITE (unit_nr, '(T2,A4,T7,A60)') &
     360          21 :                'BSE|', "Transition moments d_r^n (with r∈(x,y,z), in atomic units)"
     361             :             WRITE (unit_nr, '(T2,A4,T7,A67)') &
     362          21 :                'BSE|', "and oscillator strength f^n of excitation level n are obtained from"
     363          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     364          21 :             IF (flag_TDA) THEN
     365             :                WRITE (unit_nr, '(T2,A4,T10,A)') &
     366           8 :                   'BSE|', "d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a >  X_ia^n"
     367             :             ELSE
     368             :                WRITE (unit_nr, '(T2,A4,T10,A)') &
     369          13 :                   'BSE|', "d_r^n = sum_ia sqrt(2) < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n )"
     370             :             END IF
     371          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     372             :             WRITE (unit_nr, '(T2,A4,T14,A)') &
     373          21 :                'BSE|', "f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2"
     374          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     375             :             WRITE (unit_nr, '(T2,A4,T7,A19)') &
     376          21 :                'BSE|', "where we introduced"
     377          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     378             :             WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
     379          21 :                'BSE|', "ψ_i:", "occupied molecular orbitals,"
     380             :             WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
     381          21 :                'BSE|', "ψ_a:", "empty molecular orbitals and"
     382             :             WRITE (unit_nr, '(T2,A4,T9,A2,T14,A18)') &
     383          21 :                'BSE|', "r:", "position operator."
     384          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     385             :             WRITE (unit_nr, '(T2,A4,T7,A28)') &
     386          21 :                'BSE|', "prelim Ref.: Eqs. (23), (24)"
     387             :             WRITE (unit_nr, '(T2,A4,T7,A71)') &
     388          21 :                'BSE|', "in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290"
     389          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     390          21 :             IF (flag_TDA) THEN
     391           8 :                WRITE (unit_nr, '(T2,A4,T7,A55)') 'BSE|', &
     392          16 :                   'Optical properties from solving the BSE within the TDA:'
     393             :             ELSE
     394          13 :                WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', &
     395          26 :                   'Optical properties from solving the BSE without the TDA:'
     396             :             END IF
     397          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     398          21 :             WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T38,A5,T48,A5,T58,A5,T64,A17)') 'BSE|', &
     399          42 :                'Excitation n', "TDA/ABBA", "d_x^n", "d_y^n", "d_z^n", 'Osc. strength f^n'
     400         546 :             DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
     401             :                WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T35,F8.3,T45,F8.3,T55,F8.3,T65,F16.3)') &
     402         525 :                   'BSE|', i_exc, info_approximation, trans_mom_bse(1, 1, i_exc), trans_mom_bse(2, 1, i_exc), &
     403        1071 :                   trans_mom_bse(3, 1, i_exc), oscill_str(i_exc)
     404             :             END DO
     405          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     406          21 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     407          42 :                'Check for Thomas-Reiche-Kuhn sum rule'
     408          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     409          21 :             WRITE (unit_nr, '(T2,A4,T35,A)') 'BSE|', &
     410          42 :                'N_e = Σ_n f^n'
     411          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     412          21 :             WRITE (unit_nr, '(T2,A4,T7,A24,T65,I16)') 'BSE|', &
     413          42 :                'Number of electrons N_e:', homo_irred*2
     414          21 :             WRITE (unit_nr, '(T2,A4,T7,A,T66,F16.3)') 'BSE|', &
     415        1050 :                'Sum over oscillator strengths Σ_n f^n :', SUM(oscill_str)
     416          21 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     417          21 :             IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN
     418             :                CALL cp_warn(__LOCATION__, &
     419          21 :                             "Accuracy of TRK sum rule might suffer from cutoffs.")
     420             :             END IF
     421             :          END IF
     422             : 
     423             :          ! Compute and print the absorption spectrum to external file
     424          42 :          IF (mp2_env%bse%bse_print_spectrum) THEN
     425             :             CALL compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_ens, &
     426           4 :                                              info_approximation, unit_nr, mp2_env)
     427             :          END IF
     428             : 
     429             :       ELSE
     430           0 :          IF (unit_nr > 0) THEN
     431           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     432           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     433             :             CALL cp_warn(__LOCATION__, &
     434             :                          "Requested triplet excitation cannot couple to light. "// &
     435             :                          "Skipping calculation of transition moments, "// &
     436           0 :                          "oscillator strengths, and spectrum.")
     437             :          END IF
     438             :       END IF
     439             : 
     440          42 :       CALL timestop(handle)
     441             : 
     442          42 :    END SUBROUTINE print_optical_properties
     443             : 
     444             : ! **************************************************************************************************
     445             : !> \brief ...
     446             : !> \param fm_eigvec ...
     447             : !> \param direction_excitation ...
     448             : !> \param info_approximation ...
     449             : !> \param i_exc ...
     450             : !> \param virtual ...
     451             : !> \param homo ...
     452             : !> \param homo_irred ...
     453             : !> \param unit_nr ...
     454             : !> \param mp2_env ...
     455             : ! **************************************************************************************************
     456        1700 :    SUBROUTINE print_transition_amplitudes_core(fm_eigvec, direction_excitation, info_approximation, &
     457             :                                                i_exc, virtual, homo, homo_irred, &
     458             :                                                unit_nr, mp2_env)
     459             : 
     460             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_eigvec
     461             :       CHARACTER(LEN=2), INTENT(IN)                       :: direction_excitation
     462             :       CHARACTER(LEN=10), INTENT(IN)                      :: info_approximation
     463             :       INTEGER                                            :: i_exc, virtual, homo, homo_irred
     464             :       INTEGER, INTENT(IN)                                :: unit_nr
     465             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     466             : 
     467             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_transition_amplitudes_core'
     468             : 
     469             :       INTEGER                                            :: handle, k, num_entries
     470        1700 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: idx_homo, idx_virt
     471        1700 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigvec_entries
     472             : 
     473        1700 :       CALL timeset(routineN, handle)
     474             : 
     475             :       CALL filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
     476        1700 :                                  i_exc, virtual, num_entries, mp2_env)
     477             :       ! direction_excitation can be either => (means excitation; from fm_eigvec_X)
     478             :       ! or <= (means deexcitation; from fm_eigvec_Y)
     479        1700 :       IF (unit_nr > 0) THEN
     480        2127 :          DO k = 1, num_entries
     481             :             WRITE (unit_nr, '(T2,A4,T14,I5,T26,I5,T35,A2,T38,I5,T51,A6,T65,F16.4)') &
     482        1277 :                "BSE|", i_exc, homo_irred - homo + idx_homo(k), direction_excitation, &
     483        3404 :                homo_irred + idx_virt(k), info_approximation, ABS(eigvec_entries(k))
     484             :          END DO
     485             :       END IF
     486        1700 :       DEALLOCATE (idx_homo)
     487        1700 :       DEALLOCATE (idx_virt)
     488        1700 :       DEALLOCATE (eigvec_entries)
     489        1700 :       CALL timestop(handle)
     490             : 
     491        1700 :    END SUBROUTINE
     492             : 
     493             : ! **************************************************************************************************
     494             : !> \brief ...
     495             : !> \param exc_descr ...
     496             : !> \param ref_point_multipole ...
     497             : !> \param unit_nr ...
     498             : !> \param mp2_env ...
     499             : !> \param qs_env ...
     500             : ! **************************************************************************************************
     501           4 :    SUBROUTINE print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, mp2_env, qs_env)
     502             : 
     503             :       TYPE(exciton_descr_type), ALLOCATABLE, &
     504             :          DIMENSION(:)                                    :: exc_descr
     505             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     506             :          INTENT(IN)                                      :: ref_point_multipole
     507             :       INTEGER, INTENT(IN)                                :: unit_nr
     508             :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     509             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     510             : 
     511             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_exciton_descriptors'
     512             : 
     513             :       CHARACTER(LEN=1), DIMENSION(3)                     :: array_direction_str
     514             :       INTEGER                                            :: handle, i_dir, i_exc
     515             :       REAL(KIND=dp)                                      :: d_eh_dir, d_exc_dir, sigma_e_dir, &
     516             :                                                             sigma_h_dir
     517           4 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     518             : 
     519           4 :       CALL timeset(routineN, handle)
     520           4 :       CALL get_qs_env(qs_env, particle_set=particle_set)
     521           4 :       IF (unit_nr > 0) THEN
     522           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     523           4 :             'Exciton descriptors for excitation level n are given by'
     524           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     525           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     526           4 :             'd_eh   = | <r_h - r_e>_exc |'
     527           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     528           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     529           4 :             'σ_e    = sqrt( <r_e^2>_exc - <r_e>_exc^2 )'
     530           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     531           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     532           4 :             'σ_h    = sqrt( <r_h^2>_exc - <r_h>_exc^2 )'
     533           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     534           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     535           4 :             'COV_eh = <r_e r_h>_exc - <r_e>_exc <r_h>_exc'
     536           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     537           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     538           4 :             'd_exc  = sqrt( | < |r_h - r_e|^2 >_exc )'
     539           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     540           4 :             '       = sqrt( d_eh^2 + σ_e^2 + σ_h^2 - 2 * COV_eh )'
     541           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     542           2 :          WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     543           4 :             'R_eh   = COV_eh / (σ_e * σ_h)'
     544           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     545           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     546           4 :             'where the expectation values <.>_exc are taken with respect to the '
     547           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     548           4 :             'exciton wavefunction  of excitation n:'
     549           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     550             : 
     551           2 :          IF (exc_descr(1)%flag_TDA) THEN
     552           1 :             WRITE (unit_nr, '(T2,A4,T20,A)') 'BSE|', &
     553           2 :                '𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e)  ,'
     554             :          ELSE
     555           1 :             WRITE (unit_nr, '(T2,A4,T20,A)') 'BSE|', &
     556           2 :                '𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e)'
     557           1 :             WRITE (unit_nr, '(T2,A4,T40,A)') 'BSE|', &
     558           2 :                '+ Y_ia^n ψ_a(r_h) ψ_i(r_e) ,'
     559             :          END IF
     560           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     561           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     562           4 :             'i.e.'
     563           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     564           2 :          WRITE (unit_nr, '(T2,A4,T20,A)') 'BSE|', &
     565           4 :             '< O >_exc = < 𝚿_n | O | 𝚿_n > / < 𝚿_n | 𝚿_n >  ,'
     566           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     567           2 :          IF (exc_descr(1)%flag_TDA) THEN
     568           1 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     569           2 :                'where c_n = < 𝚿_n | 𝚿_n > = 1 within TDA.'
     570             :          ELSE
     571           1 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     572           2 :                'where c_n = < 𝚿_n | 𝚿_n > deviates from 1 without TDA.'
     573             :          END IF
     574           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     575           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     576           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     577           4 :             'Here, we introduced'
     578           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     579             :          WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
     580           2 :             'BSE|', "ψ_i:", "occupied molecular orbitals,"
     581             :          WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
     582           2 :             'BSE|', "ψ_a:", "empty molecular orbitals and"
     583             :          WRITE (unit_nr, '(T2,A4,T9,A2,T14,A)') &
     584           2 :             'BSE|', "r:", "position operator."
     585           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     586           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     587           4 :             'prelim Ref.: Eqs. (15)-(22)'
     588           2 :          WRITE (unit_nr, '(T2,A4,T7,A,A)') 'BSE|', &
     589           2 :             'JCTC 2018, 14, 710-725; ', &
     590           4 :             'http://doi.org/10.1021/acs.jctc.7b01145'
     591           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     592           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     593           2 :          IF (exc_descr(1)%flag_TDA) THEN
     594           1 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     595           2 :                'Exciton descriptors from solving the BSE within the TDA:'
     596             :          ELSE
     597           1 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     598           2 :                'Exciton descriptors from solving the BSE without the TDA:'
     599             :          END IF
     600           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     601           2 :          WRITE (unit_nr, '(T2,A4,T10,A1,6X,A3,1X,4X,A10,5X,A10,5X,A10,3X,A11,8X,A4)') 'BSE|', &
     602           4 :             'n', 'c_n', 'd_eh [Å]', 'σ_e [Å]', 'σ_h [Å]', 'd_exc [Å]', 'R_eh'
     603          52 :          DO i_exc = 1, mp2_env%bse%num_print_exc_descr
     604             :             WRITE (unit_nr, '(T2,A4,T7,I4,4X,F5.3,1X,5(2X,F10.4))') &
     605          50 :                'BSE|', i_exc, exc_descr(i_exc)%norm_XpY, &
     606          50 :                exc_descr(i_exc)%diff_r_abs*angstrom, &
     607          50 :                exc_descr(i_exc)%sigma_e*angstrom, exc_descr(i_exc)%sigma_h*angstrom, &
     608         102 :                exc_descr(i_exc)%diff_r_sqr*angstrom, exc_descr(i_exc)%corr_e_h*angstrom
     609             :          END DO
     610           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     611             :          ! For debug runs, print first d_exc separately to allow the regtests to read in
     612           2 :          IF (mp2_env%bse%bse_debug_print) THEN
     613           2 :             IF (exc_descr(1)%flag_TDA) THEN
     614           1 :                WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
     615           2 :                   'Exciton descriptor d_exc with TDA for n=1 is', exc_descr(1)%diff_r_sqr*angstrom
     616             :             ELSE
     617           1 :                WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
     618           2 :                   'Exciton descriptor d_exc without TDA for n=1 is', exc_descr(1)%diff_r_sqr*angstrom
     619             :             END IF
     620           2 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     621             :          END IF
     622           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     623             :          ! Print exciton descriptor resolved per direction
     624           2 :          IF (mp2_env%bse%print_directional_exc_descr) THEN
     625           0 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     626           0 :                'We can restrict the exciton descriptors to a specific direction,'
     627           0 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     628           0 :                'e.g. the x-components are:'
     629           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     630           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     631           0 :                'd_eh^x   = | <x_h - x_e>_exc |'
     632           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     633           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     634           0 :                'σ_e^x    = sqrt( <x_e^2>_exc - <x_e>_exc^2 )'
     635           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     636           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     637           0 :                'σ_h^x    = sqrt( <x_h^2>_exc - <x_h>_exc^2 )'
     638           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     639           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     640           0 :                'COV_eh^x = <x_e x_h>_exc - <x_e>_exc <x_h>_exc'
     641           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     642           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     643           0 :                'd_exc^x  = sqrt( | < |x_h - x_e|^2 >_exc )'
     644           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     645           0 :                '         = sqrt( (d_eh^x)^2 + (σ_e^x)^2'
     646           0 :             WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
     647           0 :                '                 + (σ_h^x)^2 - 2 * (COV_eh^x) )'
     648           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     649           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     650           0 :             IF (exc_descr(1)%flag_TDA) THEN
     651           0 :                WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     652           0 :                   'Exciton descriptors per direction from solving the BSE within the TDA:'
     653             :             ELSE
     654           0 :                WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     655           0 :                   'Exciton descriptors per direction from solving the BSE without the TDA:'
     656             :             END IF
     657           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     658           0 :             WRITE (unit_nr, '(T2,A4,T12,A1,2X,A9,5X,A12,5X,A12,5X,A12,3X,A13)') 'BSE|', &
     659           0 :                'n', 'r = x/y/z', 'd_eh^r [Å]', 'σ_e^r [Å]', 'σ_h^r [Å]', 'd_exc^r [Å]'
     660           0 :             DO i_exc = 1, mp2_env%bse%num_print_exc_descr
     661           0 :                DO i_dir = 1, 3
     662           0 :                   array_direction_str = (/"x", "y", "z"/)
     663           0 :                   d_eh_dir = ABS(exc_descr(i_exc)%r_h(i_dir) - exc_descr(i_exc)%r_e(i_dir))
     664           0 :                   sigma_e_dir = SQRT(exc_descr(i_exc)%r_e_sq(i_dir) - exc_descr(i_exc)%r_e(i_dir)**2)
     665           0 :                   sigma_h_dir = SQRT(exc_descr(i_exc)%r_h_sq(i_dir) - exc_descr(i_exc)%r_h(i_dir)**2)
     666             :                   d_exc_dir = SQRT(d_eh_dir**2 + sigma_e_dir**2 + sigma_h_dir**2 &
     667           0 :                                    - 2*exc_descr(i_exc)%cov_e_h(i_dir))
     668             :                   WRITE (unit_nr, '(T2,A4,T9,I4,10X,A1,1X,4(4X,F10.4))') &
     669           0 :                      'BSE|', i_exc, array_direction_str(i_dir), &
     670           0 :                      d_eh_dir*angstrom, &
     671           0 :                      sigma_e_dir*angstrom, sigma_h_dir*angstrom, &
     672           0 :                      d_exc_dir*angstrom
     673             :                END DO
     674           0 :                WRITE (unit_nr, '(T2,A4)') 'BSE|'
     675             :             END DO
     676           0 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     677             :          END IF
     678             :          ! Print the reference atomic geometry for the exciton descriptors
     679           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     680           4 :             'With the center of charge as reference point r_0,'
     681           2 :          WRITE (unit_nr, '(T2,A4,T15,A7,F10.4,A2,F10.4,A2,F10.4,A1)') 'BSE|', &
     682           2 :             'r_0 = (', ref_point_multipole(1)*angstrom, ', ', ref_point_multipole(2)*angstrom, ', ', &
     683           4 :             ref_point_multipole(3)*angstrom, ')'
     684           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     685           4 :             'we further obtain r_e and r_h from solving the BSE within the TDA'
     686           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     687           2 :          WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') 'BSE|', &
     688           4 :             'Excitation n', 'x_e [Å]', 'y_e [Å]', 'z_e [Å]'
     689          52 :          DO i_exc = 1, mp2_env%bse%num_print_exc_descr
     690             :             WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
     691          50 :                'BSE|', i_exc, &
     692         252 :                exc_descr(i_exc)%r_e_shift(:)*angstrom
     693             :          END DO
     694           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     695           2 :          WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') 'BSE|', &
     696           4 :             'Excitation n', 'x_h [Å]', 'y_h [Å]', 'z_h [Å]'
     697          52 :          DO i_exc = 1, mp2_env%bse%num_print_exc_descr
     698             :             WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
     699          50 :                'BSE|', i_exc, &
     700         252 :                exc_descr(i_exc)%r_h_shift(:)*angstrom
     701             :          END DO
     702           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     703           2 :          WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     704           4 :             'The reference atomic geometry for these values is given by'
     705             :       END IF
     706           4 :       CALL write_qs_particle_coordinates_bse(particle_set, unit_nr)
     707           4 :       IF (unit_nr > 0) THEN
     708           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     709             :       END IF
     710           4 :       CALL timestop(handle)
     711             : 
     712           4 :    END SUBROUTINE print_exciton_descriptors
     713             : 
     714             : ! **************************************************************************************************
     715             : !> \brief Debug function to write elements of a full matrix to file, if they are larger than a given threshold
     716             : !> \param fm ...
     717             : !> \param thresh ...
     718             : !> \param header ...
     719             : !> \param unit_nr ...
     720             : !> \param abs_vals ...
     721             : ! **************************************************************************************************
     722           0 :    SUBROUTINE fm_write_thresh(fm, thresh, header, unit_nr, abs_vals)
     723             : 
     724             :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     725             :       REAL(KIND=dp), INTENT(IN)                          :: thresh
     726             :       CHARACTER(LEN=*), INTENT(IN)                       :: header
     727             :       INTEGER, INTENT(IN)                                :: unit_nr
     728             :       LOGICAL, OPTIONAL                                  :: abs_vals
     729             : 
     730             :       CHARACTER(LEN=*), PARAMETER :: my_footer = " | ENDING WRITING OF MATRIX", &
     731             :          routineN = 'fm_write_thresh'
     732             : 
     733             :       INTEGER                                            :: handle, i, j, ncol_local, nrow_local
     734           0 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     735             :       LOGICAL                                            :: my_abs_vals
     736             : 
     737           0 :       CALL timeset(routineN, handle)
     738             : 
     739           0 :       IF (PRESENT(abs_vals)) THEN
     740           0 :          my_abs_vals = abs_vals
     741             :       ELSE
     742             :          my_abs_vals = .FALSE.
     743             :       END IF
     744             : 
     745             :       CALL cp_fm_get_info(matrix=fm, &
     746             :                           nrow_local=nrow_local, &
     747             :                           ncol_local=ncol_local, &
     748             :                           row_indices=row_indices, &
     749           0 :                           col_indices=col_indices)
     750             : 
     751           0 :       IF (unit_nr > 0) THEN
     752           0 :          WRITE (unit_nr, *) header
     753             :       END IF
     754           0 :       IF (my_abs_vals) THEN
     755           0 :          DO i = 1, nrow_local
     756           0 :             DO j = 1, ncol_local
     757           0 :                IF (ABS(fm%local_data(i, j)) > thresh) THEN
     758           0 :                   WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
     759           0 :                      ABS(fm%local_data(i, j))
     760             :                END IF
     761             :             END DO
     762             :          END DO
     763             :       ELSE
     764           0 :          DO i = 1, nrow_local
     765           0 :             DO j = 1, ncol_local
     766           0 :                IF (ABS(fm%local_data(i, j)) > thresh) THEN
     767           0 :                   WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
     768           0 :                      fm%local_data(i, j)
     769             :                END IF
     770             :             END DO
     771             :          END DO
     772             :       END IF
     773           0 :       CALL fm%matrix_struct%para_env%sync()
     774           0 :       IF (unit_nr > 0) THEN
     775           0 :          WRITE (unit_nr, *) my_footer
     776             :       END IF
     777             : 
     778           0 :       CALL timestop(handle)
     779             : 
     780           0 :    END SUBROUTINE
     781             : 
     782             : ! **************************************************************************************************
     783             : !> \brief Write the atomic coordinates to the output unit.
     784             : !> \param particle_set ...
     785             : !>       \note Adapted from particle_methods.F [MG]
     786             : !> \param unit_nr ...
     787             : ! **************************************************************************************************
     788           4 :    SUBROUTINE write_qs_particle_coordinates_bse(particle_set, unit_nr)
     789             : 
     790             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     791             :       INTEGER, INTENT(IN)                                :: unit_nr
     792             : 
     793             :       CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_particle_coordinates_bse'
     794             : 
     795             :       CHARACTER(LEN=2)                                   :: element_symbol
     796             :       INTEGER                                            :: handle, iatom, natom
     797             : 
     798           4 :       CALL timeset(routineN, handle)
     799             : 
     800           4 :       IF (unit_nr > 0) THEN
     801           2 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     802           2 :          WRITE (unit_nr, '(T2,A4,T13,A7,16X,A7,15X,A7,15X,A7)') 'BSE|', &
     803           4 :             'Element', 'x [Å]', 'y [Å]', 'z [Å]'
     804           2 :          natom = SIZE(particle_set)
     805           8 :          DO iatom = 1, natom
     806             :             CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     807           6 :                                  element_symbol=element_symbol)
     808             :             WRITE (unit_nr, '(T2,A4,T8,A12,1X,3(5X,F15.4))') &
     809          26 :                'BSE|', element_symbol, particle_set(iatom)%r(1:3)*angstrom
     810             :          END DO
     811             :       END IF
     812             : 
     813           4 :       CALL timestop(handle)
     814             : 
     815           4 :    END SUBROUTINE write_qs_particle_coordinates_bse
     816             : 
     817             : END MODULE bse_print

Generated by: LCOV version 1.15