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