LCOV - code coverage report
Current view: top level - src - bse_print.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.0 % 410 377
Test Date: 2025-07-25 12:55:17 Functions: 88.9 % 9 8

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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_and_print_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           30 :    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           30 :       CALL timeset(routineN, handle)
      62              : 
      63           30 :       IF (unit_nr > 0) THEN
      64           15 :          WRITE (unit_nr, *) ' '
      65           15 :          WRITE (unit_nr, '(T2,A79)') '*******************************************************************************'
      66           15 :          WRITE (unit_nr, '(T2,A79)') '**                                                                           **'
      67           15 :          WRITE (unit_nr, '(T2,A79)') '**           Bethe Salpeter equation (BSE) for excitation energies           **'
      68           15 :          IF (bse_tda .AND. bse_abba) THEN
      69            0 :             WRITE (unit_nr, '(T2,A79)') '**          solved with and without Tamm-Dancoff approximation (TDA)         **'
      70           15 :          ELSE IF (bse_tda) THEN
      71            6 :             WRITE (unit_nr, '(T2,A79)') '**                solved with Tamm-Dancoff approximation (TDA)               **'
      72              :          ELSE
      73            9 :             WRITE (unit_nr, '(T2,A79)') '**               solved without Tamm-Dancoff approximation (TDA)             **'
      74              :          END IF
      75              : 
      76           15 :          WRITE (unit_nr, '(T2,A79)') '**                                                                           **'
      77           15 :          WRITE (unit_nr, '(T2,A79)') '*******************************************************************************'
      78           15 :          WRITE (unit_nr, *) ' '
      79              :       END IF
      80              : 
      81           30 :       CALL timestop(handle)
      82              : 
      83           30 :    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           30 :    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           30 :       CALL timeset(routineN, handle)
     111              : 
     112           30 :       IF (unit_nr > 0) THEN
     113           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     114           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     115           15 :          IF (flag_TDA) THEN
     116            6 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     117            6 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '*   Bethe Salpeter equation (BSE) with Tamm Dancoff approximation (TDA)  *'
     118            6 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     119            6 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     120            6 :             WRITE (unit_nr, '(T2,A4,T7,A48,A23)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
     121           12 :                'the BSE within the TDA:'
     122            6 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     123            6 :             WRITE (unit_nr, '(T2,A4,T29,A16)') 'BSE|', 'A X^n = Ω^n X^n'
     124            6 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     125            6 :             WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
     126            6 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     127            6 :             WRITE (unit_nr, '(T2,A4,T7,A41)') 'BSE|', 'sum_jb ( A_ia,jb   X_jb^n ) = Ω^n X_ia^n'
     128            6 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     129            6 :             WRITE (unit_nr, '(T2,A4,T7,A30)') 'BSE|', 'prelim Ref.: Eq. (36) with B=0'
     130            6 :             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            9 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     133            9 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '*      Full ("ABBA") Bethe Salpeter equation (BSE) (i.e. without TDA)    *'
     134            9 :             WRITE (unit_nr, '(T2,A4,T7,A74)') 'BSE|', '**************************************************************************'
     135            9 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     136            9 :             WRITE (unit_nr, '(T2,A4,T7,A48,A24)') 'BSE|', 'The excitations are calculated by diagonalizing ', &
     137           18 :                'the BSE without the TDA:'
     138            9 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     139            9 :             WRITE (unit_nr, '(T2,A4,T22,A30)') 'BSE|', '|A B| |X^n|       |1  0| |X^n|'
     140            9 :             WRITE (unit_nr, '(T2,A4,T22,A31)') 'BSE|', '|B A| |Y^n| = Ω^n |0 -1| |Y^n|'
     141            9 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     142            9 :             WRITE (unit_nr, '(T2,A4,T7,A23)') 'BSE|', 'i.e. in index notation:'
     143            9 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     144            9 :             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            9 :             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           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     148           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     149           15 :          WRITE (unit_nr, '(T2,A4,T7,A4,T18,A42,T70,A1,I4,A1,I4,A1)') 'BSE|', 'i,j:', &
     150           30 :             'occupied molecular orbitals, i.e. state in', '[', homo_irred - homo + 1, ',', homo_irred, ']'
     151           15 :          WRITE (unit_nr, '(T2,A4,T7,A4,T18,A44,T70,A1,I4,A1,I4,A1)') 'BSE|', 'a,b:', &
     152           30 :             'unoccupied molecular orbitals, i.e. state in', '[', homo_irred + 1, ',', homo_irred + virtual, ']'
     153           15 :          WRITE (unit_nr, '(T2,A4,T7,A2,T18,A16)') 'BSE|', 'n:', 'Excitation index'
     154           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     155           15 :          IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
     156           12 :             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           15 :          IF (.NOT. flag_TDA) THEN
     161            9 :             IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
     162            6 :                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            9 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     167            9 :             WRITE (unit_nr, '(T2,A4,T7,A35)') 'BSE|', 'prelim Ref.: Eqs. (24-27),(30),(35)'
     168            9 :             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           15 :          IF (.NOT. flag_TDA) THEN
     171            9 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     172            9 :             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            9 :             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           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     176           15 :          WRITE (unit_nr, '(T2,A4,T7,A7,T31,A23)') 'BSE|', 'ε_...:', 'GW quasiparticle energy'
     177           15 :          WRITE (unit_nr, '(T2,A4,T7,A7,T31,A15)') 'BSE|', 'δ_...:', 'Kronecker delta'
     178           15 :          WRITE (unit_nr, '(T2,A4,T7,A3,T31,A21)') 'BSE|', 'α:', 'spin-dependent factor (Singlet/Triplet)'
     179           15 :          WRITE (unit_nr, '(T2,A4,T7,A6,T30,A34)') 'BSE|', 'v_...:', 'Electron-hole exchange interaction'
     180           15 :          IF (mp2_env%bse%screening_method == bse_screening_w0) THEN
     181           12 :             WRITE (unit_nr, '(T2,A4,T7,A,T31,A)') 'BSE|', 'W_... = 1/ϵ v_...:', &
     182           24 :                'Direct interaction screened by '
     183           12 :             WRITE (unit_nr, '(T2,A4,T30,A)') 'BSE|', &
     184           24 :                '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           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     192           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     193           15 :          WRITE (unit_nr, '(T2,A4,T7,A47,A7,A9,F3.1)') 'BSE|', &
     194           30 :             'The spin-dependent factor is for the requested ', multiplet, " is α = ", alpha
     195           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     196              :       END IF
     197              : 
     198           30 :       CALL timestop(handle)
     199              : 
     200           30 :    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           30 :    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           30 :       CALL timeset(routineN, handle)
     228              : 
     229           30 :       IF (unit_nr > 0) THEN
     230           15 :          IF (flag_TDA) THEN
     231            6 :             WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', 'Excitation energies from solving the BSE within the TDA:'
     232              :          ELSE
     233            9 :             WRITE (unit_nr, '(T2,A4,T7,A57)') 'BSE|', 'Excitation energies from solving the BSE without the TDA:'
     234              :          END IF
     235           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     236           15 :          WRITE (unit_nr, '(T2,A4,T11,A12,T26,A11,T44,A8,T55,A27)') 'BSE|', &
     237           30 :             'Excitation n', "Spin Config", 'TDA/ABBA', 'Excitation energy Ω^n (eV)'
     238              :       END IF
     239              :       !prints actual energies values
     240           30 :       IF (unit_nr > 0) THEN
     241          373 :          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          373 :                'BSE|', i_exc, multiplet, info_approximation, Exc_ens(i_exc)*evolt
     244              :          END DO
     245              :       END IF
     246              : 
     247           30 :       CALL timestop(handle)
     248              : 
     249           30 :    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           30 :    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           30 :       CALL timeset(routineN, handle)
     277              : 
     278           30 :       IF (unit_nr > 0) THEN
     279           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     280              :          WRITE (unit_nr, '(T2,A4,T7,A61)') &
     281           15 :             'BSE|', "Single-particle transitions are built up by (de-)excitations,"
     282              :          WRITE (unit_nr, '(T2,A4,T7,A18)') &
     283           15 :             'BSE|', "which we denote by"
     284           15 :          WRITE (unit_nr, '(T2,A4)') 'BSE|'
     285              :          WRITE (unit_nr, '(T2,A4,T20,A2,T30,A40)') &
     286           15 :             'BSE|', "=>", "for excitations, i.e. entries of X_ia^n,"
     287              :          WRITE (unit_nr, '(T2,A4,T20,A2,T30,A42)') &
     288           15 :             'BSE|', "<=", "for deexcitations, i.e. entries of Y_ia^n."
     289              :          WRITE (unit_nr, '(T2,A4)') &
     290           15 :             'BSE|'
     291              :          WRITE (unit_nr, '(T2,A4,T7,A73)') &
     292           15 :             'BSE|', "The following single-particle transitions have significant contributions,"
     293              :          WRITE (unit_nr, '(T2,A4,T7,A16,F5.3,A15,F5.3,A16)') &
     294           15 :             'BSE|', "i.e. |X_ia^n| > ", mp2_env%bse%eps_x, " or |Y_ia^n| > ", &
     295           30 :             mp2_env%bse%eps_x, ", respectively :"
     296              : 
     297           15 :          WRITE (unit_nr, '(T2,A4,T15,A27,I5,A13,I5,A3)') 'BSE|', '-- Quick reminder: HOMO i =', &
     298           30 :             homo_irred, ' and LUMO a =', homo_irred + 1, " --"
     299           15 :          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           15 :             "BSE|", "Excitation n", "i", "=>/<=", "a", 'TDA/ABBA', "|X_ia^n|/|Y_ia^n|"
     302              :       END IF
     303          746 :       DO i_exc = 1, MIN(homo*virtual, mp2_env%bse%num_print_exc)
     304          716 :          IF (unit_nr > 0) THEN
     305          358 :             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          716 :                                                unit_nr, mp2_env)
     311          746 :          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          450 :                                                   unit_nr, mp2_env)
     315              :          END IF
     316              :       END DO
     317              : 
     318           30 :       CALL timestop(handle)
     319              : 
     320           30 :    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           30 :    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           30 :       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           30 :       IF (mp2_env%bse%bse_spin_config == 0) THEN
     357           30 :          IF (unit_nr > 0) THEN
     358           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     359              :             WRITE (unit_nr, '(T2,A4,T7,A60)') &
     360           15 :                'BSE|', "Transition moments d_r^n (with r∈(x,y,z), in atomic units)"
     361              :             WRITE (unit_nr, '(T2,A4,T7,A67)') &
     362           15 :                'BSE|', "and oscillator strength f^n of excitation level n are obtained from"
     363           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     364           15 :             IF (flag_TDA) THEN
     365              :                WRITE (unit_nr, '(T2,A4,T10,A)') &
     366            6 :                   '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            9 :                   'BSE|', "d_r^n = sum_ia sqrt(2) < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n )"
     370              :             END IF
     371           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     372              :             WRITE (unit_nr, '(T2,A4,T14,A)') &
     373           15 :                'BSE|', "f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2"
     374           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     375              :             WRITE (unit_nr, '(T2,A4,T7,A19)') &
     376           15 :                'BSE|', "where we introduced"
     377           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     378              :             WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
     379           15 :                'BSE|', "ψ_i:", "occupied molecular orbitals,"
     380              :             WRITE (unit_nr, '(T2,A4,T7,A5,T15,A28)') &
     381           15 :                'BSE|', "ψ_a:", "empty molecular orbitals and"
     382              :             WRITE (unit_nr, '(T2,A4,T9,A2,T14,A18)') &
     383           15 :                'BSE|', "r:", "position operator."
     384           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     385              :             WRITE (unit_nr, '(T2,A4,T7,A28)') &
     386           15 :                'BSE|', "prelim Ref.: Eqs. (23), (24)"
     387              :             WRITE (unit_nr, '(T2,A4,T7,A71)') &
     388           15 :                'BSE|', "in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290"
     389           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     390           15 :             IF (flag_TDA) THEN
     391            6 :                WRITE (unit_nr, '(T2,A4,T7,A55)') 'BSE|', &
     392           12 :                   'Optical properties from solving the BSE within the TDA:'
     393              :             ELSE
     394            9 :                WRITE (unit_nr, '(T2,A4,T7,A56)') 'BSE|', &
     395           18 :                   'Optical properties from solving the BSE without the TDA:'
     396              :             END IF
     397           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     398           15 :             WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T38,A5,T48,A5,T58,A5,T64,A17)') 'BSE|', &
     399           30 :                'Excitation n', "TDA/ABBA", "d_x^n", "d_y^n", "d_z^n", 'Osc. strength f^n'
     400          373 :             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          358 :                   'BSE|', i_exc, info_approximation, trans_mom_bse(1, 1, i_exc), trans_mom_bse(2, 1, i_exc), &
     403          731 :                   trans_mom_bse(3, 1, i_exc), oscill_str(i_exc)
     404              :             END DO
     405           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     406           15 :             WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
     407           30 :                'Check for Thomas-Reiche-Kuhn sum rule'
     408           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     409           15 :             WRITE (unit_nr, '(T2,A4,T35,A15)') 'BSE|', &
     410           30 :                'N_e = Σ_n f^n'
     411           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     412           15 :             WRITE (unit_nr, '(T2,A4,T7,A24,T65,I16)') 'BSE|', &
     413           30 :                'Number of electrons N_e:', homo_irred*2
     414           15 :             WRITE (unit_nr, '(T2,A4,T7,A40,T66,F16.3)') 'BSE|', &
     415          710 :                'Sum over oscillator strengths Σ_n f^n :', SUM(oscill_str)
     416           15 :             WRITE (unit_nr, '(T2,A4)') 'BSE|'
     417           15 :             IF (mp2_env%bse%bse_cutoff_occ > 0 .OR. mp2_env%bse%bse_cutoff_empty > 0) THEN
     418              :                CALL cp_warn(__LOCATION__, &
     419           15 :                             "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           30 :          IF (mp2_env%bse%bse_print_spectrum) THEN
     425              :             CALL compute_and_print_absorption_spectrum(oscill_str, polarizability_residues, Exc_ens, &
     426            2 :                                                        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           30 :       CALL timestop(handle)
     441              : 
     442           30 :    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         1166 :    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         1166 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: idx_homo, idx_virt
     471         1166 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eigvec_entries
     472              : 
     473         1166 :       CALL timeset(routineN, handle)
     474              : 
     475              :       CALL filter_eigvec_contrib(fm_eigvec, idx_homo, idx_virt, eigvec_entries, &
     476         1166 :                                  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         1166 :       IF (unit_nr > 0) THEN
     480         1456 :          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          873 :                "BSE|", i_exc, homo_irred - homo + idx_homo(k), direction_excitation, &
     483         2329 :                homo_irred + idx_virt(k), info_approximation, ABS(eigvec_entries(k))
     484              :          END DO
     485              :       END IF
     486         1166 :       DEALLOCATE (idx_homo)
     487         1166 :       DEALLOCATE (idx_virt)
     488         1166 :       DEALLOCATE (eigvec_entries)
     489         1166 :       CALL timestop(handle)
     490              : 
     491         1166 :    END SUBROUTINE
     492              : 
     493              : ! **************************************************************************************************
     494              : !> \brief                              Prints exciton descriptors, cf. Mewes et al., JCTC 14, 710-725 (2018)
     495              : !> \param exc_descr                    Exciton descriptors with size of num_print_exc_descr
     496              : !> \param ref_point_multipole          Reference point for computation of multipole moments, e.g. center of mass
     497              : !> \param unit_nr ...
     498              : !> \param num_print_exc_descr          Number of excitation levels for which descriptors are printed
     499              : !> \param print_checkvalue             Flag, which determines if values for regtests should be printed
     500              : !> \param print_directional_exc_descr  Flag, which activates printing of directional descriptors
     501              : !> \param prefix_output                String, which is put in front of prints, i.e. "BSE|" or "" for TDDFPT
     502              : !> \param qs_env ...
     503              : ! **************************************************************************************************
     504           12 :    SUBROUTINE print_exciton_descriptors(exc_descr, ref_point_multipole, unit_nr, &
     505              :                                         num_print_exc_descr, print_checkvalue, print_directional_exc_descr, &
     506              :                                         prefix_output, qs_env)
     507              : 
     508              :       TYPE(exciton_descr_type), ALLOCATABLE, &
     509              :          DIMENSION(:)                                    :: exc_descr
     510              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
     511              :          INTENT(IN)                                      :: ref_point_multipole
     512              :       INTEGER, INTENT(IN)                                :: unit_nr, num_print_exc_descr
     513              :       LOGICAL, INTENT(IN)                                :: print_checkvalue, &
     514              :                                                             print_directional_exc_descr
     515              :       CHARACTER(LEN=4), INTENT(IN)                       :: prefix_output
     516              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     517              : 
     518              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_exciton_descriptors'
     519              : 
     520              :       CHARACTER(LEN=1), DIMENSION(3)                     :: array_direction_str
     521              :       CHARACTER(LEN=5)                                   :: method_name
     522              :       INTEGER                                            :: handle, i_dir, i_exc
     523              :       REAL(KIND=dp)                                      :: d_eh_dir, d_exc_dir, sigma_e_dir, &
     524              :                                                             sigma_h_dir
     525            6 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     526              : 
     527            6 :       IF (prefix_output == 'BSE|') THEN
     528            4 :          method_name = 'BSE'
     529              :       ELSE
     530            2 :          method_name = 'TDDFT'
     531              :       END IF
     532              : 
     533            6 :       CALL timeset(routineN, handle)
     534            6 :       CALL get_qs_env(qs_env, particle_set=particle_set)
     535            6 :       IF (unit_nr > 0) THEN
     536            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     537            6 :             'Exciton descriptors for excitation level n are given by'
     538            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     539            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     540            6 :             'd_eh   = | <r_h - r_e>_exc |'
     541            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     542            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     543            6 :             'σ_e    = sqrt( <r_e^2>_exc - <r_e>_exc^2 )'
     544            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     545            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     546            6 :             'σ_h    = sqrt( <r_h^2>_exc - <r_h>_exc^2 )'
     547            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     548            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     549            6 :             'COV_eh = <r_e r_h>_exc - <r_e>_exc <r_h>_exc'
     550            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     551            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     552            6 :             'd_exc  = sqrt( | < |r_h - r_e|^2 >_exc )'
     553            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     554            6 :             '       = sqrt( d_eh^2 + σ_e^2 + σ_h^2 - 2 * COV_eh )'
     555            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     556            3 :          WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     557            6 :             'R_eh   = COV_eh / (σ_e * σ_h)'
     558            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     559            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     560            6 :             'where the expectation values <.>_exc are taken with respect to the '
     561            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     562            6 :             'exciton wavefunction  of excitation n:'
     563            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     564              : 
     565            3 :          IF (exc_descr(1)%flag_TDA) THEN
     566            2 :             WRITE (unit_nr, '(T2,A4,T20,A)') prefix_output, &
     567            4 :                '𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e)  ,'
     568              :          ELSE
     569            1 :             WRITE (unit_nr, '(T2,A4,T20,A)') prefix_output, &
     570            2 :                '𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e)'
     571            1 :             WRITE (unit_nr, '(T2,A4,T40,A)') prefix_output, &
     572            2 :                '+ Y_ia^n ψ_a(r_h) ψ_i(r_e) ,'
     573              :          END IF
     574            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     575            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     576            6 :             'i.e.'
     577            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     578            3 :          WRITE (unit_nr, '(T2,A4,T20,A)') prefix_output, &
     579            6 :             '< O >_exc = < 𝚿_n | O | 𝚿_n > / < 𝚿_n | 𝚿_n >  ,'
     580            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     581            3 :          IF (exc_descr(1)%flag_TDA) THEN
     582            2 :             WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     583            4 :                'where c_n = < 𝚿_n | 𝚿_n > = 1 within TDA.'
     584              :          ELSE
     585            1 :             WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     586            2 :                'where c_n = < 𝚿_n | 𝚿_n > deviates from 1 without TDA.'
     587              :          END IF
     588            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     589            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     590            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     591            6 :             'Here, we introduced'
     592            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     593              :          WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
     594            3 :             prefix_output, "ψ_i:", "occupied molecular orbitals,"
     595              :          WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
     596            3 :             prefix_output, "ψ_a:", "empty molecular orbitals and"
     597              :          WRITE (unit_nr, '(T2,A4,T9,A2,T14,A)') &
     598            3 :             prefix_output, "r:", "position operator."
     599            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     600            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     601            6 :             'prelim Ref.: Eqs. (15)-(22)'
     602            3 :          WRITE (unit_nr, '(T2,A4,T7,A,A)') prefix_output, &
     603            3 :             'JCTC 2018, 14, 710-725; ', &
     604            6 :             'http://doi.org/10.1021/acs.jctc.7b01145'
     605            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     606            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     607            3 :          IF (exc_descr(1)%flag_TDA) THEN
     608            2 :             WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
     609            4 :                'Exciton descriptors from solving the ', method_name, ' within the TDA:'
     610              :          ELSE
     611            1 :             WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
     612            2 :                'Exciton descriptors from solving the ', method_name, ' without the TDA:'
     613              :          END IF
     614            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     615            3 :          WRITE (unit_nr, '(T2,A4,T10,A1,6X,A3,1X,4X,A10,5X,A10,5X,A10,3X,A11,8X,A4)') prefix_output, &
     616            6 :             'n', 'c_n', 'd_eh [Å]', 'σ_e [Å]', 'σ_h [Å]', 'd_exc [Å]', 'R_eh'
     617           58 :          DO i_exc = 1, num_print_exc_descr
     618              :             WRITE (unit_nr, '(T2,A4,T7,I4,4X,F5.3,1X,5(2X,F10.4))') &
     619           55 :                prefix_output, i_exc, exc_descr(i_exc)%norm_XpY, &
     620           55 :                exc_descr(i_exc)%diff_r_abs*angstrom, &
     621           55 :                exc_descr(i_exc)%sigma_e*angstrom, exc_descr(i_exc)%sigma_h*angstrom, &
     622          113 :                exc_descr(i_exc)%diff_r_sqr*angstrom, exc_descr(i_exc)%corr_e_h
     623              :          END DO
     624            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     625              :          ! For debug runs, print first d_exc separately to allow the regtests to read in
     626            3 :          IF (print_checkvalue) THEN
     627            3 :             WRITE (unit_nr, '(T2)')
     628            3 :             WRITE (unit_nr, '(T2,A28,T65,F16.4)') 'Checksum exciton descriptors', &
     629            6 :                exc_descr(1)%diff_r_sqr*angstrom
     630            3 :             WRITE (unit_nr, '(T2)')
     631              :          END IF
     632            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     633              :          ! Print exciton descriptor resolved per direction
     634            3 :          IF (print_directional_exc_descr) THEN
     635            1 :             WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     636            2 :                'We can restrict the exciton descriptors to a specific direction,'
     637            1 :             WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     638            2 :                'e.g. the x-components are:'
     639            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     640            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     641            2 :                'd_eh^x   = | <x_h - x_e>_exc |'
     642            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     643            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     644            2 :                'σ_e^x    = sqrt( <x_e^2>_exc - <x_e>_exc^2 )'
     645            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     646            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     647            2 :                'σ_h^x    = sqrt( <x_h^2>_exc - <x_h>_exc^2 )'
     648            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     649            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     650            2 :                'COV_eh^x = <x_e x_h>_exc - <x_e>_exc <x_h>_exc'
     651            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     652            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     653            2 :                'd_exc^x  = sqrt( | < |x_h - x_e|^2 >_exc )'
     654            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     655            2 :                '         = sqrt( (d_eh^x)^2 + (σ_e^x)^2'
     656            1 :             WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
     657            2 :                '                 + (σ_h^x)^2 - 2 * (COV_eh^x) )'
     658            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     659            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     660            1 :             IF (exc_descr(1)%flag_TDA) THEN
     661            1 :                WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
     662            2 :                   'Exciton descriptors per direction from solving the ', method_name, ' within the TDA:'
     663              :             ELSE
     664            0 :                WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
     665            0 :                   'Exciton descriptors per direction from solving the ', method_name, ' without the TDA:'
     666              :             END IF
     667            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     668            1 :             WRITE (unit_nr, '(T2,A4,T12,A1,2X,A9,5X,A12,5X,A12,5X,A12,3X,A13)') prefix_output, &
     669            2 :                'n', 'r = x/y/z', 'd_eh^r [Å]', 'σ_e^r [Å]', 'σ_h^r [Å]', 'd_exc^r [Å]'
     670            6 :             DO i_exc = 1, num_print_exc_descr
     671           20 :                DO i_dir = 1, 3
     672           60 :                   array_direction_str = (/"x", "y", "z"/)
     673           15 :                   d_eh_dir = ABS(exc_descr(i_exc)%r_h(i_dir) - exc_descr(i_exc)%r_e(i_dir))
     674           15 :                   sigma_e_dir = SQRT(exc_descr(i_exc)%r_e_sq(i_dir) - exc_descr(i_exc)%r_e(i_dir)**2)
     675           15 :                   sigma_h_dir = SQRT(exc_descr(i_exc)%r_h_sq(i_dir) - exc_descr(i_exc)%r_h(i_dir)**2)
     676              :                   d_exc_dir = SQRT(d_eh_dir**2 + sigma_e_dir**2 + sigma_h_dir**2 &
     677           15 :                                    - 2*exc_descr(i_exc)%cov_e_h(i_dir))
     678              :                   WRITE (unit_nr, '(T2,A4,T9,I4,10X,A1,1X,4(4X,F10.4))') &
     679           15 :                      prefix_output, i_exc, array_direction_str(i_dir), &
     680           15 :                      d_eh_dir*angstrom, &
     681           15 :                      sigma_e_dir*angstrom, sigma_h_dir*angstrom, &
     682           35 :                      d_exc_dir*angstrom
     683              :                END DO
     684            6 :                WRITE (unit_nr, '(T2,A4)') prefix_output
     685              :             END DO
     686            1 :             WRITE (unit_nr, '(T2,A4)') prefix_output
     687              :          END IF
     688              :          ! Print the reference atomic geometry for the exciton descriptors
     689            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     690            6 :             'With the center of charge as reference point r_0,'
     691            3 :          WRITE (unit_nr, '(T2,A4,T15,A7,F10.4,A2,F10.4,A2,F10.4,A1)') prefix_output, &
     692            3 :             'r_0 = (', ref_point_multipole(1)*angstrom, ', ', ref_point_multipole(2)*angstrom, ', ', &
     693            6 :             ref_point_multipole(3)*angstrom, ')'
     694            3 :          WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
     695            6 :             'we further obtain r_e and r_h from solving the ', method_name, ' within the TDA'
     696            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     697            3 :          WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') prefix_output, &
     698            6 :             'Excitation n', 'x_e [Å]', 'y_e [Å]', 'z_e [Å]'
     699           58 :          DO i_exc = 1, num_print_exc_descr
     700              :             WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
     701           55 :                prefix_output, i_exc, &
     702          278 :                exc_descr(i_exc)%r_e_shift(:)*angstrom
     703              :          END DO
     704            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     705            3 :          WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') prefix_output, &
     706            6 :             'Excitation n', 'x_h [Å]', 'y_h [Å]', 'z_h [Å]'
     707           58 :          DO i_exc = 1, num_print_exc_descr
     708              :             WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
     709           55 :                prefix_output, i_exc, &
     710          278 :                exc_descr(i_exc)%r_h_shift(:)*angstrom
     711              :          END DO
     712            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     713            3 :          WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
     714            6 :             'The reference atomic geometry for these values is given by'
     715              :       END IF
     716            6 :       CALL write_qs_particle_coordinates_bse(particle_set, unit_nr, prefix_output)
     717            6 :       IF (unit_nr > 0) THEN
     718            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     719              :       END IF
     720            6 :       CALL timestop(handle)
     721              : 
     722            6 :    END SUBROUTINE print_exciton_descriptors
     723              : 
     724              : ! **************************************************************************************************
     725              : !> \brief Debug function to write elements of a full matrix to file, if they are larger than a given threshold
     726              : !> \param fm ...
     727              : !> \param thresh ...
     728              : !> \param header ...
     729              : !> \param unit_nr ...
     730              : !> \param abs_vals ...
     731              : ! **************************************************************************************************
     732            0 :    SUBROUTINE fm_write_thresh(fm, thresh, header, unit_nr, abs_vals)
     733              : 
     734              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm
     735              :       REAL(KIND=dp), INTENT(IN)                          :: thresh
     736              :       CHARACTER(LEN=*), INTENT(IN)                       :: header
     737              :       INTEGER, INTENT(IN)                                :: unit_nr
     738              :       LOGICAL, OPTIONAL                                  :: abs_vals
     739              : 
     740              :       CHARACTER(LEN=*), PARAMETER :: my_footer = " | ENDING WRITING OF MATRIX", &
     741              :          routineN = 'fm_write_thresh'
     742              : 
     743              :       INTEGER                                            :: handle, i, j, ncol_local, nrow_local
     744            0 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     745              :       LOGICAL                                            :: my_abs_vals
     746              : 
     747            0 :       CALL timeset(routineN, handle)
     748              : 
     749            0 :       IF (PRESENT(abs_vals)) THEN
     750            0 :          my_abs_vals = abs_vals
     751              :       ELSE
     752              :          my_abs_vals = .FALSE.
     753              :       END IF
     754              : 
     755              :       CALL cp_fm_get_info(matrix=fm, &
     756              :                           nrow_local=nrow_local, &
     757              :                           ncol_local=ncol_local, &
     758              :                           row_indices=row_indices, &
     759            0 :                           col_indices=col_indices)
     760              : 
     761            0 :       IF (unit_nr > 0) THEN
     762            0 :          WRITE (unit_nr, *) header
     763              :       END IF
     764            0 :       IF (my_abs_vals) THEN
     765            0 :          DO i = 1, nrow_local
     766            0 :             DO j = 1, ncol_local
     767            0 :                IF (ABS(fm%local_data(i, j)) > thresh) THEN
     768            0 :                   IF (unit_nr > 0) THEN
     769            0 :                      WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
     770            0 :                         ABS(fm%local_data(i, j))
     771              :                   END IF
     772              :                END IF
     773              :             END DO
     774              :          END DO
     775              :       ELSE
     776            0 :          DO i = 1, nrow_local
     777            0 :             DO j = 1, ncol_local
     778            0 :                IF (ABS(fm%local_data(i, j)) > thresh) THEN
     779            0 :                   IF (unit_nr > 0) THEN
     780            0 :                      WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
     781            0 :                         fm%local_data(i, j)
     782              :                   END IF
     783              :                END IF
     784              :             END DO
     785              :          END DO
     786              :       END IF
     787            0 :       CALL fm%matrix_struct%para_env%sync()
     788            0 :       IF (unit_nr > 0) THEN
     789            0 :          WRITE (unit_nr, *) my_footer
     790              :       END IF
     791              : 
     792            0 :       CALL timestop(handle)
     793              : 
     794            0 :    END SUBROUTINE
     795              : 
     796              : ! **************************************************************************************************
     797              : !> \brief Write the atomic coordinates to the output unit.
     798              : !> \param particle_set ...
     799              : !>       \note Adapted from particle_methods.F [MG]
     800              : !> \param unit_nr ...
     801              : !> \param prefix_output ...
     802              : ! **************************************************************************************************
     803            6 :    SUBROUTINE write_qs_particle_coordinates_bse(particle_set, unit_nr, prefix_output)
     804              : 
     805              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     806              :       INTEGER, INTENT(IN)                                :: unit_nr
     807              :       CHARACTER(LEN=4), INTENT(IN)                       :: prefix_output
     808              : 
     809              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_particle_coordinates_bse'
     810              : 
     811              :       CHARACTER(LEN=2)                                   :: element_symbol
     812              :       INTEGER                                            :: handle, iatom, natom
     813              : 
     814            6 :       CALL timeset(routineN, handle)
     815              : 
     816            6 :       IF (unit_nr > 0) THEN
     817            3 :          WRITE (unit_nr, '(T2,A4)') prefix_output
     818            3 :          WRITE (unit_nr, '(T2,A4,T13,A7,16X,A7,15X,A7,15X,A7)') prefix_output, &
     819            6 :             'Element', 'x [Å]', 'y [Å]', 'z [Å]'
     820            3 :          natom = SIZE(particle_set)
     821           12 :          DO iatom = 1, natom
     822              :             CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     823            9 :                                  element_symbol=element_symbol)
     824              :             WRITE (unit_nr, '(T2,A4,T8,A12,1X,3(5X,F15.4))') &
     825           39 :                prefix_output, element_symbol, particle_set(iatom)%r(1:3)*angstrom
     826              :          END DO
     827              :       END IF
     828              : 
     829            6 :       CALL timestop(handle)
     830              : 
     831            6 :    END SUBROUTINE write_qs_particle_coordinates_bse
     832              : 
     833              : END MODULE bse_print
        

Generated by: LCOV version 2.0-1