(git:da6e80d)
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! **************************************************************************************************
14 USE bse_util, ONLY: fm_general_add_bse,&
18 USE cp_files, ONLY: close_file,&
22 USE cp_fm_diag, ONLY: cp_fm_svd
26 USE cp_fm_types, ONLY: cp_fm_create,&
40 USE kinds, ONLY: dp
41 USE mathconstants, ONLY: pi
42 USE mp2_types, ONLY: mp2_type
44 USE physcon, ONLY: c_light_au,&
45 evolt
48 USE qs_mo_types, ONLY: allocate_mo_set,&
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
64
65! TYPE definitions for exciton wavefunction descriptors
66
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
91CONTAINS
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 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 TYPE(cp_fm_type), DIMENSION(3) :: fm_dipole_mo_trunc_reordered, &
137 fm_dipole_per_dir, fm_trans_mom_bse
138
139 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 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 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 DO idir = 1, 3
151 CALL cp_fm_create(fm_dipole_mo_trunc_reordered(idir), matrix_struct=fm_struct_dipole_mo_trunc_reordered, &
152 name="dipoles_mo_reordered")
153 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 unit_nr, (/2, 4, 3, 1/), mp2_env)
158 CALL cp_fm_release(fm_dipole_per_dir(idir))
159 END DO
160
161 DO idir = 1, 3
162 CALL cp_fm_create(fm_trans_mom_bse(idir), matrix_struct=fm_struct_trans_mom_bse, &
163 name="excitonic_dipoles")
164 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 CALL cp_fm_create(fm_eigvec_xysum, matrix_struct=fm_eigvec_x%matrix_struct, name="excit_amplitude_sum")
169 CALL cp_fm_set_all(fm_eigvec_xysum, 0.0_dp)
170 CALL cp_fm_to_fm(fm_eigvec_x, fm_eigvec_xysum)
171 IF (PRESENT(fm_eigvec_y)) THEN
172 CALL cp_fm_scale_and_add(1.0_dp, fm_eigvec_xysum, 1.0_dp, fm_eigvec_y)
173 END IF
174 DO idir = 1, 3
175 CALL parallel_gemm('N', 'N', 1, homo_red*virtual_red, homo_red*virtual_red, sqrt(2.0_dp), &
176 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 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 ALLOCATE (trans_mom_bse(3, 1, homo_red*virtual_red))
184 ALLOCATE (polarizability_residues(3, 3, homo_red*virtual_red))
185 trans_mom_bse(:, :, :) = 0.0_dp
186
187 ! Sum over all directions
188 DO idir = 1, 3
189 CALL cp_fm_get_submatrix(fm_trans_mom_bse(idir), trans_mom_bse(idir, :, :))
190 END DO
191
192 DO n = 1, homo_red*virtual_red
193 DO idir = 1, 3
194 DO jdir = 1, 3
195 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 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 CALL cp_fm_struct_release(fm_struct_dipole_mo_trunc_reordered)
202 CALL cp_fm_struct_release(fm_struct_trans_mom_bse)
203 DO idir = 1, 3
204 CALL cp_fm_release(fm_dipole_mo_trunc_reordered(idir))
205 CALL cp_fm_release(fm_trans_mom_bse(idir))
206 CALL cp_fm_release(fm_dipole_ai_trunc(idir))
207 END DO
208 CALL cp_fm_release(fm_eigvec_xysum)
209
210 CALL timestop(handle)
211
212 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: abs_cross_section, abs_spectrum
249 REAL(kind=dp), DIMENSION(:), POINTER :: eta_list
250 TYPE(cp_logger_type), POINTER :: logger
251
252 CALL timeset(routinen, handle)
253
254 freq_step = mp2_env%bse%bse_spectrum_freq_step_size
255 freq_start = mp2_env%bse%bse_spectrum_freq_start
256 freq_end = mp2_env%bse%bse_spectrum_freq_end
257 eta_list => mp2_env%bse%bse_eta_spectrum_list
258
259 ! Calculate number of steps to fit given frequency range
260 num_steps = nint((freq_end - freq_start)/freq_step) + 1
261
262 DO k = 1, SIZE(eta_list)
263 eta = eta_list(k)
264
265 ! Some magic to get a nice formatting of the eta value in filenames
266 width_eta = max(1, int(log10(eta)) + 1) + 4
267 WRITE (width_eta_format_str, "(A2,I0,A3)") '(F', width_eta, '.3)'
268 WRITE (eta_str, width_eta_format_str) eta*evolt
269 ! Filename itself
270 file_name_spectrum = 'BSE'//trim(adjustl(info_approximation))//'eta='//trim(eta_str)//'.spectrum'
271 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 ALLOCATE (abs_spectrum(num_steps, 11))
276 abs_spectrum(:, :) = 0.0_dp
277 ! Also calculate and print the photoabsorption cross section tensor
278 ! σ_{µ µ'}(ω) = 4πω Im{α_{µ µ'}(ω)} / c
279 ALLOCATE (abs_cross_section(num_steps, 11))
280 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 DO i = 1, num_steps
295 omega = freq_start + (i - 1)*freq_step
296 abs_spectrum(i, 1) = omega
297 DO j = 1, SIZE(oscill_str)
298 abs_spectrum(i, 2) = abs_spectrum(i, 2) - oscill_str(j)* &
299 aimag(1/((omega + cmplx(0.0, eta, kind=dp))**2 - exc_ens(j)**2))
300 DO idir = 1, 3
301 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 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 DO i = 1, num_steps
314 omega = abs_spectrum(i, 1)
315 abs_cross_section(i, 1) = omega
316 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 IF (mp2_env%bse%bse_debug_print) THEN
321 IF (unit_nr > 0) THEN
322 WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
323 'Averaged dynamical dipole polarizability at 8.2 eV:', &
324 abs_spectrum(83, 2)
325 WRITE (unit_nr, '(T2,A10,T13,A,T65,F16.4)') 'BSE|DEBUG|', &
326 'Averaged photoabsorption cross section at 8.2 eV:', &
327 abs_cross_section(83, 2)
328 END IF
329 END IF
330
331 ! Print it to file
332 logger => cp_get_default_logger()
333 IF (logger%para_env%is_source()) THEN
334 unit_nr_file = cp_logger_get_default_unit_nr()
335 ELSE
336 unit_nr_file = -1
337 END IF
338
339 IF (unit_nr_file > 0) THEN
340 CALL open_file(file_name_crosssection, unit_number=unit_nr_file, &
341 file_status="UNKNOWN", file_action="WRITE")
342 WRITE (unit_nr_file, '(A,A6)') σµµωπω"# Photoabsorption cross section _{ '}() = -4/c * Im[ \sum_n "// &
343 Ωµµωη²Ω²"[2 ^n d_^n d_'^n] / [(+i)- (^n)] ] from Bethe Salpeter equation for method ", &
344 trim(adjustl(info_approximation))
345 WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", σω"_{avg}()", σω"_xx()", &
346 σω"_xy()", σω"_xz()", σω"_yx()", σω"_yy()", σω"_yz()", σω"_zx()", &
347 σω"_zy()", σω"_zz()"
348 DO i = 1, num_steps
349 WRITE (unit_nr_file, '(11(F20.8,1X))') abs_cross_section(i, 1)*evolt, abs_cross_section(i, 2:11)
350 END DO
351 CALL close_file(unit_nr_file)
352 END IF
353 DEALLOCATE (abs_cross_section)
354
355 IF (unit_nr_file > 0) THEN
356 CALL open_file(file_name_spectrum, unit_number=unit_nr_file, &
357 file_status="UNKNOWN", file_action="WRITE")
358 WRITE (unit_nr_file, '(A,A6)') αµµω"# Imaginary part of polarizability _{ '}() = -\sum_n "// &
359 Ωµµωη²Ω²"[2 ^n d_^n d_'^n] / [(+i)- (^n)] from Bethe Salpeter equation for method ", &
360 trim(adjustl(info_approximation))
361 WRITE (unit_nr_file, '(A20,1X,10(2X,A20,1X))') "# Frequency (eV)", αω"Im{_{avg}()}", αω"Im{_xx()}", &
362 αω"Im{_xy()}", αω"Im{_xz()}", αω"Im{_yx()}", αω"Im{_yy()}", αω"Im{_yz()}", αω"Im{_zx()}", &
363 αω"Im{_zy()}", αω"Im{_zz()}"
364 DO i = 1, num_steps
365 WRITE (unit_nr_file, '(11(F20.8,1X))') abs_spectrum(i, 1)*evolt, abs_spectrum(i, 2:11)
366 END DO
367 CALL close_file(unit_nr_file)
368 END IF
369 DEALLOCATE (abs_spectrum)
370 END DO
371
372 IF (unit_nr > 0) THEN
373 WRITE (unit_nr, '(T2,A4)') 'BSE|'
374 WRITE (unit_nr, '(T2,A4,T7,A,A)') &
375 'BSE|', "Printed optical absorption spectrum to local files, e.g. "
376 WRITE (unit_nr, '(T2,A4,T7,A)') &
377 'BSE|', file_name_spectrum
378 WRITE (unit_nr, '(T2,A4,T7,A,A)') &
379 'BSE|', "as well as photoabsorption cross section to, e.g. "
380 WRITE (unit_nr, '(T2,A4,T7,A)') &
381 'BSE|', file_name_crosssection
382 WRITE (unit_nr, '(T2,A4,T7,A52)') &
383 'BSE|', "using the Eq. (7.51) from C. Ullrichs Book on TDDFT:"
384 WRITE (unit_nr, '(T2,A4)') 'BSE|'
385 WRITE (unit_nr, '(T2,A4,T10,A75)') &
386 'BSE|', αωωη²Ω²"Im{_{avg}()} = -Im{\sum_{n=1}^{N_exc} \frac{f_n}{(+i) - (^n)}}"
387 WRITE (unit_nr, '(T2,A4)') 'BSE|'
388 WRITE (unit_nr, '(T2,A4,T7,A)') &
389 'BSE|', "or for the full polarizability tensor:"
390 WRITE (unit_nr, '(T2,A4,T10,A)') &
391 'BSE|', αµµωΩµµωη²Ω²"_{ '}() = -\sum_n [2 ^n d_^n d_'^n] / [(+i)- (^n)]"
392 WRITE (unit_nr, '(T2,A4)') 'BSE|'
393 WRITE (unit_nr, '(T2,A4,T7,A)') &
394 'BSE|', "as well as Eq. (7.48):"
395 WRITE (unit_nr, '(T2,A4,T10,A)') &
396 'BSE|', σµµωπωαµµω"_{ '}() = 4 Im{_{ '}()} / c"
397 WRITE (unit_nr, '(T2,A4)') 'BSE|'
398 WRITE (unit_nr, '(T2,A4,T7,A)') &
399 'BSE|', µ"with transition moments d_^n, oscillator strengths f_n,"
400 WRITE (unit_nr, '(T2,A4,T7,A)') &
401 'BSE|', Ω"excitation energies ^n and the speed of light c."
402 WRITE (unit_nr, '(T2,A4)') 'BSE|'
403 WRITE (unit_nr, '(T2,A4,T7,A)') &
404 'BSE|', "Please note that we adopt an additional minus sign for both quantities,"
405 WRITE (unit_nr, '(T2,A4,T7,A)') &
406 'BSE|', "due to the convention for charge vs particle density as done in MolGW:"
407 WRITE (unit_nr, '(T2,A4,T7,A)') &
408 'BSE|', "https://doi.org/10.1016/j.cpc.2016.06.019."
409 WRITE (unit_nr, '(T2,A4)') 'BSE|'
410 END IF
411
412 CALL timestop(handle)
413
414 END SUBROUTINE
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 SUBROUTINE calculate_ntos(fm_X, fm_Y, &
430 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 INTEGER, DIMENSION(:), POINTER :: stride
451 LOGICAL :: append_cube, cube_file
452 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigval_svd_squ
453 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 TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: nto_set
461 TYPE(section_vals_type), POINTER :: bse_section, input, nto_section
462
463 CALL timeset(routinen, handle)
464 CALL get_qs_env(qs_env=qs_env, &
465 input=input)
466 bse_section => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%GW%BSE")
467
468 nao_full = qs_env%mos(1)%nao
469 nao_trunc = homo + virtual
470 ! This is not influenced by the BSE cutoff
471 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 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 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 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 nao_full, nao_trunc)
491
492 CALL cp_fm_create(fm_mo_coeff, matrix_struct=fm_struct_mo_coeff, &
493 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 mo_coeff(1)%matrix_struct%context)
502
503 ! Print some information about the NTOs
504 IF (unit_nr > 0) THEN
505 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
506 φχ'The Natural Transition Orbital (NTO) pairs _I(r_e) and _I(r_h) for a fixed'
507 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
508 'excitation index n are obtained by singular value decomposition of T'
509 WRITE (unit_nr, '(T2,A4)') 'BSE|'
510 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
511 ' = (0 X)'
512 IF (PRESENT(fm_y)) THEN
513 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
514 'T = (Y^T 0)'
515 ELSE
516 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
517 'T = (0 0)'
518 END IF
519 WRITE (unit_nr, '(T2,A4)') 'BSE|'
520 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
521 Λ'T = U V^T'
522 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
523 φψ'_I(r_e) = \sum_p V_pI _p(r_e)'
524 WRITE (unit_nr, '(T2,A4,T15,A)') 'BSE|', &
525 χψ'_I(r_h) = \sum_p U_pI _p(r_e)'
526 WRITE (unit_nr, '(T2,A4)') 'BSE|'
527 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
528 'where we have introduced'
529 WRITE (unit_nr, '(T2,A4)') 'BSE|'
530 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
531 'BSE|', ψ"_p:", "occupied and virtual molecular orbitals,"
532 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
533 'BSE|', φ"_I(r_e):", "NTO state for the electron,"
534 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
535 'BSE|', χ"_I(r_h):", "NTO state for the hole,"
536 WRITE (unit_nr, '(T2,A4,T7,A,T20,A)') &
537 'BSE|', Λ":", λ"diagonal matrix of NTO weights _I,"
538 WRITE (unit_nr, '(T2,A4)') 'BSE|'
539 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
540 "The NTOs are calculated with the following settings:"
541 WRITE (unit_nr, '(T2,A4)') 'BSE|'
542 WRITE (unit_nr, '(T2,A4,T7,A,T71,I10)') 'BSE|', 'Number of excitations, for which NTOs are computed', &
543 mp2_env%bse%num_print_exc_ntos
544 IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp) THEN
545 WRITE (unit_nr, '(T2,A4,T7,A,T71,F10.3)') 'BSE|', 'Threshold for oscillator strength f^n', &
546 mp2_env%bse%eps_nto_osc_str
547 ELSE
548 WRITE (unit_nr, '(T2,A4,T7,A,T71,A10)') 'BSE|', 'Threshold for oscillator strength f^n', &
549 adjustl("---")
550 END IF
551 WRITE (unit_nr, '(T2,A4,T7,A,T72,F10.3)') 'BSE|', λ'Threshold for NTO weights (_I)^2', &
552 mp2_env%bse%eps_nto_eigval
553 END IF
554
555 ! Write the header of NTO info table
556 IF (unit_nr > 0) THEN
557 WRITE (unit_nr, '(T2,A4)') 'BSE|'
558 IF (.NOT. PRESENT(fm_y)) THEN
559 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
560 'NTOs from solving the BSE within the TDA:'
561 ELSE
562 WRITE (unit_nr, '(T2,A4,T7,A)') 'BSE|', &
563 'NTOs from solving the BSE without the TDA:'
564 END IF
565 WRITE (unit_nr, '(T2,A4)') 'BSE|'
566 WRITE (unit_nr, '(T2,A4,T8,A12,T22,A8,T33,A14,T62,A)') 'BSE|', &
567 'Excitation n', "TDA/ABBA", "Index of NTO I", λ'NTO weights (_I)^2'
568 END IF
569
570 DO j = 1, mp2_env%bse%num_print_exc_ntos
571 n_exc = mp2_env%bse%bse_nto_state_list_final(j)
572 ! Takes care of unallocated oscill_str array in case of Triplet
573 IF (mp2_env%bse%eps_nto_osc_str > 0.0_dp) THEN
574 ! Check actual values
575 IF (oscill_str(n_exc) < mp2_env%bse%eps_nto_osc_str) THEN
576 ! Print skipped levels to table
577 IF (unit_nr > 0) THEN
578 WRITE (unit_nr, '(T2,A4)') 'BSE|'
579 WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T42,A39)') 'BSE|', &
580 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 name="single_part_trans_dm")
588 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 name="nto_coeffs_holes")
592 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 name="nto_coeffs_particles")
596 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 .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 1, homo + 1)
607 CALL cp_fm_release(fm_x_ia)
608 ! Copy Y if present
609 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 .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 homo + 1, 1)
619
620 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 name="LEFT_SINGULAR_MATRIX")
630 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 name="RIGHT_SINGULAR_MATRIX")
634 CALL cp_fm_set_all(fm_eigvr_t, alpha=0.0_dp)
635
636 ALLOCATE (eigval_svd(nao_trunc))
637 eigval_svd(:) = 0.0_dp
638 info_svd = 0
639 CALL cp_fm_svd(fm_m, fm_eigvl, fm_eigvr_t, eigval_svd, info_svd)
640 IF (info_svd /= 0) THEN
641 IF (unit_nr > 0) THEN
642 CALL cp_warn(__location__, &
643 "SVD for computation of NTOs not successful. "// &
644 "Skipping print of NTOs.")
645 IF (info_svd > 0) THEN
646 CALL cp_warn(__location__, &
647 "PDGESVD detected heterogeneity. "// &
648 "Decreasing number of MPI ranks might solve this issue.")
649 END IF
650 END IF
651 ! Release matrices to avoid memory leaks
652 CALL cp_fm_release(fm_m)
653 CALL cp_fm_release(fm_nto_coeff_holes)
654 CALL cp_fm_release(fm_nto_coeff_particles)
655 ELSE
656 ! Rescale singular values as done in Martin2003 (10.1063/1.1558471)
657 ALLOCATE (eigval_svd_squ(nao_trunc))
658 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 IF (.NOT. PRESENT(fm_y)) THEN
661 IF (abs(sum(eigval_svd_squ) - 1) >= 1e-5) THEN
662 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 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 fm_nto_coeff_particles)
674
675 !Release intermediary work matrices
676 CALL cp_fm_release(fm_m)
677 CALL cp_fm_release(fm_eigvl)
678 CALL cp_fm_release(fm_eigvr_t)
679
680 ! Transfer NTO coefficients to sets
681 nto_name(1) = 'Hole_coord'
682 nto_name(2) = 'Particle_coord'
683 ALLOCATE (nto_set(2))
684 ! Extract number of significant NTOs
685 n_nto = 0
686 DO i_nto = 1, nao_trunc
687 IF (eigval_svd_squ(i_nto) > mp2_env%bse%eps_nto_eigval) THEN
688 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 IF (unit_nr > 0) THEN
696 WRITE (unit_nr, '(T2,A4)') 'BSE|'
697 DO i_nto = 1, n_nto
698 WRITE (unit_nr, '(T2,A4,T8,I12,T24,A6,T41,I6,T71,F10.5)') 'BSE|', &
699 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 ncol_global=n_nto)
705 CALL cp_fm_create(fm_nto_set, fm_struct_nto_set)
706 DO i = 1, 2
707 CALL allocate_mo_set(nto_set(i), nao_trunc, n_nto, 0, 0.0_dp, 2.0_dp, 0.0_dp)
708 CALL init_mo_set(nto_set(i), fm_ref=fm_nto_set, name=nto_name(i))
709 END DO
710 CALL cp_fm_release(fm_nto_set)
711 CALL cp_fm_struct_release(fm_struct_nto_set)
712
713 ! Fill NTO sets
714 CALL cp_fm_to_fm(fm_nto_coeff_holes, nto_set(1)%mo_coeff, ncol=n_nto)
715 CALL cp_fm_to_fm(fm_nto_coeff_particles, nto_set(2)%mo_coeff, ncol=n_nto)
716
717 ! Cube files
718 nto_section => section_vals_get_subs_vals(bse_section, "NTO_ANALYSIS")
719 CALL section_vals_val_get(nto_section, "CUBE_FILES", l_val=cube_file)
720 CALL section_vals_val_get(nto_section, "STRIDE", i_vals=stride)
721 CALL section_vals_val_get(nto_section, "APPEND", l_val=append_cube)
722 IF (cube_file) THEN
723 CALL print_bse_nto_cubes(qs_env, nto_set, n_exc, info_approximation, &
724 stride, append_cube, nto_section)
725 END IF
726
727 CALL cp_fm_release(fm_nto_coeff_holes)
728 CALL cp_fm_release(fm_nto_coeff_particles)
729 DEALLOCATE (eigval_svd)
730 DEALLOCATE (eigval_svd_squ)
731 DO i = 1, 2
732 CALL deallocate_mo_set(nto_set(i))
733 END DO
734 DEALLOCATE (nto_set)
735 END IF
736 END DO
737
738 CALL cp_fm_release(fm_mo_coeff)
739 CALL cp_fm_struct_release(fm_struct_m)
740 CALL cp_fm_struct_release(fm_struct_nto_holes)
741 CALL cp_fm_struct_release(fm_struct_nto_particles)
742 CALL cp_fm_struct_release(fm_struct_mo_coeff)
743
744 CALL timestop(handle)
745
746 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 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 CALL timeset(routinen, handle)
789 IF (PRESENT(fm_y_ia)) THEN
790 flag_tda = .false.
791 ELSE
792 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 mask_quadrupole = (/4, 7, 9/)
798
799 CALL cp_fm_struct_create(fm_struct_ia, &
800 context=fm_x_ia%matrix_struct%context, nrow_global=homo, ncol_global=virtual)
801 CALL cp_fm_struct_create(fm_struct_ab, &
802 context=fm_x_ia%matrix_struct%context, nrow_global=virtual, ncol_global=virtual)
803
804 r_e_x(:) = 0.0_dp
805 r_e_y(:) = 0.0_dp
806 r_h_x(:) = 0.0_dp
807 r_h_y(:) = 0.0_dp
808 r_e_sq_x(:) = 0.0_dp
809 r_h_sq_x(:) = 0.0_dp
810 r_e_sq_y(:) = 0.0_dp
811 r_h_sq_y(:) = 0.0_dp
812 r_e_h_xx(:, :) = 0.0_dp
813 r_e_h_xy(:, :) = 0.0_dp
814 r_e_h_yx(:, :) = 0.0_dp
815 r_e_h_yy(:, :) = 0.0_dp
816
817 norm_x = 0.0_dp
818 norm_y = 0.0_dp
819 norm_xpy = 0.0_dp
820
821 ! Initialize values of exciton descriptors
822 exc_descr(i_exc)%r_e(:) = 0.0_dp
823 exc_descr(i_exc)%r_h(:) = 0.0_dp
824 exc_descr(i_exc)%r_e_sq(:) = 0.0_dp
825 exc_descr(i_exc)%r_h_sq(:) = 0.0_dp
826 exc_descr(i_exc)%r_e_h(:, :) = 0.0_dp
827
828 exc_descr(i_exc)%flag_TDA = flag_tda
829 exc_descr(i_exc)%norm_XpY = 0.0_dp
830
831 ! Norm of X
832 CALL cp_fm_trace(fm_x_ia, fm_x_ia, norm_x)
833 norm_xpy = norm_x
834 ! Norm of Y
835 IF (.NOT. flag_tda) THEN
836 CALL cp_fm_trace(fm_y_ia, fm_y_ia, norm_y)
837 norm_xpy = norm_xpy + norm_y
838 END IF
839
840 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 DO i_dir = 1, 3
844 ! <r_h>_X = X_ai µ_ij X_ja + ...
845 CALL trace_exciton_descr(fm_x_ia, fm_multipole_ij_trunc(i_dir), fm_x_ia, r_h_x(i_dir))
846 r_h_x(i_dir) = r_h_x(i_dir)/norm_xpy
847 IF (.NOT. flag_tda) THEN
848 ! <r_h>_X = ... + Y_ia µ_ab Y_bi
849 CALL trace_exciton_descr(fm_y_ia, fm_y_ia, fm_multipole_ab_trunc(i_dir), r_h_y(i_dir))
850 r_h_y(i_dir) = r_h_y(i_dir)/norm_xpy
851 END IF
852 END DO
853 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 DO i_dir = 1, 3
857 ! <r_e>_X = work_ib X_bi + ... = X_ib^T work_ib + ...
858 CALL trace_exciton_descr(fm_x_ia, fm_x_ia, fm_multipole_ab_trunc(i_dir), r_e_x(i_dir))
859 r_e_x(i_dir) = r_e_x(i_dir)/norm_xpy
860 IF (.NOT. flag_tda) THEN
861 ! <r_e>_X = ... + Y_ai µ_ij Y_ja
862 CALL trace_exciton_descr(fm_y_ia, fm_multipole_ij_trunc(i_dir), fm_y_ia, r_e_y(i_dir))
863 r_e_y(i_dir) = r_e_y(i_dir)/norm_xpy
864 END IF
865 END DO
866 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 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 fm_x_ia, r_h_sq_x(i_dir))
873 r_h_sq_x(i_dir) = r_h_sq_x(i_dir)/norm_xpy
874 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 fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_h_sq_y(i_dir))
878 r_h_sq_y(i_dir) = r_h_sq_y(i_dir)/norm_xpy
879 END IF
880 END DO
881 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 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 fm_multipole_ab_trunc(mask_quadrupole(i_dir)), r_e_sq_x(i_dir))
888 r_e_sq_x(i_dir) = r_e_sq_x(i_dir)/norm_xpy
889 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 fm_y_ia, r_e_sq_y(i_dir))
893 r_e_sq_y(i_dir) = r_e_sq_y(i_dir)/norm_xpy
894 END IF
895 END DO
896 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 CALL cp_fm_create(fm_work_ia, fm_struct_ia)
903 CALL cp_fm_create(fm_work_ia_2, fm_struct_ia)
904 CALL cp_fm_create(fm_work_ba, fm_struct_ab)
905 DO i_dir = 1, 3
906 DO j_dir = 1, 3
907 ! First term - X^T µ'_ij X µ_ab
908 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
909 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 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 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 CALL cp_fm_trace(fm_x_ia, fm_work_ia_2, r_e_h_xx(i_dir, j_dir))
918 r_e_h_xx(i_dir, j_dir) = r_e_h_xx(i_dir, j_dir)/norm_xpy
919 IF (.NOT. flag_tda) THEN
920 ! Second term - Y^T µ_ij Y µ'_ab
921 CALL cp_fm_set_all(fm_work_ia, 0.0_dp)
922 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 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 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 CALL cp_fm_trace(fm_y_ia, fm_work_ia_2, r_e_h_yy(i_dir, j_dir))
931 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 CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
943 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 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 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 CALL cp_fm_trace(fm_work_ia, fm_x_ia, r_e_h_xy(i_dir, j_dir))
952 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 CALL cp_fm_set_all(fm_work_ba, 0.0_dp)
964 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 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 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 CALL cp_fm_trace(fm_work_ia, fm_y_ia, r_e_h_yx(i_dir, j_dir))
973 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 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 CALL cp_fm_release(fm_work_ia)
980 CALL cp_fm_release(fm_work_ia_2)
981 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 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 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 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 DO i_dir = 1, 3
997 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 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 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 exc_descr(i_exc)%cov_e_h_sum = 0.0_dp
1005 exc_descr(i_exc)%cov_e_h(:, :) = 0.0_dp
1006 exc_descr(i_exc)%corr_e_h_matrix(:, :) = 0.0_dp
1007 DO i_dir = 1, 3
1008 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 - 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 (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 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 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 - 2*exc_descr(i_exc)%cov_e_h_sum)
1027
1028 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 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 exc_descr(i_exc)%r_e_shift(:) = exc_descr(i_exc)%r_e(:)
1037 exc_descr(i_exc)%r_h_shift(:) = exc_descr(i_exc)%r_h(:)
1038
1039 CALL cp_fm_struct_release(fm_struct_ia)
1040 CALL cp_fm_struct_release(fm_struct_ab)
1041
1042 CALL timestop(handle)
1043
1044 END SUBROUTINE get_exciton_descriptors
1045
1046END MODULE bse_properties
Routines for computing excitonic properties, e.g. exciton diameter, from the BSE.
subroutine, public compute_and_print_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....
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...
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)
...
subroutine, public calculate_ntos(fm_x, fm_y, mo_coeff, homo, virtual, info_approximation, oscill_str, qs_env, unit_nr, mp2_env)
...
Auxiliary routines for GW + Bethe-Salpeter for computing electronic excitations.
Definition bse_util.F:13
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
subroutine, public reshuffle_eigvec(fm_eigvec, fm_eigvec_reshuffled, homo, virtual, n_exc, do_transpose, unit_nr, mp2_env)
...
Definition bse_util.F:1351
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
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
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
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:311
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:122
Basic linear algebra operations for full matrices.
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....
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
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:901
represent the structure of a full matrix
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
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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...
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
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
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,...
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
objects that represent the structure of input sections and the data contained in an input section
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
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Types needed for MP2 calculations.
Definition mp2_types.F:14
basic linear algebra operations for full matrixes
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public c_light_au
Definition physcon.F:90
real(kind=dp), parameter, public evolt
Definition physcon.F:183
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, sab_cneo, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
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...
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
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 ...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...