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