(git:ed6f26b)
Loading...
Searching...
No Matches
bse_properties.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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_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_h = 0.0_dp, &
73 r_e_shift = 0.0_dp, &
74 r_h_shift = 0.0_dp, &
75 cov_e_h = 0.0_dp
76 REAL(kind=dp) :: sigma_e = 0.0_dp, &
77 sigma_h = 0.0_dp, &
78 cov_e_h_sum = 0.0_dp, &
79 corr_e_h = 0.0_dp, &
80 diff_r_abs = 0.0_dp, &
81 diff_r_sqr = 0.0_dp, &
82 norm_xpy = 0.0_dp
83 LOGICAL :: flag_tda = .false.
84 END TYPE exciton_descr_type
85
86CONTAINS
87
88! **************************************************************************************************
89!> \brief Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n )
90!> and oscillator strengths f^n = 2/3 * Ω^n sum_r∈(x,y,z) ( d_r^n )^2
91!> Prelim Ref.: Eqs. (23), (24)
92!> in J. Chem. Phys. 152, 044105 (2020); https://doi.org/10.1063/1.5123290
93!> \param fm_eigvec_X ...
94!> \param Exc_ens ...
95!> \param fm_dipole_ai_trunc ...
96!> \param trans_mom_bse BSE dipole vectors in real space per excitation level
97!> \param oscill_str Oscillator strength per excitation level
98!> \param polarizability_residues Residues of polarizability ("tensorial oscillator strength")
99!> per excitation level
100!> \param mp2_env ...
101!> \param homo_red ...
102!> \param virtual_red ...
103!> \param unit_nr ...
104!> \param fm_eigvec_Y ...
105! **************************************************************************************************
106 SUBROUTINE get_oscillator_strengths(fm_eigvec_X, Exc_ens, fm_dipole_ai_trunc, &
107 trans_mom_bse, oscill_str, polarizability_residues, &
108 mp2_env, homo_red, virtual_red, unit_nr, &
109 fm_eigvec_Y)
110
111 TYPE(cp_fm_type), INTENT(IN) :: fm_eigvec_x
112 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
113 INTENT(IN) :: exc_ens
114 TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_ai_trunc
115 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
116 INTENT(OUT) :: trans_mom_bse
117 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
118 INTENT(OUT) :: oscill_str
119 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
120 INTENT(OUT) :: polarizability_residues
121 TYPE(mp2_type), INTENT(IN) :: mp2_env
122 INTEGER, INTENT(IN) :: homo_red, virtual_red, unit_nr
123 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_eigvec_y
124
125 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_oscillator_strengths'
126
127 INTEGER :: handle, idir, jdir, n
128 TYPE(cp_fm_struct_type), POINTER :: fm_struct_dipole_mo_trunc_reordered, &
129 fm_struct_trans_mom_bse
130 TYPE(cp_fm_type) :: fm_eigvec_xysum
131 TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_mo_trunc_reordered, &
132 fm_dipole_per_dir, fm_trans_mom_bse
133
134 CALL timeset(routinen, handle)
135
136 CALL cp_fm_struct_create(fm_struct_dipole_mo_trunc_reordered, fm_eigvec_x%matrix_struct%para_env, &
137 fm_eigvec_x%matrix_struct%context, 1, homo_red*virtual_red)
138 CALL cp_fm_struct_create(fm_struct_trans_mom_bse, fm_eigvec_x%matrix_struct%para_env, &
139 fm_eigvec_x%matrix_struct%context, 1, homo_red*virtual_red)
140
141 ! Include excitonic amplitudes in dipoles, i.e. obtain "BSE dipoles":
142 ! \vec{D}_n = sqrt(2) * sum_{i,a} \vec{D}_ai (X_{ai}^{(n)} + Y_{ai}^{(n)})
143
144 ! Reorder dipoles in order to execute the sum over i and a by parallel gemm
145 DO idir = 1, 3
146 CALL cp_fm_create(fm_dipole_mo_trunc_reordered(idir), matrix_struct=fm_struct_dipole_mo_trunc_reordered, &
147 name="dipoles_mo_reordered")
148 CALL cp_fm_set_all(fm_dipole_mo_trunc_reordered(idir), 0.0_dp)
149 CALL fm_general_add_bse(fm_dipole_mo_trunc_reordered(idir), fm_dipole_ai_trunc(idir), 1.0_dp, &
150 1, 1, &
151 1, virtual_red, &
152 unit_nr, (/2, 4, 3, 1/), mp2_env)
153 CALL cp_fm_release(fm_dipole_per_dir(idir))
154 END DO
155
156 DO idir = 1, 3
157 CALL cp_fm_create(fm_trans_mom_bse(idir), matrix_struct=fm_struct_trans_mom_bse, &
158 name="excitonic_dipoles")
159 CALL cp_fm_set_all(fm_trans_mom_bse(idir), 0.0_dp)
160 END DO
161
162 ! If TDA is invoked, Y is not present as it is simply 0
163 CALL cp_fm_create(fm_eigvec_xysum, matrix_struct=fm_eigvec_x%matrix_struct, name="excit_amplitude_sum")
164 CALL cp_fm_set_all(fm_eigvec_xysum, 0.0_dp)
165 CALL cp_fm_to_fm(fm_eigvec_x, fm_eigvec_xysum)
166 IF (PRESENT(fm_eigvec_y)) THEN
167 CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_xysum, 1.0_dp, fm_eigvec_y)
168 END IF
169 DO idir = 1, 3
170 CALL parallel_gemm('N', 'N', 1, homo_red*virtual_red, homo_red*virtual_red, sqrt(2.0_dp), &
171 fm_dipole_mo_trunc_reordered(idir), fm_eigvec_xysum, 0.0_dp, fm_trans_mom_bse(idir))
172 END DO
173
174 ! Get oscillator strengths themselves
175 ALLOCATE (oscill_str(homo_red*virtual_red))
176 ALLOCATE (trans_mom_bse(3, 1, homo_red*virtual_red))
177 ALLOCATE (polarizability_residues(3, 3, homo_red*virtual_red))
178 trans_mom_bse(:, :, :) = 0.0_dp
179
180 ! Sum over all directions
181 DO idir = 1, 3
182 CALL cp_fm_get_submatrix(fm_trans_mom_bse(idir), trans_mom_bse(idir, :, :))
183 END DO
184
185 DO n = 1, homo_red*virtual_red
186 DO idir = 1, 3
187 DO jdir = 1, 3
188 polarizability_residues(idir, jdir, n) = 2.0_dp*exc_ens(n)*trans_mom_bse(idir, 1, n)*trans_mom_bse(jdir, 1, n)
189 END DO
190 END DO
191 oscill_str(n) = 2.0_dp/3.0_dp*exc_ens(n)*sum(abs(trans_mom_bse(:, :, n))**2)
192 END DO
193
194 CALL cp_fm_struct_release(fm_struct_dipole_mo_trunc_reordered)
195 CALL cp_fm_struct_release(fm_struct_trans_mom_bse)
196 DO idir = 1, 3
197 CALL cp_fm_release(fm_dipole_mo_trunc_reordered(idir))
198 CALL cp_fm_release(fm_trans_mom_bse(idir))
199 CALL cp_fm_release(fm_dipole_ai_trunc(idir))
200 END DO
201 CALL cp_fm_release(fm_eigvec_xysum)
202
203 CALL timestop(handle)
204
205 END SUBROUTINE get_oscillator_strengths
206
207! **************************************************************************************************
208!> \brief Computes and returns absorption spectrum for the frequency range and broadening
209!> provided by the user.
210!> Prelim Ref.: C. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications
211!> (Oxford University Press, Oxford, 2012), Eq. 7.51
212!> \param oscill_str ...
213!> \param polarizability_residues ...
214!> \param Exc_ens ...
215!> \param info_approximation ...
216!> \param unit_nr ...
217!> \param mp2_env ...
218! **************************************************************************************************
219 SUBROUTINE compute_absorption_spectrum(oscill_str, polarizability_residues, Exc_ens, &
220 info_approximation, unit_nr, mp2_env)
221
222 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
223 INTENT(IN) :: oscill_str
224 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :), &
225 INTENT(IN) :: polarizability_residues
226 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
227 INTENT(IN) :: exc_ens
228 CHARACTER(LEN=10) :: info_approximation
229 INTEGER, INTENT(IN) :: unit_nr
230 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
231
232 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_absorption_spectrum'
233
234 CHARACTER(LEN=10) :: eta_str, width_eta_format_str
235 CHARACTER(LEN=40) :: file_name_crosssection, &
236 file_name_spectrum
237 INTEGER :: handle, i, idir, j, jdir, k, num_steps, &
238 unit_nr_file, width_eta
239 REAL(kind=dp) :: eta, freq_end, freq_start, freq_step, &
240 omega
241 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: abs_cross_section, abs_spectrum
242 REAL(kind=dp), DIMENSION(:), POINTER :: eta_list
243 TYPE(cp_logger_type), POINTER :: logger
244
245 CALL timeset(routinen, handle)
246
247 freq_step = mp2_env%bse%bse_spectrum_freq_step_size
248 freq_start = mp2_env%bse%bse_spectrum_freq_start
249 freq_end = mp2_env%bse%bse_spectrum_freq_end
250 eta_list => mp2_env%bse%bse_eta_spectrum_list
251
252 ! Calculate number of steps to fit given frequency range
253 num_steps = nint((freq_end - freq_start)/freq_step) + 1
254
255 DO k = 1, SIZE(eta_list)
256 eta = eta_list(k)
257
258 ! Some magic to get a nice formatting of the eta value in filenames
259 width_eta = max(1, int(log10(eta)) + 1) + 4
260 WRITE (width_eta_format_str, "(A2,I0,A3)") '(F', width_eta, '.3)'
261 WRITE (eta_str, width_eta_format_str) eta*evolt
262 ! Filename itself
263 file_name_spectrum = 'BSE'//trim(adjustl(info_approximation))//'eta='//trim(eta_str)//'.spectrum'
264 file_name_crosssection = 'BSE'//trim(adjustl(info_approximation))//'eta='//trim(eta_str)//'.crosssection'
265
266 ! First column is frequency in eV, second column is imaginary part of the trace of the polarizability
267 ! The following 9 columns are the entries of the polarizability tensor
268 ALLOCATE (abs_spectrum(num_steps, 11))
269 abs_spectrum(:, :) = 0.0_dp
270 ! Also calculate and print the photoabsorption cross section tensor
271 ! σ_{µ µ'}(ω) = 4πω Im{α_{µ µ'}(ω)} / c
272 ALLOCATE (abs_cross_section(num_steps, 11))
273 abs_cross_section(:, :) = 0.0_dp
274
275 ! Calculate the imaginary part of the mean dipole polarizability α_{avg}(ω)
276 ! which is given by (cf. C. Ullrichs Book on TDDFT, Eq. 7.51)
277 ! We introduce an additional - due to his convention for charge vs particle density, see also:
278 ! Computer Physics Communications, 208:149–161, November 2016
279 ! https://doi.org/10.1016/j.cpc.2016.06.019
280 ! α_{avg}(ω) = - \sum_{n=1}^{N_exc} \frac{f_n}{(ω+iη)² - (Ω^n)²}
281 ! and then the imaginary part is (in the limit η -> 0)
282 ! Im[α_{avg}(ω)] = - \sum_{n=1}^{N_exc} f_n * η / ((ω - Ω^n)² + η²)
283 ! where f_n are the oscillator strengths and E_exc the excitation energies
284 ! For the full polarizability tensor, we have
285 ! α_{µ µ'}(ω) = - \sum_n [2 Ω^n d^n_µ d^n_µ'] / [(ω+iη)^2- (Ω^n)^2]
286 ! = - \sum_n "polarizability_residues" / [(ω+iη)^2- (Ω^n)^2]
287 DO i = 1, num_steps
288 omega = freq_start + (i - 1)*freq_step
289 abs_spectrum(i, 1) = omega
290 DO j = 1, SIZE(oscill_str)
291 abs_spectrum(i, 2) = abs_spectrum(i, 2) - oscill_str(j)* &
292 aimag(1/((omega + cmplx(0.0, eta, kind=dp))**2 - exc_ens(j)**2))
293 DO idir = 1, 3
294 DO jdir = 1, 3
295 ! Factor 2 from formula for tensor is already in the polarizability_residues
296 ! to follow the same convention as the oscillator strengths
297 abs_spectrum(i, 2 + (idir - 1)*3 + jdir) = abs_spectrum(i, 2 + (idir - 1)*3 + jdir) &
298 - polarizability_residues(idir, jdir, j)* &
299 aimag(1/((omega + cmplx(0.0, eta, kind=dp))**2 - exc_ens(j)**2))
300 END DO
301 END DO
302 END DO
303 END DO
304
305 ! Extract cross section σ from polarizability tensor
306 DO i = 1, num_steps
307 omega = abs_spectrum(i, 1)
308 abs_cross_section(i, 1) = omega
309 abs_cross_section(i, 2:) = 4.0_dp*pi*abs_spectrum(i, 2:)*omega/c_light_au
310 END DO
311
312 !For debug runs: Export an entry of the two tensors to allow regtests on spectra
313 IF (mp2_env%bse%bse_debug_print) THEN
314 IF (unit_nr > 0) THEN
315 WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
316 'Averaged dynamical dipole polarizability at 8.2 eV:', &
317 abs_spectrum(83, 2)
318 WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
319 'Averaged photoabsorption cross section at 8.2 eV:', &
320 abs_cross_section(83, 2)
321 END IF
322 END IF
323
324 ! Print it to file
325 logger => cp_get_default_logger()
326 IF (logger%para_env%is_source()) THEN
327 unit_nr_file = cp_logger_get_default_unit_nr()
328 ELSE
329 unit_nr_file = -1
330 END IF
331
332 IF (unit_nr_file > 0) THEN
333 CALL open_file(file_name_crosssection, unit_number=unit_nr_file, &
334 file_status="UNKNOWN", file_action="WRITE")
335 WRITE (unit_nr_file, '(A,A6)') σµµωπω"# Photoabsorption cross section _{ '}() = -4/c * Im[ \sum_n "// &
336 Ωµµωη²Ω²"[2 ^n d_^n d_'^n] / [(+i)- (^n)] ] from Bethe Salpeter equation for method ", &
337 trim(adjustl(info_approximation))
338 WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", σω"_{avg}()", σω"_xx()", &
339 σω"_xy()", σω"_xz()", σω"_yx()", σω"_yy()", σω"_yz()", σω"_zx()", &
340 σω"_zy()", σω"_zz()"
341 DO i = 1, num_steps
342 WRITE (unit_nr_file, '(11(F20.8,1X))') abs_cross_section(i, 1)*evolt, abs_cross_section(i, 2:11)
343 END DO
344 CALL close_file(unit_nr_file)
345 END IF
346 DEALLOCATE (abs_cross_section)
347
348 IF (unit_nr_file > 0) THEN
349 CALL open_file(file_name_spectrum, unit_number=unit_nr_file, &
350 file_status="UNKNOWN", file_action="WRITE")
351 WRITE (unit_nr_file, '(A,A6)') αµµω"# Imaginary part of polarizability _{ '}() = -\sum_n "// &
352 Ωµµωη²Ω²"[2 ^n d_^n d_'^n] / [(+i)- (^n)] from Bethe Salpeter equation for method ", &
353 trim(adjustl(info_approximation))
354 WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", αω"Im{_{avg}()}", αω"Im{_xx()}", &
355 αω"Im{_xy()}", αω"Im{_xz()}", αω"Im{_yx()}", αω"Im{_yy()}", αω"Im{_yz()}", αω"Im{_zx()}", &
356 αω"Im{_zy()}", αω"Im{_zz()}"
357 DO i = 1, num_steps
358 WRITE (unit_nr_file, '(11(F20.8,1X))') abs_spectrum(i, 1)*evolt, abs_spectrum(i, 2:11)
359 END DO
360 CALL close_file(unit_nr_file)
361 END IF
362 DEALLOCATE (abs_spectrum)
363 END DO
364
365 IF (unit_nr > 0) THEN
366 WRITE (unit_nr, '(T2,A4)') 'BSE|'
367 WRITE (unit_nr, '(T2,A4,T7,A,A)') &
368 'BSE|', "Printed optical absorption spectrum to local files, e.g. "
369 WRITE (unit_nr, '(T2,A4,T7,A)') &
370 'BSE|', file_name_spectrum
371 WRITE (unit_nr, '(T2,A4,T7,A,A)') &
372 'BSE|', "as well as photoabsorption cross section to, e.g. "
373 WRITE (unit_nr, '(T2,A4,T7,A)') &
374 'BSE|', file_name_crosssection
375 WRITE (unit_nr, '(T2,A4,T7,A52)') &
376 'BSE|', "using the Eq. (7.51) from C. Ullrichs Book on TDDFT:"
377 WRITE (unit_nr, '(T2,A4)') 'BSE|'
378 WRITE (unit_nr, '(T2,A4,T10,A75)') &
379 'BSE|', αωωη²Ω²"Im{_{avg}()} = -Im{\sum_{n=1}^{N_exc} \frac{f_n}{(+i) - (^n)}}"
380 WRITE (unit_nr, '(T2,A4)') 'BSE|'
381 WRITE (unit_nr, '(T2,A4,T7,A)') &
382 'BSE|', "or for the full polarizability tensor:"
383 WRITE (unit_nr, '(T2,A4,T10,A)') &
384 'BSE|', αµµωΩµµωη²Ω²"_{ '}() = -\sum_n [2 ^n d_^n d_'^n] / [(+i)- (^n)]"
385 WRITE (unit_nr, '(T2,A4)') 'BSE|'
386 WRITE (unit_nr, '(T2,A4,T7,A)') &
387 'BSE|', "as well as Eq. (7.48):"
388 WRITE (unit_nr, '(T2,A4,T10,A)') &
389 'BSE|', σµµωπωαµµω"_{ '}() = 4 Im{_{ '}()} / c"
390 WRITE (unit_nr, '(T2,A4)') 'BSE|'
391 WRITE (unit_nr, '(T2,A4,T7,A)') &
392 'BSE|', µ"with transition moments d_^n, oscillator strengths f_n,"
393 WRITE (unit_nr, '(T2,A4,T7,A)') &
394 'BSE|', Ω"excitation energies ^n and the speed of light c."
395 WRITE (unit_nr, '(T2,A4)') 'BSE|'
396 WRITE (unit_nr, '(T2,A4,T7,A)') &
397 'BSE|', "Please note that we adopt an additional minus sign for both quantities,"
398 WRITE (unit_nr, '(T2,A4,T7,A)') &
399 'BSE|', "due to the convention for charge vs particle density as done in MolGW:"
400 WRITE (unit_nr, '(T2,A4,T7,A)') &
401 'BSE|', "https://doi.org/10.1016/j.cpc.2016.06.019."
402 WRITE (unit_nr, '(T2,A4)') 'BSE|'
403 END IF
404
405 CALL timestop(handle)
406
407 END SUBROUTINE
408
409! **************************************************************************************************
410!> \brief ...
411!> \param fm_X ...
412!> \param fm_Y ...
413!> \param mo_coeff ...
414!> \param homo ...
415!> \param virtual ...
416!> \param info_approximation ...
417!> \param oscill_str ...
418!> \param qs_env ...
419!> \param unit_nr ...
420!> \param mp2_env ...
421! **************************************************************************************************
422 SUBROUTINE calculate_ntos(fm_X, fm_Y, &
423 mo_coeff, homo, virtual, &
424 info_approximation, &
425 oscill_str, &
426 qs_env, unit_nr, mp2_env)
427
428 TYPE(cp_fm_type), INTENT(IN) :: fm_x
429 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_y
430 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
431 INTEGER, INTENT(IN) :: homo, virtual
432 CHARACTER(LEN=10) :: info_approximation
433 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: oscill_str
434 TYPE(qs_environment_type), POINTER :: qs_env
435 INTEGER, INTENT(IN) :: unit_nr
436 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
437
438 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_NTOs'
439
440 CHARACTER(LEN=20), DIMENSION(2) :: nto_name
441 INTEGER :: handle, homo_irred, i, i_nto, info_svd, &
442 j, n_exc, n_nto, nao_full, nao_trunc
443 INTEGER, DIMENSION(:), POINTER :: stride
444 LOGICAL :: append_cube, cube_file
445 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval_svd_squ
446 REAL(kind=dp), DIMENSION(:), POINTER :: eigval_svd
447 TYPE(cp_fm_struct_type), POINTER :: fm_struct_m, fm_struct_mo_coeff, &
448 fm_struct_nto_holes, &
449 fm_struct_nto_particles, &
450 fm_struct_nto_set
451 TYPE(cp_fm_type) :: fm_eigvl, fm_eigvr_t, fm_m, fm_mo_coeff, fm_nto_coeff_holes, &
452 fm_nto_coeff_particles, fm_nto_set, fm_x_ia, fm_y_ai
453 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
454 TYPE(section_vals_type), POINTER :: bse_section, input, nto_section
455
456 CALL timeset(routinen, handle)
457 CALL get_qs_env(qs_env=qs_env, &
458 input=input)
459 bse_section => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%GW%BSE")
460
461 nao_full = qs_env%mos(1)%nao
462 nao_trunc = homo + virtual
463 ! This is not influenced by the BSE cutoff
464 homo_irred = qs_env%mos(1)%homo
465 ! M will have a block structure and is quadratic in homo+virtual, i.e.
466 ! occ virt
467 ! | 0 X_i,a | occ = homo
468 ! M = | Y_a,i 0 | virt = virtual
469 !
470 ! X and Y are here not the eigenvectors X_ia,n - instead we fix n and reshape the combined ia index
471 ! Notice the index structure of the lower block, i.e. X is transposed
472 CALL cp_fm_struct_create(fm_struct_m, &
473 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
474 nao_trunc, nao_trunc)
475 CALL cp_fm_struct_create(fm_struct_mo_coeff, &
476 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
477 nao_full, nao_trunc)
478 CALL cp_fm_struct_create(fm_struct_nto_holes, &
479 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
480 nao_full, nao_trunc)
481 CALL cp_fm_struct_create(fm_struct_nto_particles, &
482 fm_x%matrix_struct%para_env, fm_x%matrix_struct%context, &
483 nao_full, nao_trunc)
484
485 CALL cp_fm_create(fm_mo_coeff, matrix_struct=fm_struct_mo_coeff, &
486 name="mo_coeff")
487 ! Here, we take care of possible cutoffs
488 ! Simply truncating the matrix causes problems with the print function
489 ! Therefore, we keep the dimension, but set the coefficients of truncated indices to 0
490 CALL cp_fm_to_fm_submat_general(mo_coeff(1), fm_mo_coeff, &
491 nao_full, nao_trunc, &
492 1, homo_irred - homo + 1, &
493 1, 1, &
494 mo_coeff(1)%matrix_struct%context)
495
496 ! Print some information about the NTOs
497 IF (unit_nr > 0) THEN
498 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
499 φχ'The Natural Transition Orbital (NTO) pairs _I(r_e) and _I(r_h) for a fixed'
500 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
501 'excitation index n are obtained by singular value decomposition of T'
502 WRITE (unit_nr, '(T2,A4)') 'BSE|'
503 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
504 ' = (0 X)'
505 IF (PRESENT(fm_y)) THEN
506 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
507 'T = (Y^T 0)'
508 ELSE
509 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
510 'T = (0 0)'
511 END IF
512 WRITE (unit_nr, '(T2,A4)') 'BSE|'
513 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
514 Λ'T = U V^T'
515 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
516 φψ'_I(r_e) = \sum_p V_pI _p(r_e)'
517 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
518 χψ'_I(r_h) = \sum_p U_pI _p(r_e)'
519 WRITE (unit_nr, '(T2,A4)') 'BSE|'
520 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
521 'where we have introduced'
522 WRITE (unit_nr, '(T2,A4)') 'BSE|'
523 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
524 'BSE|', ψ"_p:", "occupied and virtual molecular orbitals,"
525 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
526 'BSE|', φ"_I(r_e):", "NTO state for the electron,"
527 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
528 'BSE|', χ"_I(r_h):", "NTO state for the hole,"
529 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
530 'BSE|', Λ":", λ"diagonal matrix of NTO weights _I,"
531 WRITE (unit_nr, '(T2,A4)') 'BSE|'
532 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
533 "The NTOs are calculated with the following settings:"
534 WRITE (unit_nr, '(T2,A4)') 'BSE|'
535 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Number of excitations, for which NTOs are computed', &
536 mp2_env%bse%num_print_exc_ntos
537 IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp) THEN
538 WRITE (unit_nr, '(T2,A4,T7,A,T71,F10.3)') 'BSE|', 'Threshold for oscillator strength f^n', &
539 mp2_env%bse%eps_nto_osc_str
540 ELSE
541 WRITE (unit_nr, '(T2,A4,T7,A,T71,A10)') 'BSE|', 'Threshold for oscillator strength f^n', &
542 adjustl("---")
543 END IF
544 WRITE (unit_nr, '(T2,A4,T7,A,T72,F10.3)') 'BSE|', λ'Threshold for NTO weights (_I)^2', &
545 mp2_env%bse%eps_nto_eigval
546 END IF
547
548 ! Write the header of NTO info table
549 IF (unit_nr > 0) THEN
550 WRITE (unit_nr, '(T2,A4)') 'BSE|'
551 IF (.NOT. PRESENT(fm_y)) THEN
552 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
553 'NTOs from solving the BSE within the TDA:'
554 ELSE
555 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
556 'NTOs from solving the BSE without the TDA:'
557 END IF
558 WRITE (unit_nr, '(T2,A4)') 'BSE|'
559 WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T33,A14,T62,A)') 'BSE|', &
560 'Excitation n', "TDA/ABBA", "Index of NTO I", λ'NTO weights (_I)^2'
561 END IF
562
563 DO j = 1, mp2_env%bse%num_print_exc_ntos
564 n_exc = mp2_env%bse%bse_nto_state_list_final(j)
565 ! Takes care of unallocated oscill_str array in case of Triplet
566 IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp) THEN
567 ! Check actual values
568 IF (oscill_str(n_exc) < mp2_env%bse%eps_nto_osc_str) THEN
569 ! Print skipped levels to table
570 IF (unit_nr > 0) THEN
571 WRITE (unit_nr, '(T2,A4)') 'BSE|'
572 WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T42,A39)') 'BSE|', &
573 n_exc, info_approximation, "Skipped (Oscillator strength too small)"
574 END IF
575 cycle
576 END IF
577 END IF
578
579 CALL cp_fm_create(fm_m, matrix_struct=fm_struct_m, &
580 name="single_part_trans_dm")
581 CALL cp_fm_set_all(fm_m, 0.0_dp)
582
583 CALL cp_fm_create(fm_nto_coeff_holes, matrix_struct=fm_struct_nto_holes, &
584 name="nto_coeffs_holes")
585 CALL cp_fm_set_all(fm_nto_coeff_holes, 0.0_dp)
586
587 CALL cp_fm_create(fm_nto_coeff_particles, matrix_struct=fm_struct_nto_particles, &
588 name="nto_coeffs_particles")
589 CALL cp_fm_set_all(fm_nto_coeff_particles, 0.0_dp)
590
591 ! Reshuffle from X_ia,n_exc to X_i,a
592 CALL reshuffle_eigvec(fm_x, fm_x_ia, homo, virtual, n_exc, &
593 .false., unit_nr, mp2_env)
594
595 ! Copy X to upper block in M, i.e. starting from column homo+1
596 CALL cp_fm_to_fm_submat(fm_x_ia, fm_m, &
597 homo, virtual, &
598 1, 1, &
599 1, homo + 1)
600 CALL cp_fm_release(fm_x_ia)
601 ! Copy Y if present
602 IF (PRESENT(fm_y)) THEN
603 ! Reshuffle from Y_ia,n_exc to Y_a,i
604 CALL reshuffle_eigvec(fm_y, fm_y_ai, homo, virtual, n_exc, &
605 .true., unit_nr, mp2_env)
606
607 ! Copy Y^T to lower block in M, i.e. starting from row homo+1
608 CALL cp_fm_to_fm_submat(fm_y_ai, fm_m, &
609 virtual, homo, &
610 1, 1, &
611 homo + 1, 1)
612
613 CALL cp_fm_release(fm_y_ai)
614
615 END IF
616
617 ! Now we compute the SVD of M_{occ+virt,occ+virt}, which yields
618 ! M = U * Lambda * V^T
619 ! Initialize matrices and arrays to store left/right eigenvectors and singular values
620 CALL cp_fm_create(matrix=fm_eigvl, &
621 matrix_struct=fm_m%matrix_struct, &
622 name="LEFT_SINGULAR_MATRIX")
623 CALL cp_fm_set_all(fm_eigvl, alpha=0.0_dp)
624 CALL cp_fm_create(matrix=fm_eigvr_t, &
625 matrix_struct=fm_m%matrix_struct, &
626 name="RIGHT_SINGULAR_MATRIX")
627 CALL cp_fm_set_all(fm_eigvr_t, alpha=0.0_dp)
628
629 ALLOCATE (eigval_svd(nao_trunc))
630 eigval_svd(:) = 0.0_dp
631 info_svd = 0
632 CALL cp_fm_svd(fm_m, fm_eigvl, fm_eigvr_t, eigval_svd, info_svd)
633 IF (info_svd /= 0) THEN
634 IF (unit_nr > 0) THEN
635 CALL cp_warn(__location__, &
636 "SVD for computation of NTOs not successful. "// &
637 "Skipping print of NTOs.")
638 IF (info_svd > 0) THEN
639 CALL cp_warn(__location__, &
640 "PDGESVD detected heterogeneity. "// &
641 "Decreasing number of MPI ranks might solve this issue.")
642 END IF
643 END IF
644 ! Release matrices to avoid memory leaks
645 CALL cp_fm_release(fm_m)
646 CALL cp_fm_release(fm_nto_coeff_holes)
647 CALL cp_fm_release(fm_nto_coeff_particles)
648 ELSE
649 ! Rescale singular values as done in Martin2003 (10.1063/1.1558471)
650 ALLOCATE (eigval_svd_squ(nao_trunc))
651 eigval_svd_squ(:) = eigval_svd(:)**2
652 ! Sanity check for TDA: In case of TDA, the sum should be \sum_ia |X_ia|^2 = 1
653 IF (.NOT. PRESENT(fm_y)) THEN
654 IF (abs(sum(eigval_svd_squ) - 1) >= 1e-5) THEN
655 cpwarn("Sum of NTO coefficients deviates from 1!")
656 END IF
657 END IF
658
659 ! Create NTO coefficients for later print to grid via TDDFT routine
660 ! Apply U = fm_eigvl to MO coeffs, which yields hole states
661 CALL parallel_gemm("N", "N", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvl, 0.0_dp, &
662 fm_nto_coeff_holes)
663
664 ! Apply V^T = fm_eigvr_t to MO coeffs, which yields particle states
665 CALL parallel_gemm("N", "T", nao_full, nao_trunc, nao_trunc, 1.0_dp, fm_mo_coeff, fm_eigvr_t, 0.0_dp, &
666 fm_nto_coeff_particles)
667
668 !Release intermediary work matrices
669 CALL cp_fm_release(fm_m)
670 CALL cp_fm_release(fm_eigvl)
671 CALL cp_fm_release(fm_eigvr_t)
672
673 ! Transfer NTO coefficients to sets
674 nto_name(1) = 'Hole_coord'
675 nto_name(2) = 'Particle_coord'
676 ALLOCATE (nto_set(2))
677 ! Extract number of significant NTOs
678 n_nto = 0
679 DO i_nto = 1, nao_trunc
680 IF (eigval_svd_squ(i_nto) > mp2_env%bse%eps_nto_eigval) THEN
681 n_nto = n_nto + 1
682 ELSE
683 ! Since svd orders in descending order, we can exit the loop if smaller
684 EXIT
685 END IF
686 END DO
687
688 IF (unit_nr > 0) THEN
689 WRITE (unit_nr, '(T2,A4)') 'BSE|'
690 DO i_nto = 1, n_nto
691 WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T41,I6,T71,F10.5)') 'BSE|', &
692 n_exc, info_approximation, i_nto, eigval_svd_squ(i_nto)
693 END DO
694 END IF
695
696 CALL cp_fm_struct_create(fm_struct_nto_set, template_fmstruct=fm_struct_nto_holes, &
697 ncol_global=n_nto)
698 CALL cp_fm_create(fm_nto_set, fm_struct_nto_set)
699 DO i = 1, 2
700 CALL allocate_mo_set(nto_set(i), nao_trunc, n_nto, 0, 0.0_dp, 2.0_dp, 0.0_dp)
701 CALL init_mo_set(nto_set(i), fm_ref=fm_nto_set, name=nto_name(i))
702 END DO
703 CALL cp_fm_release(fm_nto_set)
704 CALL cp_fm_struct_release(fm_struct_nto_set)
705
706 ! Fill NTO sets
707 CALL cp_fm_to_fm(fm_nto_coeff_holes, nto_set(1)%mo_coeff, ncol=n_nto)
708 CALL cp_fm_to_fm(fm_nto_coeff_particles, nto_set(2)%mo_coeff, ncol=n_nto)
709
710 ! Cube files
711 nto_section => section_vals_get_subs_vals(bse_section, "NTO_ANALYSIS")
712 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
713 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
714 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
715 IF (cube_file) THEN
716 CALL print_bse_nto_cubes(qs_env, nto_set, n_exc, info_approximation, &
717 stride, append_cube, nto_section)
718 END IF
719
720 CALL cp_fm_release(fm_nto_coeff_holes)
721 CALL cp_fm_release(fm_nto_coeff_particles)
722 DEALLOCATE (eigval_svd)
723 DEALLOCATE (eigval_svd_squ)
724 DO i = 1, 2
725 CALL deallocate_mo_set(nto_set(i))
726 END DO
727 DEALLOCATE (nto_set)
728 END IF
729 END DO
730
731 CALL cp_fm_release(fm_mo_coeff)
732 CALL cp_fm_struct_release(fm_struct_m)
733 CALL cp_fm_struct_release(fm_struct_nto_holes)
734 CALL cp_fm_struct_release(fm_struct_nto_particles)
735 CALL cp_fm_struct_release(fm_struct_mo_coeff)
736
737 CALL timestop(handle)
738
739 END SUBROUTINE calculate_ntos
740
741! **************************************************************************************************
742!> \brief ...
743!> \param exc_descr Allocated and initialized on exit
744!> \param fm_X_ia ...
745!> \param fm_multipole_ij_trunc ...
746!> \param fm_multipole_ab_trunc ...
747!> \param fm_multipole_ai_trunc ...
748!> \param i_exc ...
749!> \param homo ...
750!> \param virtual ...
751!> \param fm_Y_ia ...
752! **************************************************************************************************
753 SUBROUTINE get_exciton_descriptors(exc_descr, fm_X_ia, &
754 fm_multipole_ij_trunc, fm_multipole_ab_trunc, &
755 fm_multipole_ai_trunc, &
756 i_exc, homo, virtual, &
757 fm_Y_ia)
758
759 TYPE(exciton_descr_type), ALLOCATABLE, &
760 DIMENSION(:) :: exc_descr
761 TYPE(cp_fm_type), INTENT(IN) :: fm_x_ia
762 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
763 INTENT(IN) :: fm_multipole_ij_trunc, &
764 fm_multipole_ab_trunc, &
765 fm_multipole_ai_trunc
766 INTEGER, INTENT(IN) :: i_exc, homo, virtual
767 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_y_ia
768
769 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_exciton_descriptors'
770
771 INTEGER :: handle, i_dir
772 INTEGER, DIMENSION(3) :: mask_quadrupole
773 LOGICAL :: flag_tda
774 REAL(kind=dp) :: norm_x, norm_xpy, norm_y
775 REAL(kind=dp), DIMENSION(3) :: r_e_h_xx, r_e_h_xy, r_e_h_yx, r_e_h_yy, &
776 r_e_sq_x, r_e_sq_y, r_e_x, r_e_y, &
777 r_h_sq_x, r_h_sq_y, r_h_x, r_h_y
778 TYPE(cp_fm_struct_type), POINTER :: fm_struct_ab, fm_struct_ia
779 TYPE(cp_fm_type) :: fm_work_ba, fm_work_ia, fm_work_ia_2
780
781 CALL timeset(routinen, handle)
782 IF (PRESENT(fm_y_ia)) THEN
783 flag_tda = .false.
784 ELSE
785 flag_tda = .true.
786 END IF
787
788 ! translates 1,2,3 to diagonal entries of quadrupoles xx, yy, zz
789 ! Ordering in quadrupole moments is x, y, z, xx, xy, xz, yy, yz, zz
790 mask_quadrupole = (/4, 7, 9/)
791
792 CALL cp_fm_struct_create(fm_struct_ia, &
793 context=fm_x_ia%matrix_struct%context, nrow_global=homo, ncol_global=virtual)
794 CALL cp_fm_struct_create(fm_struct_ab, &
795 context=fm_x_ia%matrix_struct%context, nrow_global=virtual, ncol_global=virtual)
796
797 r_e_x(:) = 0.0_dp
798 r_e_y(:) = 0.0_dp
799 r_h_x(:) = 0.0_dp
800 r_h_y(:) = 0.0_dp
801 r_e_sq_x(:) = 0.0_dp
802 r_h_sq_x(:) = 0.0_dp
803 r_e_sq_y(:) = 0.0_dp
804 r_h_sq_y(:) = 0.0_dp
805 r_e_h_xx(:) = 0.0_dp
806 r_e_h_xy(:) = 0.0_dp
807 r_e_h_yx(:) = 0.0_dp
808 r_e_h_yy(:) = 0.0_dp
809
810 norm_x = 0.0_dp
811 norm_y = 0.0_dp
812 norm_xpy = 0.0_dp
813
814 ! Initialize values of exciton descriptors
815 exc_descr(i_exc)%r_e(:) = 0.0_dp
816 exc_descr(i_exc)%r_h(:) = 0.0_dp
817 exc_descr(i_exc)%r_e_sq(:) = 0.0_dp
818 exc_descr(i_exc)%r_h_sq(:) = 0.0_dp
819 exc_descr(i_exc)%r_e_h(:) = 0.0_dp
820
821 exc_descr(i_exc)%flag_TDA = flag_tda
822 exc_descr(i_exc)%norm_XpY = 0.0_dp
823
824 ! Norm of X
825 CALL cp_fm_trace(fm_x_ia, fm_x_ia, norm_x)
826 norm_xpy = norm_x
827 ! Norm of Y
828 IF (.NOT. flag_tda) THEN
829 CALL cp_fm_trace(fm_y_ia, fm_y_ia, norm_y)
830 norm_xpy = norm_xpy + norm_y
831 END IF
832
833 exc_descr(i_exc)%norm_XpY = norm_xpy
834
835 ! <r_h>_X = Tr[ X^T µ_ij X + Y µ_ab Y^T ] = X_ai µ_ij X_ja + Y_ia µ_ab Y_bi
836 DO i_dir = 1, 3
837 ! <r_h>_X = X_ai µ_ij X_ja + ...
838 CALL trace_exciton_descr(fm_x_ia, fm_multipole_ij_trunc(i_dir), fm_x_ia, r_h_x(i_dir))
839 r_h_x(i_dir) = r_h_x(i_dir)/norm_xpy
840 IF (.NOT. flag_tda) THEN
841 ! <r_h>_X = ... + Y_ia µ_ab Y_bi
842 CALL trace_exciton_descr(fm_y_ia, fm_y_ia, fm_multipole_ab_trunc(i_dir), r_h_y(i_dir))
843 r_h_y(i_dir) = r_h_y(i_dir)/norm_xpy
844 END IF
845 END DO
846 exc_descr(i_exc)%r_h(:) = r_h_x(:) + r_h_y(:)
847
848 ! <r_e>_X = Tr[ X µ_ab X^T + Y^T µ_ij Y ] = X_ia µ_ab X_bi + Y_ai µ_ij Y_ja
849 DO i_dir = 1, 3
850 ! <r_e>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
851 CALL trace_exciton_descr(fm_x_ia, fm_x_ia, fm_multipole_ab_trunc(i_dir), r_e_x(i_dir))
852 r_e_x(i_dir) = r_e_x(i_dir)/norm_xpy
853 IF (.NOT. flag_tda) THEN
854 ! <r_e>_X = ... + Y_ai µ_ij Y_ja
855 CALL trace_exciton_descr(fm_y_ia, fm_multipole_ij_trunc(i_dir), fm_y_ia, r_e_y(i_dir))
856 r_e_y(i_dir) = r_e_y(i_dir)/norm_xpy
857 END IF
858 END DO
859 exc_descr(i_exc)%r_e(:) = r_e_x(:) + r_e_y(:)
860
861 ! <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
862 DO i_dir = 1, 3
863 ! <r_h^2>_X = X_ai M_ij X_ja + ...
864 CALL trace_exciton_descr(fm_x_ia, fm_multipole_ij_trunc(mask_quadrupole(i_dir)), &
865 fm_x_ia, r_h_sq_x(i_dir))
866 r_h_sq_x(i_dir) = r_h_sq_x(i_dir)/norm_xpy
867 IF (.NOT. flag_tda) THEN
868 ! <r_h^2>_X = ... + Y_ia M_ab Y_bi
869 CALL trace_exciton_descr(fm_y_ia, fm_y_ia, &
870 fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_h_sq_y(i_dir))
871 r_h_sq_y(i_dir) = r_h_sq_y(i_dir)/norm_xpy
872 END IF
873 END DO
874 exc_descr(i_exc)%r_h_sq(:) = r_h_sq_x(:) + r_h_sq_y(:)
875
876 ! <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
877 DO i_dir = 1, 3
878 ! <r_e^2>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
879 CALL trace_exciton_descr(fm_x_ia, fm_x_ia, &
880 fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_e_sq_x(i_dir))
881 r_e_sq_x(i_dir) = r_e_sq_x(i_dir)/norm_xpy
882 IF (.NOT. flag_tda) THEN
883 ! <r_e^2>_X = ... + Y_ai M_ij Y_ja
884 CALL trace_exciton_descr(fm_y_ia, fm_multipole_ij_trunc(mask_quadrupole(i_dir)), &
885 fm_y_ia, r_e_sq_y(i_dir))
886 r_e_sq_y(i_dir) = r_e_sq_y(i_dir)/norm_xpy
887 END IF
888 END DO
889 exc_descr(i_exc)%r_e_sq(:) = r_e_sq_x(:) + r_e_sq_y(:)
890
891 ! <r_h r_e>_X
892 ! = Tr[ X^T µ_ij X µ_ab + Y^T µ_ij Y µ_ab + X µ_ai Y µ_ai + Y µ_ai X µ_ai]
893 ! = X_bj µ_ji X_ia µ_ab + Y_bj µ_ji Y_ia µ_ab + X_ia µ_aj Y_jb µ_bi + Y_ia µ_aj X_jb µ_bi
894 CALL cp_fm_create(fm_work_ia, fm_struct_ia)
895 CALL cp_fm_create(fm_work_ia_2, fm_struct_ia)
896 CALL cp_fm_create(fm_work_ba, fm_struct_ab)
897 DO i_dir = 1, 3
898 ! First term - X^T µ_ij X µ_ab
899 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
900 CALL cp_fm_set_all(fm_work_ia_2, 0.0_dp)
901 ! work_ib = X_ia µ_ab
902 CALL parallel_gemm("N", "N", homo, virtual, virtual, 1.0_dp, &
903 fm_x_ia, fm_multipole_ab_trunc(i_dir), 0.0_dp, fm_work_ia)
904 ! work_ja_2 = µ_ji work_ia
905 CALL parallel_gemm("N", "N", homo, virtual, homo, 1.0_dp, &
906 fm_multipole_ij_trunc(i_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
907 ! <r_h r_e>_X = work_ia_2 X_bj + ... = X^T work_ia_2 + ...
908 CALL cp_fm_trace(fm_x_ia, fm_work_ia_2, r_e_h_xx(i_dir))
909 r_e_h_xx(i_dir) = r_e_h_xx(i_dir)/norm_xpy
910 IF (.NOT. flag_tda) THEN
911 ! Second term - Y^T µ_ij Y µ_ab
912 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
913 CALL cp_fm_set_all(fm_work_ia_2, 0.0_dp)
914 ! work_ib = Y_ia µ_ab
915 CALL parallel_gemm("N", "N", homo, virtual, virtual, 1.0_dp, &
916 fm_y_ia, fm_multipole_ab_trunc(i_dir), 0.0_dp, fm_work_ia)
917 ! work_ja_2 = µ_ji work_ia
918 CALL parallel_gemm("N", "N", homo, virtual, homo, 1.0_dp, &
919 fm_multipole_ij_trunc(i_dir), fm_work_ia, 0.0_dp, fm_work_ia_2)
920 ! <r_h r_e>_X = work_ia_2 Y_bj + ... = Y^T work_ia_2 + ...
921 CALL cp_fm_trace(fm_y_ia, fm_work_ia_2, r_e_h_yy(i_dir))
922 r_e_h_yy(i_dir) = r_e_h_yy(i_dir)/norm_xpy
923
924 ! Third term - X µ_ai Y µ_ai = X_ia µ_aj Y_jb µ_bi
925 ! Reshuffle for usage of trace (where first argument is transposed)
926 ! = µ_aj Y_jb µ_bi X_ia =
927 ! \__________/
928 ! fm_work_ai
929 ! fm_work_ai = µ_aj Y_jb µ_bi
930 ! fm_work_ia = µ_ib Y_bj µ_ja
931 ! \_____/
932 ! fm_work_ba
933 CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
934 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
935 ! fm_work_ba = Y_bj µ_ja
936 CALL parallel_gemm("T", "T", virtual, virtual, homo, 1.0_dp, &
937 fm_y_ia, fm_multipole_ai_trunc(i_dir), 0.0_dp, fm_work_ba)
938 ! fm_work_ia = µ_ib fm_work_ba
939 CALL parallel_gemm("T", "N", homo, virtual, virtual, 1.0_dp, &
940 fm_multipole_ai_trunc(i_dir), fm_work_ba, 0.0_dp, fm_work_ia)
941 ! <r_h r_e>_X = ... + X_ia µ_aj Y_jb µ_bi
942 CALL cp_fm_trace(fm_work_ia, fm_x_ia, r_e_h_xy(i_dir))
943 r_e_h_xy(i_dir) = r_e_h_xy(i_dir)/norm_xpy
944
945 ! Fourth term - Y µ_ai X µ_ai = Y_ia µ_aj X_jb µ_bi
946 ! Reshuffle for usage of trace (where first argument is transposed)
947 ! = µ_aj X_jb µ_bi Y_ia =
948 ! \__________/
949 ! fm_work_ai
950 ! fm_work_ai = µ_aj X_jb µ_bi
951 ! fm_work_ia = µ_ib X_bj µ_ja
952 ! \_____/
953 ! fm_work_ba
954 CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
955 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
956 ! fm_work_ba = Y_bj µ_ja
957 CALL parallel_gemm("T", "T", virtual, virtual, homo, 1.0_dp, &
958 fm_x_ia, fm_multipole_ai_trunc(i_dir), 0.0_dp, fm_work_ba)
959 ! fm_work_ia = µ_ib fm_work_ba
960 CALL parallel_gemm("T", "N", homo, virtual, virtual, 1.0_dp, &
961 fm_multipole_ai_trunc(i_dir), fm_work_ba, 0.0_dp, fm_work_ia)
962 ! <r_h r_e>_X = ... + X_ia µ_aj Y_jb µ_bi
963 CALL cp_fm_trace(fm_work_ia, fm_y_ia, r_e_h_yx(i_dir))
964 r_e_h_yx(i_dir) = r_e_h_yx(i_dir)/norm_xpy
965 END IF
966 END DO
967 exc_descr(i_exc)%r_e_h(:) = r_e_h_xx(:) + r_e_h_xy(:) + r_e_h_yx(:) + r_e_h_yy(:)
968
969 CALL cp_fm_release(fm_work_ia)
970 CALL cp_fm_release(fm_work_ia_2)
971 CALL cp_fm_release(fm_work_ba)
972
973 ! diff_r_abs = |<r_h>_X - <r_e>_X|
974 exc_descr(i_exc)%diff_r_abs = sqrt(sum((exc_descr(i_exc)%r_h(:) - exc_descr(i_exc)%r_e(:))**2))
975
976 ! σ_e = sqrt( <r_e^2>_X - <r_e>_X^2 )
977 exc_descr(i_exc)%sigma_e = sqrt(sum(exc_descr(i_exc)%r_e_sq(:)) - sum(exc_descr(i_exc)%r_e(:)**2))
978
979 ! σ_h = sqrt( <r_h^2>_X - <r_h>_X^2 )
980 exc_descr(i_exc)%sigma_h = sqrt(sum(exc_descr(i_exc)%r_h_sq(:)) - sum(exc_descr(i_exc)%r_h(:)**2))
981
982 ! COV(r_e, r_h) = < r_e r_h >_X - < r_e >_X < r_h >_X
983 exc_descr(i_exc)%cov_e_h_sum = 0.0_dp
984 exc_descr(i_exc)%cov_e_h(:) = 0.0_dp
985 DO i_dir = 1, 3
986 exc_descr(i_exc)%cov_e_h(i_dir) = exc_descr(i_exc)%r_e_h(i_dir) &
987 - exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(i_dir)
988 exc_descr(i_exc)%cov_e_h_sum = exc_descr(i_exc)%cov_e_h_sum + &
989 exc_descr(i_exc)%r_e_h(i_dir) - exc_descr(i_exc)%r_e(i_dir)*exc_descr(i_exc)%r_h(i_dir)
990 END DO
991
992 ! root-mean-square e-h separation
993 exc_descr(i_exc)%diff_r_sqr = sqrt(exc_descr(i_exc)%diff_r_abs**2 + &
994 exc_descr(i_exc)%sigma_e**2 + exc_descr(i_exc)%sigma_h**2 &
995 - 2*exc_descr(i_exc)%cov_e_h_sum)
996
997 ! e-h-correlation coefficient R_eh = COV(r_e, r_h) / ( σ_e σ_h )
998 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)
999
1000 ! Expectation values of r_e and r_h
1001 exc_descr(i_exc)%r_e_shift(:) = exc_descr(i_exc)%r_e(:)
1002 exc_descr(i_exc)%r_h_shift(:) = exc_descr(i_exc)%r_h(:)
1003
1004 CALL cp_fm_struct_release(fm_struct_ia)
1005 CALL cp_fm_struct_release(fm_struct_ab)
1006
1007 CALL timestop(handle)
1008
1009 END SUBROUTINE get_exciton_descriptors
1010
1011END MODULE bse_properties
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:75
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:82
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
bse_properties
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
Definition bse_properties.F:13
bse_properties::compute_absorption_spectrum
subroutine, public compute_absorption_spectrum(oscill_str, polarizability_residues, exc_ens, info_approximation, unit_nr, mp2_env)
Computes and returns absorption spectrum for the frequency range and broadening provided by the user....
Definition bse_properties.F:221
bse_properties::get_oscillator_strengths
subroutine, public get_oscillator_strengths(fm_eigvec_x, exc_ens, fm_dipole_ai_trunc, trans_mom_bse, oscill_str, polarizability_residues, mp2_env, homo_red, virtual_red, unit_nr, fm_eigvec_y)
Compute and return BSE dipoles d_r^n = sqrt(2) sum_ia < ψ_i | r | ψ_a > ( X_ia^n + Y_ia^n ) and oscil...
Definition bse_properties.F:110
bse_properties::get_exciton_descriptors
subroutine, public get_exciton_descriptors(exc_descr, fm_x_ia, fm_multipole_ij_trunc, fm_multipole_ab_trunc, fm_multipole_ai_trunc, i_exc, homo, virtual, fm_y_ia)
...
Definition bse_properties.F:758
bse_properties::calculate_ntos
subroutine, public calculate_ntos(fm_x, fm_y, mo_coeff, homo, virtual, info_approximation, oscill_str, qs_env, unit_nr, mp2_env)
...
Definition bse_properties.F:427
bse_util
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
bse_util::trace_exciton_descr
subroutine, public trace_exciton_descr(fm_a, fm_b, fm_c, alpha)
Computes trace of form Tr{A^T B C} for exciton descriptors.
Definition bse_util.F:1893
bse_util::reshuffle_eigvec
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, mp2_env)
...
Definition bse_util.F:1351
bse_util::print_bse_nto_cubes
subroutine, public print_bse_nto_cubes(qs_env, mos, istate, info_approximation, stride, append_cube, print_section)
Borrowed from the tddfpt module with slight adaptions.
Definition bse_util.F:1425
bse_util::fm_general_add_bse
subroutine, public fm_general_add_bse(fm_out, fm_in, beta, nrow_secidx_in, ncol_secidx_in, nrow_secidx_out, ncol_secidx_out, unit_nr, reordering, mp2_env)
Adds and reorders full matrices with a combined index structure, e.g. adding W_ij,...
Definition bse_util.F:199
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_svd
subroutine, public cp_fm_svd(matrix_a, matrix_eigvl, matrix_eigvr_t, eigval, info)
decomposes a quadratic matrix into its singular value decomposition
Definition cp_fm_diag.F:916
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_to_fm_submat_general
subroutine, public cp_fm_to_fm_submat_general(source, destination, nrows, ncols, s_firstrow, s_firstcol, d_firstrow, d_firstcol, global_context)
General copy of a submatrix of fm matrix to a submatrix of another fm matrix. The two matrices can ha...
Definition cp_fm_types.F:2002
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1472
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:894
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::c_light_au
real(kind=dp), parameter, public c_light_au
Definition physcon.F:90
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:206
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition qs_mo_types.F:352
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:245
bse_properties::exciton_descr_type
Definition bse_properties.F:67
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
mp2_types::mp2_type
Definition mp2_types.F:345
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46