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