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 computing excitonic properties, e.g. exciton diameter, from the BSE
10 : !> \par History
11 : !> 10.2024 created [Maximilian Graml]
12 : ! **************************************************************************************************
13 : MODULE bse_properties
14 : USE bse_types, ONLY: bse_env_type
15 : USE bse_util, ONLY: print_bse_nto_cubes,&
16 : reshuffle_eigvec,&
17 : trace_exciton_descr
18 : USE cp_files, ONLY: close_file,&
19 : open_file
20 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
21 : cp_fm_trace
22 : USE cp_fm_diag, ONLY: cp_fm_svd
23 : USE cp_fm_redistribution, ONLY: cp_fm_permute
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_release,&
26 : cp_fm_struct_type
27 : USE cp_fm_types, ONLY: &
28 : cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_release, cp_fm_set_all, &
29 : cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_to_fm_submat_general, cp_fm_type
30 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 : section_vals_type,&
32 : section_vals_val_get
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: pi
35 : USE parallel_gemm_api, ONLY: parallel_gemm
36 : USE physcon, ONLY: c_light_au,&
37 : evolt
38 : USE qs_mo_types, ONLY: allocate_mo_set,&
39 : deallocate_mo_set,&
40 : init_mo_set,&
41 : mo_set_type
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'bse_properties'
49 :
50 : PUBLIC :: exciton_descr_type
51 :
52 : PUBLIC :: get_exciton_descriptors, get_oscillator_strengths, compute_and_print_absorption_spectrum, &
53 : calculate_NTOs
54 :
55 : ! TYPE definitions for exciton wavefunction descriptors
56 :
57 : TYPE exciton_descr_type
58 : REAL(KIND=dp), DIMENSION(3) :: r_e = 0.0_dp, &
59 : r_h = 0.0_dp, &
60 : r_e_sq = 0.0_dp, &
61 : r_h_sq = 0.0_dp, &
62 : r_e_shift = 0.0_dp, &
63 : r_h_shift = 0.0_dp, &
64 : d_eh_dir = 0.0_dp, &
65 : sigma_e_dir = 0.0_dp, &
66 : sigma_h_dir = 0.0_dp, &
67 : d_exc_dir = 0.0_dp
68 : REAL(KIND=dp), DIMENSION(3, 3) :: r_e_h = 0.0_dp, &
69 : cov_e_h = 0.0_dp, &
70 : corr_e_h_matrix = 0.0_dp
71 : REAL(KIND=dp) :: sigma_e = 0.0_dp, &
72 : sigma_h = 0.0_dp, &
73 : cov_e_h_sum = 0.0_dp, &
74 : corr_e_h = 0.0_dp, &
75 : diff_r_abs = 0.0_dp, &
76 : diff_r_sqr = 0.0_dp, &
77 : norm_XpY = 0.0_dp
78 : LOGICAL :: flag_TDA = .FALSE.
79 : END TYPE exciton_descr_type
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n )
85 : !> and oscillator strengths f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2
86 : !> Prelim Ref.: Eqs. (23), (24)
87 : !> in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290
88 : !> \param fm_eigvec_X ...
89 : !> \param Exc_ens ...
90 : !> \param fm_dipole_ai_trunc ...
91 : !> \param trans_mom_bse BSE dipole vectors in real space per excitation level
92 : !> \param oscill_str Oscillator strength per excitation level
93 : !> \param polarizability_residues Residues of polarizability ("tensorial oscillator strength")
94 : !> per excitation level
95 : !> \param homo_red ...
96 : !> \param virtual_red ...
97 : !> \param fm_eigvec_Y ...
98 : ! **************************************************************************************************
99 46 : SUBROUTINE get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
100 : trans_mom_bse, oscill_str, polarizability_residues, &
101 : homo_red, virtual_red, fm_eigvec_Y)
102 :
103 : TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_X
104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
105 : INTENT(IN) :: Exc_ens
106 : TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_ai_trunc
107 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
108 : INTENT(OUT) :: trans_mom_bse
109 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
110 : INTENT(OUT) :: oscill_str
111 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
112 : INTENT(OUT) :: polarizability_residues
113 : INTEGER, INTENT(IN) :: homo_red, virtual_red
114 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_Y
115 :
116 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_oscillator_strengths'
117 :
118 : INTEGER :: handle, idir, jdir, n, n_exc
119 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_dipole_MO_trunc_reordered, &
120 : fm_struct_trans_mom_bse
121 : TYPE(cp_fm_type) :: fm_eigvec_XYsum
122 184 : TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_MO_trunc_reordered, &
123 368 : fm_dipole_per_dir, fm_trans_mom_bse
124 :
125 46 : CALL timeset(routineN, handle)
126 :
127 : ! an iterative solver provides fewer excitations than transitions
128 46 : CALL cp_fm_get_info(fm_eigvec_X, ncol_global=n_exc)
129 :
130 : CALL cp_fm_struct_create(fm_struct_dipole_MO_trunc_reordered, fm_eigvec_X%matrix_struct%para_env, &
131 46 : fm_eigvec_X%matrix_struct%context, 1, homo_red*virtual_red)
132 : CALL cp_fm_struct_create(fm_struct_trans_mom_bse, fm_eigvec_X%matrix_struct%para_env, &
133 46 : fm_eigvec_X%matrix_struct%context, 1, n_exc)
134 :
135 : ! Include excitonic amplitudes in dipoles, i.e. obtain "BSE dipoles":
136 : ! \vec{D}_n = sqrt(2) * sum_{i,a} \vec{D}_ai (X_{ai}^{(n)} + Y_{ai}^{(n)})
137 :
138 : ! Reorder dipoles in order to execute the sum over i and a by parallel gemm
139 184 : DO idir = 1, 3
140 : CALL cp_fm_create(fm_dipole_MO_trunc_reordered(idir), matrix_struct=fm_struct_dipole_MO_trunc_reordered, &
141 138 : name="dipoles_mo_reordered")
142 138 : CALL cp_fm_set_all(fm_dipole_MO_trunc_reordered(idir), 0.0_dp)
143 : CALL cp_fm_permute(fm_dipole_MO_trunc_reordered(idir), fm_dipole_ai_trunc(idir), 1.0_dp, 1, 1, 1, &
144 138 : virtual_red, [2, 4, 3, 1])
145 184 : CALL cp_fm_release(fm_dipole_per_dir(idir))
146 : END DO
147 :
148 184 : DO idir = 1, 3
149 : CALL cp_fm_create(fm_trans_mom_bse(idir), matrix_struct=fm_struct_trans_mom_bse, &
150 138 : name="excitonic_dipoles")
151 184 : CALL cp_fm_set_all(fm_trans_mom_bse(idir), 0.0_dp)
152 : END DO
153 :
154 : ! If TDA is invoked, Y is not present as it is simply 0
155 46 : CALL cp_fm_create(fm_eigvec_XYsum, matrix_struct=fm_eigvec_X%matrix_struct, name="excit_amplitude_sum")
156 46 : CALL cp_fm_set_all(fm_eigvec_XYsum, 0.0_dp)
157 46 : CALL cp_fm_to_fm(fm_eigvec_X, fm_eigvec_XYsum)
158 46 : IF (PRESENT(fm_eigvec_Y)) THEN
159 26 : CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_XYsum, 1.0_dp, fm_eigvec_Y)
160 : END IF
161 184 : DO idir = 1, 3
162 : CALL parallel_gemm('N', 'N', 1, n_exc, homo_red*virtual_red, SQRT(2.0_dp), &
163 184 : fm_dipole_MO_trunc_reordered(idir), fm_eigvec_XYsum, 0.0_dp, fm_trans_mom_bse(idir))
164 : END DO
165 :
166 : ! Get oscillator strengths themselves
167 138 : ALLOCATE (oscill_str(n_exc))
168 : ! trans_mom_bse needs to be a 2D array per direction idir, such that cp_fm_get_submatrix can directly
169 : ! write to it
170 138 : ALLOCATE (trans_mom_bse(3, 1, n_exc))
171 138 : ALLOCATE (polarizability_residues(3, 3, n_exc))
172 46 : trans_mom_bse(:, :, :) = 0.0_dp
173 :
174 : ! Sum over all directions
175 184 : DO idir = 1, 3
176 184 : CALL cp_fm_get_submatrix(fm_trans_mom_bse(idir), trans_mom_bse(idir, :, :))
177 : END DO
178 :
179 1984 : DO n = 1, n_exc
180 7752 : DO idir = 1, 3
181 25194 : DO jdir = 1, 3
182 23256 : polarizability_residues(idir, jdir, n) = 2.0_dp*Exc_ens(n)*trans_mom_bse(idir, 1, n)*trans_mom_bse(jdir, 1, n)
183 : END DO
184 : END DO
185 7798 : oscill_str(n) = 2.0_dp/3.0_dp*Exc_ens(n)*SUM(ABS(trans_mom_bse(:, 1, n))**2)
186 : END DO
187 :
188 46 : CALL cp_fm_struct_release(fm_struct_dipole_MO_trunc_reordered)
189 46 : CALL cp_fm_struct_release(fm_struct_trans_mom_bse)
190 184 : DO idir = 1, 3
191 138 : CALL cp_fm_release(fm_dipole_MO_trunc_reordered(idir))
192 138 : CALL cp_fm_release(fm_trans_mom_bse(idir))
193 184 : CALL cp_fm_release(fm_dipole_ai_trunc(idir))
194 : END DO
195 46 : CALL cp_fm_release(fm_eigvec_XYsum)
196 :
197 46 : CALL timestop(handle)
198 :
199 138 : END SUBROUTINE get_oscillator_strengths
200 :
201 : ! **************************************************************************************************
202 : !> \brief Computes and returns absorption spectrum for the frequency range and broadening
203 : !> provided by the user.
204 : !> Prelim Ref.: C. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications
205 : !> (Oxford University Press, Oxford, 2012), Eq. 7.51
206 : !> \param oscill_str ...
207 : !> \param polarizability_residues ...
208 : !> \param Exc_ens ...
209 : !> \param info_approximation ...
210 : !> \param unit_nr ...
211 : !> \param bse_env the BSE environment (settings; para_env, whose ionode writes the files)
212 : ! **************************************************************************************************
213 2 : SUBROUTINE compute_and_print_absorption_spectrum(oscill_str, polarizability_residues, Exc_ens, &
214 : info_approximation, unit_nr, bse_env)
215 :
216 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
217 : INTENT(IN) :: oscill_str
218 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :), &
219 : INTENT(IN) :: polarizability_residues
220 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
221 : INTENT(IN) :: Exc_ens
222 : CHARACTER(LEN=10) :: info_approximation
223 : INTEGER, INTENT(IN) :: unit_nr
224 : TYPE(bse_env_type), INTENT(IN) :: bse_env
225 :
226 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_and_print_absorption_spectrum'
227 :
228 : CHARACTER(LEN=10) :: eta_str, width_eta_format_str
229 : CHARACTER(LEN=40) :: file_name_crosssection, &
230 : file_name_spectrum
231 : INTEGER :: handle, i, idir, j, jdir, k, num_steps, &
232 : unit_nr_file, width_eta
233 : REAL(KIND=dp) :: eta, freq_end, freq_start, freq_step, &
234 : omega
235 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: abs_cross_section, abs_spectrum
236 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: eta_list
237 :
238 2 : CALL timeset(routineN, handle)
239 :
240 2 : freq_step = bse_env%bse_spectrum_freq_step_size
241 2 : freq_start = bse_env%bse_spectrum_freq_start
242 2 : freq_end = bse_env%bse_spectrum_freq_end
243 2 : eta_list => bse_env%bse_eta_spectrum_list
244 :
245 : ! Calculate number of steps to fit given frequency range
246 2 : num_steps = NINT((freq_end - freq_start)/freq_step) + 1
247 :
248 4 : DO k = 1, SIZE(eta_list)
249 2 : eta = eta_list(k)
250 :
251 : ! Some magic to get a nice formatting of the eta value in filenames
252 2 : width_eta = MAX(1, INT(LOG10(eta)) + 1) + 4
253 2 : WRITE (width_eta_format_str, "(A2,I0,A3)") '(F', width_eta, '.3)'
254 2 : WRITE (eta_str, width_eta_format_str) eta*evolt
255 : ! Filename itself
256 2 : file_name_spectrum = 'BSE'//TRIM(ADJUSTL(info_approximation))//'eta='//TRIM(eta_str)//'.spectrum'
257 2 : file_name_crosssection = 'BSE'//TRIM(ADJUSTL(info_approximation))//'eta='//TRIM(eta_str)//'.crosssection'
258 :
259 : ! First column is frequency in eV, second column is imaginary part of the trace of the polarizability
260 : ! The following 9 columns are the entries of the polarizability tensor
261 6 : ALLOCATE (abs_spectrum(num_steps, 11))
262 2 : abs_spectrum(:, :) = 0.0_dp
263 : ! Also calculate and print the photoabsorption cross section tensor
264 : ! σ_{µ µ'}(ω) = 4πω Im{α_{µ µ'}(ω)} / c
265 4 : ALLOCATE (abs_cross_section(num_steps, 11))
266 2 : abs_cross_section(:, :) = 0.0_dp
267 :
268 : ! Calculate the imaginary part of the mean dipole polarizability α_{avg}(ω)
269 : ! which is given by (cf. C. Ullrichs Book on TDDFT, Eq. 7.51)
270 : ! We introduce an additional - due to his convention for charge vs particle density, see also:
271 : ! Computer Physics Communications, 208:149–161, November 2016
272 : ! https://doi.org/10.1016/j.cpc.2016.06.019
273 : ! α_{avg}(ω) = - \sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}
274 : ! and then the imaginary part is (in the limit η -> 0)
275 : ! Im[α_{avg}(ω)] = - \sum_{n=1}^{N_exc} f_n * η / ((ω - Ω^n)² + η²)
276 : ! where f_n are the oscillator strengths and E_exc the excitation energies
277 : ! For the full polarizability tensor, we have
278 : ! α_{µ µ'}(ω) = - \sum_n [2 Ω^n d^n_µ d^n_µ'] / [(ω+iη)^2- (Ω^n)^2]
279 : ! = - \sum_n "polarizability_residues" / [(ω+iη)^2- (Ω^n)^2]
280 2004 : DO i = 1, num_steps
281 2002 : omega = freq_start + (i - 1)*freq_step
282 2002 : abs_spectrum(i, 1) = omega
283 98100 : DO j = 1, SIZE(oscill_str)
284 : abs_spectrum(i, 2) = abs_spectrum(i, 2) - oscill_str(j)* &
285 96096 : AIMAG(1/((omega + CMPLX(0.0, eta, kind=dp))**2 - Exc_ens(j)**2))
286 386386 : DO idir = 1, 3
287 1249248 : DO jdir = 1, 3
288 : ! Factor 2 from formula for tensor is already in the polarizability_residues
289 : ! to follow the same convention as the oscillator strengths
290 : abs_spectrum(i, 2 + (idir - 1)*3 + jdir) = abs_spectrum(i, 2 + (idir - 1)*3 + jdir) &
291 : - polarizability_residues(idir, jdir, j)* &
292 1153152 : AIMAG(1/((omega + CMPLX(0.0, eta, kind=dp))**2 - Exc_ens(j)**2))
293 : END DO
294 : END DO
295 : END DO
296 : END DO
297 :
298 : ! Extract cross section σ from polarizability tensor
299 2004 : DO i = 1, num_steps
300 2002 : omega = abs_spectrum(i, 1)
301 2002 : abs_cross_section(i, 1) = omega
302 22024 : abs_cross_section(i, 2:) = 4.0_dp*pi*abs_spectrum(i, 2:)*omega/c_light_au
303 : END DO
304 :
305 : !For debug runs: Export an entry of the two tensors to allow regtests on spectra
306 2 : IF (bse_env%bse_debug_print) THEN
307 2 : IF (unit_nr > 0) THEN
308 1 : WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
309 1 : 'Averaged dynamical dipole polarizability at 8.2 eV:', &
310 2 : abs_spectrum(83, 2)
311 1 : WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
312 1 : 'Averaged photoabsorption cross section at 8.2 eV:', &
313 2 : abs_cross_section(83, 2)
314 : END IF
315 : END IF
316 :
317 : ! Print it to file from the ionode of the BSE communicator, at every print level; the
318 : ! default logger inside the MP2/GW driver is the subgroup's and would make every group's root write
319 2 : IF (bse_env%para_env%is_source()) THEN
320 : CALL open_file(file_name_crosssection, unit_number=unit_nr_file, &
321 1 : file_status="UNKNOWN", file_action="WRITE")
322 : WRITE (unit_nr_file, '(A,A6)') "# Photoabsorption cross section σ_{µ µ'}(ω) = -4πω/c * Im[ \sum_n "// &
323 1 : "[2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] ] from Bethe Salpeter equation for method ", &
324 2 : TRIM(ADJUSTL(info_approximation))
325 1 : WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", "σ_{avg}(ω)", "σ_xx(ω)", &
326 1 : "σ_xy(ω)", "σ_xz(ω)", "σ_yx(ω)", "σ_yy(ω)", "σ_yz(ω)", "σ_zx(ω)", &
327 2 : "σ_zy(ω)", "σ_zz(ω)"
328 1002 : DO i = 1, num_steps
329 1002 : WRITE (unit_nr_file, '(11(F20.8,1X))') abs_cross_section(i, 1)*evolt, abs_cross_section(i, 2:11)
330 : END DO
331 1 : CALL close_file(unit_nr_file)
332 : END IF
333 2 : DEALLOCATE (abs_cross_section)
334 :
335 2 : IF (bse_env%para_env%is_source()) THEN
336 : CALL open_file(file_name_spectrum, unit_number=unit_nr_file, &
337 1 : file_status="UNKNOWN", file_action="WRITE")
338 : WRITE (unit_nr_file, '(A,A6)') "# Imaginary part of polarizability α_{µ µ'}(ω) = -\sum_n "// &
339 1 : "[2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²] from Bethe Salpeter equation for method ", &
340 2 : TRIM(ADJUSTL(info_approximation))
341 1 : WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", "Im{α_{avg}(ω)}", "Im{α_xx(ω)}", &
342 1 : "Im{α_xy(ω)}", "Im{α_xz(ω)}", "Im{α_yx(ω)}", "Im{α_yy(ω)}", "Im{α_yz(ω)}", "Im{α_zx(ω)}", &
343 2 : "Im{α_zy(ω)}", "Im{α_zz(ω)}"
344 1002 : DO i = 1, num_steps
345 1002 : WRITE (unit_nr_file, '(11(F20.8,1X))') abs_spectrum(i, 1)*evolt, abs_spectrum(i, 2:11)
346 : END DO
347 1 : CALL close_file(unit_nr_file)
348 : END IF
349 4 : DEALLOCATE (abs_spectrum)
350 : END DO
351 :
352 2 : IF (unit_nr > 0) THEN
353 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
354 : WRITE (unit_nr, '(T2,A4,T7,A,A)') &
355 1 : 'BSE|', "Printed optical absorption spectrum to local files, e.g. "
356 : WRITE (unit_nr, '(T2,A4,T7,A)') &
357 1 : 'BSE|', file_name_spectrum
358 : WRITE (unit_nr, '(T2,A4,T7,A,A)') &
359 1 : 'BSE|', "as well as photoabsorption cross section to, e.g. "
360 : WRITE (unit_nr, '(T2,A4,T7,A)') &
361 1 : 'BSE|', file_name_crosssection
362 : WRITE (unit_nr, '(T2,A4,T7,A52)') &
363 1 : 'BSE|', "using the Eq. (7.51) from C. Ullrichs Book on TDDFT:"
364 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
365 : WRITE (unit_nr, '(T2,A4,T10,A75)') &
366 1 : 'BSE|', "Im{α_{avg}(ω)} = -Im{\sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}}"
367 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
368 : WRITE (unit_nr, '(T2,A4,T7,A)') &
369 1 : 'BSE|', "or for the full polarizability tensor:"
370 : WRITE (unit_nr, '(T2,A4,T10,A)') &
371 1 : 'BSE|', "α_{µ µ'}(ω) = -\sum_n [2 Ω^n d_µ^n d_µ'^n] / [(ω+iη)²- (Ω^n)²]"
372 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
373 : WRITE (unit_nr, '(T2,A4,T7,A)') &
374 1 : 'BSE|', "as well as Eq. (7.48):"
375 : WRITE (unit_nr, '(T2,A4,T10,A)') &
376 1 : 'BSE|', "σ_{µ µ'}(ω) = 4πω Im{α_{µ µ'}(ω)} / c"
377 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
378 : WRITE (unit_nr, '(T2,A4,T7,A)') &
379 1 : 'BSE|', "with transition moments d_µ^n, oscillator strengths f_n,"
380 : WRITE (unit_nr, '(T2,A4,T7,A)') &
381 1 : 'BSE|', "excitation energies Ω^n and the speed of light c."
382 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
383 : WRITE (unit_nr, '(T2,A4,T7,A)') &
384 1 : 'BSE|', "Please note that we adopt an additional minus sign for both quantities,"
385 : WRITE (unit_nr, '(T2,A4,T7,A)') &
386 1 : 'BSE|', "due to the convention for charge vs particle density as done in MolGW:"
387 : WRITE (unit_nr, '(T2,A4,T7,A)') &
388 1 : 'BSE|', "https://doi.org/10.1016/j.cpc.2016.06.019."
389 1 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
390 : END IF
391 :
392 2 : CALL timestop(handle)
393 :
394 4 : END SUBROUTINE compute_and_print_absorption_spectrum
395 :
396 : ! **************************************************************************************************
397 : !> \brief ...
398 : !> \param fm_X ...
399 : !> \param fm_Y ...
400 : !> \param mo_coeff ...
401 : !> \param homo ...
402 : !> \param virtual ...
403 : !> \param info_approximation ...
404 : !> \param oscill_str ...
405 : !> \param bse_env the BSE environment (settings, the MO sets and the section of the NTO print keys)
406 : !> \param unit_nr ...
407 : ! **************************************************************************************************
408 4 : SUBROUTINE calculate_NTOs(fm_X, fm_Y, &
409 4 : mo_coeff, homo, virtual, &
410 : info_approximation, &
411 : oscill_str, &
412 : bse_env, unit_nr)
413 :
414 : TYPE(cp_fm_type), INTENT(IN) :: fm_X
415 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_Y
416 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
417 : INTEGER, INTENT(IN) :: homo, virtual
418 : CHARACTER(LEN=10) :: info_approximation
419 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str
420 : TYPE(bse_env_type), INTENT(IN) :: bse_env
421 : INTEGER, INTENT(IN) :: unit_nr
422 :
423 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_NTOs'
424 : REAL(KIND=dp), PARAMETER :: coeff_err = 1.0E-5_dp
425 :
426 : CHARACTER(LEN=20), DIMENSION(2) :: nto_name
427 : INTEGER :: handle, homo_irred, i, i_nto, info_svd, &
428 : j, n_exc, n_nto, nao_full, nao_trunc
429 4 : INTEGER, DIMENSION(:), POINTER :: stride
430 : LOGICAL :: append_cube, cube_file
431 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigval_svd_squ
432 4 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigval_svd
433 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_m, fm_struct_mo_coeff, &
434 : fm_struct_nto_holes, &
435 : fm_struct_nto_particles, &
436 : fm_struct_nto_set
437 : TYPE(cp_fm_type) :: fm_eigvl, fm_eigvr_t, fm_m, fm_mo_coeff, fm_nto_coeff_holes, &
438 : fm_nto_coeff_particles, fm_nto_set, fm_X_ia, fm_Y_ai
439 4 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
440 : TYPE(section_vals_type), POINTER :: bse_section, nto_section
441 :
442 4 : CALL timeset(routineN, handle)
443 4 : bse_section => bse_env%bse_section
444 :
445 4 : nao_full = bse_env%mos(1)%nao
446 4 : nao_trunc = homo + virtual
447 : ! This is not influenced by the BSE cutoff
448 4 : homo_irred = bse_env%mos(1)%homo
449 : ! M will have a block structure and is quadratic in homo+virtual, i.e.
450 : ! occ virt
451 : ! | 0 X_i,a | occ = homo
452 : ! M = | Y_a,i 0 | virt = virtual
453 : !
454 : ! X and Y are here not the eigenvectors X_ia,n - instead we fix n and reshape the combined ia index
455 : ! Notice the index structure of the lower block, i.e. X is transposed
456 : CALL cp_fm_struct_create(fm_struct_m, &
457 : fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
458 4 : nao_trunc, nao_trunc)
459 : CALL cp_fm_struct_create(fm_struct_mo_coeff, &
460 : fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
461 4 : nao_full, nao_trunc)
462 : CALL cp_fm_struct_create(fm_struct_nto_holes, &
463 : fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
464 4 : nao_full, nao_trunc)
465 : CALL cp_fm_struct_create(fm_struct_nto_particles, &
466 : fm_X%matrix_struct%para_env, fm_X%matrix_struct%context, &
467 4 : nao_full, nao_trunc)
468 :
469 : CALL cp_fm_create(fm_mo_coeff, matrix_struct=fm_struct_mo_coeff, &
470 4 : name="mo_coeff")
471 : ! Here, we take care of possible cutoffs
472 : ! Simply truncating the matrix causes problems with the print function
473 : ! Therefore, we keep the dimension, but set the coefficients of truncated indices to 0
474 : CALL cp_fm_to_fm_submat_general(mo_coeff(1), fm_mo_coeff, &
475 : nao_full, nao_trunc, &
476 : 1, homo_irred - homo + 1, &
477 : 1, 1, &
478 4 : mo_coeff(1)%matrix_struct%context)
479 :
480 : ! Print some information about the NTOs
481 4 : IF (unit_nr > 0) THEN
482 2 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
483 4 : 'The Natural Transition Orbital (NTO) pairs φ_I(r_e) and χ_I(r_h) for a fixed'
484 2 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
485 4 : 'excitation index n are obtained by singular value decomposition of T'
486 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
487 2 : WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
488 4 : ' = (0 X)'
489 2 : IF (PRESENT(fm_Y)) THEN
490 1 : WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
491 2 : 'T = (Y^T 0)'
492 : ELSE
493 1 : WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
494 2 : 'T = (0 0)'
495 : END IF
496 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
497 2 : WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
498 4 : 'T = U Λ V^T'
499 2 : WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
500 4 : 'φ_I(r_e) = \sum_p V_pI ψ_p(r_e)'
501 2 : WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
502 4 : 'χ_I(r_h) = \sum_p U_pI ψ_p(r_e)'
503 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
504 2 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
505 4 : 'where we have introduced'
506 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
507 : WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
508 2 : 'BSE|', "ψ_p:", "occupied and virtual molecular orbitals,"
509 : WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
510 2 : 'BSE|', "φ_I(r_e):", "NTO state for the electron,"
511 : WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
512 2 : 'BSE|', "χ_I(r_h):", "NTO state for the hole,"
513 : WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
514 2 : 'BSE|', "Λ:", "diagonal matrix of NTO weights λ_I,"
515 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
516 2 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
517 4 : "The NTOs are calculated with the following settings:"
518 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
519 2 : WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Number of excitations, for which NTOs are computed', &
520 4 : bse_env%num_print_exc_ntos
521 2 : IF (bse_env%eps_nto_osc_str > 0.0_dp) THEN
522 0 : WRITE (unit_nr, '(T2,A4,T7,A,T71,F10.3)') 'BSE|', 'Threshold for oscillator strength f^n', &
523 0 : bse_env%eps_nto_osc_str
524 : ELSE
525 2 : WRITE (unit_nr, '(T2,A4,T7,A,T71,A10)') 'BSE|', 'Threshold for oscillator strength f^n', &
526 4 : ADJUSTL("---")
527 : END IF
528 2 : WRITE (unit_nr, '(T2,A4,T7,A,T72,F10.3)') 'BSE|', 'Threshold for NTO weights (λ_I)^2', &
529 4 : bse_env%eps_nto_eigval
530 : END IF
531 :
532 : ! Write the header of NTO info table
533 4 : IF (unit_nr > 0) THEN
534 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
535 2 : IF (.NOT. PRESENT(fm_Y)) THEN
536 1 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
537 2 : 'NTOs from solving the BSE within the TDA:'
538 : ELSE
539 1 : WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
540 2 : 'NTOs from solving the BSE without the TDA:'
541 : END IF
542 2 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
543 2 : WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T33,A14,T62,A)') 'BSE|', &
544 4 : 'Excitation n', "TDA/ABBA", "Index of NTO I", 'NTO weights (λ_I)^2'
545 : END IF
546 :
547 104 : DO j = 1, bse_env%num_print_exc_ntos
548 100 : n_exc = bse_env%bse_nto_state_list_final(j)
549 : ! Takes care of unallocated oscill_str array in case of Triplet
550 100 : IF (bse_env%eps_nto_osc_str > 0.0_dp) THEN
551 : ! Check actual values
552 0 : IF (oscill_str(n_exc) < bse_env%eps_nto_osc_str) THEN
553 : ! Print skipped levels to table
554 0 : IF (unit_nr > 0) THEN
555 0 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
556 0 : WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T42,A39)') 'BSE|', &
557 0 : n_exc, info_approximation, "Skipped (Oscillator strength too small)"
558 : END IF
559 : CYCLE
560 : END IF
561 : END IF
562 :
563 : CALL cp_fm_create(fm_m, matrix_struct=fm_struct_m, &
564 100 : name="single_part_trans_dm")
565 100 : CALL cp_fm_set_all(fm_m, 0.0_dp)
566 :
567 : CALL cp_fm_create(fm_nto_coeff_holes, matrix_struct=fm_struct_nto_holes, &
568 100 : name="nto_coeffs_holes")
569 100 : CALL cp_fm_set_all(fm_nto_coeff_holes, 0.0_dp)
570 :
571 : CALL cp_fm_create(fm_nto_coeff_particles, matrix_struct=fm_struct_nto_particles, &
572 100 : name="nto_coeffs_particles")
573 100 : CALL cp_fm_set_all(fm_nto_coeff_particles, 0.0_dp)
574 :
575 : ! Reshuffle from X_ia,n_exc to X_i,a
576 100 : CALL reshuffle_eigvec(fm_X, fm_X_ia, homo, virtual, n_exc, .FALSE.)
577 :
578 : ! Copy X to upper block in M, i.e. starting from column homo+1
579 : CALL cp_fm_to_fm_submat(fm_X_ia, fm_m, &
580 : homo, virtual, &
581 : 1, 1, &
582 100 : 1, homo + 1)
583 100 : CALL cp_fm_release(fm_X_ia)
584 : ! Copy Y if present
585 100 : IF (PRESENT(fm_Y)) THEN
586 : ! Reshuffle from Y_ia,n_exc to Y_a,i
587 50 : CALL reshuffle_eigvec(fm_Y, fm_Y_ai, homo, virtual, n_exc, .TRUE.)
588 :
589 : ! Copy Y^T to lower block in M, i.e. starting from row homo+1
590 : CALL cp_fm_to_fm_submat(fm_Y_ai, fm_m, &
591 : virtual, homo, &
592 : 1, 1, &
593 50 : homo + 1, 1)
594 :
595 50 : CALL cp_fm_release(fm_Y_ai)
596 :
597 : END IF
598 :
599 : ! Now we compute the SVD of M_{occ+virt,occ+virt}, which yields
600 : ! M = U * Lambda * V^T
601 : ! Initialize matrices and arrays to store left/right eigenvectors and singular values
602 : CALL cp_fm_create(matrix=fm_eigvl, &
603 : matrix_struct=fm_m%matrix_struct, &
604 100 : name="LEFT_SINGULAR_MATRIX")
605 100 : CALL cp_fm_set_all(fm_eigvl, alpha=0.0_dp)
606 : CALL cp_fm_create(matrix=fm_eigvr_t, &
607 : matrix_struct=fm_m%matrix_struct, &
608 100 : name="RIGHT_SINGULAR_MATRIX")
609 100 : CALL cp_fm_set_all(fm_eigvr_t, alpha=0.0_dp)
610 :
611 300 : ALLOCATE (eigval_svd(nao_trunc))
612 1700 : eigval_svd(:) = 0.0_dp
613 : info_svd = 0
614 100 : CALL cp_fm_svd(fm_m, fm_eigvl, fm_eigvr_t, eigval_svd, info_svd)
615 604 : IF (info_svd /= 0) THEN
616 0 : IF (unit_nr > 0) THEN
617 : CALL cp_warn(__LOCATION__, &
618 : "SVD for computation of NTOs not successful. "// &
619 0 : "Skipping print of NTOs.")
620 0 : IF (info_svd > 0) THEN
621 : CALL cp_warn(__LOCATION__, &
622 : "PDGESVD detected heterogeneity. "// &
623 0 : "Decreasing number of MPI ranks might solve this issue.")
624 : END IF
625 : END IF
626 : ! Release matrices to avoid memory leaks
627 0 : CALL cp_fm_release(fm_m)
628 0 : CALL cp_fm_release(fm_nto_coeff_holes)
629 0 : CALL cp_fm_release(fm_nto_coeff_particles)
630 : ELSE
631 : ! Rescale singular values as done in Martin2003 (10.1063/1.1558471)
632 200 : ALLOCATE (eigval_svd_squ(nao_trunc))
633 1700 : eigval_svd_squ(:) = eigval_svd(:)**2
634 : ! Sanity check for TDA: In case of TDA, the sum should be \sum_ia |X_ia|^2 = 1
635 100 : IF (.NOT. PRESENT(fm_Y)) THEN
636 850 : IF (ABS(SUM(eigval_svd_squ) - 1) >= coeff_err) THEN
637 0 : CPWARN("Sum of NTO coefficients deviates from 1!")
638 : END IF
639 : END IF
640 :
641 : ! Create NTO coefficients for later print to grid via TDDFT routine
642 : ! Apply U = fm_eigvl to MO coeffs, which yields hole states
643 : CALL parallel_gemm("N", "N", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvl, 0.0_dp, &
644 100 : fm_nto_coeff_holes)
645 :
646 : ! Apply V^T = fm_eigvr_t to MO coeffs, which yields particle states
647 : CALL parallel_gemm("N", "T", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvr_t, 0.0_dp, &
648 100 : fm_nto_coeff_particles)
649 :
650 : !Release intermediary work matrices
651 100 : CALL cp_fm_release(fm_m)
652 100 : CALL cp_fm_release(fm_eigvl)
653 100 : CALL cp_fm_release(fm_eigvr_t)
654 :
655 : ! Transfer NTO coefficients to sets
656 100 : nto_name(1) = 'Hole_coord'
657 100 : nto_name(2) = 'Particle_coord'
658 300 : ALLOCATE (nto_set(2))
659 : ! Extract number of significant NTOs
660 100 : n_nto = 0
661 320 : DO i_nto = 1, nao_trunc
662 320 : IF (eigval_svd_squ(i_nto) > bse_env%eps_nto_eigval) THEN
663 220 : n_nto = n_nto + 1
664 : ELSE
665 : ! Since svd orders in descending order, we can exit the loop if smaller
666 : EXIT
667 : END IF
668 : END DO
669 :
670 100 : IF (unit_nr > 0) THEN
671 50 : WRITE (unit_nr, '(T2,A4)') 'BSE|'
672 160 : DO i_nto = 1, n_nto
673 110 : WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T41,I6,T71,F10.5)') 'BSE|', &
674 270 : n_exc, info_approximation, i_nto, eigval_svd_squ(i_nto)
675 : END DO
676 : END IF
677 :
678 : CALL cp_fm_struct_create(fm_struct_nto_set, template_fmstruct=fm_struct_nto_holes, &
679 100 : ncol_global=n_nto)
680 100 : CALL cp_fm_create(fm_nto_set, fm_struct_nto_set)
681 300 : DO i = 1, 2
682 200 : CALL allocate_mo_set(nto_set(i), nao_trunc, n_nto, 0, 0.0_dp, 2.0_dp, 0.0_dp)
683 300 : CALL init_mo_set(nto_set(i), fm_ref=fm_nto_set, name=nto_name(i))
684 : END DO
685 100 : CALL cp_fm_release(fm_nto_set)
686 100 : CALL cp_fm_struct_release(fm_struct_nto_set)
687 :
688 : ! Fill NTO sets
689 100 : CALL cp_fm_to_fm(fm_nto_coeff_holes, nto_set(1)%mo_coeff, ncol=n_nto)
690 100 : CALL cp_fm_to_fm(fm_nto_coeff_particles, nto_set(2)%mo_coeff, ncol=n_nto)
691 :
692 : ! Cube files
693 100 : nto_section => section_vals_get_subs_vals(bse_section, "NTO_ANALYSIS")
694 100 : CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
695 100 : CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
696 100 : CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
697 100 : IF (cube_file) THEN
698 : CALL print_bse_nto_cubes(bse_env, nto_set, n_exc, info_approximation, &
699 0 : stride, append_cube, nto_section)
700 : END IF
701 :
702 100 : CALL cp_fm_release(fm_nto_coeff_holes)
703 100 : CALL cp_fm_release(fm_nto_coeff_particles)
704 100 : DEALLOCATE (eigval_svd)
705 100 : DEALLOCATE (eigval_svd_squ)
706 300 : DO i = 1, 2
707 300 : CALL deallocate_mo_set(nto_set(i))
708 : END DO
709 400 : DEALLOCATE (nto_set)
710 : END IF
711 : END DO
712 :
713 4 : CALL cp_fm_release(fm_mo_coeff)
714 4 : CALL cp_fm_struct_release(fm_struct_m)
715 4 : CALL cp_fm_struct_release(fm_struct_nto_holes)
716 4 : CALL cp_fm_struct_release(fm_struct_nto_particles)
717 4 : CALL cp_fm_struct_release(fm_struct_mo_coeff)
718 :
719 4 : CALL timestop(handle)
720 :
721 8 : END SUBROUTINE calculate_NTOs
722 :
723 : ! **************************************************************************************************
724 : !> \brief ...
725 : !> \param exc_descr Allocated and initialized on exit
726 : !> \param fm_X_ia ...
727 : !> \param fm_multipole_ij_trunc ...
728 : !> \param fm_multipole_ab_trunc ...
729 : !> \param fm_multipole_ai_trunc ...
730 : !> \param i_exc ...
731 : !> \param homo ...
732 : !> \param virtual ...
733 : !> \param fm_Y_ia ...
734 : ! **************************************************************************************************
735 110 : SUBROUTINE get_exciton_descriptors(exc_descr, fm_X_ia, &
736 : fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
737 : fm_multipole_ai_trunc, &
738 : i_exc, homo, virtual, &
739 : fm_Y_ia)
740 :
741 : TYPE(exciton_descr_type), ALLOCATABLE, &
742 : DIMENSION(:) :: exc_descr
743 : TYPE(cp_fm_type), INTENT(IN) :: fm_X_ia
744 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
745 : INTENT(IN) :: fm_multipole_ij_trunc, &
746 : fm_multipole_ab_trunc, &
747 : fm_multipole_ai_trunc
748 : INTEGER, INTENT(IN) :: i_exc, homo, virtual
749 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_Y_ia
750 :
751 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_exciton_descriptors'
752 :
753 : INTEGER :: handle, i_dir, j_dir
754 : INTEGER, DIMENSION(3) :: mask_quadrupole
755 : LOGICAL :: flag_TDA
756 : REAL(KIND=dp) :: norm_X, norm_XpY, norm_Y
757 : REAL(KIND=dp), DIMENSION(3) :: r_e_sq_X, r_e_sq_Y, r_e_X, r_e_Y, &
758 : r_h_sq_X, r_h_sq_Y, r_h_X, r_h_Y
759 : REAL(KIND=dp), DIMENSION(3, 3) :: r_e_h_XX, r_e_h_XY, r_e_h_YY
760 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia
761 : TYPE(cp_fm_type) :: fm_work_ba, fm_work_ia, fm_work_ia_2
762 :
763 110 : CALL timeset(routineN, handle)
764 110 : IF (PRESENT(fm_Y_ia)) THEN
765 : flag_TDA = .FALSE.
766 : ELSE
767 60 : flag_TDA = .TRUE.
768 : END IF
769 :
770 : ! translates 1,2,3 to diagonal entries of quadrupoles xx, yy, zz
771 : ! Ordering in quadrupole moments is x, y, z, xx, xy, xz, yy, yz, zz
772 110 : mask_quadrupole = [4, 7, 9]
773 :
774 : CALL cp_fm_struct_create(fm_struct_ia, &
775 110 : context=fm_X_ia%matrix_struct%context, nrow_global=homo, ncol_global=virtual)
776 : CALL cp_fm_struct_create(fm_struct_ab, &
777 110 : context=fm_X_ia%matrix_struct%context, nrow_global=virtual, ncol_global=virtual)
778 :
779 110 : r_e_X(:) = 0.0_dp
780 110 : r_e_Y(:) = 0.0_dp
781 110 : r_h_X(:) = 0.0_dp
782 110 : r_h_Y(:) = 0.0_dp
783 110 : r_e_sq_X(:) = 0.0_dp
784 110 : r_h_sq_X(:) = 0.0_dp
785 110 : r_e_sq_Y(:) = 0.0_dp
786 110 : r_h_sq_Y(:) = 0.0_dp
787 110 : r_e_h_XX(:, :) = 0.0_dp
788 110 : r_e_h_XY(:, :) = 0.0_dp
789 110 : r_e_h_YY(:, :) = 0.0_dp
790 :
791 : norm_X = 0.0_dp
792 : norm_Y = 0.0_dp
793 110 : norm_XpY = 0.0_dp
794 :
795 : ! Initialize values of exciton descriptors
796 440 : exc_descr(i_exc)%r_e(:) = 0.0_dp
797 440 : exc_descr(i_exc)%r_h(:) = 0.0_dp
798 440 : exc_descr(i_exc)%r_e_sq(:) = 0.0_dp
799 440 : exc_descr(i_exc)%r_h_sq(:) = 0.0_dp
800 1430 : exc_descr(i_exc)%r_e_h(:, :) = 0.0_dp
801 :
802 110 : exc_descr(i_exc)%flag_TDA = flag_TDA
803 110 : exc_descr(i_exc)%norm_XpY = 0.0_dp
804 :
805 : ! Norm of X
806 110 : CALL cp_fm_trace(fm_X_ia, fm_X_ia, norm_X)
807 110 : norm_XpY = norm_X
808 : ! Norm of Y
809 110 : IF (.NOT. flag_TDA) THEN
810 50 : CALL cp_fm_trace(fm_Y_ia, fm_Y_ia, norm_Y)
811 50 : norm_XpY = norm_XpY + norm_Y
812 : END IF
813 :
814 110 : exc_descr(i_exc)%norm_XpY = norm_XpY
815 :
816 : ! <r_h>_X = Tr[ X^T µ_ij X + Y µ_ab Y^T ] = X_ai µ_ij X_ja + Y_ia µ_ab Y_bi
817 440 : DO i_dir = 1, 3
818 : ! <r_h>_X = X_ai µ_ij X_ja + ...
819 330 : CALL trace_exciton_descr(fm_X_ia, fm_multipole_ij_trunc(i_dir), fm_X_ia, r_h_X(i_dir))
820 330 : r_h_X(i_dir) = r_h_X(i_dir)/norm_XpY
821 440 : IF (.NOT. flag_TDA) THEN
822 : ! <r_h>_X = ... + Y_ia µ_ab Y_bi
823 150 : CALL trace_exciton_descr(fm_Y_ia, fm_Y_ia, fm_multipole_ab_trunc(i_dir), r_h_Y(i_dir))
824 150 : r_h_Y(i_dir) = r_h_Y(i_dir)/norm_XpY
825 : END IF
826 : END DO
827 440 : exc_descr(i_exc)%r_h(:) = r_h_X(:) + r_h_Y(:)
828 :
829 : ! <r_e>_X = Tr[ X µ_ab X^T + Y^T µ_ij Y ] = X_ia µ_ab X_bi + Y_ai µ_ij Y_ja
830 440 : DO i_dir = 1, 3
831 : ! <r_e>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
832 330 : CALL trace_exciton_descr(fm_X_ia, fm_X_ia, fm_multipole_ab_trunc(i_dir), r_e_X(i_dir))
833 330 : r_e_X(i_dir) = r_e_X(i_dir)/norm_XpY
834 440 : IF (.NOT. flag_TDA) THEN
835 : ! <r_e>_X = ... + Y_ai µ_ij Y_ja
836 150 : CALL trace_exciton_descr(fm_Y_ia, fm_multipole_ij_trunc(i_dir), fm_Y_ia, r_e_Y(i_dir))
837 150 : r_e_Y(i_dir) = r_e_Y(i_dir)/norm_XpY
838 : END IF
839 : END DO
840 440 : exc_descr(i_exc)%r_e(:) = r_e_X(:) + r_e_Y(:)
841 :
842 : ! <r_h^2>_X = Tr[ X^T M_ij X + Y M_ab Y^T ] = X_ai M_ij X_ja + Y_ia M_ab Y_bi
843 440 : DO i_dir = 1, 3
844 : ! <r_h^2>_X = X_ai M_ij X_ja + ...
845 : CALL trace_exciton_descr(fm_X_ia, fm_multipole_ij_trunc(mask_quadrupole(i_dir)), &
846 330 : fm_X_ia, r_h_sq_X(i_dir))
847 330 : r_h_sq_X(i_dir) = r_h_sq_X(i_dir)/norm_XpY
848 440 : IF (.NOT. flag_TDA) THEN
849 : ! <r_h^2>_X = ... + Y_ia M_ab Y_bi
850 : CALL trace_exciton_descr(fm_Y_ia, fm_Y_ia, &
851 150 : fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_h_sq_Y(i_dir))
852 150 : r_h_sq_Y(i_dir) = r_h_sq_Y(i_dir)/norm_XpY
853 : END IF
854 : END DO
855 440 : exc_descr(i_exc)%r_h_sq(:) = r_h_sq_X(:) + r_h_sq_Y(:)
856 :
857 : ! <r_e^2>_X = Tr[ X M_ab X^T + Y^T M_ij Y ] = X_ia M_ab X_bi + Y_ai M_ij Y_ja
858 440 : DO i_dir = 1, 3
859 : ! <r_e^2>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
860 : CALL trace_exciton_descr(fm_X_ia, fm_X_ia, &
861 330 : fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_e_sq_X(i_dir))
862 330 : r_e_sq_X(i_dir) = r_e_sq_X(i_dir)/norm_XpY
863 440 : IF (.NOT. flag_TDA) THEN
864 : ! <r_e^2>_X = ... + Y_ai M_ij Y_ja
865 : CALL trace_exciton_descr(fm_Y_ia, fm_multipole_ij_trunc(mask_quadrupole(i_dir)), &
866 150 : fm_Y_ia, r_e_sq_Y(i_dir))
867 150 : r_e_sq_Y(i_dir) = r_e_sq_Y(i_dir)/norm_XpY
868 : END IF
869 : END DO
870 440 : exc_descr(i_exc)%r_e_sq(:) = r_e_sq_X(:) + r_e_sq_Y(:)
871 :
872 : ! <r_e^\mu r_h^\mu'>_X
873 : ! = Tr[ X^T µ'_ij X µ_ab + Y^T µ_ij Y µ'_ab + 2 X µ_ai Y µ'_ai ]
874 : ! = X_bj µ'_ji X_ia µ_ab + Y_bj µ_ji Y_ia µ'_ab + 2 X_ia µ_aj Y_jb µ'_bi
875 : ! The i_dir and j_dir convert to mu and mu'. µ (electron) sits between a (from X) and j (from Y),
876 : ! µ' (hole) between i (from X) and b (from Y); Tr[ Y µ'_ai X µ_ai ] = Tr[ X µ_ai Y µ'_ai ] by cyclicity, hence the 2.
877 110 : CALL cp_fm_create(fm_work_ia, fm_struct_ia)
878 110 : CALL cp_fm_create(fm_work_ia_2, fm_struct_ia)
879 110 : CALL cp_fm_create(fm_work_ba, fm_struct_ab)
880 440 : DO i_dir = 1, 3
881 1430 : DO j_dir = 1, 3
882 : ! First term - X^T µ'_ij X µ_ab
883 990 : CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
884 990 : CALL cp_fm_set_all(fm_work_ia_2, 0.0_dp)
885 : ! work_ib = X_ia µ_ab
886 : CALL parallel_gemm("N", "N", homo, virtual, virtual, 1.0_dp, &
887 990 : fm_X_ia, fm_multipole_ab_trunc(i_dir), 0.0_dp, fm_work_ia)
888 : ! work_ja_2 = µ'_ji work_ia
889 : CALL parallel_gemm("N", "N", homo, virtual, homo, 1.0_dp, &
890 990 : fm_multipole_ij_trunc(j_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
891 : ! <r_e^\mu r_h^\mu'>_X = work_ia_2 X_bj + ... = X^T work_ia_2 + ...
892 990 : CALL cp_fm_trace(fm_X_ia, fm_work_ia_2, r_e_h_XX(i_dir, j_dir))
893 990 : r_e_h_XX(i_dir, j_dir) = r_e_h_XX(i_dir, j_dir)/norm_XpY
894 1320 : IF (.NOT. flag_TDA) THEN
895 : ! Second term - Y^T µ_ij Y µ'_ab
896 450 : CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
897 450 : CALL cp_fm_set_all(fm_work_ia_2, 0.0_dp)
898 : ! work_ib = Y_ia µ'_ab
899 : CALL parallel_gemm("N", "N", homo, virtual, virtual, 1.0_dp, &
900 450 : fm_Y_ia, fm_multipole_ab_trunc(j_dir), 0.0_dp, fm_work_ia)
901 : ! work_ja_2 = µ_ji work_ia
902 : CALL parallel_gemm("N", "N", homo, virtual, homo, 1.0_dp, &
903 450 : fm_multipole_ij_trunc(i_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
904 : ! <r_h r_e>_X = work_ia_2 Y_bj + ... = Y^T work_ia_2 + ...
905 450 : CALL cp_fm_trace(fm_Y_ia, fm_work_ia_2, r_e_h_YY(i_dir, j_dir))
906 450 : r_e_h_YY(i_dir, j_dir) = r_e_h_YY(i_dir, j_dir)/norm_XpY
907 :
908 : ! Third term (counted twice) - X µ_ai Y µ'_ai = X_ia µ_aj Y_jb µ'_bi
909 : ! Reshuffle for usage of trace (where first argument is transposed)
910 : ! = µ_aj Y_jb µ'_bi X_ia =
911 : ! \___________/
912 : ! fm_work_ai
913 : ! fm_work_ai = µ_aj Y_jb µ'_bi
914 : ! fm_work_ia = µ'_ib Y_bj µ_ja
915 : ! \_____/
916 : ! fm_work_ba
917 450 : CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
918 450 : CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
919 : ! fm_work_ba = Y_bj µ_ja
920 : CALL parallel_gemm("T", "T", virtual, virtual, homo, 1.0_dp, &
921 450 : fm_Y_ia, fm_multipole_ai_trunc(i_dir), 0.0_dp, fm_work_ba)
922 : ! fm_work_ia = µ'_ib fm_work_ba
923 : CALL parallel_gemm("T", "N", homo, virtual, virtual, 1.0_dp, &
924 450 : fm_multipole_ai_trunc(j_dir), fm_work_ba, 0.0_dp, fm_work_ia)
925 : ! <r_e r_h>_X = ... + 2 X_ia µ_aj Y_jb µ'_bi
926 450 : CALL cp_fm_trace(fm_work_ia, fm_X_ia, r_e_h_XY(i_dir, j_dir))
927 450 : r_e_h_XY(i_dir, j_dir) = 2.0_dp*r_e_h_XY(i_dir, j_dir)/norm_XpY
928 : END IF
929 : END DO
930 : END DO
931 1430 : exc_descr(i_exc)%r_e_h(:, :) = r_e_h_XX(:, :) + r_e_h_XY(:, :) + r_e_h_YY(:, :)
932 :
933 110 : CALL cp_fm_release(fm_work_ia)
934 110 : CALL cp_fm_release(fm_work_ia_2)
935 110 : CALL cp_fm_release(fm_work_ba)
936 :
937 : ! Now we compute all the descriptors and correlation coefficients
938 : ! Order is: Directional ones, then covariances and correlation coefficients and
939 :
940 : ! diff_r_abs = |<r_h>_X - <r_e>_X|
941 440 : exc_descr(i_exc)%diff_r_abs = SQRT(SUM((exc_descr(i_exc)%r_h(:) - exc_descr(i_exc)%r_e(:))**2))
942 :
943 : ! σ_e = sqrt( <r_e^2>_X - <r_e>_X^2 )
944 770 : exc_descr(i_exc)%sigma_e = SQRT(SUM(exc_descr(i_exc)%r_e_sq(:)) - SUM(exc_descr(i_exc)%r_e(:)**2))
945 :
946 : ! σ_h = sqrt( <r_h^2>_X - <r_h>_X^2 )
947 770 : exc_descr(i_exc)%sigma_h = SQRT(SUM(exc_descr(i_exc)%r_h_sq(:)) - SUM(exc_descr(i_exc)%r_h(:)**2))
948 :
949 : ! Now directed ones
950 440 : DO i_dir = 1, 3
951 330 : exc_descr(i_exc)%d_eh_dir(i_dir) = ABS(exc_descr(i_exc)%r_h(i_dir) - exc_descr(i_exc)%r_e(i_dir))
952 330 : exc_descr(i_exc)%sigma_e_dir(i_dir) = SQRT(exc_descr(i_exc)%r_e_sq(i_dir) - exc_descr(i_exc)%r_e(i_dir)**2)
953 440 : exc_descr(i_exc)%sigma_h_dir(i_dir) = SQRT(exc_descr(i_exc)%r_h_sq(i_dir) - exc_descr(i_exc)%r_h(i_dir)**2)
954 : END DO
955 :
956 : ! Covariance and correlation coefficient (as well as crosscorrelation matrices)
957 : ! COV(r_e, r_h) = < r_e r_h >_X - < r_e >_X < r_h >_X
958 110 : exc_descr(i_exc)%cov_e_h_sum = 0.0_dp
959 1430 : exc_descr(i_exc)%cov_e_h(:, :) = 0.0_dp
960 1430 : exc_descr(i_exc)%corr_e_h_matrix(:, :) = 0.0_dp
961 440 : DO i_dir = 1, 3
962 1320 : DO j_dir = 1, 3
963 : exc_descr(i_exc)%cov_e_h(i_dir, j_dir) = exc_descr(i_exc)%r_e_h(i_dir, j_dir) &
964 990 : - exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(j_dir)
965 : exc_descr(i_exc)%corr_e_h_matrix(i_dir, j_dir) = &
966 : exc_descr(i_exc)%cov_e_h(i_dir, j_dir)/ &
967 1320 : (exc_descr(i_exc)%sigma_e_dir(i_dir)*exc_descr(i_exc)%sigma_h_dir(j_dir))
968 : END DO
969 : exc_descr(i_exc)%cov_e_h_sum = exc_descr(i_exc)%cov_e_h_sum + &
970 : exc_descr(i_exc)%r_e_h(i_dir, i_dir) - &
971 440 : exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(i_dir)
972 : END DO
973 :
974 : ! e-h-correlation coefficient R_eh = COV(r_e, r_h) / ( σ_e σ_h )
975 110 : exc_descr(i_exc)%corr_e_h = exc_descr(i_exc)%cov_e_h_sum/(exc_descr(i_exc)%sigma_e*exc_descr(i_exc)%sigma_h)
976 :
977 : ! root-mean-square e-h separation
978 : exc_descr(i_exc)%diff_r_sqr = SQRT(exc_descr(i_exc)%diff_r_abs**2 + &
979 : exc_descr(i_exc)%sigma_e**2 + exc_descr(i_exc)%sigma_h**2 &
980 110 : - 2*exc_descr(i_exc)%cov_e_h_sum)
981 :
982 440 : DO i_dir = 1, 3
983 : exc_descr(i_exc)%d_exc_dir(i_dir) = SQRT(exc_descr(i_exc)%d_eh_dir(i_dir)**2 + &
984 : exc_descr(i_exc)%sigma_e_dir(i_dir)**2 + &
985 : exc_descr(i_exc)%sigma_h_dir(i_dir)**2 - &
986 440 : 2*exc_descr(i_exc)%cov_e_h(i_dir, i_dir))
987 : END DO
988 :
989 : ! Expectation values of r_e and r_h
990 440 : exc_descr(i_exc)%r_e_shift(:) = exc_descr(i_exc)%r_e(:)
991 440 : exc_descr(i_exc)%r_h_shift(:) = exc_descr(i_exc)%r_h(:)
992 :
993 110 : CALL cp_fm_struct_release(fm_struct_ia)
994 110 : CALL cp_fm_struct_release(fm_struct_ab)
995 :
996 110 : CALL timestop(handle)
997 :
998 330 : END SUBROUTINE get_exciton_descriptors
999 :
1000 0 : END MODULE bse_properties
|
