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