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