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 print_BSE_start_flag
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 print_output_header
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 1640 : 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 1057 : "BSE|", i_exc, homo_irred - homo + idx_homo(k), direction_excitation, &
483 2697 : 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 print_transition_amplitudes_core
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 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
524 :
525 6 : IF (prefix_output == 'BSE|') THEN
526 4 : method_name = 'BSE'
527 : ELSE
528 2 : method_name = 'TDDFT'
529 : END IF
530 :
531 6 : CALL timeset(routineN, handle)
532 6 : CALL get_qs_env(qs_env, particle_set=particle_set)
533 6 : IF (unit_nr > 0) THEN
534 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
535 6 : 'Exciton descriptors for excitation level n are given by'
536 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
537 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
538 6 : 'd_eh = | <r_h - r_e>_exc |'
539 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
540 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
541 6 : 'σ_e = sqrt( <r_e^2>_exc - <r_e>_exc^2 )'
542 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
543 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
544 6 : 'σ_h = sqrt( <r_h^2>_exc - <r_h>_exc^2 )'
545 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
546 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
547 6 : 'COV_eh = <r_e r_h>_exc - <r_e>_exc <r_h>_exc'
548 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
549 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
550 6 : 'd_exc = sqrt( | < |r_h - r_e|^2 >_exc )'
551 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
552 6 : ' = sqrt( d_eh^2 + σ_e^2 + σ_h^2 - 2 * COV_eh )'
553 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
554 3 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
555 6 : 'R_eh = COV_eh / (σ_e * σ_h)'
556 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
557 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
558 6 : 'where the expectation values <.>_exc are taken with respect to the '
559 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
560 6 : 'exciton wavefunction of excitation n:'
561 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
562 :
563 3 : IF (exc_descr(1)%flag_TDA) THEN
564 2 : WRITE (unit_nr, '(T2,A4,T20,A)') prefix_output, &
565 4 : '𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e) ,'
566 : ELSE
567 1 : WRITE (unit_nr, '(T2,A4,T20,A)') prefix_output, &
568 2 : '𝚿_n(r_e,r_h) = Σ_{i,a} X_ia^n ψ_i(r_h) ψ_a(r_e)'
569 1 : WRITE (unit_nr, '(T2,A4,T40,A)') prefix_output, &
570 2 : '+ Y_ia^n ψ_a(r_h) ψ_i(r_e) ,'
571 : END IF
572 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
573 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
574 6 : 'i.e.'
575 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
576 3 : WRITE (unit_nr, '(T2,A4,T20,A)') prefix_output, &
577 6 : '< O >_exc = < 𝚿_n | O | 𝚿_n > / < 𝚿_n | 𝚿_n > ,'
578 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
579 3 : IF (exc_descr(1)%flag_TDA) THEN
580 2 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
581 4 : 'where c_n = < 𝚿_n | 𝚿_n > = 1 within TDA.'
582 : ELSE
583 1 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
584 2 : 'where c_n = < 𝚿_n | 𝚿_n > deviates from 1 without TDA.'
585 : END IF
586 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
587 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
588 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
589 6 : 'Here, we introduced'
590 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
591 : WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
592 3 : prefix_output, "ψ_i:", "occupied molecular orbitals,"
593 : WRITE (unit_nr, '(T2,A4,T7,A5,T15,A)') &
594 3 : prefix_output, "ψ_a:", "empty molecular orbitals and"
595 : WRITE (unit_nr, '(T2,A4,T9,A2,T14,A)') &
596 3 : prefix_output, "r:", "position operator."
597 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
598 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
599 6 : 'prelim Ref.: Eqs. (15)-(22)'
600 3 : WRITE (unit_nr, '(T2,A4,T7,A,A)') prefix_output, &
601 3 : 'JCTC 2018, 14, 710-725; ', &
602 6 : 'http://doi.org/10.1021/acs.jctc.7b01145'
603 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
604 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
605 3 : IF (exc_descr(1)%flag_TDA) THEN
606 2 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
607 4 : 'Exciton descriptors from solving the ', method_name, ' within the TDA:'
608 : ELSE
609 1 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
610 2 : 'Exciton descriptors from solving the ', method_name, ' without the TDA:'
611 : END IF
612 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
613 3 : WRITE (unit_nr, '(T2,A4,T10,A1,6X,A3,1X,4X,A10,5X,A10,5X,A10,3X,A11,8X,A4)') prefix_output, &
614 6 : 'n', 'c_n', 'd_eh [Å]', 'σ_e [Å]', 'σ_h [Å]', 'd_exc [Å]', 'R_eh'
615 58 : DO i_exc = 1, num_print_exc_descr
616 : WRITE (unit_nr, '(T2,A4,T7,I4,4X,F5.3,1X,5(2X,F10.4))') &
617 55 : prefix_output, i_exc, exc_descr(i_exc)%norm_XpY, &
618 55 : exc_descr(i_exc)%diff_r_abs*angstrom, &
619 55 : exc_descr(i_exc)%sigma_e*angstrom, exc_descr(i_exc)%sigma_h*angstrom, &
620 113 : exc_descr(i_exc)%diff_r_sqr*angstrom, exc_descr(i_exc)%corr_e_h
621 : END DO
622 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
623 : ! For debug runs, print first d_exc separately to allow the regtests to read in
624 3 : IF (print_checkvalue) THEN
625 3 : WRITE (unit_nr, '(T2)')
626 3 : WRITE (unit_nr, '(T2,A28,T65,F16.4)') 'Checksum exciton descriptors', &
627 6 : exc_descr(1)%diff_r_sqr*angstrom
628 3 : WRITE (unit_nr, '(T2)')
629 : END IF
630 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
631 : ! Print exciton descriptor resolved per direction
632 3 : IF (print_directional_exc_descr) THEN
633 1 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
634 2 : 'We can restrict the exciton descriptors to a specific direction,'
635 1 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
636 2 : 'e.g. the x-components are:'
637 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
638 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
639 2 : 'd_eh^x = | <x_h - x_e>_exc |'
640 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
641 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
642 2 : 'σ_e^x = sqrt( <x_e^2>_exc - <x_e>_exc^2 )'
643 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
644 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
645 2 : 'σ_h^x = sqrt( <x_h^2>_exc - <x_h>_exc^2 )'
646 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
647 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
648 2 : "COV_eh^{μμ'} = <r^μ_e r^μ'_h>_exc - <r^μ_e>_exc <r^μ'_h>_exc"
649 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
650 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
651 2 : 'd_exc^x = sqrt( | < |x_h - x_e|^2 >_exc )'
652 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
653 2 : ' = sqrt( (d_eh^x)^2 + (σ_e^x)^2'
654 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
655 2 : " + (σ_h^x)^2 - 2 * (COV_eh^{xx}) )"
656 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
657 1 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
658 2 : "Subsequently, the cross-correlation matrix R_eh^{μμ'} is printed"
659 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
660 1 : WRITE (unit_nr, '(T2,A4,T15,A)') prefix_output, &
661 2 : "R_eh^{μμ'} = COV_eh^{μμ'}/(σ^μ_e σ^μ_h) "
662 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
663 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
664 1 : IF (exc_descr(1)%flag_TDA) THEN
665 1 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
666 2 : 'Exciton descriptors per direction from solving the ', method_name, ' within the TDA:'
667 : ELSE
668 0 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
669 0 : 'Exciton descriptors per direction from solving the ', method_name, ' without the TDA:'
670 : END IF
671 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
672 1 : WRITE (unit_nr, '(T2,A4,T12,A1,2X,A9,5X,A12,5X,A12,5X,A12,3X,A13)') prefix_output, &
673 2 : 'n', 'r = x/y/z', 'd_eh^r [Å]', 'σ_e^r [Å]', 'σ_h^r [Å]', 'd_exc^r [Å]'
674 6 : DO i_exc = 1, num_print_exc_descr
675 20 : DO i_dir = 1, 3
676 60 : array_direction_str = ["x", "y", "z"]
677 : WRITE (unit_nr, '(T2,A4,T9,I4,10X,A1,1X,4(4X,F10.4))') &
678 15 : prefix_output, i_exc, array_direction_str(i_dir), &
679 15 : exc_descr(i_exc)%d_eh_dir(i_dir)*angstrom, &
680 15 : exc_descr(i_exc)%sigma_e_dir(i_dir)*angstrom, &
681 15 : exc_descr(i_exc)%sigma_h_dir(i_dir)*angstrom, &
682 35 : exc_descr(i_exc)%d_exc_dir(i_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 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
688 1 : IF (exc_descr(1)%flag_TDA) THEN
689 1 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
690 2 : 'Crosscorrelation matrix from solving the ', method_name, ' within the TDA:'
691 : ELSE
692 0 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
693 0 : 'Crosscorrelation matrix from solving the ', method_name, ' without the TDA:'
694 : END IF
695 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
696 1 : WRITE (unit_nr, '(T2,A4,T12,A1,8X,6(8X,A2))') prefix_output, &
697 2 : 'n', 'xx', 'yy', 'zz', 'xy', 'xz', 'yz'
698 6 : DO i_exc = 1, num_print_exc_descr
699 : WRITE (unit_nr, '(T2,A4,T9,I4,8X,6(3X,F7.4),3X,F7.4)') &
700 5 : prefix_output, i_exc, &
701 5 : exc_descr(i_exc)%corr_e_h_matrix(1, 1), &
702 5 : exc_descr(i_exc)%corr_e_h_matrix(2, 2), &
703 5 : exc_descr(i_exc)%corr_e_h_matrix(3, 3), &
704 5 : exc_descr(i_exc)%corr_e_h_matrix(1, 2), &
705 5 : exc_descr(i_exc)%corr_e_h_matrix(1, 3), &
706 11 : exc_descr(i_exc)%corr_e_h_matrix(2, 3)
707 : END DO
708 1 : WRITE (unit_nr, '(T2,A4)') prefix_output
709 : END IF
710 : ! Print the reference atomic geometry for the exciton descriptors
711 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
712 6 : 'With the center of charge as reference point r_0,'
713 3 : WRITE (unit_nr, '(T2,A4,T15,A7,F10.4,A2,F10.4,A2,F10.4,A1)') prefix_output, &
714 3 : 'r_0 = (', ref_point_multipole(1)*angstrom, ', ', ref_point_multipole(2)*angstrom, ', ', &
715 6 : ref_point_multipole(3)*angstrom, ')'
716 3 : IF (exc_descr(1)%flag_TDA) THEN
717 2 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
718 4 : 'we further obtain r_e and r_h from solving the ', method_name, ' within the TDA'
719 : ELSE
720 1 : WRITE (unit_nr, '(T2,A4,T7,A,A,A)') prefix_output, &
721 2 : 'we further obtain r_e and r_h from solving the ', method_name, ' without the TDA'
722 : END IF
723 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
724 3 : WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') prefix_output, &
725 6 : 'Excitation n', 'x_e [Å]', 'y_e [Å]', 'z_e [Å]'
726 58 : DO i_exc = 1, num_print_exc_descr
727 : WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
728 55 : prefix_output, i_exc, &
729 278 : exc_descr(i_exc)%r_e_shift(:)*angstrom
730 : END DO
731 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
732 3 : WRITE (unit_nr, '(T2,A4,T8,A12,1X,13X,A9,13X,A9,13X,A9)') prefix_output, &
733 6 : 'Excitation n', 'x_h [Å]', 'y_h [Å]', 'z_h [Å]'
734 58 : DO i_exc = 1, num_print_exc_descr
735 : WRITE (unit_nr, '(T2,A4,T8,I12,1X,3(5X,F15.4))') &
736 55 : prefix_output, i_exc, &
737 278 : exc_descr(i_exc)%r_h_shift(:)*angstrom
738 : END DO
739 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
740 3 : WRITE (unit_nr, '(T2,A4,T7,A)') prefix_output, &
741 6 : 'The reference atomic geometry for these values is given by'
742 : END IF
743 6 : CALL write_qs_particle_coordinates_bse(particle_set, unit_nr, prefix_output)
744 6 : IF (unit_nr > 0) THEN
745 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
746 : END IF
747 6 : CALL timestop(handle)
748 :
749 6 : END SUBROUTINE print_exciton_descriptors
750 :
751 : ! **************************************************************************************************
752 : !> \brief Debug function to write elements of a full matrix to file, if they are larger than a given threshold
753 : !> \param fm ...
754 : !> \param thresh ...
755 : !> \param header ...
756 : !> \param unit_nr ...
757 : !> \param abs_vals ...
758 : ! **************************************************************************************************
759 0 : SUBROUTINE fm_write_thresh(fm, thresh, header, unit_nr, abs_vals)
760 :
761 : TYPE(cp_fm_type), INTENT(IN) :: fm
762 : REAL(KIND=dp), INTENT(IN) :: thresh
763 : CHARACTER(LEN=*), INTENT(IN) :: header
764 : INTEGER, INTENT(IN) :: unit_nr
765 : LOGICAL, OPTIONAL :: abs_vals
766 :
767 : CHARACTER(LEN=*), PARAMETER :: my_footer = " | ENDING WRITING OF MATRIX", &
768 : routineN = 'fm_write_thresh'
769 :
770 : INTEGER :: handle, i, j, ncol_local, nrow_local
771 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
772 : LOGICAL :: my_abs_vals
773 :
774 0 : CALL timeset(routineN, handle)
775 :
776 0 : IF (PRESENT(abs_vals)) THEN
777 0 : my_abs_vals = abs_vals
778 : ELSE
779 : my_abs_vals = .FALSE.
780 : END IF
781 :
782 : CALL cp_fm_get_info(matrix=fm, &
783 : nrow_local=nrow_local, &
784 : ncol_local=ncol_local, &
785 : row_indices=row_indices, &
786 0 : col_indices=col_indices)
787 :
788 0 : IF (unit_nr > 0) THEN
789 0 : WRITE (unit_nr, *) header
790 : END IF
791 0 : IF (my_abs_vals) THEN
792 0 : DO i = 1, nrow_local
793 0 : DO j = 1, ncol_local
794 0 : IF (ABS(fm%local_data(i, j)) > thresh) THEN
795 0 : IF (unit_nr > 0) THEN
796 0 : WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
797 0 : ABS(fm%local_data(i, j))
798 : END IF
799 : END IF
800 : END DO
801 : END DO
802 : ELSE
803 0 : DO i = 1, nrow_local
804 0 : DO j = 1, ncol_local
805 0 : IF (ABS(fm%local_data(i, j)) > thresh) THEN
806 0 : IF (unit_nr > 0) THEN
807 0 : WRITE (unit_nr, "(A7,T10,I5,T20,I5,T30,F13.5)") header, row_indices(i), col_indices(j), &
808 0 : fm%local_data(i, j)
809 : END IF
810 : END IF
811 : END DO
812 : END DO
813 : END IF
814 0 : CALL fm%matrix_struct%para_env%sync()
815 0 : IF (unit_nr > 0) THEN
816 0 : WRITE (unit_nr, *) my_footer
817 : END IF
818 :
819 0 : CALL timestop(handle)
820 :
821 0 : END SUBROUTINE fm_write_thresh
822 :
823 : ! **************************************************************************************************
824 : !> \brief Write the atomic coordinates to the output unit.
825 : !> \param particle_set ...
826 : !> \note Adapted from particle_methods.F [MG]
827 : !> \param unit_nr ...
828 : !> \param prefix_output ...
829 : ! **************************************************************************************************
830 6 : SUBROUTINE write_qs_particle_coordinates_bse(particle_set, unit_nr, prefix_output)
831 :
832 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
833 : INTEGER, INTENT(IN) :: unit_nr
834 : CHARACTER(LEN=4), INTENT(IN) :: prefix_output
835 :
836 : CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_particle_coordinates_bse'
837 :
838 : CHARACTER(LEN=2) :: element_symbol
839 : INTEGER :: handle, iatom, natom
840 :
841 6 : CALL timeset(routineN, handle)
842 :
843 6 : IF (unit_nr > 0) THEN
844 3 : WRITE (unit_nr, '(T2,A4)') prefix_output
845 3 : WRITE (unit_nr, '(T2,A4,T13,A7,16X,A7,15X,A7,15X,A7)') prefix_output, &
846 6 : 'Element', 'x [Å]', 'y [Å]', 'z [Å]'
847 3 : natom = SIZE(particle_set)
848 12 : DO iatom = 1, natom
849 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
850 9 : element_symbol=element_symbol)
851 : WRITE (unit_nr, '(T2,A4,T8,A12,1X,3(5X,F15.4))') &
852 39 : prefix_output, element_symbol, particle_set(iatom)%r(1:3)*angstrom
853 : END DO
854 : END IF
855 :
856 6 : CALL timestop(handle)
857 :
858 6 : END SUBROUTINE write_qs_particle_coordinates_bse
859 :
860 : END MODULE bse_print
|
