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