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