Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the propagation of the linearized RT-BSE equations of motion.
10 : !> Propagates the first-order density matrix response Δρ within the active MO window,
11 : !> in the Tamm-Dancoff approximation or with the full (A, B) coupling, instead of the
12 : !> lesser Green's function propagated by rt_bse. Also provides the Liouvillian eigenvalue
13 : !> diagnostic, which builds the Liouvillian by probing the kernel with canonical basis
14 : !> vectors and diagonalizes it.
15 : !> \note The control is handed directly from cp2k_runs
16 : !> The initialization and delta-kick routines are adapted from the full RT-BSE
17 : !> propagator in rt_bse.F.
18 : !> \author Maximilian Graml (03.26)
19 : !> \author Stepan Marek (09.24) - original RT-BSE routines adapted here
20 : ! **************************************************************************************************
21 :
22 : MODULE rt_bse_linearized
23 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm,&
24 : cp_cfm_norm,&
25 : cp_cfm_scale,&
26 : cp_cfm_scale_and_add,&
27 : cp_cfm_transpose
28 : USE cp_cfm_diag, ONLY: cp_cfm_heevd
29 : USE cp_cfm_types, ONLY: &
30 : cp_cfm_get_info, cp_cfm_get_submatrix, cp_cfm_set_all, cp_cfm_set_element, &
31 : cp_cfm_set_submatrix, cp_cfm_to_cfm, cp_cfm_to_fm, cp_cfm_type, cp_fm_to_cfm
32 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
33 : dbcsr_copy,&
34 : dbcsr_get_info,&
35 : dbcsr_p_type,&
36 : dbcsr_release,&
37 : dbcsr_set
38 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
39 : copy_fm_to_dbcsr
40 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
41 : cp_fm_scale_and_add,&
42 : cp_fm_transpose
43 : USE cp_fm_types, ONLY: cp_fm_create,&
44 : cp_fm_get_diag,&
45 : cp_fm_get_info,&
46 : cp_fm_release,&
47 : cp_fm_set_all,&
48 : cp_fm_to_fm_submat_general,&
49 : cp_fm_type
50 : USE cp_log_handling, ONLY: cp_get_default_logger,&
51 : cp_logger_type
52 : USE cp_output_handling, ONLY: cp_add_iter_level,&
53 : cp_iterate,&
54 : cp_print_key_finished_output,&
55 : cp_print_key_unit_nr,&
56 : cp_rm_iter_level
57 : USE dbt_api, ONLY: dbt_copy_matrix_to_tensor
58 : USE force_env_methods, ONLY: force_env_calc_energy_force
59 : USE force_env_types, ONLY: force_env_type
60 : USE input_constants, ONLY: evGW0,&
61 : rtp_bse_ham_gw,&
62 : use_mom_ref_zero,&
63 : use_rt_restart
64 : USE kinds, ONLY: dp
65 : USE machine, ONLY: m_walltime
66 : USE mathconstants, ONLY: twopi
67 : USE moments_utils, ONLY: get_reference_point
68 : USE parallel_gemm_api, ONLY: parallel_gemm
69 : USE physcon, ONLY: evolt,&
70 : seconds
71 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
72 : USE qs_environment_types, ONLY: get_qs_env
73 : USE qs_moments, ONLY: build_local_moment_matrix
74 : USE rpa_gw_kpoints_util, ONLY: cp_cfm_power
75 : USE rt_bse, ONLY: get_hartree,&
76 : get_sigma,&
77 : init_hartree,&
78 : initialize_rtbse_env,&
79 : propagate_density,&
80 : rho_metric
81 : USE rt_bse_io, ONLY: &
82 : check_restart_eps_consistency, output_field, output_moments, output_mos_contravariant, &
83 : output_restart_linearized, print_timestep_info, read_restart_C, read_restart_density, &
84 : read_restart_info, read_restart_trace
85 : USE rt_bse_ri_rs, ONLY: compute_hartree_ri_rs,&
86 : compute_hartree_ri_rs_complex,&
87 : compute_hartree_ri_rs_from_diag,&
88 : rt_bse_ri_rs_ensure_V_grid,&
89 : rt_bse_ri_rs_ensure_W0_grid
90 : USE rt_bse_types, ONLY: create_rtbse_env,&
91 : release_rtbse_env,&
92 : rtbse_env_type
93 : USE rt_propagation_output, ONLY: print_ft
94 : #include "../base/base_uses.f90"
95 :
96 : IMPLICIT NONE
97 :
98 : PRIVATE
99 :
100 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_bse_linearized'
101 :
102 : ! build_shared_sex_and_hartree input-convention selector. The mask choice also fixes the input
103 : ! Hermiticity, which gates the Hartree imaginary channel (Im computed iff non-Hermitian = OV only).
104 : INTEGER, PARAMETER, PRIVATE :: kernel_input_ov = 1, kernel_input_ovvo = 2, kernel_input_full = 3
105 :
106 : PUBLIC :: run_propagation_linearized_bse
107 :
108 : CONTAINS
109 :
110 : ! **************************************************************************************************
111 : !> \brief Runs the electron-only real time propagation of the linearized BSE
112 : !> \param force_env Force environment data, entry point of the calculation
113 : ! **************************************************************************************************
114 58 : SUBROUTINE run_propagation_linearized_bse(force_env)
115 : TYPE(force_env_type), POINTER :: force_env
116 :
117 : CHARACTER(len=*), PARAMETER :: routineN = 'run_propagation_linearized_bse'
118 :
119 : INTEGER :: handle, i, j
120 : REAL(kind=dp) :: t_phys, t_start, timestep_walltime, &
121 : timestep_walltime_start
122 : REAL(kind=dp), DIMENSION(2) :: enum_im, enum_re
123 58 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_lab
124 : TYPE(cp_logger_type), POINTER :: logger
125 : TYPE(rtbse_env_type), POINTER :: rtbse_env
126 :
127 : ! Per-spin (alpha/beta) electron numbers; only 1:n_spin entries are used.
128 :
129 58 : CALL timeset(routineN, handle)
130 :
131 : CALL cp_warn(__LOCATION__, &
132 : "Linearized RT-BSE is under active development. Make sure you understand "// &
133 58 : "the method and validate results before using it for production calculations.")
134 :
135 : ! To Do: Bibliography information
136 :
137 58 : logger => cp_get_default_logger()
138 :
139 : ! Run the initial SCF calculation / read SCF restart information
140 58 : CALL force_env_calc_energy_force(force_env, calc_force=.FALSE., consistent_energies=.FALSE.)
141 :
142 : ! Allocate all persistant storage and read input that does not need further processing
143 58 : CALL create_rtbse_env(rtbse_env, force_env, linearized=.TRUE.)
144 :
145 : ! Restart phase 1a: read sim_start + the original run's dt (from the trace header) BEFORE
146 : ! ENFORCE_MAX_DT, so the continuation inherits that dt instead of a window-dependent one.
147 58 : IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
148 14 : CALL read_restart_info(rtbse_env)
149 : END IF
150 :
151 58 : CALL initialize_maximum_timestep(rtbse_env)
152 :
153 : ! Restart phase 1b: load the trace prefix now that ENFORCE_MAX_DT has sized the trace arrays.
154 58 : IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
155 14 : IF (rtbse_env%sim_start >= rtbse_env%sim_nsteps) THEN
156 0 : CPABORT("RT_RESTART: restart step >= STEPS - increase MOTION%MD%STEPS")
157 : END IF
158 14 : CALL read_restart_trace(rtbse_env)
159 : END IF
160 :
161 58 : CALL print_linrtbse_header_info(rtbse_env)
162 :
163 : ! Build the truncated MO coefficient slabs C_active(:, first_active_mo..last_active_mo)
164 : ! used by all AO<->MO transforms in the linearized path.
165 58 : CALL populate_C_active(rtbse_env)
166 :
167 : ! Initiate iteration level "MD" in order to copy the structure of other RTP codes
168 58 : CALL cp_add_iter_level(logger%iter_info, "MD")
169 : ! Initialize non-trivial values
170 : ! - calculates the moment operators
171 58 : CALL initialize_moments(rtbse_env)
172 : ! - populates overlap and inverse overlap matrices
173 58 : CALL initialize_rtbse_env(rtbse_env)
174 :
175 : ! - populates the fresh SCF density matrix rho^0 (and the rho_orig reference for delta rho)
176 58 : CALL initialize_density_matrix(rtbse_env)
177 :
178 : ! Restart phase 2: overwrite rho from the lab-frame restart files, bridge into this run's MO
179 : ! gauge, then enter this run's rotating frame (rotate_rho_phase is a no-op when omega_shift=0)
180 58 : IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart) THEN
181 14 : CALL read_restart_density(rtbse_env)
182 14 : IF (rtbse_env%restart_extracted) THEN
183 14 : CALL apply_restart_basis_bridge(rtbse_env)
184 14 : t_start = REAL(rtbse_env%sim_start, dp)*rtbse_env%sim_dt
185 30 : DO i = 1, rtbse_env%n_spin
186 30 : CALL rotate_rho_phase(rtbse_env, rtbse_env%rho(i), i, -rtbse_env%omega_shift*t_start)
187 : END DO
188 : END IF
189 : END IF
190 : ! - calculates/populates the G0W0/KS Hamiltonian, respectively
191 58 : CALL initialize_singleparticle_hamiltonian(rtbse_env)
192 : ! Restart Hamiltonian-consistency heads-up: eps_active exists only now, so compare here
193 58 : CALL check_restart_eps_consistency(rtbse_env)
194 : ! Transform initial density matrix to AO basis for use in Hartree and self-energy calculations
195 124 : DO i = 1, rtbse_env%n_spin
196 124 : CALL transform_mo_to_ao_contravariant_cfm(rtbse_env, rtbse_env%rho_orig(i), rtbse_env%rho_ao_scratch(i), i)
197 : END DO
198 : ! - calculates the Hartree reference potential
199 124 : DO i = 1, rtbse_env%n_spin
200 124 : CALL cp_cfm_set_all(rtbse_env%ham_reference(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
201 : END DO
202 58 : CALL initialize_hartree_potential(rtbse_env)
203 : ! - calculates the SEX reference self-energy
204 58 : CALL initialize_sex_selfenergy(rtbse_env)
205 :
206 : ! Liouvillian eigenvalue diagnostic (one-shot at init, TDA or ABBA via dispatcher).
207 : ! Detached from the propagator; safe to call after the reference-init routines.
208 58 : IF (rtbse_env%diagnose_liouvillian_eig) THEN
209 58 : CALL diagnose_liouvillian_eigenvalues(rtbse_env)
210 : END IF
211 :
212 : ! Setup the time based on the starting step
213 : ! Assumes identical dt between two runs
214 58 : rtbse_env%sim_time = REAL(rtbse_env%sim_start, dp)*rtbse_env%sim_dt
215 58 : NULLIFY (rho_lab)
216 : ! Output 0 time moments and field
217 58 : IF (.NOT. rtbse_env%restart_extracted) THEN
218 44 : CALL output_field(rtbse_env)
219 44 : CALL build_rho_lab(rtbse_env, rtbse_env%rho, rtbse_env%sim_time, rho_lab)
220 44 : CALL output_moments(rtbse_env, rho_lab)
221 : END IF
222 :
223 : ! Do not apply the delta kick if we are doing a restart calculation
224 58 : IF (rtbse_env%dft_control%rtp_control%apply_delta_pulse .AND. (.NOT. rtbse_env%restart_extracted)) THEN
225 44 : CALL apply_delta_pulse_MO(rtbse_env)
226 : END IF
227 :
228 : ! ********************** Start the time loop **********************
229 : ! NOTE : Time-loop starts at index sim_start = 0, unless restarted or configured otherwise
230 888 : DO i = rtbse_env%sim_start, rtbse_env%sim_nsteps - 1
231 830 : timestep_walltime_start = m_walltime()
232 :
233 : ! Update the simulation time
234 830 : rtbse_env%sim_time = REAL(i, dp)*rtbse_env%sim_dt
235 830 : rtbse_env%sim_step = i
236 :
237 830 : CALL solve_rk4_timestep(rtbse_env, rtbse_env%rho, rtbse_env%rho_new)
238 : CALL get_electron_number_MO(rtbse_env, rtbse_env%rho_new, &
239 830 : enum_re(1:rtbse_env%n_spin), enum_im(1:rtbse_env%n_spin))
240 830 : timestep_walltime = m_walltime() - timestep_walltime_start
241 830 : CALL print_timestep_info(rtbse_env, i, enum_re(1:rtbse_env%n_spin), step_walltime=timestep_walltime)
242 830 : CALL cp_iterate(logger%iter_info, iter_nr=i, last=(i == rtbse_env%sim_nsteps - 1))
243 :
244 : ! Update rho
245 1780 : DO j = 1, rtbse_env%n_spin
246 1780 : CALL cp_cfm_to_cfm(rtbse_env%rho_new(j), rtbse_env%rho(j))
247 : END DO
248 : ! Print the updated field
249 830 : CALL output_field(rtbse_env)
250 : ! rho is the rotating-frame density at physical time t_phys = (i+1)*dt.
251 : ! Build a lab-frame copy once and feed it to all observable/restart sinks.
252 830 : t_phys = REAL(i + 1, dp)*rtbse_env%sim_dt
253 830 : CALL build_rho_lab(rtbse_env, rtbse_env%rho, t_phys, rho_lab)
254 : ! If needed, print out the density matrix in MO basis
255 830 : CALL output_mos_contravariant(rtbse_env, rho_lab, rtbse_env%rho_section)
256 : ! Also handles outputting to memory
257 830 : CALL output_moments(rtbse_env, rho_lab)
258 : ! Output restart files, so that the restart resumes at the step recorded in .info
259 888 : CALL output_restart_linearized(rtbse_env, rho_lab)
260 : END DO
261 : ! ********************** End the time loop **********************
262 :
263 58 : CALL cp_rm_iter_level(logger%iter_info, "MD")
264 :
265 : ! Carry out the FT
266 : CALL print_ft(rtbse_env%rtp_section, &
267 : rtbse_env%moments_trace, &
268 : rtbse_env%time_trace, &
269 : rtbse_env%field_trace, &
270 : rtbse_env%dft_control%rtp_control, &
271 58 : info_opt=rtbse_env%unit_nr)
272 :
273 : ! Deallocate everything
274 58 : CALL release_rtbse_env(rtbse_env)
275 :
276 58 : CALL timestop(handle)
277 58 : END SUBROUTINE run_propagation_linearized_bse
278 :
279 : ! **************************************************************************************************
280 : !> \brief Computes the analytic RK4 stability bound t* = 2√2 / Ω_max (a.u.) from the largest active
281 : !> GW/KS gap, and (TDA + first-peak) the rotating-frame shift Ω_0 = ε^ai_min.
282 : !> Writes the timestep diagnostics to stdout; optionally rewrites TIMESTEP/STEPS under
283 : !> ENFORCE_MAX_DT; on restart it inherits the original dt from the trace and only rescales STEPS.
284 : !> \param rtbse_env Entry point - rtbse environment
285 : ! **************************************************************************************************
286 58 : SUBROUTINE initialize_maximum_timestep(rtbse_env)
287 : TYPE(rtbse_env_type), POINTER :: rtbse_env
288 :
289 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_maximum_timestep'
290 :
291 : CHARACTER(len=256) :: hint_msg
292 : INTEGER :: handle, i_first, i_last, ispin, &
293 : n_steps_new, n_steps_old
294 : REAL(kind=dp) :: eps_max_ai, eps_min_ai, eps_occ_max, eps_occ_min, eps_virt_max, &
295 : eps_virt_min, ev_tmp, grace_factor, omega_max, sim_dt_as, total_time
296 :
297 58 : CALL timeset(routineN, handle)
298 :
299 58 : i_first = rtbse_env%first_active_mo
300 58 : i_last = rtbse_env%last_active_mo
301 :
302 58 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
303 : omega_max = MAXVAL(rtbse_env%bs_env%eigenval_GW(i_first:i_last, :, :)) - &
304 1940 : MINVAL(rtbse_env%bs_env%eigenval_GW(i_first:i_last, :, :))
305 : ELSE
306 : omega_max = MAXVAL(rtbse_env%bs_env%eigenval_scf_Gamma(i_first:i_last, :)) - &
307 0 : MINVAL(rtbse_env%bs_env%eigenval_scf_Gamma(i_first:i_last, :))
308 : END IF
309 :
310 : ! First-peak shift (TDA only): set Ω_0 = eps_min_ai so the lowest
311 : ! active OV mode oscillates at zero frequency in the rotating frame
312 : ! (RK4-exact for peak 1). omega_max is the full active OV width
313 : ! Delta = eps_max_ai - eps_min_ai, where eps_ai = eps_a - eps_i runs
314 : ! over the active OV pairs only (i in active occupied, a in active
315 : ! virtual).
316 58 : rtbse_env%omega_shift = 0.0_dp
317 58 : IF (rtbse_env%tda_active .AND. rtbse_env%tda_shift_to_first_peak) THEN
318 6 : eps_occ_min = HUGE(0.0_dp)
319 6 : eps_occ_max = -HUGE(0.0_dp)
320 6 : eps_virt_min = HUGE(0.0_dp)
321 6 : eps_virt_max = -HUGE(0.0_dp)
322 12 : DO ispin = 1, rtbse_env%n_spin
323 : ! Active occupied window: first_active_mo .. n_occ(ispin)
324 6 : IF (rtbse_env%n_occ(ispin) >= i_first) THEN
325 6 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
326 18 : ev_tmp = MINVAL(rtbse_env%bs_env%eigenval_GW(i_first:rtbse_env%n_occ(ispin), :, ispin))
327 6 : eps_occ_min = MIN(eps_occ_min, ev_tmp)
328 18 : ev_tmp = MAXVAL(rtbse_env%bs_env%eigenval_GW(i_first:rtbse_env%n_occ(ispin), :, ispin))
329 6 : eps_occ_max = MAX(eps_occ_max, ev_tmp)
330 : ELSE
331 0 : ev_tmp = MINVAL(rtbse_env%bs_env%eigenval_scf_Gamma(i_first:rtbse_env%n_occ(ispin), ispin))
332 0 : eps_occ_min = MIN(eps_occ_min, ev_tmp)
333 0 : ev_tmp = MAXVAL(rtbse_env%bs_env%eigenval_scf_Gamma(i_first:rtbse_env%n_occ(ispin), ispin))
334 0 : eps_occ_max = MAX(eps_occ_max, ev_tmp)
335 : END IF
336 : END IF
337 : ! Active virtual window: n_occ(ispin)+1 .. last_active_mo
338 12 : IF (rtbse_env%n_occ(ispin) < i_last) THEN
339 6 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
340 78 : ev_tmp = MINVAL(rtbse_env%bs_env%eigenval_GW(rtbse_env%n_occ(ispin) + 1:i_last, :, ispin))
341 6 : eps_virt_min = MIN(eps_virt_min, ev_tmp)
342 78 : ev_tmp = MAXVAL(rtbse_env%bs_env%eigenval_GW(rtbse_env%n_occ(ispin) + 1:i_last, :, ispin))
343 6 : eps_virt_max = MAX(eps_virt_max, ev_tmp)
344 : ELSE
345 0 : ev_tmp = MINVAL(rtbse_env%bs_env%eigenval_scf_Gamma(rtbse_env%n_occ(ispin) + 1:i_last, ispin))
346 0 : eps_virt_min = MIN(eps_virt_min, ev_tmp)
347 0 : ev_tmp = MAXVAL(rtbse_env%bs_env%eigenval_scf_Gamma(rtbse_env%n_occ(ispin) + 1:i_last, ispin))
348 0 : eps_virt_max = MAX(eps_virt_max, ev_tmp)
349 : END IF
350 : END IF
351 : END DO
352 :
353 6 : IF (eps_occ_max > -HUGE(0.0_dp) .AND. eps_virt_min < HUGE(0.0_dp)) THEN
354 6 : eps_min_ai = eps_virt_min - eps_occ_max
355 6 : eps_max_ai = eps_virt_max - eps_occ_min
356 6 : rtbse_env%omega_shift = eps_min_ai
357 6 : omega_max = eps_max_ai - eps_min_ai
358 6 : IF (rtbse_env%unit_nr > 0) THEN
359 : WRITE (rtbse_env%unit_nr, '(A)') &
360 3 : " RTBSE| ---------- First-peak shift diagnostics (TDA, active OV pairs) ----------"
361 : WRITE (rtbse_env%unit_nr, '(A,F14.6,A,F14.6)') &
362 3 : " RTBSE| eps_occ [eV] min / max =", eps_occ_min*evolt, &
363 6 : " /", eps_occ_max*evolt
364 : WRITE (rtbse_env%unit_nr, '(A,F14.6,A,F14.6)') &
365 3 : " RTBSE| eps_virt [eV] min / max =", eps_virt_min*evolt, &
366 6 : " /", eps_virt_max*evolt
367 : WRITE (rtbse_env%unit_nr, '(A,F14.6,A,F14.6)') &
368 3 : " RTBSE| eps_ai [eV] min / max =", eps_min_ai*evolt, &
369 6 : " /", eps_max_ai*evolt
370 : WRITE (rtbse_env%unit_nr, '(A,F14.6)') &
371 3 : " RTBSE| omega_shift [eV] =", rtbse_env%omega_shift*evolt
372 : WRITE (rtbse_env%unit_nr, '(A,F14.6)') &
373 3 : " RTBSE| omega_max [eV] (full) =", omega_max*evolt
374 : WRITE (rtbse_env%unit_nr, '(A)') &
375 3 : " RTBSE| ------------------------------------------------------------------------"
376 : END IF
377 : ELSE
378 : ! Active window has no genuine OV pair - fall back to no shift
379 0 : rtbse_env%omega_shift = 0.0_dp
380 0 : rtbse_env%tda_shift_to_first_peak = .FALSE.
381 : END IF
382 : END IF
383 :
384 58 : rtbse_env%omega_max = omega_max
385 :
386 58 : IF (omega_max > 0.0_dp) THEN
387 : ! t* = 2√2 / Ω_max (a.u.; imaginary-axis RK4 bound |R(iy)| ≤ 1 at y = 2√2)
388 58 : rtbse_env%maximum_timestep = 2.0_dp*SQRT(2.0_dp)/omega_max
389 : ELSE
390 : CALL cp_abort(__LOCATION__, &
391 : "Error in estimating maximum timestep: largest KS/GW gap is "// &
392 : "non-positive. Check the active MO window (cutoffs) and the "// &
393 0 : "eigenvalues.")
394 : END IF
395 :
396 58 : IF (rtbse_env%sim_dt <= 0.0_dp) THEN
397 : CALL cp_abort(__LOCATION__, &
398 : "TIMESTEP must be positive for linearized RT-BSE. Use RTBSE%ENFORCE_MAX_DT "// &
399 0 : "with a positive TIMESTEP to automatically rewrite TIMESTEP and STEPS.")
400 : END IF
401 :
402 58 : n_steps_old = MAX(0, rtbse_env%sim_nsteps)
403 :
404 58 : IF (rtbse_env%enforce_max_dt) THEN
405 8 : total_time = REAL(n_steps_old, dp)*rtbse_env%sim_dt
406 : ! Full code needs a grace factor of 4
407 8 : IF (rtbse_env%tda_active) THEN
408 : grace_factor = 1.0_dp
409 : ELSE
410 0 : grace_factor = 4.0_dp
411 : END IF
412 8 : IF (rtbse_env%dft_control%rtp_control%initial_wfn == use_rt_restart .AND. &
413 : rtbse_env%sim_dt_restart > 0.0_dp) THEN
414 : ! Continuation: dt is frozen in the trace, so inherit it and only rescale the step count
415 : ! to the requested window. Recomputing dt from the (longer) window would desync the trace
416 : ! time-grid and trip the continuation guard in read_restart_trace.
417 2 : rtbse_env%sim_dt = rtbse_env%sim_dt_restart
418 2 : n_steps_new = MAX(1, NINT(total_time/rtbse_env%sim_dt))
419 2 : rtbse_env%sim_nsteps = n_steps_new
420 2 : sim_dt_as = rtbse_env%sim_dt*seconds*1e18_dp
421 : WRITE (hint_msg, '(A,F16.4,A,I0,A)') &
422 2 : 'ENFORCE_MAX_DT on restart: inheriting original TIMESTEP ', sim_dt_as, &
423 4 : ' as and setting STEPS to ', n_steps_new, '.'
424 2 : CALL cp_hint(__LOCATION__, TRIM(hint_msg))
425 : ! The inherited dt was stable in the original run; warn only if this run's stability
426 : ! window shrank below it (e.g. the recomputed GW eigenvalues shifted the gap).
427 2 : IF (rtbse_env%sim_dt > rtbse_env%maximum_timestep/grace_factor) THEN
428 : CALL cp_warn(__LOCATION__, &
429 : "ENFORCE_MAX_DT restart: inherited dt exceeds this run's stability "// &
430 0 : "limit - the recomputed Hamiltonian may make the propagation unstable.")
431 : END IF
432 : ELSE
433 6 : n_steps_new = MAX(1, CEILING(total_time/(rtbse_env%maximum_timestep/grace_factor)))
434 6 : rtbse_env%sim_dt = total_time/REAL(n_steps_new, dp)
435 6 : rtbse_env%sim_nsteps = n_steps_new
436 6 : sim_dt_as = rtbse_env%sim_dt*seconds*1e18_dp
437 : WRITE (hint_msg, '(A,F16.4,A,I0,A)') &
438 6 : 'ENFORCE_MAX_DT enabled. Resetting TIMESTEP to ', sim_dt_as, &
439 12 : ' as and STEPS to ', n_steps_new, '.'
440 6 : CALL cp_hint(__LOCATION__, TRIM(hint_msg))
441 : END IF
442 : END IF
443 :
444 58 : IF (rtbse_env%sim_nsteps /= n_steps_old) THEN
445 8 : CALL reallocate_ft_traces(rtbse_env)
446 : END IF
447 :
448 58 : CALL timestop(handle)
449 58 : END SUBROUTINE initialize_maximum_timestep
450 :
451 : ! **************************************************************************************************
452 : !> \brief Reallocate FT trace buffers after ENFORCE_MAX_DT rewrites the step count.
453 : !> \param rtbse_env Entry point - rtbse environment
454 : ! **************************************************************************************************
455 8 : SUBROUTINE reallocate_ft_traces(rtbse_env)
456 : TYPE(rtbse_env_type), POINTER :: rtbse_env
457 :
458 : CHARACTER(len=*), PARAMETER :: routineN = 'reallocate_ft_traces'
459 :
460 : INTEGER :: handle
461 :
462 8 : CALL timeset(routineN, handle)
463 :
464 8 : IF (ASSOCIATED(rtbse_env%moments_trace)) DEALLOCATE (rtbse_env%moments_trace)
465 8 : IF (ASSOCIATED(rtbse_env%field_trace)) DEALLOCATE (rtbse_env%field_trace)
466 8 : IF (ASSOCIATED(rtbse_env%time_trace)) DEALLOCATE (rtbse_env%time_trace)
467 :
468 : ALLOCATE (rtbse_env%moments_trace(rtbse_env%n_spin, 3, rtbse_env%sim_nsteps + 1), &
469 536 : source=CMPLX(0.0_dp, 0.0_dp, kind=dp))
470 : ALLOCATE (rtbse_env%field_trace(3, rtbse_env%sim_nsteps + 1), &
471 312 : source=CMPLX(0.0_dp, 0.0_dp, kind=dp))
472 96 : ALLOCATE (rtbse_env%time_trace(rtbse_env%sim_nsteps + 1), source=0.0_dp)
473 :
474 8 : CALL timestop(handle)
475 8 : END SUBROUTINE reallocate_ft_traces
476 :
477 : ! **************************************************************************************************
478 : !> \brief Prints the linRTBSE run header to stdout: active-MO window (first/last/count) and the
479 : !> occupied/virtual energy cutoffs.
480 : !> \param rtbse_env Entry point - rtbse environment
481 : ! **************************************************************************************************
482 58 : SUBROUTINE print_linrtbse_header_info(rtbse_env)
483 : TYPE(rtbse_env_type) :: rtbse_env
484 :
485 : INTEGER :: ispin, n_steps
486 : REAL(kind=dp) :: e_first, e_last, fft_resolution, &
487 : nyquist_frequency, total_time
488 : TYPE(cp_logger_type), POINTER :: logger
489 :
490 58 : logger => cp_get_default_logger()
491 58 : n_steps = MAX(0, rtbse_env%sim_nsteps)
492 58 : total_time = REAL(n_steps, dp)*rtbse_env%sim_dt
493 58 : fft_resolution = 0.0_dp
494 58 : nyquist_frequency = 0.0_dp
495 58 : IF (total_time > 0.0_dp) fft_resolution = twopi/total_time
496 58 : IF (rtbse_env%sim_dt > 0.0_dp) nyquist_frequency = twopi/(2.0_dp*rtbse_env%sim_dt)
497 :
498 58 : IF (rtbse_env%unit_nr > 0) THEN
499 29 : WRITE (rtbse_env%unit_nr, *) ''
500 : WRITE (rtbse_env%unit_nr, '(A)') ' /-----------------------------------------------'// &
501 29 : '------------------------------\'
502 : WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
503 29 : ' |'
504 : WRITE (rtbse_env%unit_nr, '(A)') ' | Linearized Real Time Bethe-Salpeter Propagation'// &
505 29 : ' |'
506 : WRITE (rtbse_env%unit_nr, '(A)') ' | '// &
507 29 : ' |'
508 : WRITE (rtbse_env%unit_nr, '(A)') ' \-----------------------------------------------'// &
509 29 : '------------------------------/'
510 29 : WRITE (rtbse_env%unit_nr, *) ''
511 :
512 29 : WRITE (rtbse_env%unit_nr, '(A18,L62)') ' Apply delta pulse', &
513 58 : rtbse_env%dft_control%rtp_control%apply_delta_pulse
514 29 : WRITE (rtbse_env%unit_nr, '(A)') ''
515 29 : WRITE (rtbse_env%unit_nr, '(A18,L62)') ' Use Tamm-Dancoff approximation', &
516 58 : rtbse_env%tda_active
517 29 : IF (rtbse_env%tda_active) THEN
518 20 : WRITE (rtbse_env%unit_nr, '(A,T71,L10)') ' TDA first-peak shift active', &
519 40 : rtbse_env%tda_shift_to_first_peak
520 20 : IF (rtbse_env%tda_shift_to_first_peak) THEN
521 3 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.6)') ' TDA first-peak shift Omega_0 [eV]:', &
522 6 : rtbse_env%omega_shift*evolt
523 : END IF
524 : END IF
525 :
526 29 : WRITE (rtbse_env%unit_nr, '(A)') ''
527 :
528 29 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.4)') ' Estimated maximum timestep within stability region [as]:', &
529 58 : rtbse_env%maximum_timestep*seconds*1e18_dp
530 29 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.4)') ' Applied timestep [as]:', &
531 58 : rtbse_env%sim_dt*seconds*1e18_dp
532 29 : WRITE (rtbse_env%unit_nr, '(A,T71,I10)') ' Number of propagation steps:', n_steps
533 29 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.4)') ' Total propagation time [as]:', &
534 58 : total_time*seconds*1e18_dp
535 29 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.6)') ' Estimated FFT frequency resolution without interpolation [eV]:', &
536 58 : fft_resolution*evolt
537 29 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.6)') ' Nyquist frequency [eV]:', &
538 58 : nyquist_frequency*evolt
539 29 : WRITE (rtbse_env%unit_nr, '(A,T65,F16.6)') ' Estimated maximum oscillation frequency (gap-based) [eV]:', &
540 58 : rtbse_env%omega_max*evolt
541 :
542 : ! Which single-particle eigenvalues the propagation uses. bs_env%eigenval_GW carries the
543 : ! result of the highest GW flavour requested.
544 29 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
545 29 : IF (rtbse_env%bs_env%gw_flavour == evGW0) THEN
546 : WRITE (rtbse_env%unit_nr, '(A,T75,A6)') &
547 0 : ' GW flavor for computing GW eigenvalues used in RT-BSE:', ' evGW0'
548 : ELSE
549 : WRITE (rtbse_env%unit_nr, '(A,T75,A6)') &
550 29 : ' GW flavor for computing GW eigenvalues used in RT-BSE:', ' G0W0'
551 : END IF
552 : ELSE
553 : WRITE (rtbse_env%unit_nr, '(A,T75,A6)') &
554 0 : ' Single-particle eigenvalues used in RT-BSE:', ' KS'
555 : END IF
556 :
557 : ! Active MO window (energy-cutoff truncation) for linearized RT-BSE
558 29 : IF (rtbse_env%rtbse_energy_cutoff_occ > 0.0_dp) THEN
559 1 : WRITE (rtbse_env%unit_nr, '(A,T71,F10.3)') ' Active-window occupied energy cutoff [eV]:', &
560 2 : rtbse_env%rtbse_energy_cutoff_occ*evolt
561 : ELSE
562 28 : WRITE (rtbse_env%unit_nr, '(A,T71,A10)') ' Active-window occupied energy cutoff [eV]:', ' disabled'
563 : END IF
564 29 : IF (rtbse_env%rtbse_energy_cutoff_empty > 0.0_dp) THEN
565 1 : WRITE (rtbse_env%unit_nr, '(A,T71,F10.3)') ' Active-window virtual energy cutoff [eV]:', &
566 2 : rtbse_env%rtbse_energy_cutoff_empty*evolt
567 : ELSE
568 28 : WRITE (rtbse_env%unit_nr, '(A,T71,A10)') ' Active-window virtual energy cutoff [eV]:', ' disabled'
569 : END IF
570 29 : WRITE (rtbse_env%unit_nr, '(A,T71,I10)') ' First active occupied MO index:', rtbse_env%first_active_mo
571 29 : WRITE (rtbse_env%unit_nr, '(A,T71,I10)') ' Last active virtual MO index:', rtbse_env%last_active_mo
572 29 : WRITE (rtbse_env%unit_nr, '(A,T71,I10)') ' Number of active MOs:', rtbse_env%mo_active
573 29 : IF (rtbse_env%active_mo_truncation) THEN
574 : ! The window is cut on the DFT axis but propagated on the QP axis, so the QP edges may
575 : ! exceed the nominal cutoff. Print both so the window can be checked against the input.
576 2 : DO ispin = 1, rtbse_env%n_spin
577 : e_first = (rtbse_env%bs_env%eigenval_scf_Gamma(rtbse_env%first_active_mo, ispin) - &
578 1 : rtbse_env%bs_env%eigenval_scf_Gamma(rtbse_env%n_occ(ispin), ispin))*evolt
579 : e_last = (rtbse_env%bs_env%eigenval_scf_Gamma(rtbse_env%last_active_mo, ispin) - &
580 1 : rtbse_env%bs_env%eigenval_scf_Gamma(rtbse_env%n_occ(ispin) + 1, ispin))*evolt
581 1 : WRITE (rtbse_env%unit_nr, '(A,I1,A,T71,F10.3)') ' Spin ', ispin, &
582 2 : ' first active MO, E - E_HOMO (KS) [eV]:', e_first
583 1 : WRITE (rtbse_env%unit_nr, '(A,I1,A,T71,F10.3)') ' Spin ', ispin, &
584 2 : ' last active MO, E - E_LUMO (KS) [eV]:', e_last
585 2 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
586 : e_first = (rtbse_env%bs_env%eigenval_GW(rtbse_env%first_active_mo, 1, ispin) - &
587 1 : rtbse_env%bs_env%eigenval_GW(rtbse_env%n_occ(ispin), 1, ispin))*evolt
588 : e_last = (rtbse_env%bs_env%eigenval_GW(rtbse_env%last_active_mo, 1, ispin) - &
589 1 : rtbse_env%bs_env%eigenval_GW(rtbse_env%n_occ(ispin) + 1, 1, ispin))*evolt
590 1 : WRITE (rtbse_env%unit_nr, '(A,I1,A,T71,F10.3)') ' Spin ', ispin, &
591 2 : ' first active MO, E - E_HOMO (QP) [eV]:', e_first
592 1 : WRITE (rtbse_env%unit_nr, '(A,I1,A,T71,F10.3)') ' Spin ', ispin, &
593 2 : ' last active MO, E - E_LUMO (QP) [eV]:', e_last
594 : END IF
595 : END DO
596 : END IF
597 : END IF
598 :
599 58 : END SUBROUTINE print_linrtbse_header_info
600 :
601 : ! **************************************************************************************************
602 : !> \brief Populates rtbse_env%C_active(i_spin) (n_ao x mo_active) by extracting columns
603 : !> first_active_mo..last_active_mo from bs_env%fm_mo_coeff_Gamma(i_spin).
604 : !> \param rtbse_env RT-BSE environment
605 : !> \author Maximilian Graml (05.26)
606 : ! **************************************************************************************************
607 58 : SUBROUTINE populate_C_active(rtbse_env)
608 : TYPE(rtbse_env_type), POINTER :: rtbse_env
609 :
610 : CHARACTER(len=*), PARAMETER :: routineN = 'populate_C_active'
611 :
612 : INTEGER :: handle, i
613 :
614 58 : CALL timeset(routineN, handle)
615 :
616 124 : DO i = 1, rtbse_env%n_spin
617 66 : CALL cp_fm_set_all(rtbse_env%C_active(i), 0.0_dp)
618 : CALL cp_fm_to_fm_submat_general( &
619 : rtbse_env%bs_env%fm_mo_coeff_Gamma(i), rtbse_env%C_active(i), &
620 : rtbse_env%n_ao, rtbse_env%mo_active, &
621 : 1, rtbse_env%first_active_mo, &
622 : 1, 1, &
623 124 : rtbse_env%bs_env%fm_mo_coeff_Gamma(i)%matrix_struct%context)
624 : END DO
625 :
626 58 : CALL timestop(handle)
627 58 : END SUBROUTINE populate_C_active
628 :
629 : ! **************************************************************************************************
630 : !> \brief Builds the dipole moment operators r_mn in the active-MO basis (per axis, per spin) from
631 : !> the AO moment matrices: moments(k,σ) at the reference point, moments_field(k,σ) at origin.
632 : !> \param rtbse_env RT-BSE environment
633 : !> \author Stepan Marek (09.24)
634 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
635 : ! **************************************************************************************************
636 58 : SUBROUTINE initialize_moments(rtbse_env)
637 : TYPE(rtbse_env_type), POINTER :: rtbse_env
638 :
639 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_moments'
640 :
641 : INTEGER :: handle, i_spin, k
642 : REAL(kind=dp), DIMENSION(3) :: rpoint
643 : TYPE(cp_fm_type) :: tmp_ao
644 58 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments_dbcsr_p
645 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
646 :
647 58 : CALL timeset(routineN, handle)
648 : ! Get pointers to parameters from qs_env
649 58 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env, matrix_s=matrix_s)
650 :
651 : ! AO-sized scratch buffer for moment matrices before transform to MO basis
652 58 : CALL cp_fm_create(tmp_ao, bs_env%fm_s_Gamma%matrix_struct)
653 :
654 : ! ****** START MOMENTS OPERATOR CALCULATION
655 : ! Construct moments from dbcsr
656 : NULLIFY (moments_dbcsr_p)
657 232 : ALLOCATE (moments_dbcsr_p(3))
658 232 : DO k = 1, 3
659 : ! Make sure the pointer is empty
660 174 : NULLIFY (moments_dbcsr_p(k)%matrix)
661 : ! Allocate a new matrix that the pointer points to
662 174 : ALLOCATE (moments_dbcsr_p(k)%matrix)
663 : ! Create the matrix storage - matrix copies the structure of overlap matrix
664 232 : CALL dbcsr_copy(moments_dbcsr_p(k)%matrix, matrix_s(1)%matrix)
665 : END DO
666 : ! Run the moment calculation
667 : ! check for presence to prevent memory errors
668 58 : rpoint(:) = 0.0_dp
669 : CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, &
670 58 : reference=rtbse_env%moment_ref_type, ref_point=rtbse_env%user_moment_ref_point)
671 58 : CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint)
672 : ! Copy to AO scratch then transform to MO-active
673 232 : DO k = 1, 3
674 174 : CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, tmp_ao)
675 430 : DO i_spin = 1, rtbse_env%n_spin
676 372 : CALL transform_ao_to_mo_covariant_fm(rtbse_env, tmp_ao, rtbse_env%moments(k, i_spin), i_spin)
677 : END DO
678 : END DO
679 : ! TODO: remove moments_field (only needed for the TDDFT comparison)
680 : ! Now, repeat without reference point to get the moments for field
681 : CALL get_reference_point(rpoint, qs_env=rtbse_env%qs_env, &
682 58 : reference=use_mom_ref_zero)
683 58 : CALL build_local_moment_matrix(rtbse_env%qs_env, moments_dbcsr_p, 1, rpoint)
684 232 : DO k = 1, 3
685 174 : CALL copy_dbcsr_to_fm(moments_dbcsr_p(k)%matrix, tmp_ao)
686 430 : DO i_spin = 1, rtbse_env%n_spin
687 372 : CALL transform_ao_to_mo_covariant_fm(rtbse_env, tmp_ao, rtbse_env%moments_field(k, i_spin), i_spin)
688 : END DO
689 : END DO
690 :
691 : ! Now can deallocate dbcsr matrices
692 232 : DO k = 1, 3
693 174 : CALL dbcsr_release(moments_dbcsr_p(k)%matrix)
694 232 : DEALLOCATE (moments_dbcsr_p(k)%matrix)
695 : END DO
696 58 : DEALLOCATE (moments_dbcsr_p)
697 58 : CALL cp_fm_release(tmp_ao)
698 : ! ****** END MOMENTS OPERATOR CALCULATION
699 :
700 58 : CALL timestop(handle)
701 116 : END SUBROUTINE initialize_moments
702 :
703 : ! **************************************************************************************************
704 : !> \brief Initial MO density ρ^0_mn = f_m δ_mn (f_m = 1 on active occupied, 0 on virtual), copied to
705 : !> rho_orig as the reference for the δ-kick. Imaginary part zero.
706 : !> \param rtbse_env RT-BSE environment
707 : !> \author Stepan Marek (09.24)
708 : !> \author Maximilian Graml (03.26) - adapted to the linearized active-MO path
709 : ! **************************************************************************************************
710 116 : SUBROUTINE initialize_density_matrix(rtbse_env)
711 : TYPE(rtbse_env_type), POINTER :: rtbse_env
712 :
713 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_density_matrix'
714 :
715 : INTEGER :: handle, i, i_row_global, ii, &
716 : j_col_global, jj, ncol_global, &
717 : ncol_local, nrow_global, nrow_local
718 58 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
719 :
720 58 : CALL timeset(routineN, handle)
721 :
722 : ! Get distribution of MO-active workspace
723 : CALL cp_fm_get_info(rtbse_env%real_workspace_mo(1), &
724 : nrow_global=nrow_global, ncol_global=ncol_global, &
725 : nrow_local=nrow_local, ncol_local=ncol_local, &
726 58 : row_indices=row_indices, col_indices=col_indices)
727 :
728 : ! Iterate over both spins
729 124 : DO i = 1, rtbse_env%n_spin
730 : !Ensure that workspace is set to 0
731 66 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(1), 0.0_dp)
732 456 : DO ii = 1, nrow_local
733 390 : i_row_global = row_indices(ii)
734 5100 : DO jj = 1, ncol_local
735 4644 : j_col_global = col_indices(jj)
736 4644 : IF (i_row_global == j_col_global .AND. &
737 390 : (i_row_global + rtbse_env%first_active_mo - 1) <= rtbse_env%n_occ(i)) THEN
738 33 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = 1.0_dp
739 : END IF
740 : END DO
741 : END DO
742 : ! Sets imaginary part to zero
743 66 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), mtarget=rtbse_env%rho(i))
744 : ! Save the reference value for the case of delta kick
745 124 : CALL cp_cfm_to_cfm(rtbse_env%rho(i), rtbse_env%rho_orig(i))
746 : END DO
747 : ! rho_orig stays the SCF reference for delta rho; a restart overwrites only rho, downstream in
748 : ! the driver (read_restart_density + apply_restart_basis_bridge + rotate_rho_phase).
749 :
750 58 : CALL timestop(handle)
751 58 : END SUBROUTINE initialize_density_matrix
752 :
753 : ! **************************************************************************************************
754 : !> \brief Single-particle reference Hamiltonian in the active-MO basis: H^0_mn = ε^GW_m δ_mn (or KS
755 : !> ε^scf), diagonal. TDA first-peak adds +Ω_0/2 on occupied, -Ω_0/2 on virtual diagonals.
756 : !> \param rtbse_env RT-BSE environment
757 : !> \author Stepan Marek (09.24)
758 : !> \author Maximilian Graml (03.26) - refactor in prep. of linearized propagation
759 : ! **************************************************************************************************
760 116 : SUBROUTINE initialize_singleparticle_hamiltonian(rtbse_env)
761 : TYPE(rtbse_env_type), POINTER :: rtbse_env
762 :
763 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_singleparticle_hamiltonian'
764 :
765 : INTEGER :: abs_mo_idx, handle, i, i_row_global, ii, &
766 : j_col_global, jj, ncol_global, &
767 : ncol_local, nrow_global, nrow_local
768 58 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
769 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
770 :
771 58 : CALL timeset(routineN, handle)
772 : ! Get pointers to parameters from qs_env
773 58 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
774 :
775 : ! Get distribution of MO-active workspace
776 : CALL cp_fm_get_info(rtbse_env%real_workspace_mo(1), &
777 : nrow_global=nrow_global, ncol_global=ncol_global, &
778 : nrow_local=nrow_local, ncol_local=ncol_local, &
779 58 : row_indices=row_indices, col_indices=col_indices)
780 :
781 : !Ensure that workspace is set to 0
782 58 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(1), 0.0_dp)
783 904 : rtbse_env%eps_active(:, :) = 0.0_dp
784 : ! ****** START SINGLE PARTICLE HAMILTONIAN CALCULATION
785 124 : DO i = 1, rtbse_env%n_spin
786 456 : DO ii = 1, nrow_local
787 390 : i_row_global = row_indices(ii)
788 5100 : DO jj = 1, ncol_local
789 4644 : j_col_global = col_indices(jj)
790 5034 : IF (i_row_global == j_col_global) THEN
791 390 : abs_mo_idx = i_row_global + rtbse_env%first_active_mo - 1
792 390 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
793 : ! GW Hamiltonian
794 390 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = bs_env%eigenval_GW(abs_mo_idx, 1, i)
795 : ELSE
796 : ! KS Hamiltonian
797 0 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = bs_env%eigenval_scf_Gamma(abs_mo_idx, i)
798 : END IF
799 : ! First-peak shift (TDA only): rotate the single-particle Hamiltonian
800 : ! into a frame where the lowest active OV mode oscillates at zero, with
801 : ! Ω_0 = eps_min_ai. Adds +Ω_0/2 on active occupied diagonals and
802 : ! -Ω_0/2 on active virtual diagonals so that [h', rho]_OV picks up an
803 : ! overall (- eps_ai + Ω_0) and OO/VV blocks remain commutator-free.
804 390 : IF (rtbse_env%tda_active .AND. rtbse_env%tda_shift_to_first_peak) THEN
805 36 : IF (abs_mo_idx <= rtbse_env%n_occ(i)) THEN
806 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = &
807 3 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) + 0.5_dp*rtbse_env%omega_shift
808 : ELSE
809 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = &
810 33 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) - 0.5_dp*rtbse_env%omega_shift
811 : END IF
812 : END IF
813 : ! Mirror the finalized diagonal into the replicated active-energy array
814 390 : rtbse_env%eps_active(i_row_global, i) = rtbse_env%real_workspace_mo(1)%local_data(ii, jj)
815 : END IF
816 : END DO
817 : END DO
818 124 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), mtarget=rtbse_env%ham_reference_singleparticle(i))
819 : END DO
820 : ! Each diagonal element was set on its single owner rank; sum to replicate eps_active.
821 1750 : CALL rtbse_env%real_workspace_mo(1)%matrix_struct%para_env%sum(rtbse_env%eps_active)
822 : ! ****** END SINGLE PARTICLE HAMILTONIAN CALCULATION
823 :
824 58 : CALL timestop(handle)
825 58 : END SUBROUTINE initialize_singleparticle_hamiltonian
826 :
827 : ! **************************************************************************************************
828 : !> \brief Reference Hartree subtraction: builds V^H[ρ^0] (AO-RI or RI-RS) and subtracts it into
829 : !> ham_reference, realizing H_eff = ... + V_H[ρ] - V_H[ρ^0]. Only for non-TDA n_spin=1 (in
830 : !> TDA the OV/VO projection of the OO-diagonal ρ^0 vanishes, so the reference is zero).
831 : !> \param rtbse_env RT-BSE environment
832 : !> \author Stepan Marek (09.24)
833 : !> \author Maximilian Graml (03.26) - add transform to MO
834 : ! **************************************************************************************************
835 58 : SUBROUTINE initialize_hartree_potential(rtbse_env)
836 : TYPE(rtbse_env_type), POINTER :: rtbse_env
837 :
838 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_hartree_potential'
839 :
840 : INTEGER :: handle, i, n_grid
841 : LOGICAL :: use_hartree_reference, use_rirs_kernel
842 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
843 :
844 58 : CALL timeset(routineN, handle)
845 : ! Get pointers to parameters from qs_env
846 58 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
847 : use_hartree_reference = (.NOT. rtbse_env%tda_active) .AND. (rtbse_env%n_spin == 1) .AND. &
848 58 : (.NOT. rtbse_env%debug_disable_hartree)
849 58 : use_rirs_kernel = rtbse_env%rirs_kernel
850 :
851 : ! Make sure the RI-RS V_grid kernel is available if we need it for Hartree.
852 58 : IF (use_rirs_kernel) CALL rt_bse_ri_rs_ensure_V_grid(bs_env, rtbse_env%qs_env)
853 :
854 : ! Spin-summed grid-density accumulators for the RI-RS Hartree reuse (diag(φρφ^T) harvested in SEX;
855 : ! mat_phi_mu_l is grid x AO, so its row count is n_grid). Allocated ONLY for RI-RS + Hartree, so the
856 : ! get_sigma harvest calls (passed unconditionally below) see an absent optional and self-disable on
857 : ! AO-RI / Hartree-off (F2008 unallocated-allocatable -> absent optional). Sole allocator, runs once.
858 22 : IF (use_rirs_kernel .AND. (.NOT. rtbse_env%debug_disable_hartree)) THEN
859 22 : CALL dbcsr_get_info(bs_env%ri_rs%mat_phi_mu_l, nfullrows_total=n_grid)
860 88 : ALLOCATE (rtbse_env%hartree_diag_re(n_grid), rtbse_env%hartree_diag_im(n_grid))
861 : END IF
862 :
863 : ! The RI-RS Hartree normally reuses the grid density harvested by the SEX kernel. With SEX
864 : ! disabled but Hartree on, that diagonal is never produced, so the Hartree must rebuild the full
865 : ! real-space density grid itself every RK4 stage - much slower. Warn once (this is a debug-only
866 : ! configuration); the rebuild fallback lives in the use_sex branches of the Hartree kernels.
867 58 : IF (use_rirs_kernel .AND. (.NOT. rtbse_env%debug_disable_hartree) .AND. rtbse_env%debug_disable_sex) THEN
868 : CALL cp_warn(__LOCATION__, &
869 : "RI-RS Hartree rebuilds the full density grid every RK4 stage because SEX is "// &
870 : "disabled (DEBUG_DISABLE_SEX) and no SEX-harvested diagonal is available to reuse. "// &
871 0 : "This slows down the Hartree computation considerably.")
872 : END IF
873 :
874 : ! ****** START HARTREE POTENTIAL REFERENCE CALCULATION
875 : ! v_dbcsr is needed by either AO-RI Hartree (here) or by AO-RI SEX (W = V + W^c assembly).
876 : IF (.NOT. use_rirs_kernel) THEN
877 36 : CALL init_hartree(rtbse_env, rtbse_env%v_dbcsr)
878 : END IF
879 : ! Always zero ham_reference here (this routine is the first to touch it).
880 : ! The Hartree reference subtraction is then conditionally added below.
881 124 : DO i = 1, rtbse_env%n_spin
882 124 : CALL cp_cfm_set_all(rtbse_env%ham_reference(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
883 : END DO
884 : ! Calculate the original Hartree potential
885 : ! Uses rho_orig - same as rho for initial run but different for continued run
886 : ! In TDA the propagator evaluates separate OV/VO-projected kernel passes.
887 : ! rho_orig is OO-diagonal in the MO basis, so its OV/VO projections vanish
888 : ! and the corresponding Hartree reference is identically zero.
889 58 : IF (use_hartree_reference) THEN
890 28 : DO i = 1, rtbse_env%n_spin
891 14 : IF (use_rirs_kernel) THEN
892 : ! V^H_λσ = sum_l φ_λ(r_l) v_l φ_σ(r_l), v_l = sum_l' V_ll' n_l' (RI-RS)
893 : ! AO-RI get_hartree uses only Re(rho); mirror that on the RI-RS path.
894 : CALL cp_cfm_to_fm(msource=rtbse_env%rho_ao_scratch(i), &
895 6 : mtargetr=rtbse_env%real_workspace(1))
896 : CALL compute_hartree_ri_rs(bs_env, rtbse_env%real_workspace(1), &
897 6 : rtbse_env%hartree_curr_ao(i))
898 : ELSE
899 : ! V^H_λσ = sum_PQ (λσ|P) V_PQ [sum_µν (µν|Q) ρ^0_µν] (AO-RI; reference density)
900 8 : CALL get_hartree(rtbse_env, rtbse_env%rho_ao_scratch(i), rtbse_env%hartree_curr_ao(i))
901 : END IF
902 : ! Scaling by spin degeneracy
903 14 : CALL cp_fm_scale(rtbse_env%spin_degeneracy, rtbse_env%hartree_curr_ao(i))
904 : ! Transform to MO basis (AO scratch -> MO-active result)
905 14 : CALL transform_ao_to_mo_covariant_fm(rtbse_env, rtbse_env%hartree_curr_ao(i), rtbse_env%hartree_curr(i), i)
906 : ! Apply occupation factor f_n-f_m
907 14 : CALL transform_mo_occupation_factor_diff_fm(rtbse_env, rtbse_env%hartree_curr(i), i)
908 : ! Subtract the reference from the reference Hamiltonian
909 : ! following H_eff = ... + V_Hartree(rho) - V_Hartree(rho_0),
910 14 : CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr(i), mtarget=rtbse_env%ham_workspace(1))
911 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
912 28 : CMPLX(-1.0, 0.0, kind=dp), rtbse_env%ham_workspace(1))
913 : END DO
914 : END IF
915 : ! ****** END HARTREE POTENTIAL REFERENCE CALCULATION
916 :
917 58 : CALL timestop(handle)
918 58 : END SUBROUTINE initialize_hartree_potential
919 :
920 : ! **************************************************************************************************
921 : !> \brief Reference SEX self-energy subtraction into ham_reference (non-TDA, n_spin=1). Assembles
922 : !> W = V + W^c, then Σ^SX = -W ρ^0, (f_n - f_m)-weighted, subtracted.
923 : !> \param rtbse_env RT-BSE environment
924 : !> \author Stepan Marek (09.24)
925 : !> \author Maximilian Graml (03.26) - add transform to MO
926 : ! **************************************************************************************************
927 58 : SUBROUTINE initialize_sex_selfenergy(rtbse_env)
928 : TYPE(rtbse_env_type), POINTER :: rtbse_env
929 :
930 : CHARACTER(len=*), PARAMETER :: routineN = 'initialize_sex_selfenergy'
931 :
932 : INTEGER :: handle, i
933 : LOGICAL :: use_rirs_kernel, use_sex_reference
934 :
935 58 : CALL timeset(routineN, handle)
936 : use_sex_reference = (.NOT. rtbse_env%tda_active) .AND. (rtbse_env%n_spin == 1) .AND. &
937 58 : (.NOT. rtbse_env%debug_disable_sex)
938 58 : use_rirs_kernel = rtbse_env%rirs_kernel
939 :
940 : ! Make sure the RI-RS W0_grid kernel is available if we need it for SEX.
941 58 : IF (use_rirs_kernel) CALL rt_bse_ri_rs_ensure_W0_grid(rtbse_env%bs_env, rtbse_env%qs_env)
942 :
943 : ! ****** START SEX REFERENCE CALCULATION
944 : ! w_dbcsr and screened_dbt are needed for get_sigma routines (AO-RI path only).
945 : ! For RI-RS the W = V + W^c kernel is on the real-space grid in mat_W0_grid_rtbse.
946 : IF (.NOT. use_rirs_kernel) THEN
947 36 : IF (rtbse_env%ham_reference_type == rtp_bse_ham_gw) THEN
948 : ! W(w=0) is built by the GW step only under its RTBSE rtp_method gate; reaching this
949 : ! consumer without it means gate and consumer disagree. The HF branch below needs no W.
950 36 : IF (.NOT. ASSOCIATED(rtbse_env%bs_env%fm_W_MIC_freq_zero%matrix_struct)) THEN
951 : CALL cp_abort(__LOCATION__, &
952 : "RT-BSE AO-RI kernel needs the screened interaction W(w=0), which the "// &
953 : "GW step did not build. Select the RT-BSE propagator with '&RTBSE' or "// &
954 0 : "'&RTBSE RTBSE', not '&RTBSE TDDFT'.")
955 : END IF
956 : ! In a non-HF calculation, copy the actual correlation part of the interaction
957 36 : CALL copy_fm_to_dbcsr(rtbse_env%bs_env%fm_W_MIC_freq_zero, rtbse_env%w_dbcsr)
958 : ELSE
959 : ! In HF, correlation is set to zero
960 0 : CALL dbcsr_set(rtbse_env%w_dbcsr, 0.0_dp)
961 : END IF
962 : ! Add the Hartree to the screened_dbt tensor - now W = V + W^c
963 36 : CALL dbcsr_add(rtbse_env%w_dbcsr, rtbse_env%v_dbcsr, 1.0_dp, 1.0_dp)
964 36 : CALL dbt_copy_matrix_to_tensor(rtbse_env%w_dbcsr, rtbse_env%screened_dbt)
965 : END IF
966 : ! Calculate the SEX starting energies
967 124 : DO i = 1, rtbse_env%n_spin
968 : ! Calculate the exchange (SEX) part for this spin channel
969 : ! Uses rho_orig - same as rho for initial run but different for continued run
970 : ! For KS reference this is the time-dependent Fock exchange (w_dbcsr = v only).
971 : ! In TDA the propagator evaluates separate OV/VO-projected kernel passes.
972 : ! rho_orig is OO-diagonal in the MO basis, so its OV/VO projections
973 : ! vanish and the SEX reference must remain zero in TDA.
974 124 : IF (use_sex_reference) THEN
975 : ! Σ^SX_λσ = -sum_νQ [sum_µ (λµ|Q) ρ^0_µν][sum_P (νσ|P) W_PQ] (reference)
976 14 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX_ao(i), -1.0_dp, rtbse_env%rho_ao_scratch(i))
977 : ! Transform to MO basis (AO scratch -> MO-active result)
978 14 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%sigma_SEX_ao(i), rtbse_env%sigma_SEX(i), i)
979 : ! Apply occupation factor f_n-f_m
980 14 : CALL transform_mo_occupation_factor_diff_cfm(rtbse_env, rtbse_env%sigma_SEX(i), i)
981 : ! Subtract from the complex reference Hamiltonian
982 : CALL cp_cfm_scale_and_add(CMPLX(1.0, 0.0, kind=dp), rtbse_env%ham_reference(i), &
983 14 : CMPLX(-1.0, 0.0, kind=dp), rtbse_env%sigma_SEX(i))
984 : END IF
985 : END DO
986 : ! ****** END SEX REFERENCE CALCULATION
987 :
988 58 : CALL timestop(handle)
989 58 : END SUBROUTINE initialize_sex_selfenergy
990 :
991 : ! **************************************************************************************************
992 : !> \brief Propagates the density one timestep by RK4 for ∂_t ρ = f(t,ρ) =
993 : !> -i( Δε Δρ + (f_n - f_m) V_Hartree(Δρ) + ΔΣ(Δρ) ); see body for the 4-stage scheme.
994 : !> Spin loop is inner to each stage (cross-spin Hartree coupling).
995 : !> \param rtbse_env Entry point - rtbse environment
996 : !> \param rho_start Initial density matrix
997 : !> \param rho_end Final density matrix
998 : ! **************************************************************************************************
999 830 : SUBROUTINE solve_rk4_timestep(rtbse_env, rho_start, rho_end)
1000 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1001 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_start, rho_end
1002 :
1003 : CHARACTER(len=*), PARAMETER :: routineN = 'solve_rk4_timestep'
1004 :
1005 : INTEGER :: handle, i
1006 :
1007 830 : CALL timeset(routineN, handle)
1008 : ! RK4 follows the typical scheme
1009 : ! d/dt ρ = f(t, ρ)
1010 : ! f(t,ρ) = -i( Δε Δρ + (f_n-f_m) V_Hartree(Δρ) + ΔΣ(Δρ) )
1011 : ! i.e. RK4 reads
1012 : ! k_1 = f(t, ρ_start)
1013 : ! k_2 = f(t + dt/2, ρ_start + dt/2 * k_1)
1014 : ! k_3 = f(t + dt/2, ρ_start + dt/2 * k_2)
1015 : ! k_4 = f(t + dt, ρ_start + dt * k_3)
1016 : ! ρ_end = ρ_start + dt/6 * (k_1 + 2*k_2 + 2*k_3 + k_4)
1017 : ! Note that the effective Hamiltonian needs to be updated for each evaluation of f,
1018 : ! as it depends on the density matrix at the respective time
1019 :
1020 : ! Spin loop is INNER to each RK4 stage (inside do_rk4_stage): cross-spin-coupled kernels (Hartree in
1021 : ! open shell) need every spin's stage density before any spin advances. rk4_coefficients(i) holds
1022 : ! spin i's CURRENT-stage k (indexed by spin, not stage - see its allocation in create_rtbse_env),
1023 : ! reused across stages; rho_workspace(i) holds spin i's stage density. Bit-identical for n_spin=1.
1024 1780 : DO i = 1, rtbse_env%n_spin
1025 1780 : CALL cp_cfm_to_cfm(rho_start(i), rho_end(i))
1026 : END DO
1027 : ! Each stage: do_rk4_stage evaluates k = f(stage density) for all spins, accumulates
1028 : ! rho_end += result_weight*dt*k (Butcher b = 1/6, 1/3, 1/3, 1/6) and forms the next stage density
1029 : ! rho_workspace = rho_start + advance_weight*dt*k (node c = 1/2, 1/2, 1; omitted on k_4, which only
1030 : ! accumulates). The dt factor is applied inside do_rk4_stage, so the calls show the bare weights.
1031 : ! k_1 = f(t, rho_start)
1032 : CALL do_rk4_stage(rtbse_env, rho_start, rho_start, rho_end, &
1033 830 : result_weight=1.0_dp/6.0_dp, advance_weight=0.5_dp)
1034 : ! k_2 = f(t + dt/2, rho_start + dt/2 * k_1)
1035 : CALL do_rk4_stage(rtbse_env, rtbse_env%rho_workspace, rho_start, rho_end, &
1036 830 : result_weight=1.0_dp/3.0_dp, advance_weight=0.5_dp)
1037 : ! k_3 = f(t + dt/2, rho_start + dt/2 * k_2)
1038 : CALL do_rk4_stage(rtbse_env, rtbse_env%rho_workspace, rho_start, rho_end, &
1039 830 : result_weight=1.0_dp/3.0_dp, advance_weight=1.0_dp)
1040 : ! k_4 = f(t + dt, rho_start + dt * k_3)
1041 : CALL do_rk4_stage(rtbse_env, rtbse_env%rho_workspace, rho_start, rho_end, &
1042 830 : result_weight=1.0_dp/6.0_dp)
1043 :
1044 : ! Update bookkeeping to the next timestep similar to logic of etrs_scf_loop
1045 830 : rtbse_env%sim_step = rtbse_env%sim_step + 1
1046 830 : rtbse_env%sim_time = rtbse_env%sim_time + rtbse_env%sim_dt
1047 :
1048 830 : CALL timestop(handle)
1049 830 : END SUBROUTINE solve_rk4_timestep
1050 :
1051 : ! **************************************************************************************************
1052 : !> \brief Takes one RK4 stage. Evaluates k = f(t, rho_eval) for every spin into
1053 : !> rtbse_env%rk4_coefficients, accumulates it into the running result (rho_end += result_weight*dt*k)
1054 : !> and - unless this is the last stage - forms the next stage density
1055 : !> (rho_workspace = rho_base + advance_weight*dt*k). Builds the stage's shared kernels first: the
1056 : !> cross-spin Hartree is built once per stage and consumed by every spin in update_effective_ham_MO.
1057 : !> All shell/kernel combinations go through build_shared_sex_and_hartree (mask_mode selects the
1058 : !> input convention); update_effective_ham_MO is then a pure consumer.
1059 : !> No timeset/timestop: the callees are individually timed and this runs 4x per RK4 timestep.
1060 : !> \param rtbse_env RT-BSE environment
1061 : !> \param rho_eval Per-spin MO density f is evaluated at (rho_start for k_1, rho_workspace otherwise)
1062 : !> \param rho_base Per-spin MO density the next stage advances from (the step's rho_start)
1063 : !> \param rho_end Per-spin RK4 result accumulator (= rho_start + dt/6*(k1+2k2+2k3+k4) after all 4 stages)
1064 : !> \param result_weight RK4 Butcher weight b (dt factor applied inside) for accumulating k into rho_end
1065 : !> \param advance_weight RK4 node c (dt factor applied inside) for the next stage density; ABSENT on the last stage
1066 : !> \author Maximilian Graml
1067 : ! **************************************************************************************************
1068 3320 : SUBROUTINE do_rk4_stage(rtbse_env, rho_eval, rho_base, rho_end, result_weight, advance_weight)
1069 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1070 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_eval, rho_base, rho_end
1071 : REAL(kind=dp), INTENT(IN) :: result_weight
1072 : REAL(kind=dp), INTENT(IN), OPTIONAL :: advance_weight
1073 :
1074 : INTEGER :: i, mask_mode
1075 :
1076 : ! Input convention for this stage's shared kernel build:
1077 3320 : IF (rtbse_env%tda_active) THEN
1078 2040 : mask_mode = kernel_input_ov ! TDA (any shell): OV only, non-Hermitian
1079 1280 : ELSE IF (rtbse_env%n_spin > 1) THEN
1080 160 : mask_mode = kernel_input_ovvo ! open-shell ABBA: OV+VO, Hermitian
1081 : ELSE
1082 1120 : mask_mode = kernel_input_full ! closed-shell ABBA: full ρ, Hermitian
1083 : END IF
1084 : ! Build this stage's shared SEX + bare Hartree; consumers are per-spin in update_effective_ham_MO.
1085 3320 : CALL build_shared_sex_and_hartree(rtbse_env, rho_eval, mask_mode)
1086 7120 : DO i = 1, rtbse_env%n_spin
1087 3800 : CALL update_effective_ham_MO(rtbse_env, rho_eval(i), rtbse_env%rk4_coefficients(i), i)
1088 7120 : IF (rtbse_env%tda_active .OR. rtbse_env%n_spin > 1) THEN
1089 2680 : CALL project_drho_to_ov(rtbse_env, rtbse_env%rk4_coefficients(i), i)
1090 : END IF
1091 : END DO
1092 : ! Fold each spin's k into the RK4 result and (unless last stage) form the next stage density.
1093 : ! The dt factor lives here so the call sites carry the bare RK4 weights.
1094 7120 : DO i = 1, rtbse_env%n_spin
1095 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rho_end(i), &
1096 3800 : CMPLX(result_weight*rtbse_env%sim_dt, 0.0_dp, kind=dp), rtbse_env%rk4_coefficients(i))
1097 7120 : IF (PRESENT(advance_weight)) THEN
1098 2850 : CALL cp_cfm_to_cfm(rho_base(i), rtbse_env%rho_workspace(i))
1099 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%rho_workspace(i), &
1100 2850 : CMPLX(advance_weight*rtbse_env%sim_dt, 0.0_dp, kind=dp), rtbse_env%rk4_coefficients(i))
1101 : END IF
1102 : END DO
1103 3320 : END SUBROUTINE do_rk4_stage
1104 :
1105 : ! **************************************************************************************************
1106 : !> \brief Builds the spin-summed complex Hartree potential in AO, once per RK4 stage (AO-RI or RI-RS).
1107 : !> rho_total = spin_degeneracy * sum_sigma (OV-masked rho^sigma -> AO); hartree_total_ao =
1108 : !> V_H[rho_total]. update_effective_ham_MO back-transforms it with each spin's C, so this
1109 : !> single call feeds every spin block - the cross-spin Hartree coupling of open shell.
1110 : !> Side effect: leaves rho_ao_scratch(sigma) = OV-masked AO density (recomputed per spin in
1111 : !> update_effective_ham_MO; built here only to form the sum).
1112 : !> \param rtbse_env RT-BSE environment
1113 : !> \param rho_stage Per-spin MO density at the current RK4 stage
1114 : !> \param keep_ovvo .FALSE. = OV source mask (TDA); .TRUE. = OV+VO (open-shell ABBA).
1115 : !> \author Maximilian Graml
1116 : ! **************************************************************************************************
1117 176 : SUBROUTINE build_shared_hartree_ao(rtbse_env, rho_stage, keep_ovvo)
1118 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1119 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_stage
1120 : LOGICAL, INTENT(IN) :: keep_ovvo
1121 :
1122 : CHARACTER(len=*), PARAMETER :: routineN = 'build_shared_hartree_ao'
1123 :
1124 : INTEGER :: handle, isp
1125 :
1126 176 : CALL timeset(routineN, handle)
1127 : ! Per spin: mask the stage density and project to AO.
1128 528 : DO isp = 1, rtbse_env%n_spin
1129 352 : CALL cp_cfm_to_cfm(rho_stage(isp), rtbse_env%rho_delta_mo(isp))
1130 : CALL mask_mo_block_cfm(rtbse_env, rtbse_env%rho_delta_mo(isp), isp, &
1131 352 : keep_OV=.TRUE., keep_ovvo=keep_ovvo)
1132 : CALL transform_mo_to_ao_contravariant_cfm(rtbse_env, rtbse_env%rho_delta_mo(isp), &
1133 528 : rtbse_env%rho_ao_scratch(isp), isp)
1134 : END DO
1135 : ! Sum: rho_total = spin_degeneracy * sum_isp rho_ao_scratch(isp).
1136 176 : CALL cp_cfm_set_all(rtbse_env%rho_total_ao_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
1137 528 : DO isp = 1, rtbse_env%n_spin
1138 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%rho_total_ao_scratch, &
1139 : CMPLX(rtbse_env%spin_degeneracy, 0.0_dp, kind=dp), &
1140 528 : rtbse_env%rho_ao_scratch(isp))
1141 : END DO
1142 : ! One complex Hartree contraction on the summed density. RI backend = kernel choice: both
1143 : ! get_hartree_complex and compute_hartree_ri_rs_complex take AO x AO in/out and are spin-blind,
1144 : ! so the spin-summed cross-spin density feeds whichever kernel is active.
1145 176 : IF (rtbse_env%rirs_kernel) THEN
1146 : ! V^H_λσ = sum_l φ_λ(r_l) v_l φ_σ(r_l), v_l = sum_l' V_ll' n_l'[ρ^total] (RI-RS)
1147 : CALL compute_hartree_ri_rs_complex(rtbse_env%bs_env, rtbse_env%rho_total_ao_scratch, &
1148 88 : rtbse_env%hartree_total_ao)
1149 : ELSE
1150 : ! V^H_λσ = sum_PQ (λσ|P) V_PQ [sum_µν (µν|Q) ρ^total_µν] (AO-RI)
1151 : CALL get_hartree_complex(rtbse_env, rtbse_env%rho_total_ao_scratch, &
1152 88 : rtbse_env%hartree_total_ao, 1)
1153 : END IF
1154 176 : CALL timestop(handle)
1155 176 : END SUBROUTINE build_shared_hartree_ao
1156 :
1157 : ! **************************************************************************************************
1158 : !> \brief Unified cross-spin kernel builder, once per RK4 stage for every shell. Computes the
1159 : !> per-spin SEX self-energy (stashed in sigma_SEX_ao(σ)) and the single shared bare Hartree
1160 : !> (hartree_total_ao), so update_effective_ham_MO only consumes them.
1161 : !> mask_mode selects the source-density convention (OV / OV+VO / full-ρ); it also determines
1162 : !> the input Hermiticity, which gates the Hartree imaginary channel: Im computed only when
1163 : !> mask_mode = kernel_input_ov (non-Hermitian TDA input); for Hermitian input Im ≡ 0 analytically
1164 : !> and is skipped. Hartree emitted BARE (no spin_degeneracy); consumer scales g on the MO output.
1165 : !> \param rtbse_env RT-BSE environment
1166 : !> \param rho_stage Per-spin MO density at the current RK4 stage
1167 : !> \param mask_mode Input-convention selector: kernel_input_ov / kernel_input_ovvo / kernel_input_full
1168 : !> \author Maximilian Graml
1169 : ! **************************************************************************************************
1170 3320 : SUBROUTINE build_shared_sex_and_hartree(rtbse_env, rho_stage, mask_mode)
1171 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1172 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_stage
1173 : INTEGER, INTENT(IN) :: mask_mode
1174 :
1175 : CHARACTER(len=*), PARAMETER :: routineN = 'build_shared_sex_and_hartree'
1176 :
1177 : INTEGER :: handle, isp
1178 : LOGICAL :: harvest_im, use_hartree, &
1179 : use_rirs_kernel, use_sex
1180 :
1181 3320 : CALL timeset(routineN, handle)
1182 3320 : use_hartree = .NOT. rtbse_env%debug_disable_hartree
1183 3320 : use_sex = .NOT. rtbse_env%debug_disable_sex
1184 3320 : use_rirs_kernel = rtbse_env%rirs_kernel
1185 : ! Im channel iff the masked input is non-Hermitian (OV-only, TDA). For Hermitian input (OV+VO
1186 : ! or full ρ) Im(ρ) is antisymmetric and Coulomb factors are symmetric, so V_H[Im] ≡ 0.
1187 3320 : harvest_im = (mask_mode == kernel_input_ov)
1188 :
1189 : ! Pre-zero the RI-RS spin-summed grid-diagonal accumulators (the SEX harvest target).
1190 3320 : IF (use_rirs_kernel .AND. use_sex .AND. use_hartree) THEN
1191 456640 : rtbse_env%hartree_diag_re(:) = 0.0_dp
1192 272640 : IF (harvest_im) rtbse_env%hartree_diag_im(:) = 0.0_dp
1193 : END IF
1194 : ! Per spin: stage ρ into AO (masked per mask_mode); SEX (stash sigma_SEX_ao(σ)); RI-RS harvests
1195 : ! the diagonal.
1196 7120 : DO isp = 1, rtbse_env%n_spin
1197 3800 : CALL cp_cfm_to_cfm(rho_stage(isp), rtbse_env%rho_delta_mo(isp))
1198 2360 : SELECT CASE (mask_mode)
1199 : CASE (kernel_input_ov)
1200 : CALL mask_mo_block_cfm(rtbse_env, rtbse_env%rho_delta_mo(isp), isp, &
1201 2360 : keep_OV=.TRUE., keep_ovvo=.FALSE.)
1202 : CASE (kernel_input_ovvo)
1203 : CALL mask_mo_block_cfm(rtbse_env, rtbse_env%rho_delta_mo(isp), isp, &
1204 320 : keep_OV=.TRUE., keep_ovvo=.TRUE.)
1205 : CASE (kernel_input_full)
1206 : ! no mask: full ρ (closed-shell ABBA; reference subtracted later via ham_reference)
1207 : CASE DEFAULT
1208 3800 : CPABORT("Unknown mask_mode in build_shared_sex_and_hartree")
1209 : END SELECT
1210 : CALL transform_mo_to_ao_contravariant_cfm(rtbse_env, rtbse_env%rho_delta_mo(isp), &
1211 3800 : rtbse_env%rho_ao_scratch(isp), isp)
1212 7120 : IF (use_sex) THEN
1213 : ! Σ^SX_λσ = -sum_νQ [sum_µ (λµ|Q) Δρ_µν][sum_P (νσ|P) W_PQ]. Im accumulator passed only
1214 : ! when harvesting; absent (unallocated optional) on AO-RI / Hartree-off / Hermitian.
1215 3480 : IF (harvest_im) THEN
1216 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX_ao(isp), -1.0_dp, rtbse_env%rho_ao_scratch(isp), &
1217 : grid_diag_re_accum=rtbse_env%hartree_diag_re, &
1218 2200 : grid_diag_im_accum=rtbse_env%hartree_diag_im)
1219 : ELSE
1220 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX_ao(isp), -1.0_dp, rtbse_env%rho_ao_scratch(isp), &
1221 1280 : grid_diag_re_accum=rtbse_env%hartree_diag_re)
1222 : END IF
1223 : END IF
1224 : END DO
1225 : ! Shared cross-spin Hartree, emitted BARE (no spin_degeneracy — consumer scales g on output).
1226 : ! RI-RS reuses the spin-summed grid diagonal; AO-RI contracts rho_total = sum_σ ρ_σ.
1227 : ! Hermitian input: real Hartree path (Im ≡ 0); non-Hermitian: complex.
1228 3320 : IF (use_hartree) THEN
1229 3000 : IF (use_rirs_kernel .AND. use_sex) THEN
1230 1600 : IF (harvest_im) THEN
1231 : CALL compute_hartree_ri_rs_from_diag(rtbse_env%bs_env, rtbse_env%hartree_diag_re, &
1232 960 : rtbse_env%hartree_total_ao, n_im=rtbse_env%hartree_diag_im)
1233 : ELSE
1234 : CALL compute_hartree_ri_rs_from_diag(rtbse_env%bs_env, rtbse_env%hartree_diag_re, &
1235 640 : rtbse_env%hartree_total_ao)
1236 : END IF
1237 : ELSE
1238 1400 : CALL cp_cfm_set_all(rtbse_env%rho_total_ao_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
1239 2960 : DO isp = 1, rtbse_env%n_spin
1240 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%rho_total_ao_scratch, &
1241 2960 : CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%rho_ao_scratch(isp))
1242 : END DO
1243 1400 : IF (use_rirs_kernel) THEN
1244 0 : IF (harvest_im) THEN
1245 : CALL compute_hartree_ri_rs_complex(rtbse_env%bs_env, rtbse_env%rho_total_ao_scratch, &
1246 0 : rtbse_env%hartree_total_ao)
1247 : ELSE
1248 : ! Hermitian input: real RI-RS on Re(ρ_total) -> cfm with Im ≡ 0.
1249 : CALL cp_cfm_to_fm(msource=rtbse_env%rho_total_ao_scratch, &
1250 0 : mtargetr=rtbse_env%real_workspace(1))
1251 : CALL compute_hartree_ri_rs(rtbse_env%bs_env, rtbse_env%real_workspace(1), &
1252 0 : rtbse_env%hartree_curr_ao(1))
1253 0 : CALL cp_cfm_set_all(rtbse_env%hartree_total_ao, CMPLX(0.0_dp, 0.0_dp, kind=dp))
1254 : CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr_ao(1), &
1255 0 : mtarget=rtbse_env%hartree_total_ao)
1256 : END IF
1257 : ELSE
1258 1400 : IF (harvest_im) THEN
1259 : CALL get_hartree_complex(rtbse_env, rtbse_env%rho_total_ao_scratch, &
1260 920 : rtbse_env%hartree_total_ao, 1)
1261 : ELSE
1262 : ! Hermitian input: real AO-RI on Re(ρ_total) -> cfm with Im ≡ 0.
1263 : CALL get_hartree(rtbse_env, rtbse_env%rho_total_ao_scratch, &
1264 480 : rtbse_env%hartree_curr_ao(1))
1265 480 : CALL cp_cfm_set_all(rtbse_env%hartree_total_ao, CMPLX(0.0_dp, 0.0_dp, kind=dp))
1266 : CALL cp_fm_to_cfm(msourcer=rtbse_env%hartree_curr_ao(1), &
1267 480 : mtarget=rtbse_env%hartree_total_ao)
1268 : END IF
1269 : END IF
1270 : END IF
1271 : END IF
1272 3320 : CALL timestop(handle)
1273 3320 : END SUBROUTINE build_shared_sex_and_hartree
1274 :
1275 : ! **************************************************************************************************
1276 : !> \brief Assembles the linearized RT-BSE right-hand side in the MO basis and returns it scaled by -i:
1277 : !> ham_effective <- -i( [H^0, ρ] + (f_n - f_m)(ΔΣ^SX[Δρ] + ΔV^H[Δρ]) ), i.e. the
1278 : !> f(t,ρ) of ∂_t Δρ_mn = -i( (ε_m - ε_n)Δρ_mn + (f_n - f_m)(V^H_mn + Σ^SX_mn) ).
1279 : !> ham_reference already carries KS+G0W0 minus the reference SEX/Hartree, so the kernels
1280 : !> enter as differences vs the reference. Forks: tda_active (drop B-coupling: one OV kernel
1281 : !> pass + VO conjugate) vs full ABBA (OV+VO); n_spin and the KERNEL_RI (rirs_kernel) flag
1282 : !> select the AO-RI or RI-RS backend per term.
1283 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
1284 : !> \param rho Real and imaginary parts ( + spin) of the density at current time
1285 : !> \param ham_effective Effective Hamiltonian in the MO basis that is updated in this routine
1286 : !> \param ispin Spin channel σ being assembled
1287 : ! **************************************************************************************************
1288 7600 : SUBROUTINE update_effective_ham_MO(rtbse_env, rho, ham_effective, ispin)
1289 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1290 : TYPE(cp_cfm_type) :: rho, ham_effective
1291 : INTEGER :: ispin
1292 :
1293 : CHARACTER(len=*), PARAMETER :: routineN = 'update_effective_ham_MO'
1294 :
1295 : INTEGER :: handle, i_global, i_loc, j_global, &
1296 : j_loc, ncl, nrl
1297 3800 : INTEGER, DIMENSION(:), POINTER :: c_idx, r_idx
1298 : LOGICAL :: use_hartree, use_sex
1299 :
1300 3800 : CALL timeset(routineN, handle)
1301 3800 : use_hartree = .NOT. rtbse_env%debug_disable_hartree
1302 3800 : use_sex = .NOT. rtbse_env%debug_disable_sex
1303 :
1304 : ! Reset the effective Hamiltonian to KS Hamiltonian + G0W0 - reference SEX - reference Hartree
1305 : ! Sets the imaginary part to zero
1306 3800 : CALL cp_cfm_to_cfm(rtbse_env%ham_reference(ispin), ham_effective)
1307 : ! [H^0, ρ]_mn = (ε_m - ε_n) ρ_mn exactly (H^0 diagonal in active-MO basis). Element-wise
1308 : ! local pass replacing the two gemms; ham_effective and rho share fm_struct_mo_active so
1309 : ! their local layouts coincide. Covers all (m,n) (OO/VV retained for closed-shell ABBA).
1310 : CALL cp_cfm_get_info(matrix=ham_effective, nrow_local=nrl, ncol_local=ncl, &
1311 3800 : row_indices=r_idx, col_indices=c_idx)
1312 48440 : DO j_loc = 1, ncl
1313 44640 : j_global = c_idx(j_loc)
1314 313400 : DO i_loc = 1, nrl
1315 264960 : i_global = r_idx(i_loc)
1316 : ham_effective%local_data(i_loc, j_loc) = ham_effective%local_data(i_loc, j_loc) &
1317 : + CMPLX(rtbse_env%eps_active(i_global, ispin) - rtbse_env%eps_active(j_global, ispin), &
1318 309600 : 0.0_dp, kind=dp)*rho%local_data(i_loc, j_loc)
1319 : END DO
1320 : END DO
1321 : ! Determine the field at current time
1322 3800 : IF (rtbse_env%dft_control%apply_efield_field) THEN
1323 : CALL cp_abort(__LOCATION__, &
1324 : "Continuous/pulsed E(t) field coupling is not implemented for linearized "// &
1325 : "RT-BSE. Only the delta-kick (impulsive) absorption spectrum is supported; "// &
1326 0 : "use APPLY_DELTA_PULSE.")
1327 : ELSE
1328 : ! No field
1329 15200 : rtbse_env%field(:) = 0.0_dp
1330 : END IF
1331 3800 : IF (.NOT. rtbse_env%tda_active) THEN
1332 : ! ===== ABBA: consume the prebuilt per-spin SEX (sigma_SEX_ao(σ)) and shared bare Hartree
1333 : ! (hartree_total_ao). (f_n-f_m) zeros OO/VV and sets OV/VO signs. Closed shell uses full-ρ
1334 : ! input (no mask, reference subtracted via ham_reference); open shell uses OV+VO mask. =====
1335 1440 : IF (use_sex) THEN
1336 : ! Σ^SX_λσ = -sum_νQ [sum_µ (λµ|Q) Δρ_µν][sum_P (νσ|P) W_PQ]
1337 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%sigma_SEX_ao(ispin), &
1338 1280 : rtbse_env%sigma_SEX(ispin), ispin)
1339 1280 : CALL transform_mo_occupation_factor_diff_cfm(rtbse_env, rtbse_env%sigma_SEX(ispin), ispin)
1340 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), ham_effective, &
1341 1280 : CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%sigma_SEX(ispin))
1342 : END IF
1343 1440 : IF (use_hartree) THEN
1344 : ! Builder emits bare V_H (no spin_degeneracy); fold g here (post-occ-factor).
1345 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%hartree_total_ao, &
1346 1280 : rtbse_env%ham_workspace(1), ispin)
1347 1280 : CALL transform_mo_occupation_factor_diff_cfm(rtbse_env, rtbse_env%ham_workspace(1), ispin)
1348 1280 : CALL cp_cfm_scale(CMPLX(rtbse_env%spin_degeneracy, 0.0_dp, kind=dp), rtbse_env%ham_workspace(1))
1349 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), ham_effective, &
1350 1280 : CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%ham_workspace(1))
1351 : END IF
1352 : ELSE
1353 : ! ----- TDA: drop B-coupling. Consume the OV-input kernels prebuilt by
1354 : ! build_shared_sex_and_hartree. K_MO[Δρ_VO] = (K_MO[Δρ_OV])^C (real C + symmetric AO kernels),
1355 : ! so evaluate on OV, mask MO to OV, add VO as conjugate transpose.
1356 : ! Signs: (f_n - f_m) = -1 on OV, +1 on VO; applied explicitly. -----
1357 :
1358 : ! SEX: AO->MO, mask OV, stash (assembled after Hartree).
1359 2360 : IF (use_sex) THEN
1360 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%sigma_SEX_ao(ispin), &
1361 2200 : rtbse_env%sigma_SEX(ispin), ispin)
1362 2200 : CALL mask_mo_block_cfm(rtbse_env, rtbse_env%sigma_SEX(ispin), ispin, keep_OV=.TRUE.)
1363 : END IF
1364 :
1365 : ! Hartree: hartree_total_ao is bare (no spin_degeneracy); fold g here (post-mask).
1366 : ! VO = (OV)^C; rho_delta_mo(ispin) is idle scratch.
1367 2360 : IF (use_hartree) THEN
1368 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%hartree_total_ao, &
1369 2200 : rtbse_env%ham_workspace(1), ispin)
1370 2200 : CALL mask_mo_block_cfm(rtbse_env, rtbse_env%ham_workspace(1), ispin, keep_OV=.TRUE.)
1371 2200 : CALL cp_cfm_scale(CMPLX(rtbse_env%spin_degeneracy, 0.0_dp, kind=dp), rtbse_env%ham_workspace(1))
1372 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), ham_effective, &
1373 2200 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%ham_workspace(1))
1374 2200 : CALL cp_cfm_transpose(rtbse_env%ham_workspace(1), 'C', rtbse_env%rho_delta_mo(ispin))
1375 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), ham_effective, &
1376 2200 : CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%rho_delta_mo(ispin))
1377 : END IF
1378 :
1379 : ! SEX assembled (stashed sigma_SEX, OV-masked). OV sign -1; VO = (OV)^C sign +1.
1380 2360 : IF (use_sex) THEN
1381 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), ham_effective, &
1382 2200 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%sigma_SEX(ispin))
1383 2200 : CALL cp_cfm_transpose(rtbse_env%sigma_SEX(ispin), 'C', rtbse_env%rho_delta_mo(ispin))
1384 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), ham_effective, &
1385 2200 : CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%rho_delta_mo(ispin))
1386 : END IF
1387 :
1388 : ! Restore rho_ao_scratch to the AO image of the current full rho, so any
1389 : ! post-routine consumer (output_mos) sees the same invariant.
1390 2360 : CALL transform_mo_to_ao_contravariant_cfm(rtbse_env, rho, rtbse_env%rho_ao_scratch(ispin), ispin)
1391 : END IF
1392 : ! Return the actual RHS f(t,rho) = -i * A(rho) for RK4
1393 3800 : CALL cp_cfm_scale(CMPLX(0.0_dp, -1.0_dp, kind=dp), ham_effective)
1394 :
1395 3800 : CALL timestop(handle)
1396 3800 : END SUBROUTINE update_effective_ham_MO
1397 :
1398 : ! **************************************************************************************************
1399 : !> \brief Single-spin Liouvillian matvec: apply L^sigma_sigma to one spin's Delta rho_MO and return
1400 : !> L * Delta rho_MO in MO basis, OV+VO blocks populated. Drives the n_spin=1 TDA diagnostic
1401 : !> and the (n_spin=1) ABBA diagnostic; the open-shell TDA path uses the array routine
1402 : !> apply_liouvillian_to_drho instead. Detached from the propagator (does NOT touch
1403 : !> rho / rho_orig / ham_effective / ham_reference); only env scratches drho_probe(ispin),
1404 : !> rho_ao_scratch, sigma_SEX_ao, hartree_total_ao, ham_workspace(1), sigma_SEX,
1405 : !> rho_delta_mo(ispin), real_workspace_mo(1) are used.
1406 : !>
1407 : !> Matrix elements correspond to the Casida-A matrix:
1408 : !> L_{ia,jb} = (eps_a - eps_i) delta_{ij} delta_{ab} + (ia|jb) - W_{ij,ab}
1409 : !> with no (f_n - f_m) factor (the propagator path applies -1 on OV; we want raw +K).
1410 : !> Honors rtbse_env%rirs_kernel to dispatch Hartree to the RI-RS grid kernel
1411 : !> (compute_hartree_ri_rs_complex); the SX RI-RS dispatch also reads rirs_kernel inside get_sigma.
1412 : !> AO-RI path: get_hartree_complex / get_sigma.
1413 : !>
1414 : !> \param rtbse_env RT-BSE environment (TDA, n_spin=1).
1415 : !> \param drho_in Input Delta rho_MO (mo_active x mo_active complex).
1416 : !> \param L_drho_out Output L * Delta rho_MO; OV+VO blocks populated; OO/VV zero.
1417 : !> \param ispin Spin index.
1418 : !> \author Maximilian Graml (05.26)
1419 : ! **************************************************************************************************
1420 1076 : SUBROUTINE apply_liouvillian_to_drho_spin(rtbse_env, drho_in, L_drho_out, ispin)
1421 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1422 : TYPE(cp_cfm_type), INTENT(IN) :: drho_in
1423 : TYPE(cp_cfm_type), INTENT(INOUT) :: L_drho_out
1424 : INTEGER, INTENT(IN) :: ispin
1425 :
1426 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_liouvillian_to_drho_spin'
1427 :
1428 : INTEGER :: abs_mo_idx, handle, i_row_global, ii, &
1429 : j_col_global, jj, ncol_local, &
1430 : nrow_local
1431 538 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1432 : LOGICAL :: use_hartree, use_sex
1433 :
1434 538 : CALL timeset(routineN, handle)
1435 :
1436 : ! Mirror update_effective_ham_MO kernel gating: only Hartree + SX are gated
1437 : ! (no COH kernel in linRTBSE).
1438 538 : use_hartree = .NOT. rtbse_env%debug_disable_hartree
1439 538 : use_sex = .NOT. rtbse_env%debug_disable_sex
1440 : ! RI-RS dispatch reads rtbse_env%rirs_kernel directly below.
1441 : ! V_grid / W0_grid are populated by initialize_hartree_potential /
1442 : ! initialize_sex_selfenergy, which run before this diagnostic.
1443 :
1444 : ! 1. Stage drho_in into drho_probe. TDA: mask to OV (propagator carries only OV
1445 : ! and the VO contribution is recovered later via the (.)^C shortcut). ABBA: keep
1446 : ! the full OV+VO content - drho_in carries both blocks independently.
1447 538 : CALL cp_cfm_to_cfm(drho_in, rtbse_env%drho_probe(ispin))
1448 538 : IF (rtbse_env%tda_active) THEN
1449 374 : CALL mask_mo_block_cfm(rtbse_env, rtbse_env%drho_probe(ispin), ispin, keep_OV=.TRUE.)
1450 : END IF
1451 :
1452 : ! 2. Project Delta rho_OV (MO -> AO, contravariant).
1453 : CALL transform_mo_to_ao_contravariant_cfm(rtbse_env, rtbse_env%drho_probe(ispin), &
1454 538 : rtbse_env%rho_ao_scratch(ispin), ispin)
1455 :
1456 : ! 3. Initialize the output accumulator.
1457 538 : CALL cp_cfm_set_all(L_drho_out, CMPLX(0.0_dp, 0.0_dp, kind=dp))
1458 :
1459 : ! 4. eps_{ai} diagonal contribution, kept before the kernel steps (5/6). drho_in aliases
1460 : ! rtbse_env%drho_probe(ispin) at the caller; the kernels' VO Hermitian-conjugate scratch
1461 : ! is rho_delta_mo(ispin) (NOT drho_probe), so they no longer corrupt drho_in - but eps
1462 : ! first is the clean ordering. Build H_eps as a real diagonal fm in real_workspace_mo(1),
1463 : ! convert to cfm in ham_workspace(1), accumulate [drho_in, H_eps] = drho * H - H * drho.
1464 : ! On OV: ([drho, H])_{ia} = (eps_a - eps_i) * drho_{ia} = +eps_{ai} * drho_{ia}.
1465 : ! On VO: -eps_{ai} * drho_{ai} (sign flips); irrelevant - driver only reads OV.
1466 : ! Bare GW eigenvalues (lab frame) so eigenvalues compare 1:1 to bse_full.F.
1467 : CALL cp_fm_get_info(rtbse_env%real_workspace_mo(1), &
1468 : nrow_local=nrow_local, ncol_local=ncol_local, &
1469 538 : row_indices=row_indices, col_indices=col_indices)
1470 538 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(1), 0.0_dp)
1471 3736 : DO ii = 1, nrow_local
1472 3198 : i_row_global = row_indices(ii)
1473 41932 : DO jj = 1, ncol_local
1474 38196 : j_col_global = col_indices(jj)
1475 41394 : IF (i_row_global == j_col_global) THEN
1476 3198 : abs_mo_idx = i_row_global + rtbse_env%first_active_mo - 1
1477 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = &
1478 3198 : rtbse_env%bs_env%eigenval_GW(abs_mo_idx, 1, ispin)
1479 : END IF
1480 : END DO
1481 : END DO
1482 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), &
1483 538 : mtarget=rtbse_env%ham_workspace(1))
1484 : ! drho * H_eps
1485 : CALL cp_cfm_gemm('N', 'N', rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
1486 : CMPLX(1.0_dp, 0.0_dp, kind=dp), drho_in, rtbse_env%ham_workspace(1), &
1487 538 : CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out)
1488 : ! -H_eps * drho
1489 : CALL cp_cfm_gemm('N', 'N', rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
1490 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%ham_workspace(1), drho_in, &
1491 538 : CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out)
1492 :
1493 : ! 5. Hartree contribution: get_hartree_complex (Re/Im split) on AO Delta rho, then
1494 : ! AO->MO covariant, mask OV+VO, accumulate +spin_degeneracy on OV and on VO.
1495 : ! No (f_n - f_m) factor: we want raw +K (Casida convention), not the propagator's -K_OV.
1496 : ! VO contribution comes from the Hermitian conjugate of the OV result (P2 in
1497 : ! rt_bse_pitfalls_physics.md): for real C_active + AO-pair-symmetric kernel,
1498 : ! K_MO[Delta rho_VO] = (K_MO[Delta rho_OV])^C.
1499 538 : IF (use_hartree) THEN
1500 494 : IF (rtbse_env%rirs_kernel) THEN
1501 : ! V^H_λσ = sum_l φ_λ(r_l) v_l φ_σ(r_l), v_l = sum_l' V_ll' n_l' (RI-RS)
1502 : CALL compute_hartree_ri_rs_complex(rtbse_env%bs_env, rtbse_env%rho_ao_scratch(ispin), &
1503 198 : rtbse_env%hartree_total_ao)
1504 : ELSE
1505 : ! V^H_λσ = sum_PQ (λσ|P) V_PQ [sum_µν (µν|Q) Δρ_µν] (AO-RI)
1506 : CALL get_hartree_complex(rtbse_env, rtbse_env%rho_ao_scratch(ispin), &
1507 296 : rtbse_env%hartree_total_ao, ispin)
1508 : END IF
1509 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%hartree_total_ao, &
1510 494 : rtbse_env%ham_workspace(1), ispin)
1511 : CALL add_K_MO_to_L_drho(rtbse_env, rtbse_env%ham_workspace(1), L_drho_out, &
1512 494 : CMPLX(rtbse_env%spin_degeneracy, 0.0_dp, kind=dp), ispin)
1513 : END IF
1514 :
1515 : ! 6. Screened-exchange contribution: get_sigma is already complex-aware. The
1516 : ! -1.0_dp factor passed to get_sigma builds the -W contribution; we then add
1517 : ! +1.0 here (no occupation-factor flip), giving raw K^SX = -W as required by
1518 : ! K = (ia|jb) - W_{ij,ab} (Casida convention with +K).
1519 538 : IF (use_sex) THEN
1520 : ! Σ^SX_λσ = -sum_νQ [sum_µ (λµ|Q) Δρ_µν][sum_P (νσ|P) W_PQ]
1521 494 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX_ao(ispin), -1.0_dp, rtbse_env%rho_ao_scratch(ispin))
1522 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%sigma_SEX_ao(ispin), &
1523 494 : rtbse_env%sigma_SEX(ispin), ispin)
1524 : CALL add_K_MO_to_L_drho(rtbse_env, rtbse_env%sigma_SEX(ispin), L_drho_out, &
1525 494 : CMPLX(1.0_dp, 0.0_dp, kind=dp), ispin)
1526 : END IF
1527 :
1528 538 : CALL timestop(handle)
1529 538 : END SUBROUTINE apply_liouvillian_to_drho_spin
1530 :
1531 : ! **************************************************************************************************
1532 : !> \brief Adds the MO-domain kernel contribution K_MO * Delta rho to L_drho_out, branching on
1533 : !> rtbse_env%tda_active. TDA: masks K_MO to OV in place, adds scale*K_MO on OV, then
1534 : !> builds the VO contribution as (K_MO_OV)^C and adds it (valid only under the TDA
1535 : !> assumption Delta rho_VO = (Delta rho_OV)^C). ABBA: adds the full mo_active x mo_active
1536 : !> K_MO directly - OV+VO blocks are computed naturally from the full OV+VO input.
1537 : !> K_MO is INTENT(INOUT); in the TDA branch it is masked in place (treat as scratch
1538 : !> after this call). rho_delta_mo(ispin) is used as the VO-transpose scratch in the TDA
1539 : !> branch (NOT drho_probe, which the open-shell array driver keeps as the live probe).
1540 : !> \param rtbse_env RT-BSE environment.
1541 : !> \param K_MO Kernel contribution in MO basis (mo_active x mo_active). Scratched in TDA branch.
1542 : !> \param L_drho_out Accumulator (mo_active x mo_active).
1543 : !> \param scale Complex scale factor (spin_degeneracy for Hartree, 1.0 for SX).
1544 : !> \param ispin Spin index.
1545 : !> \author Maximilian Graml (05.26)
1546 : ! **************************************************************************************************
1547 1692 : SUBROUTINE add_K_MO_to_L_drho(rtbse_env, K_MO, L_drho_out, scale, ispin)
1548 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1549 : TYPE(cp_cfm_type), INTENT(INOUT) :: K_MO, L_drho_out
1550 : COMPLEX(kind=dp), INTENT(IN) :: scale
1551 : INTEGER, INTENT(IN) :: ispin
1552 :
1553 : CHARACTER(len=*), PARAMETER :: routineN = 'add_K_MO_to_L_drho'
1554 :
1555 : INTEGER :: handle
1556 :
1557 1692 : CALL timeset(routineN, handle)
1558 :
1559 1692 : IF (rtbse_env%tda_active) THEN
1560 1232 : CALL mask_mo_block_cfm(rtbse_env, K_MO, ispin, keep_OV=.TRUE.)
1561 1232 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out, scale, K_MO)
1562 1232 : CALL cp_cfm_transpose(K_MO, 'C', rtbse_env%rho_delta_mo(ispin))
1563 1232 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out, scale, rtbse_env%rho_delta_mo(ispin))
1564 : ELSE
1565 460 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out, scale, K_MO)
1566 : END IF
1567 :
1568 1692 : CALL timestop(handle)
1569 1692 : END SUBROUTINE add_K_MO_to_L_drho
1570 :
1571 : ! **************************************************************************************************
1572 : !> \brief Open-shell TDA Liouvillian matvec: apply the joint spin-block L_TDA to a per-spin probe
1573 : !> and return L * Delta rho on every spin block. A probe on spin sigma feeds the diagonal
1574 : !> block A^{sigma,sigma} (eps + Coulomb + SX) AND the off-diagonal Coulomb block
1575 : !> A^{sigma',sigma} for every other output spin sigma' (cross-spin Hartree - the term that
1576 : !> produces the singlet/triplet split). AO-RI only: the Phase-D guard forbids RIRS
1577 : !> for n_spin>1, so there is no RIRS branch here; the n_spin=1 diagnostic routes
1578 : !> through apply_liouvillian_to_drho_spin instead (which keeps the RIRS path).
1579 : !> \param rtbse_env RT-BSE environment (TDA).
1580 : !> \param drho_in Per-spin probe Delta rho_MO (mo_active x mo_active each); zero on non-probed spins.
1581 : !> \param L_drho_out Per-spin output L * Delta rho; OV+VO blocks populated, OO/VV zero.
1582 : !> \author Maximilian Graml
1583 : ! **************************************************************************************************
1584 352 : SUBROUTINE apply_liouvillian_to_drho(rtbse_env, drho_in, L_drho_out)
1585 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1586 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: drho_in, L_drho_out
1587 :
1588 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_liouvillian_to_drho'
1589 :
1590 : INTEGER :: abs_mo_idx, handle, i_row_global, ii, &
1591 : isp, isp_out, j_col_global, jj, &
1592 : ncol_local, nrow_local
1593 176 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1594 : LOGICAL :: use_hartree, use_sex
1595 :
1596 176 : CALL timeset(routineN, handle)
1597 :
1598 176 : use_hartree = .NOT. rtbse_env%debug_disable_hartree
1599 176 : use_sex = .NOT. rtbse_env%debug_disable_sex
1600 :
1601 : ! Zero every output spin block before accumulating.
1602 528 : DO isp = 1, rtbse_env%n_spin
1603 528 : CALL cp_cfm_set_all(L_drho_out(isp), CMPLX(0.0_dp, 0.0_dp, kind=dp))
1604 : END DO
1605 :
1606 : ! eps^sigma commutator -> diagonal block only (eps is spin-diagonal). MUST precede the
1607 : ! kernel steps (add_K_MO_to_L_drho scratches MO buffers). Lab-frame bare GW eigenvalues
1608 : ! so they compare 1:1 to bse_full_diag.F.
1609 : CALL cp_fm_get_info(rtbse_env%real_workspace_mo(1), &
1610 : nrow_local=nrow_local, ncol_local=ncol_local, &
1611 176 : row_indices=row_indices, col_indices=col_indices)
1612 528 : DO isp = 1, rtbse_env%n_spin
1613 352 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(1), 0.0_dp)
1614 2464 : DO ii = 1, nrow_local
1615 2112 : i_row_global = row_indices(ii)
1616 27808 : DO jj = 1, ncol_local
1617 25344 : j_col_global = col_indices(jj)
1618 27456 : IF (i_row_global == j_col_global) THEN
1619 2112 : abs_mo_idx = i_row_global + rtbse_env%first_active_mo - 1
1620 : rtbse_env%real_workspace_mo(1)%local_data(ii, jj) = &
1621 2112 : rtbse_env%bs_env%eigenval_GW(abs_mo_idx, 1, isp)
1622 : END IF
1623 : END DO
1624 : END DO
1625 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), &
1626 352 : mtarget=rtbse_env%ham_workspace(1))
1627 : ! [drho, H_eps] = drho * H_eps - H_eps * drho
1628 : CALL cp_cfm_gemm('N', 'N', rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
1629 : CMPLX(1.0_dp, 0.0_dp, kind=dp), drho_in(isp), rtbse_env%ham_workspace(1), &
1630 352 : CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out(isp))
1631 : CALL cp_cfm_gemm('N', 'N', rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
1632 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%ham_workspace(1), drho_in(isp), &
1633 528 : CMPLX(1.0_dp, 0.0_dp, kind=dp), L_drho_out(isp))
1634 : END DO
1635 :
1636 : ! Project per-spin OV-masked rho_ao(sigma) (consumed by SX) and build the single cross-spin
1637 : ! Hartree V_H[spin_degeneracy * sum_sigma rho_ao(sigma)] (Coulomb is spin-blind).
1638 : ! build_shared_hartree_ao does both; skip only when neither kernel is active.
1639 176 : IF (use_hartree .OR. use_sex) THEN
1640 176 : CALL build_shared_hartree_ao(rtbse_env, drho_in, keep_ovvo=.FALSE.)
1641 : END IF
1642 :
1643 : ! Hartree: the one shared V_H read back with each output spin's C fills the diagonal
1644 : ! A^{sigma,sigma} AND the off-diagonal A^{sigma',sigma}. Coeff 1.0 (spin_degeneracy lives in
1645 : ! the summed density); no (f_n - f_m) factor (Casida +K convention).
1646 176 : IF (use_hartree) THEN
1647 528 : DO isp_out = 1, rtbse_env%n_spin
1648 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%hartree_total_ao, &
1649 352 : rtbse_env%ham_workspace(1), isp_out)
1650 : CALL add_K_MO_to_L_drho(rtbse_env, rtbse_env%ham_workspace(1), L_drho_out(isp_out), &
1651 528 : CMPLX(1.0_dp, 0.0_dp, kind=dp), isp_out)
1652 : END DO
1653 : END IF
1654 :
1655 : ! Screened exchange: spin-diagonal (W^sigma acts only within spin sigma). get_sigma builds
1656 : ! -W; add raw +1 (no occupation-factor flip), giving K^SX = -W.
1657 176 : IF (use_sex) THEN
1658 528 : DO isp = 1, rtbse_env%n_spin
1659 : ! Σ^SX_λσ = -sum_νQ [sum_µ (λµ|Q) Δρ_µν][sum_P (νσ|P) W_PQ] (per spin)
1660 352 : CALL get_sigma(rtbse_env, rtbse_env%sigma_SEX_ao(isp), -1.0_dp, rtbse_env%rho_ao_scratch(isp))
1661 : CALL transform_ao_to_mo_covariant_cfm(rtbse_env, rtbse_env%sigma_SEX_ao(isp), &
1662 352 : rtbse_env%sigma_SEX(isp), isp)
1663 : CALL add_K_MO_to_L_drho(rtbse_env, rtbse_env%sigma_SEX(isp), L_drho_out(isp), &
1664 528 : CMPLX(1.0_dp, 0.0_dp, kind=dp), isp)
1665 : END DO
1666 : END IF
1667 :
1668 176 : CALL timestop(handle)
1669 176 : END SUBROUTINE apply_liouvillian_to_drho
1670 :
1671 : ! **************************************************************************************************
1672 : !> \brief Public entry for the Liouvillian eigenvalue diagnostic.
1673 : !> Dispatches to the TDA branch (Casida-A via cp_cfm_heevd) or the ABBA branch
1674 : !> (Furche reduction via cp_cfm_power) based on rtbse_env%tda_active.
1675 : !> Called once at job init from run_propagation_linearized_bse, gated on
1676 : !> rtbse_env%diagnose_liouvillian_eig (n_spin = 1 enforced at env creation).
1677 : !> \param rtbse_env RT-BSE environment with diagnostic scratch already allocated.
1678 : !> \author Maximilian Graml (05.26)
1679 : ! **************************************************************************************************
1680 58 : SUBROUTINE diagnose_liouvillian_eigenvalues(rtbse_env)
1681 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1682 :
1683 : CHARACTER(len=*), PARAMETER :: routineN = 'diagnose_liouvillian_eigenvalues'
1684 :
1685 : INTEGER :: handle
1686 :
1687 58 : CALL timeset(routineN, handle)
1688 :
1689 58 : IF (rtbse_env%tda_active) THEN
1690 40 : CALL diagnose_TDA_liouvillian(rtbse_env)
1691 : ELSE
1692 18 : CALL diagnose_ABBA_liouvillian(rtbse_env)
1693 : END IF
1694 :
1695 58 : CALL timestop(handle)
1696 58 : END SUBROUTINE diagnose_liouvillian_eigenvalues
1697 :
1698 : ! **************************************************************************************************
1699 : !> \brief TDA branch of the Liouvillian eigenvalue diagnostic. Probes the Liouvillian with
1700 : !> canonical OV unit vectors and assembles the joint spin-block Casida-A matrix as L_pairs
1701 : !> (N_OV_joint x N_OV_joint, spin blocks stacked), then diagonalizes via cp_cfm_heevd.
1702 : !> A probe on spin sigma fills its column block and the response on every output spin lands
1703 : !> in that spin's row block (the off-diagonal blocks carry the cross-spin Coulomb that
1704 : !> splits singlet/triplet). The matvec is dispatched on n_spin: n_spin=1 uses the
1705 : !> single-spin apply_liouvillian_to_drho_spin (keeps the RIRS path, bit-identical to the
1706 : !> closed-shell baseline); n_spin=2 uses the AO-RI array apply_liouvillian_to_drho.
1707 : !> Eigenvalues go to stdout (RTBSE|) and to the LIOUVILLIAN_EIG .dat file. Detached from
1708 : !> RK4 state. Called from the dispatcher when tda_active=.TRUE.
1709 : !> \param rtbse_env RT-BSE environment with TDA diagnostic scratch already allocated.
1710 : !> \author Maximilian Graml (05.26)
1711 : ! **************************************************************************************************
1712 40 : SUBROUTINE diagnose_TDA_liouvillian(rtbse_env)
1713 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1714 :
1715 : CHARACTER(len=*), PARAMETER :: routineN = 'diagnose_TDA_liouvillian'
1716 :
1717 40 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ov_block
1718 : INTEGER :: b, eig_unit, handle, j, k_col, k_local, &
1719 : n, n_ov_joint, sigma_out, sigma_probe
1720 40 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_act_occ, n_act_virt, n_ov, off
1721 : REAL(kind=dp) :: residual_max
1722 : TYPE(cp_logger_type), POINTER :: logger
1723 :
1724 40 : CALL timeset(routineN, handle)
1725 40 : logger => cp_get_default_logger()
1726 :
1727 : ! Per-spin OV counts + offsets into the stacked joint Liouvillian. off(1)=0,
1728 : ! off(2)=n_ov(1); n_ov_joint = sum_sigma n_ov(sigma). n_spin=1 -> single block.
1729 : ALLOCATE (n_act_occ(rtbse_env%n_spin), n_act_virt(rtbse_env%n_spin), &
1730 240 : n_ov(rtbse_env%n_spin), off(rtbse_env%n_spin))
1731 40 : n_ov_joint = 0
1732 86 : DO sigma_probe = 1, rtbse_env%n_spin
1733 46 : n_act_occ(sigma_probe) = rtbse_env%n_occ(sigma_probe) - rtbse_env%first_active_mo + 1
1734 46 : n_act_virt(sigma_probe) = rtbse_env%last_active_mo - rtbse_env%n_occ(sigma_probe)
1735 46 : n_ov(sigma_probe) = n_act_occ(sigma_probe)*n_act_virt(sigma_probe)
1736 46 : off(sigma_probe) = n_ov_joint
1737 86 : n_ov_joint = n_ov_joint + n_ov(sigma_probe)
1738 : END DO
1739 :
1740 40 : IF (rtbse_env%unit_nr > 0) THEN
1741 20 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE| ----- TDA Liouvillian diagnostic -----'
1742 : WRITE (rtbse_env%unit_nr, '(A,I0,A,I0)') &
1743 20 : ' RTBSE| n_spin = ', rtbse_env%n_spin, ', joint N_OV = ', n_ov_joint
1744 20 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE|'
1745 : END IF
1746 :
1747 : ! Single joint block to the .dat file (REWIND). ignore_should_output=.TRUE. so this fires
1748 : ! at init regardless of MD-iteration cadence.
1749 : eig_unit = cp_print_key_unit_nr(logger, rtbse_env%eig_section, &
1750 : extension=".dat", &
1751 : file_form="FORMATTED", &
1752 : file_position="REWIND", &
1753 40 : ignore_should_output=.TRUE.)
1754 40 : IF (eig_unit > 0) THEN
1755 20 : WRITE (eig_unit, '(A)') '# Joint spin-block TDA Liouvillian eigenvalues'
1756 20 : IF (rtbse_env%n_spin == 1) THEN
1757 17 : WRITE (eig_unit, '(A,I0,A,I0)') '# n_spin = ', rtbse_env%n_spin, ', N_OV = ', n_ov(1)
1758 : ELSE
1759 3 : WRITE (eig_unit, '(A,I0,A,I0,A,I0)') '# n_spin = ', rtbse_env%n_spin, &
1760 6 : ', N_OV(1) = ', n_ov(1), ', N_OV(2) = ', n_ov(2)
1761 : END IF
1762 : END IF
1763 :
1764 : ! Assemble the joint Casida-A matrix. A probe on spin sigma_probe with canonical OV unit
1765 : ! vector e_{(j,b)} fills column off(sigma_probe)+k_local; the response on each output spin
1766 : ! sigma_out lands in its row block [off(sigma_out)+1 .. +n_ov(sigma_out)]. Block-level
1767 : ! transfers keep this at O(N_OV_joint) collective ops. Column-major OV index
1768 : ! k_local = (b_local-1)*n_act_occ + j_local matches the Fortran layout of ov_block, so
1769 : ! RESHAPE without padding gives the right (n_ov, 1) column.
1770 86 : DO sigma_probe = 1, rtbse_env%n_spin
1771 592 : DO b = rtbse_env%n_occ(sigma_probe) + 1, rtbse_env%last_active_mo
1772 1058 : DO j = rtbse_env%first_active_mo, rtbse_env%n_occ(sigma_probe)
1773 : k_local = (b - rtbse_env%n_occ(sigma_probe) - 1)*n_act_occ(sigma_probe) + &
1774 506 : (j - rtbse_env%first_active_mo + 1)
1775 506 : k_col = off(sigma_probe) + k_local
1776 :
1777 1144 : DO sigma_out = 1, rtbse_env%n_spin
1778 1144 : CALL cp_cfm_set_all(rtbse_env%drho_probe(sigma_out), CMPLX(0.0_dp, 0.0_dp, kind=dp))
1779 : END DO
1780 : CALL cp_cfm_set_element(rtbse_env%drho_probe(sigma_probe), &
1781 : j - rtbse_env%first_active_mo + 1, &
1782 : b - rtbse_env%first_active_mo + 1, &
1783 506 : CMPLX(1.0_dp, 0.0_dp, kind=dp))
1784 :
1785 : ! n_spin=1 keeps the single-spin matvec (RIRS-capable, bit-identical baseline);
1786 : ! n_spin=2 is AO-RI-only (Phase-D guard) -> cross-spin array matvec.
1787 506 : IF (rtbse_env%n_spin == 1) THEN
1788 : CALL apply_liouvillian_to_drho_spin(rtbse_env, rtbse_env%drho_probe(1), &
1789 374 : rtbse_env%L_drho(1), 1)
1790 : ELSE
1791 132 : CALL apply_liouvillian_to_drho(rtbse_env, rtbse_env%drho_probe, rtbse_env%L_drho)
1792 : END IF
1793 :
1794 : ! Stack each output spin's OV response into the joint column.
1795 1650 : DO sigma_out = 1, rtbse_env%n_spin
1796 2552 : ALLOCATE (ov_block(n_act_occ(sigma_out), n_act_virt(sigma_out)))
1797 : CALL cp_cfm_get_submatrix(rtbse_env%L_drho(sigma_out), ov_block, &
1798 : start_row=1, start_col=n_act_occ(sigma_out) + 1, &
1799 638 : n_rows=n_act_occ(sigma_out), n_cols=n_act_virt(sigma_out))
1800 : CALL cp_cfm_set_submatrix(rtbse_env%L_pairs, &
1801 : RESHAPE(ov_block, [n_ov(sigma_out), 1]), &
1802 : start_row=off(sigma_out) + 1, start_col=k_col, &
1803 1914 : n_rows=n_ov(sigma_out), n_cols=1)
1804 1144 : DEALLOCATE (ov_block)
1805 : END DO
1806 : END DO
1807 : END DO
1808 : END DO
1809 :
1810 : ! Hermitian residual on the joint matrix: real-orbital BSE => L real symmetric, so
1811 : ! ||L - L^H||_max should be at the FP floor. D = L^H - L into eigvecs_pairs (idle here;
1812 : ! heevd overwrites it), then its max-element norm via the BLACS-native pzlange path.
1813 40 : CALL cp_cfm_transpose(rtbse_env%L_pairs, 'C', rtbse_env%eigvecs_pairs)
1814 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%eigvecs_pairs, &
1815 40 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%L_pairs)
1816 40 : residual_max = cp_cfm_norm(rtbse_env%eigvecs_pairs, 'M')
1817 40 : IF (rtbse_env%unit_nr > 0) THEN
1818 : WRITE (rtbse_env%unit_nr, '(A,ES16.6)') &
1819 20 : ' RTBSE| Hermitian residual ||L - L^H||_max = ', residual_max
1820 20 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE|'
1821 : END IF
1822 40 : IF (residual_max > 1.0E-6_dp) THEN
1823 0 : CPABORT("Liouvillian Hermitian residual > 1e-6 - check kernel signs / symmetry.")
1824 : END IF
1825 :
1826 : ! Diagonalize the joint matrix. cp_cfm_heevd returns ascending real eigenvalues.
1827 : CALL cp_cfm_heevd(rtbse_env%L_pairs, rtbse_env%eigvecs_pairs, &
1828 40 : rtbse_env%eigenvalues_liouvillian)
1829 :
1830 : ! Stdout table (eV, F12.4 right-aligned to col 80, L-7) + .dat (a.u. + eV).
1831 40 : IF (rtbse_env%unit_nr > 0) THEN
1832 : WRITE (rtbse_env%unit_nr, '(A,T26,A,T59,A)') &
1833 20 : ' RTBSE|', "Excitation index n", "Excitation energy (eV)"
1834 273 : DO n = 1, n_ov_joint
1835 : WRITE (rtbse_env%unit_nr, '(A,T40,I4,T69,F12.4)') &
1836 273 : ' RTBSE|', n, rtbse_env%eigenvalues_liouvillian(n)*evolt
1837 : END DO
1838 20 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE|'
1839 : END IF
1840 40 : IF (eig_unit > 0) THEN
1841 20 : WRITE (eig_unit, '(A)') '# n Omega [a.u.] Omega [eV]'
1842 273 : DO n = 1, n_ov_joint
1843 253 : WRITE (eig_unit, '(I5,4X,ES24.14E3,4X,ES24.14E3)') n, &
1844 253 : rtbse_env%eigenvalues_liouvillian(n), &
1845 526 : rtbse_env%eigenvalues_liouvillian(n)*evolt
1846 : END DO
1847 : END IF
1848 :
1849 40 : CALL cp_print_key_finished_output(eig_unit, logger, rtbse_env%eig_section)
1850 :
1851 40 : DEALLOCATE (n_act_occ, n_act_virt, n_ov, off)
1852 :
1853 40 : CALL timestop(handle)
1854 80 : END SUBROUTINE diagnose_TDA_liouvillian
1855 :
1856 : ! **************************************************************************************************
1857 : !> \brief ABBA branch of the Liouvillian eigenvalue diagnostic. Assembles the joint spin-block
1858 : !> A and B by probing apply_liouvillian_to_drho with canonical OV unit vectors: a probe on
1859 : !> spin sigma_probe fills joint column off(sigma_probe)+k_local; each output spin's OV
1860 : !> response -> A, VO response -> -B^* (recovered by sign-flip + conjugation). One joint
1861 : !> Furche reduction follows: (A-B)>0 gate (independent cp_cfm_heevd), (A-B)^{1/2} via
1862 : !> cp_cfm_power, C = (A-B)^{1/2}(A+B)(A-B)^{1/2}, cp_cfm_heevd, Ω_n = √(C). The
1863 : !> matvec is dispatched on n_spin: n_spin=1 -> apply_liouvillian_to_drho_spin (RIRS-capable,
1864 : !> bit-identical to the closed-shell baseline); n_spin=2 -> the AO-RI cross-spin array
1865 : !> apply_liouvillian_to_drho. For n_spin=1 the routine reduces to the single-block path.
1866 : !> Output: a single joint spectrum (stdout RTBSE| + LIOUVILLIAN_EIG .dat). All eigenvalues
1867 : !> are retained, including optically dark triplet modes (the kernel-correctness gate).
1868 : !> \param rtbse_env RT-BSE environment with ABBA diagnostic scratch allocated.
1869 : !> \author Maximilian Graml (05.26)
1870 : ! **************************************************************************************************
1871 18 : SUBROUTINE diagnose_ABBA_liouvillian(rtbse_env)
1872 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1873 :
1874 : CHARACTER(len=*), PARAMETER :: routineN = 'diagnose_ABBA_liouvillian'
1875 :
1876 18 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ov_block, vo_block
1877 18 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: B_local
1878 : INTEGER :: b, eig_unit, handle, j, k_col, k_local, &
1879 : n, n_ov_joint, sigma_out, sigma_probe
1880 18 : INTEGER, ALLOCATABLE, DIMENSION(:) :: n_act_occ, n_act_virt, n_ov, off
1881 : REAL(kind=dp) :: lambda_min_AmB, residual_A, residual_B
1882 : TYPE(cp_logger_type), POINTER :: logger
1883 :
1884 18 : CALL timeset(routineN, handle)
1885 18 : logger => cp_get_default_logger()
1886 :
1887 : ! Per-spin OV counts + offsets into the stacked joint A/B (same layout as the TDA
1888 : ! diagnostic). off(1)=0, off(2)=n_ov(1); n_ov_joint = sum_sigma n_ov(sigma). For
1889 : ! n_spin=1 this is the single block, bit-identical to the closed-shell ABBA path.
1890 : ALLOCATE (n_act_occ(rtbse_env%n_spin), n_act_virt(rtbse_env%n_spin), &
1891 108 : n_ov(rtbse_env%n_spin), off(rtbse_env%n_spin))
1892 18 : n_ov_joint = 0
1893 38 : DO sigma_probe = 1, rtbse_env%n_spin
1894 20 : n_act_occ(sigma_probe) = rtbse_env%n_occ(sigma_probe) - rtbse_env%first_active_mo + 1
1895 20 : n_act_virt(sigma_probe) = rtbse_env%last_active_mo - rtbse_env%n_occ(sigma_probe)
1896 20 : n_ov(sigma_probe) = n_act_occ(sigma_probe)*n_act_virt(sigma_probe)
1897 20 : off(sigma_probe) = n_ov_joint
1898 38 : n_ov_joint = n_ov_joint + n_ov(sigma_probe)
1899 : END DO
1900 :
1901 18 : IF (rtbse_env%unit_nr > 0) THEN
1902 9 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE| ----- ABBA Liouvillian diagnostic -----'
1903 : WRITE (rtbse_env%unit_nr, '(A,I0,A,I0)') &
1904 9 : ' RTBSE| n_spin = ', rtbse_env%n_spin, ', joint N_OV = ', n_ov_joint
1905 9 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE|'
1906 : END IF
1907 :
1908 : eig_unit = cp_print_key_unit_nr(logger, rtbse_env%eig_section, &
1909 : extension=".dat", &
1910 : file_form="FORMATTED", &
1911 : file_position="REWIND", &
1912 18 : ignore_should_output=.TRUE.)
1913 18 : IF (eig_unit > 0) THEN
1914 9 : WRITE (eig_unit, '(A)') '# Joint spin-block ABBA Liouvillian eigenvalues'
1915 9 : IF (rtbse_env%n_spin == 1) THEN
1916 8 : WRITE (eig_unit, '(A,I0,A,I0)') '# n_spin = ', rtbse_env%n_spin, ', N_OV = ', n_ov(1)
1917 : ELSE
1918 1 : WRITE (eig_unit, '(A,I0,A,I0,A,I0)') '# n_spin = ', rtbse_env%n_spin, &
1919 2 : ', N_OV(1) = ', n_ov(1), ', N_OV(2) = ', n_ov(2)
1920 : END IF
1921 : END IF
1922 :
1923 : ! Assemble joint A and B. Probe on spin sigma_probe, OV pair (j,b) -> joint column
1924 : ! k_col = off(sigma_probe)+k_local (column-major k_local, as TDA). Each output spin's
1925 : ! OV response -> A rows [off(sigma_out)+1 ..]; VO response -> B rows (TRANSPOSE to OV
1926 : ! layout), recovered as -B^* below. The n_spin=2 array matvec produces the full OV+VO
1927 : ! readout (add_K_MO_to_L_drho .NOT.tda_active branch); one cross-spin Hartree fills
1928 : ! both A^{s,s'} and B^{s,s'} Coulomb; SX stays spin-diagonal.
1929 38 : DO sigma_probe = 1, rtbse_env%n_spin
1930 246 : DO b = rtbse_env%n_occ(sigma_probe) + 1, rtbse_env%last_active_mo
1931 436 : DO j = rtbse_env%first_active_mo, rtbse_env%n_occ(sigma_probe)
1932 : k_local = (b - rtbse_env%n_occ(sigma_probe) - 1)*n_act_occ(sigma_probe) + &
1933 208 : (j - rtbse_env%first_active_mo + 1)
1934 208 : k_col = off(sigma_probe) + k_local
1935 :
1936 460 : DO sigma_out = 1, rtbse_env%n_spin
1937 460 : CALL cp_cfm_set_all(rtbse_env%drho_probe(sigma_out), CMPLX(0.0_dp, 0.0_dp, kind=dp))
1938 : END DO
1939 : CALL cp_cfm_set_element(rtbse_env%drho_probe(sigma_probe), &
1940 : j - rtbse_env%first_active_mo + 1, &
1941 : b - rtbse_env%first_active_mo + 1, &
1942 208 : CMPLX(1.0_dp, 0.0_dp, kind=dp))
1943 :
1944 208 : IF (rtbse_env%n_spin == 1) THEN
1945 : CALL apply_liouvillian_to_drho_spin(rtbse_env, rtbse_env%drho_probe(1), &
1946 164 : rtbse_env%L_drho(1), 1)
1947 : ELSE
1948 44 : CALL apply_liouvillian_to_drho(rtbse_env, rtbse_env%drho_probe, rtbse_env%L_drho)
1949 : END IF
1950 :
1951 668 : DO sigma_out = 1, rtbse_env%n_spin
1952 1008 : ALLOCATE (ov_block(n_act_occ(sigma_out), n_act_virt(sigma_out)))
1953 756 : ALLOCATE (vo_block(n_act_virt(sigma_out), n_act_occ(sigma_out)))
1954 : CALL cp_cfm_get_submatrix(rtbse_env%L_drho(sigma_out), ov_block, &
1955 : start_row=1, start_col=n_act_occ(sigma_out) + 1, &
1956 252 : n_rows=n_act_occ(sigma_out), n_cols=n_act_virt(sigma_out))
1957 : CALL cp_cfm_set_submatrix(rtbse_env%A_mat, &
1958 : RESHAPE(ov_block, [n_ov(sigma_out), 1]), &
1959 : start_row=off(sigma_out) + 1, start_col=k_col, &
1960 756 : n_rows=n_ov(sigma_out), n_cols=1)
1961 : ! TRANSPOSE puts vo_block in (n_act_occ, n_act_virt) layout, matching the OV pack.
1962 : CALL cp_cfm_get_submatrix(rtbse_env%L_drho(sigma_out), vo_block, &
1963 : start_row=n_act_occ(sigma_out) + 1, start_col=1, &
1964 252 : n_rows=n_act_virt(sigma_out), n_cols=n_act_occ(sigma_out))
1965 : CALL cp_cfm_set_submatrix(rtbse_env%B_mat, &
1966 : RESHAPE(TRANSPOSE(vo_block), [n_ov(sigma_out), 1]), &
1967 : start_row=off(sigma_out) + 1, start_col=k_col, &
1968 756 : n_rows=n_ov(sigma_out), n_cols=1)
1969 460 : DEALLOCATE (ov_block, vo_block)
1970 : END DO
1971 : END DO
1972 : END DO
1973 : END DO
1974 :
1975 : ! Recover B from -B^* via rank-local pass on the cfm's MPI-local data (documented
1976 : ! exception to the fm/cfm-routines-only rule); B_recovered = -CONJG(stored).
1977 18 : B_local => rtbse_env%B_mat%local_data
1978 1582 : B_local = -CONJG(B_local)
1979 :
1980 : ! Block-symmetry residuals on the JOINT matrices: A Hermitian, B real-symmetric.
1981 18 : CALL cp_cfm_transpose(rtbse_env%A_mat, 'C', rtbse_env%eigvecs_pairs)
1982 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%eigvecs_pairs, &
1983 18 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%A_mat)
1984 18 : residual_A = cp_cfm_norm(rtbse_env%eigvecs_pairs, 'M')
1985 :
1986 18 : CALL cp_cfm_transpose(rtbse_env%B_mat, 'T', rtbse_env%eigvecs_pairs)
1987 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%eigvecs_pairs, &
1988 18 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%B_mat)
1989 18 : residual_B = cp_cfm_norm(rtbse_env%eigvecs_pairs, 'M')
1990 :
1991 18 : IF (rtbse_env%unit_nr > 0) THEN
1992 : WRITE (rtbse_env%unit_nr, '(A,ES16.6)') &
1993 9 : ' RTBSE| Hermitian residual ||A - A^H||_max = ', residual_A
1994 : WRITE (rtbse_env%unit_nr, '(A,ES16.6)') &
1995 9 : ' RTBSE| Symmetry residual ||B - B^T||_max = ', residual_B
1996 : END IF
1997 18 : IF (residual_A > 1.0E-6_dp) THEN
1998 0 : CPABORT("A is not Hermitian within 1e-6 - check kernel signs / symmetry.")
1999 : END IF
2000 18 : IF (residual_B > 1.0E-6_dp) THEN
2001 0 : CPABORT("B is not symmetric within 1e-6 - check kernel signs / symmetry.")
2002 : END IF
2003 :
2004 : ! A +/- B in scratches (both Hermitian).
2005 18 : CALL cp_cfm_to_cfm(rtbse_env%A_mat, rtbse_env%AmB_scratch)
2006 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%AmB_scratch, &
2007 18 : CMPLX(-1.0_dp, 0.0_dp, kind=dp), rtbse_env%B_mat)
2008 18 : CALL cp_cfm_to_cfm(rtbse_env%A_mat, rtbse_env%ApB_scratch)
2009 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%ApB_scratch, &
2010 18 : CMPLX(1.0_dp, 0.0_dp, kind=dp), rtbse_env%B_mat)
2011 :
2012 : ! (A-B) positivity gate on the joint matrix: heevd on a copy; abort if lambda_min < 0.
2013 18 : CALL cp_cfm_to_cfm(rtbse_env%AmB_scratch, rtbse_env%L_pairs)
2014 18 : CALL cp_cfm_heevd(rtbse_env%L_pairs, rtbse_env%eigvecs_pairs, rtbse_env%eigenvalues_liouvillian)
2015 18 : lambda_min_AmB = rtbse_env%eigenvalues_liouvillian(1)
2016 18 : IF (rtbse_env%unit_nr > 0) THEN
2017 : WRITE (rtbse_env%unit_nr, '(A,ES16.6,A,F12.6,A)') &
2018 9 : ' RTBSE| lambda_min(A - B) = ', &
2019 18 : lambda_min_AmB, ' a.u. (', lambda_min_AmB*evolt, ' eV)'
2020 : END IF
2021 : ! Hard abort by design: a non-positive (A-B) breaks the Furche reduction.
2022 18 : IF (lambda_min_AmB < 0.0_dp) THEN
2023 : CALL cp_abort(__LOCATION__, &
2024 : "(A - B) not positive definite - this may hint at a triplet or "// &
2025 0 : "charge-transfer instability of the reference state.")
2026 : END IF
2027 :
2028 : ! In-place AmB_scratch -> (A-B)^{1/2} via cp_cfm_power. threshold=0 substitution
2029 : ! codepath is unreachable here (lambda_min_AmB > 0 already enforced above).
2030 18 : CALL cp_cfm_power(rtbse_env%AmB_scratch, threshold=0.0_dp, exponent=0.5_dp)
2031 :
2032 : ! C = (A-B)^{1/2}(A+B)(A-B)^{1/2}. B_mat free after step "A+/-B" -> reuse as T scratch.
2033 : CALL cp_cfm_gemm('N', 'N', n_ov_joint, n_ov_joint, n_ov_joint, &
2034 : CMPLX(1.0_dp, 0.0_dp, kind=dp), &
2035 : rtbse_env%AmB_scratch, rtbse_env%ApB_scratch, &
2036 : CMPLX(0.0_dp, 0.0_dp, kind=dp), &
2037 18 : rtbse_env%B_mat)
2038 : CALL cp_cfm_gemm('N', 'N', n_ov_joint, n_ov_joint, n_ov_joint, &
2039 : CMPLX(1.0_dp, 0.0_dp, kind=dp), &
2040 : rtbse_env%B_mat, rtbse_env%AmB_scratch, &
2041 : CMPLX(0.0_dp, 0.0_dp, kind=dp), &
2042 18 : rtbse_env%L_pairs)
2043 :
2044 : ! Diagonalize C -> Ω_n^2; take +√. Safety clamp on tiny-negative noise.
2045 18 : CALL cp_cfm_heevd(rtbse_env%L_pairs, rtbse_env%eigvecs_pairs, rtbse_env%eigenvalues_liouvillian)
2046 226 : DO n = 1, n_ov_joint
2047 226 : IF (rtbse_env%eigenvalues_liouvillian(n) < 0.0_dp) THEN
2048 0 : rtbse_env%eigenvalues_liouvillian(n) = 0.0_dp
2049 : ELSE
2050 208 : rtbse_env%eigenvalues_liouvillian(n) = SQRT(rtbse_env%eigenvalues_liouvillian(n))
2051 : END IF
2052 : END DO
2053 :
2054 18 : IF (rtbse_env%unit_nr > 0) THEN
2055 9 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE|'
2056 : WRITE (rtbse_env%unit_nr, '(A,T26,A,T59,A)') &
2057 9 : ' RTBSE|', "Excitation index n", "Excitation energy (eV)"
2058 113 : DO n = 1, n_ov_joint
2059 : WRITE (rtbse_env%unit_nr, '(A,T40,I4,T69,F12.4)') &
2060 113 : ' RTBSE|', n, rtbse_env%eigenvalues_liouvillian(n)*evolt
2061 : END DO
2062 9 : WRITE (rtbse_env%unit_nr, '(A)') ' RTBSE|'
2063 : END IF
2064 18 : IF (eig_unit > 0) THEN
2065 9 : WRITE (eig_unit, '(A)') '# n Omega [a.u.] Omega [eV]'
2066 113 : DO n = 1, n_ov_joint
2067 104 : WRITE (eig_unit, '(I5,4X,ES24.14E3,4X,ES24.14E3)') n, &
2068 104 : rtbse_env%eigenvalues_liouvillian(n), &
2069 217 : rtbse_env%eigenvalues_liouvillian(n)*evolt
2070 : END DO
2071 : END IF
2072 :
2073 18 : CALL cp_print_key_finished_output(eig_unit, logger, rtbse_env%eig_section)
2074 :
2075 18 : DEALLOCATE (n_act_occ, n_act_virt, n_ov, off)
2076 :
2077 18 : CALL timestop(handle)
2078 36 : END SUBROUTINE diagnose_ABBA_liouvillian
2079 :
2080 : ! **************************************************************************************************
2081 : !> \brief Covariant AO->MO transform of an operator (real): M^MO_mn = sum_µν C_µm M^AO_µν C_νn.
2082 : !> For operator-like kernels (Σ^SX, V^H); the density uses the contravariant routine.
2083 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
2084 : !> \param fm_ao operator in the AO basis (n_ao x n_ao), input
2085 : !> \param fm_mo operator in the active-MO basis (mo_active x mo_active), output
2086 : !> \param i_spin spin channel σ; selects C_active(σ)
2087 : ! **************************************************************************************************
2088 410 : SUBROUTINE transform_ao_to_mo_covariant_fm(rtbse_env, fm_ao, fm_mo, i_spin)
2089 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2090 : TYPE(cp_fm_type) :: fm_ao, fm_mo
2091 : INTEGER, INTENT(IN) :: i_spin
2092 :
2093 : CHARACTER(len=*), PARAMETER :: routineN = 'transform_ao_to_mo_covariant_fm'
2094 :
2095 : INTEGER :: handle
2096 :
2097 410 : CALL timeset(routineN, handle)
2098 :
2099 : ! step 1: T_µn = sum_ν M^AO_µν C_νn (n_ao x mo_active)
2100 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%mo_active, rtbse_env%n_ao, &
2101 : 1.0_dp, fm_ao, rtbse_env%C_active(i_spin), &
2102 410 : 0.0_dp, rtbse_env%ao_mo_workspace(1))
2103 : ! step 2: M^MO_mn = sum_µ C_µm T_µn (mo_active x mo_active)
2104 : CALL parallel_gemm("T", "N", rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%n_ao, &
2105 : 1.0_dp, rtbse_env%C_active(i_spin), rtbse_env%ao_mo_workspace(1), &
2106 410 : 0.0_dp, fm_mo)
2107 :
2108 410 : CALL timestop(handle)
2109 410 : END SUBROUTINE transform_ao_to_mo_covariant_fm
2110 :
2111 : ! **************************************************************************************************
2112 : !> \brief Covariant AO->MO transform (complex): M^MO_mn = sum_µν C_µm M^AO_µν C_νn, applied to the
2113 : !> real and imaginary AO parts separately (2x the real cost).
2114 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
2115 : !> \param fm_ao operator in the AO basis (n_ao x n_ao) cfm, input
2116 : !> \param fm_mo operator in the active-MO basis (mo_active x mo_active) cfm, output
2117 : !> \param i_spin spin channel σ; selects C_active(σ)
2118 : ! **************************************************************************************************
2119 8666 : SUBROUTINE transform_ao_to_mo_covariant_cfm(rtbse_env, fm_ao, fm_mo, i_spin)
2120 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2121 : TYPE(cp_cfm_type) :: fm_ao, fm_mo
2122 : INTEGER, INTENT(IN) :: i_spin
2123 :
2124 : CHARACTER(len=*), PARAMETER :: routineN = 'transform_ao_to_mo_covariant_cfm'
2125 :
2126 : INTEGER :: handle
2127 :
2128 8666 : CALL timeset(routineN, handle)
2129 :
2130 : ! Decompose into real/imag AO-sized parts
2131 : CALL cp_cfm_to_fm(msource=fm_ao, mtargetr=rtbse_env%real_workspace(1), &
2132 8666 : mtargeti=rtbse_env%real_workspace(2))
2133 : ! Re(M^MO)_mn = sum_µν C_µm Re(M^AO)_µν C_νn (two gemms via ao_mo_workspace)
2134 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%mo_active, rtbse_env%n_ao, &
2135 : 1.0_dp, rtbse_env%real_workspace(1), rtbse_env%C_active(i_spin), &
2136 8666 : 0.0_dp, rtbse_env%ao_mo_workspace(1))
2137 : CALL parallel_gemm("T", "N", rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%n_ao, &
2138 : 1.0_dp, rtbse_env%C_active(i_spin), rtbse_env%ao_mo_workspace(1), &
2139 8666 : 0.0_dp, rtbse_env%real_workspace_mo(1))
2140 : ! Im(M^MO)_mn = sum_µν C_µm Im(M^AO)_µν C_νn
2141 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%mo_active, rtbse_env%n_ao, &
2142 : 1.0_dp, rtbse_env%real_workspace(2), rtbse_env%C_active(i_spin), &
2143 8666 : 0.0_dp, rtbse_env%ao_mo_workspace(1))
2144 : CALL parallel_gemm("T", "N", rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%n_ao, &
2145 : 1.0_dp, rtbse_env%C_active(i_spin), rtbse_env%ao_mo_workspace(1), &
2146 8666 : 0.0_dp, rtbse_env%real_workspace_mo(2))
2147 : ! Reassemble into MO-sized cfm
2148 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), &
2149 : msourcei=rtbse_env%real_workspace_mo(2), &
2150 8666 : mtarget=fm_mo)
2151 :
2152 8666 : CALL timestop(handle)
2153 8666 : END SUBROUTINE transform_ao_to_mo_covariant_cfm
2154 :
2155 : ! **************************************************************************************************
2156 : !> \brief Contravariant MO->AO transform of the density (complex): Δρ^AO_µν = sum_mn C_µm Δρ^MO_mn C_νn.
2157 : !> Density-like (expands MO indices), unlike the covariant operator transform.
2158 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
2159 : !> \param fm_mo density in the active-MO basis (mo_active x mo_active) cfm, input
2160 : !> \param fm_ao density in the AO basis (n_ao x n_ao) cfm, output
2161 : !> \param i_spin spin channel σ; selects C_active(σ)
2162 : ! **************************************************************************************************
2163 7116 : SUBROUTINE transform_mo_to_ao_contravariant_cfm(rtbse_env, fm_mo, fm_ao, i_spin)
2164 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2165 : TYPE(cp_cfm_type) :: fm_mo, fm_ao
2166 : INTEGER, INTENT(IN) :: i_spin
2167 :
2168 : CHARACTER(len=*), PARAMETER :: routineN = 'transform_mo_to_ao_contravariant_cfm'
2169 :
2170 : INTEGER :: handle
2171 :
2172 7116 : CALL timeset(routineN, handle)
2173 :
2174 : ! Re/Im split of Δρ^MO (mo_active x mo_active) into persistent MO-sized scratch
2175 : CALL cp_cfm_to_fm(msource=fm_mo, mtargetr=rtbse_env%real_workspace_mo(1), &
2176 7116 : mtargeti=rtbse_env%real_workspace_mo(2))
2177 : ! Re(Δρ^AO)_µν = sum_mn C_µm Re(Δρ^MO)_mn C_νn (C·ρ via ao_mo_workspace, then ·C^T)
2178 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%mo_active, rtbse_env%mo_active, &
2179 : 1.0_dp, rtbse_env%C_active(i_spin), rtbse_env%real_workspace_mo(1), &
2180 7116 : 0.0_dp, rtbse_env%ao_mo_workspace(1))
2181 : CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%mo_active, &
2182 : 1.0_dp, rtbse_env%ao_mo_workspace(1), rtbse_env%C_active(i_spin), &
2183 7116 : 0.0_dp, rtbse_env%real_workspace(1))
2184 : ! Im(Δρ^AO)_µν = sum_mn C_µm Im(Δρ^MO)_mn C_νn
2185 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%mo_active, rtbse_env%mo_active, &
2186 : 1.0_dp, rtbse_env%C_active(i_spin), rtbse_env%real_workspace_mo(2), &
2187 7116 : 0.0_dp, rtbse_env%ao_mo_workspace(1))
2188 : CALL parallel_gemm("N", "T", rtbse_env%n_ao, rtbse_env%n_ao, rtbse_env%mo_active, &
2189 : 1.0_dp, rtbse_env%ao_mo_workspace(1), rtbse_env%C_active(i_spin), &
2190 7116 : 0.0_dp, rtbse_env%real_workspace(2))
2191 : ! Reassemble into AO-sized cfm
2192 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(1), &
2193 7116 : msourcei=rtbse_env%real_workspace(2), mtarget=fm_ao)
2194 :
2195 7116 : CALL timestop(handle)
2196 7116 : END SUBROUTINE transform_mo_to_ao_contravariant_cfm
2197 :
2198 : ! **************************************************************************************************
2199 : !> \brief Scales each MO-active element (m,n) by the occupation prefactor (f_n - f_m), f in {0,1}
2200 : !> (complex): OV -> -1, VO -> +1, OO/VV -> 0. The only trace of "ρ^0 diagonal in MO".
2201 : !> Applied to whatever kernel the caller passes (Σ^SX, V^H in MO) - bound by the caller.
2202 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
2203 : !> \param cfm MO-active kernel matrix (mo_active x mo_active) cfm, scaled in place
2204 : !> \param i_spin spin channel σ; OV/VO boundary set by n_occ(σ)
2205 : ! **************************************************************************************************
2206 2574 : SUBROUTINE transform_mo_occupation_factor_diff_cfm(rtbse_env, cfm, i_spin)
2207 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2208 : TYPE(cp_cfm_type) :: cfm
2209 : INTEGER :: i_spin
2210 :
2211 : CHARACTER(len=*), PARAMETER :: routineN = 'transform_mo_occupation_factor_diff_cfm'
2212 :
2213 : INTEGER :: handle
2214 :
2215 2574 : CALL timeset(routineN, handle)
2216 :
2217 : CALL cp_cfm_to_fm(msource=cfm, mtargetr=rtbse_env%real_workspace_mo(1), &
2218 2574 : mtargeti=rtbse_env%real_workspace_mo(2))
2219 : ! (f_n - f_m) applied to the real part
2220 2574 : CALL transform_mo_occupation_factor_diff_fm(rtbse_env, rtbse_env%real_workspace_mo(1), i_spin)
2221 : ! (f_n - f_m) applied to the imaginary part
2222 2574 : CALL transform_mo_occupation_factor_diff_fm(rtbse_env, rtbse_env%real_workspace_mo(2), i_spin)
2223 : ! Copy back to cfm
2224 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), &
2225 : msourcei=rtbse_env%real_workspace_mo(2), &
2226 2574 : mtarget=cfm)
2227 :
2228 2574 : CALL timestop(handle)
2229 2574 : END SUBROUTINE transform_mo_occupation_factor_diff_cfm
2230 :
2231 : ! **************************************************************************************************
2232 : !> \brief Scales each MO-active element (m,n) by the occupation prefactor (f_n - f_m), f in {0,1}
2233 : !> (real): OV -> -1, VO -> +1, OO/VV -> 0. Real-input worker for the cfm variant; the only
2234 : !> trace of "ρ^0 diagonal in MO". Bound by the caller to the kernel being scaled.
2235 : !> \param rtbse_env Entry point of the calculation - contains current state of variables
2236 : !> \param fm MO-active kernel matrix (mo_active x mo_active) fm, scaled in place
2237 : !> \param i_spin spin channel σ; OV/VO boundary set by n_occ(σ)
2238 : ! **************************************************************************************************
2239 10324 : SUBROUTINE transform_mo_occupation_factor_diff_fm(rtbse_env, fm, i_spin)
2240 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2241 : TYPE(cp_fm_type) :: fm
2242 : INTEGER :: i_spin
2243 :
2244 : CHARACTER(len=*), PARAMETER :: routineN = 'transform_mo_occupation_factor_diff_fm'
2245 :
2246 : INTEGER :: handle, i_global, i_global_mo, i_local, &
2247 : j_global, j_global_mo, j_local, n_occ, &
2248 : ncol_local, nrow_local, shift
2249 5162 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2250 : REAL(kind=dp) :: occ_factor
2251 5162 : REAL(kind=dp), DIMENSION(:, :), POINTER :: local_data
2252 :
2253 5162 : CALL timeset(routineN, handle)
2254 :
2255 5162 : n_occ = rtbse_env%n_occ(i_spin)
2256 : ! Shift mapping local active-window index to absolute MO index
2257 5162 : shift = rtbse_env%first_active_mo - 1
2258 :
2259 : CALL cp_fm_get_info(matrix=fm, &
2260 : nrow_local=nrow_local, ncol_local=ncol_local, &
2261 5162 : row_indices=row_indices, col_indices=col_indices)
2262 :
2263 5162 : local_data => fm%local_data
2264 :
2265 34196 : DO i_local = 1, nrow_local
2266 29034 : i_global = row_indices(i_local)
2267 29034 : i_global_mo = i_global + shift
2268 370976 : DO j_local = 1, ncol_local
2269 336780 : j_global = col_indices(j_local)
2270 336780 : j_global_mo = j_global + shift
2271 :
2272 336780 : IF (i_global_mo <= n_occ .AND. j_global_mo > n_occ) THEN
2273 : occ_factor = -1.0_dp
2274 310327 : ELSE IF (i_global_mo > n_occ .AND. j_global_mo <= n_occ) THEN
2275 : occ_factor = 1.0_dp
2276 : ELSE
2277 283874 : occ_factor = 0.0_dp
2278 : END IF
2279 :
2280 365814 : local_data(i_local, j_local) = occ_factor*local_data(i_local, j_local)
2281 : END DO
2282 : END DO
2283 :
2284 5162 : CALL timestop(handle)
2285 5162 : END SUBROUTINE transform_mo_occupation_factor_diff_fm
2286 :
2287 : ! **************************************************************************************************
2288 : !> \brief Mask an MO-active cfm: keep either OV or VO block, zero everything else.
2289 : !> \param rtbse_env RT-BSE environment
2290 : !> \param cfm MO-active cfm to mask in place
2291 : !> \param i_spin Spin index
2292 : !> \param keep_OV .TRUE. keeps the (occ row, virt col) block; .FALSE. keeps (virt row, occ col)
2293 : !> \param keep_ovvo if present and .TRUE., keep BOTH off-diagonal blocks (OV and VO) and zero
2294 : !> OO/VV; overrides keep_OV. Absent/false reproduces the keep_OV behaviour.
2295 : ! **************************************************************************************************
2296 18076 : SUBROUTINE mask_mo_block_cfm(rtbse_env, cfm, i_spin, keep_OV, keep_ovvo)
2297 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2298 : TYPE(cp_cfm_type) :: cfm
2299 : INTEGER, INTENT(IN) :: i_spin
2300 : LOGICAL, INTENT(IN) :: keep_OV
2301 : LOGICAL, INTENT(IN), OPTIONAL :: keep_ovvo
2302 :
2303 : CHARACTER(len=*), PARAMETER :: routineN = 'mask_mo_block_cfm'
2304 :
2305 9038 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: local_data
2306 : INTEGER :: handle, i_global, i_global_mo, i_local, &
2307 : j_global, j_global_mo, j_local, n_occ, &
2308 : ncol_local, nrow_local, shift
2309 9038 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2310 : LOGICAL :: keep, l_keep_ovvo
2311 :
2312 9038 : CALL timeset(routineN, handle)
2313 :
2314 9038 : l_keep_ovvo = .FALSE.
2315 9038 : IF (PRESENT(keep_ovvo)) l_keep_ovvo = keep_ovvo
2316 :
2317 9038 : n_occ = rtbse_env%n_occ(i_spin)
2318 9038 : shift = rtbse_env%first_active_mo - 1
2319 :
2320 : CALL cp_cfm_get_info(matrix=cfm, &
2321 : nrow_local=nrow_local, ncol_local=ncol_local, &
2322 9038 : row_indices=row_indices, col_indices=col_indices)
2323 :
2324 9038 : local_data => cfm%local_data
2325 :
2326 63266 : DO i_local = 1, nrow_local
2327 54228 : i_global = row_indices(i_local)
2328 54228 : i_global_mo = i_global + shift
2329 714002 : DO j_local = 1, ncol_local
2330 650736 : j_global = col_indices(j_local)
2331 650736 : j_global_mo = j_global + shift
2332 650736 : IF (l_keep_ovvo) THEN
2333 : ! keep both off-diagonal blocks (OV and VO); drop OO/VV
2334 23040 : keep = ((i_global_mo <= n_occ) .NEQV. (j_global_mo <= n_occ))
2335 627696 : ELSE IF (keep_OV) THEN
2336 627696 : keep = (i_global_mo <= n_occ .AND. j_global_mo > n_occ)
2337 : ELSE
2338 0 : keep = (i_global_mo > n_occ .AND. j_global_mo <= n_occ)
2339 : END IF
2340 704964 : IF (.NOT. keep) local_data(i_local, j_local) = CMPLX(0.0_dp, 0.0_dp, kind=dp)
2341 : END DO
2342 : END DO
2343 :
2344 9038 : CALL timestop(handle)
2345 9038 : END SUBROUTINE mask_mo_block_cfm
2346 :
2347 : ! **************************************************************************************************
2348 : !> \brief Mask an MO-active fm: keep either OV or VO block, zero everything else.
2349 : !> \param rtbse_env RT-BSE environment
2350 : !> \param fm MO-active fm to mask in place
2351 : !> \param i_spin Spin index
2352 : !> \param keep_OV .TRUE. keeps the (occ row, virt col) block; .FALSE. keeps (virt row, occ col)
2353 : ! **************************************************************************************************
2354 0 : SUBROUTINE mask_mo_block_fm(rtbse_env, fm, i_spin, keep_OV)
2355 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2356 : TYPE(cp_fm_type) :: fm
2357 : INTEGER, INTENT(IN) :: i_spin
2358 : LOGICAL, INTENT(IN) :: keep_OV
2359 :
2360 : CHARACTER(len=*), PARAMETER :: routineN = 'mask_mo_block_fm'
2361 :
2362 : INTEGER :: handle, i_global, i_global_mo, i_local, &
2363 : j_global, j_global_mo, j_local, n_occ, &
2364 : ncol_local, nrow_local, shift
2365 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2366 : LOGICAL :: keep
2367 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: local_data
2368 :
2369 0 : CALL timeset(routineN, handle)
2370 :
2371 0 : n_occ = rtbse_env%n_occ(i_spin)
2372 0 : shift = rtbse_env%first_active_mo - 1
2373 :
2374 : CALL cp_fm_get_info(matrix=fm, &
2375 : nrow_local=nrow_local, ncol_local=ncol_local, &
2376 0 : row_indices=row_indices, col_indices=col_indices)
2377 :
2378 0 : local_data => fm%local_data
2379 :
2380 0 : DO i_local = 1, nrow_local
2381 0 : i_global = row_indices(i_local)
2382 0 : i_global_mo = i_global + shift
2383 0 : DO j_local = 1, ncol_local
2384 0 : j_global = col_indices(j_local)
2385 0 : j_global_mo = j_global + shift
2386 0 : IF (keep_OV) THEN
2387 0 : keep = (i_global_mo <= n_occ .AND. j_global_mo > n_occ)
2388 : ELSE
2389 0 : keep = (i_global_mo > n_occ .AND. j_global_mo <= n_occ)
2390 : END IF
2391 0 : IF (.NOT. keep) local_data(i_local, j_local) = 0.0_dp
2392 : END DO
2393 : END DO
2394 :
2395 0 : CALL timestop(handle)
2396 0 : END SUBROUTINE mask_mo_block_fm
2397 :
2398 : ! **************************************************************************************************
2399 : !> \brief Multiply a MO-active cfm rho by the TDA symmetric-shift rotation:
2400 : !> rho_OV *= exp(i * phase), rho_VO *= exp(-i * phase).
2401 : !> OO/VV blocks are left untouched. For direction='to_lab' pass phase = -Ω_0*t;
2402 : !> for direction='to_rotating' pass phase = +Ω_0*t. No-op when omega_shift = 0.
2403 : !> Preserves Hermiticity since the two phases are complex conjugates of each other.
2404 : !> \param rtbse_env RT-BSE environment
2405 : !> \param rho MO-active cfm rotated in place
2406 : !> \param i_spin Spin index
2407 : !> \param phase Real phase argument (radians); typically +/- Ω_0 * t
2408 : ! **************************************************************************************************
2409 164 : SUBROUTINE rotate_rho_phase(rtbse_env, rho, i_spin, phase)
2410 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2411 : TYPE(cp_cfm_type) :: rho
2412 : INTEGER, INTENT(IN) :: i_spin
2413 : REAL(kind=dp), INTENT(IN) :: phase
2414 :
2415 : CHARACTER(len=*), PARAMETER :: routineN = 'rotate_rho_phase'
2416 :
2417 : COMPLEX(kind=dp) :: phase_ov, phase_vo
2418 64 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: local_data
2419 : INTEGER :: handle, i_global, i_global_mo, i_local, &
2420 : j_global, j_global_mo, j_local, n_occ, &
2421 : ncol_local, nrow_local, shift
2422 64 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2423 :
2424 14 : IF (rtbse_env%omega_shift == 0.0_dp) RETURN
2425 :
2426 50 : CALL timeset(routineN, handle)
2427 :
2428 50 : n_occ = rtbse_env%n_occ(i_spin)
2429 50 : shift = rtbse_env%first_active_mo - 1
2430 50 : phase_ov = CMPLX(COS(phase), SIN(phase), kind=dp)
2431 50 : phase_vo = CONJG(phase_ov)
2432 :
2433 : CALL cp_cfm_get_info(matrix=rho, &
2434 : nrow_local=nrow_local, ncol_local=ncol_local, &
2435 50 : row_indices=row_indices, col_indices=col_indices)
2436 50 : local_data => rho%local_data
2437 :
2438 350 : DO i_local = 1, nrow_local
2439 300 : i_global = row_indices(i_local)
2440 300 : i_global_mo = i_global + shift
2441 3950 : DO j_local = 1, ncol_local
2442 3600 : j_global = col_indices(j_local)
2443 3600 : j_global_mo = j_global + shift
2444 3900 : IF (i_global_mo <= n_occ .AND. j_global_mo > n_occ) THEN
2445 : ! ρ_OV *= e^{+iφ}
2446 275 : local_data(i_local, j_local) = phase_ov*local_data(i_local, j_local)
2447 3325 : ELSE IF (i_global_mo > n_occ .AND. j_global_mo <= n_occ) THEN
2448 : ! ρ_VO *= e^{-iφ}
2449 275 : local_data(i_local, j_local) = phase_vo*local_data(i_local, j_local)
2450 : END IF
2451 : END DO
2452 : END DO
2453 :
2454 50 : CALL timestop(handle)
2455 64 : END SUBROUTINE rotate_rho_phase
2456 :
2457 : ! **************************************************************************************************
2458 : !> \brief Build a lab-frame copy of the (possibly rotating-frame) density rho for I/O.
2459 : !> On the TDA + symmetric-shift path returns rho_lab(t) by multiplying OV/VO by
2460 : !> exp(-/+ i Ω_0 t). Otherwise returns a plain copy. Writes into rho_new_last
2461 : !> (mo_struct-sized, idle outside ETRS) and returns a pointer to it; falls back to
2462 : !> the input rho when no scratch is available.
2463 : !> \param rtbse_env RT-BSE environment
2464 : !> \param rho_in Rotating-frame density (per spin)
2465 : !> \param t_phys Physical time associated with rho_in
2466 : !> \param rho_lab On exit, points to a per-spin cfm array holding rho in the lab frame.
2467 : ! **************************************************************************************************
2468 922 : SUBROUTINE build_rho_lab(rtbse_env, rho_in, t_phys, rho_lab)
2469 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2470 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_in
2471 : REAL(kind=dp), INTENT(IN) :: t_phys
2472 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_lab
2473 :
2474 : CHARACTER(len=*), PARAMETER :: routineN = 'build_rho_lab'
2475 :
2476 : INTEGER :: handle, i
2477 :
2478 874 : IF (rtbse_env%omega_shift == 0.0_dp .OR. .NOT. ASSOCIATED(rtbse_env%rho_new_last)) THEN
2479 826 : rho_lab => rho_in
2480 826 : RETURN
2481 : END IF
2482 :
2483 48 : CALL timeset(routineN, handle)
2484 :
2485 96 : DO i = 1, rtbse_env%n_spin
2486 48 : CALL cp_cfm_to_cfm(rho_in(i), rtbse_env%rho_new_last(i))
2487 : ! Lab-frame rho_OV(t) = exp(+i*Ω_0*t) * rho_tilde_OV(t)
2488 : ! (sign derived from [h_shifted, rho_tilde]_OV = (-eps_ai + Ω_0) rho_tilde_OV).
2489 96 : CALL rotate_rho_phase(rtbse_env, rtbse_env%rho_new_last(i), i, rtbse_env%omega_shift*t_phys)
2490 : END DO
2491 48 : rho_lab => rtbse_env%rho_new_last
2492 :
2493 48 : CALL timestop(handle)
2494 : END SUBROUTINE build_rho_lab
2495 :
2496 : ! **************************************************************************************************
2497 : !> \brief Bridge the restart density from the previous run's active-MO gauge into this run's:
2498 : !> overlap-metric basis change U_mn = sum_µν C2_µm S_µν C1_νn (mo_active × mo_active),
2499 : !> then ρ_mn ← sum_pq U_mp ρ_pq U_nq (ρ ← U ρ U^T, U real orthogonal up to FP).
2500 : !> Exact under per-MO sign flips and degenerate-subspace rotations of the SCF solution.
2501 : !> Diagnostics per spin: max|U−1| (total gauge correction), sign-flip count, max off-diag
2502 : !> (degenerate rotation), max|U^T U−1| (representability loss; warn ≥1e-10, abort ≥1e-3),
2503 : !> max|U_OV| (occ/virt mixing; warn ≥1e-6). No-op (U=1) when the two gauges agree.
2504 : !> \param rtbse_env RT-BSE environment
2505 : ! **************************************************************************************************
2506 14 : SUBROUTINE apply_restart_basis_bridge(rtbse_env)
2507 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2508 :
2509 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_restart_basis_bridge'
2510 : COMPLEX(kind=dp), PARAMETER :: c_one = CMPLX(1.0_dp, 0.0_dp, kind=dp), &
2511 : c_zero = CMPLX(0.0_dp, 0.0_dp, kind=dp)
2512 :
2513 : INTEGER :: handle, i, i_glob, i_mo, ii, j_glob, &
2514 : j_mo, jj, n_flip, ncol_local, &
2515 : nrow_local
2516 14 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2517 : LOGICAL :: found
2518 : REAL(kind=dp) :: dev_ident, dev_offdiag, dev_ov, &
2519 : dev_unitary
2520 14 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: u_diag
2521 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
2522 14 : POINTER :: u_data
2523 : TYPE(cp_fm_type) :: SC_old
2524 14 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: C_old
2525 :
2526 14 : CALL timeset(routineN, handle)
2527 :
2528 14 : NULLIFY (C_old)
2529 58 : ALLOCATE (C_old(rtbse_env%n_spin))
2530 30 : DO i = 1, rtbse_env%n_spin
2531 30 : CALL cp_fm_create(C_old(i), rtbse_env%fm_struct_ao_mo_active)
2532 : END DO
2533 14 : CALL read_restart_C(rtbse_env, C_old, found)
2534 14 : IF (.NOT. found) THEN
2535 0 : DO i = 1, rtbse_env%n_spin
2536 0 : CALL cp_fm_release(C_old(i))
2537 : END DO
2538 0 : DEALLOCATE (C_old)
2539 0 : CALL timestop(handle)
2540 : RETURN
2541 : END IF
2542 :
2543 14 : CALL cp_fm_create(SC_old, rtbse_env%fm_struct_ao_mo_active)
2544 42 : ALLOCATE (u_diag(rtbse_env%mo_active))
2545 :
2546 30 : DO i = 1, rtbse_env%n_spin
2547 : ! S C1 : [S C1]_µn = sum_ν S_µν C1_νn
2548 : CALL parallel_gemm("N", "N", rtbse_env%n_ao, rtbse_env%mo_active, rtbse_env%n_ao, &
2549 16 : 1.0_dp, rtbse_env%S_fm, C_old(i), 0.0_dp, SC_old)
2550 : ! U = C2^T (S C1) : U_mn = sum_µ C2_µm [S C1]_µn -> real_workspace_mo(1)
2551 : CALL parallel_gemm("T", "N", rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%n_ao, &
2552 16 : 1.0_dp, rtbse_env%C_active(i), SC_old, 0.0_dp, rtbse_env%real_workspace_mo(1))
2553 :
2554 : ! Diagnostics on U (local blocks + global MAX reduction; diagonal is gathered globally)
2555 16 : CALL cp_fm_get_diag(rtbse_env%real_workspace_mo(1), u_diag)
2556 208 : n_flip = COUNT(u_diag < 0.0_dp)
2557 : CALL cp_fm_get_info(rtbse_env%real_workspace_mo(1), nrow_local=nrow_local, ncol_local=ncol_local, &
2558 16 : row_indices=row_indices, col_indices=col_indices, local_data=u_data)
2559 16 : dev_ident = 0.0_dp; dev_offdiag = 0.0_dp; dev_ov = 0.0_dp
2560 112 : DO ii = 1, nrow_local
2561 96 : i_glob = row_indices(ii)
2562 96 : i_mo = i_glob + rtbse_env%first_active_mo - 1
2563 1264 : DO jj = 1, ncol_local
2564 1152 : j_glob = col_indices(jj)
2565 1152 : j_mo = j_glob + rtbse_env%first_active_mo - 1
2566 1152 : IF (i_glob == j_glob) THEN
2567 96 : dev_ident = MAX(dev_ident, ABS(u_data(ii, jj) - 1.0_dp))
2568 : ELSE
2569 1056 : dev_ident = MAX(dev_ident, ABS(u_data(ii, jj)))
2570 1056 : dev_offdiag = MAX(dev_offdiag, ABS(u_data(ii, jj)))
2571 : END IF
2572 1248 : IF ((i_mo <= rtbse_env%n_occ(i)) .NEQV. (j_mo <= rtbse_env%n_occ(i))) THEN
2573 176 : dev_ov = MAX(dev_ov, ABS(u_data(ii, jj)))
2574 : END IF
2575 : END DO
2576 : END DO
2577 : ! U^T U − 1 : [U^T U]_mn = sum_p U_pm U_pn -> real_workspace_mo(2)
2578 : CALL parallel_gemm("T", "N", rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
2579 : 1.0_dp, rtbse_env%real_workspace_mo(1), rtbse_env%real_workspace_mo(1), &
2580 16 : 0.0_dp, rtbse_env%real_workspace_mo(2))
2581 : CALL cp_fm_get_info(rtbse_env%real_workspace_mo(2), nrow_local=nrow_local, ncol_local=ncol_local, &
2582 16 : row_indices=row_indices, col_indices=col_indices, local_data=u_data)
2583 16 : dev_unitary = 0.0_dp
2584 112 : DO ii = 1, nrow_local
2585 96 : i_glob = row_indices(ii)
2586 1264 : DO jj = 1, ncol_local
2587 1152 : j_glob = col_indices(jj)
2588 2304 : dev_unitary = MAX(dev_unitary, ABS(u_data(ii, jj) - MERGE(1.0_dp, 0.0_dp, i_glob == j_glob)))
2589 : END DO
2590 : END DO
2591 16 : CALL rtbse_env%real_workspace_mo(1)%matrix_struct%para_env%max(dev_ident)
2592 16 : CALL rtbse_env%real_workspace_mo(1)%matrix_struct%para_env%max(dev_offdiag)
2593 16 : CALL rtbse_env%real_workspace_mo(1)%matrix_struct%para_env%max(dev_ov)
2594 16 : CALL rtbse_env%real_workspace_mo(1)%matrix_struct%para_env%max(dev_unitary)
2595 :
2596 16 : IF (rtbse_env%unit_nr > 0) THEN
2597 8 : WRITE (rtbse_env%unit_nr, '(A,I3,A)') " RTBSE| Restart basis bridge U = C2^T S C1 (spin ", i, "):"
2598 8 : WRITE (rtbse_env%unit_nr, '(A,ES12.3,A)') " RTBSE| max |U - 1| ", dev_ident, &
2599 16 : " (total gauge correction)"
2600 8 : WRITE (rtbse_env%unit_nr, '(A,I12)') " RTBSE| sign flips (U_ii<0)", n_flip
2601 8 : WRITE (rtbse_env%unit_nr, '(A,ES12.3,A)') " RTBSE| max offdiag |U_ij| ", dev_offdiag, &
2602 16 : " (degenerate-subspace rotation)"
2603 8 : WRITE (rtbse_env%unit_nr, '(A,ES12.3,A)') " RTBSE| max |U^T U - 1| ", dev_unitary, &
2604 16 : " (representability loss)"
2605 8 : WRITE (rtbse_env%unit_nr, '(A,ES12.3,A)') " RTBSE| max |U_OV| ", dev_ov, &
2606 16 : " (occ/virt structure change)"
2607 : END IF
2608 16 : IF (dev_unitary >= 1.0e-3_dp) THEN
2609 : CALL cp_abort(__LOCATION__, &
2610 0 : "Restart basis bridge: active spaces of the two runs differ severely (|U^T U - 1| >= 1e-3)")
2611 : END IF
2612 16 : IF (dev_unitary >= 1.0e-10_dp .AND. dev_unitary < 1.0e-3_dp) THEN
2613 : CALL cp_warn(__LOCATION__, &
2614 0 : "Restart basis bridge: representability loss above 1e-10 - active windows differ slightly.")
2615 : END IF
2616 : ! 1e-6 floor: benign SCF reconvergence gives ~1e-9 occ/virt gauge noise (the bridge maps it
2617 : ! correctly either way); only a genuine occupation-structure change reaches this threshold.
2618 16 : IF (dev_ov >= 1.0e-6_dp) THEN
2619 : CALL cp_warn(__LOCATION__, &
2620 0 : "Restart basis bridge: occupied/virtual mixing above 1e-6 - occupation structure changed.")
2621 : END IF
2622 :
2623 : ! ρ ← U ρ U^T : lift U to complex, [Uρ]_mn = sum_p U_mp ρ_pn, then ρ_mn = sum_q [Uρ]_mq U_nq
2624 16 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace_mo(1), mtarget=rtbse_env%rho_workspace(1))
2625 : CALL cp_cfm_gemm('N', 'N', rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
2626 16 : c_one, rtbse_env%rho_workspace(1), rtbse_env%rho(i), c_zero, rtbse_env%rho_workspace(2))
2627 : CALL cp_cfm_gemm('N', 'C', rtbse_env%mo_active, rtbse_env%mo_active, rtbse_env%mo_active, &
2628 62 : c_one, rtbse_env%rho_workspace(2), rtbse_env%rho_workspace(1), c_zero, rtbse_env%rho(i))
2629 : END DO
2630 :
2631 14 : DEALLOCATE (u_diag)
2632 14 : CALL cp_fm_release(SC_old)
2633 30 : DO i = 1, rtbse_env%n_spin
2634 30 : CALL cp_fm_release(C_old(i))
2635 : END DO
2636 14 : DEALLOCATE (C_old)
2637 14 : CALL timestop(handle)
2638 42 : END SUBROUTINE apply_restart_basis_bridge
2639 :
2640 : ! **************************************************************************************************
2641 : !> \brief Complex-linear Hartree contraction. Calls the real-input get_hartree on
2642 : !> Re(rho_AO) and on Im(rho_AO) separately and assembles
2643 : !> v_AO = V_H[Re(rho_AO)] + i * V_H[Im(rho_AO)] .
2644 : !> Required by the TDA propagator where the per-pass input Delta rho_OV (or
2645 : !> Delta rho_VO) is non-Hermitian, so the imaginary part must be carried.
2646 : !> The real kernel get_hartree realises V^H_λσ = sum_PQ (λσ|P) V_PQ [sum_µν (µν|Q) Δρ_µν].
2647 : !> \param rtbse_env RT-BSE environment
2648 : !> \param rho_cfm AO complex input density
2649 : !> \param v_cfm AO complex Hartree output (overwritten)
2650 : !> \param ispin Spin index (selects scratch slots in rtbse_env)
2651 : ! **************************************************************************************************
2652 1304 : SUBROUTINE get_hartree_complex(rtbse_env, rho_cfm, v_cfm, ispin)
2653 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2654 : TYPE(cp_cfm_type), INTENT(IN) :: rho_cfm
2655 : TYPE(cp_cfm_type) :: v_cfm
2656 : INTEGER, INTENT(IN) :: ispin
2657 :
2658 : CHARACTER(len=*), PARAMETER :: routineN = 'get_hartree_complex'
2659 :
2660 : INTEGER :: handle
2661 :
2662 : MARK_USED(ispin)
2663 :
2664 1304 : CALL timeset(routineN, handle)
2665 :
2666 : ! Mirrors get_sigma_complex: split rho_cfm into real and imaginary fm parts,
2667 : ! call the real-input Hartree contraction on each, then assemble
2668 : ! v_cfm = V_H[Re(rho_cfm)] + i*V_H[Im(rho_cfm)].
2669 : ! Scratch usage:
2670 : ! real_workspace(1) - holds Re(rho) then V_H[Re]
2671 : ! real_workspace(2) - holds Im(rho) then V_H[Im]
2672 : ! sigma_complex_workspace(1) - cfm wrapper feeding the real-input slot of get_hartree
2673 :
2674 : ! V^H[Re(Δρ^AO)] -> real_workspace(1)
2675 1304 : CALL cp_cfm_to_fm(msource=rho_cfm, mtargetr=rtbse_env%real_workspace(1))
2676 1304 : CALL cp_cfm_set_all(rtbse_env%sigma_complex_workspace(1), CMPLX(0.0_dp, 0.0_dp, kind=dp))
2677 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(1), &
2678 1304 : mtarget=rtbse_env%sigma_complex_workspace(1))
2679 : CALL get_hartree(rtbse_env, rtbse_env%sigma_complex_workspace(1), &
2680 1304 : rtbse_env%real_workspace(1))
2681 :
2682 : ! Imaginary part: extract Im(rho_cfm) into real_workspace(2)
2683 1304 : CALL cp_cfm_to_fm(msource=rho_cfm, mtargeti=rtbse_env%real_workspace(2))
2684 1304 : CALL cp_cfm_set_all(rtbse_env%sigma_complex_workspace(1), CMPLX(0.0_dp, 0.0_dp, kind=dp))
2685 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(2), &
2686 1304 : mtarget=rtbse_env%sigma_complex_workspace(1))
2687 : CALL get_hartree(rtbse_env, rtbse_env%sigma_complex_workspace(1), &
2688 1304 : rtbse_env%real_workspace(2))
2689 :
2690 : ! Assemble v_cfm = real_workspace(1) + i * real_workspace(2)
2691 : CALL cp_fm_to_cfm(msourcer=rtbse_env%real_workspace(1), &
2692 : msourcei=rtbse_env%real_workspace(2), &
2693 1304 : mtarget=v_cfm)
2694 :
2695 1304 : CALL timestop(handle)
2696 1304 : END SUBROUTINE get_hartree_complex
2697 :
2698 : ! **************************************************************************************************
2699 : !> \brief δ-kick (Marek2025) seeding the linearized EOM: builds the MO-active dipole operator
2700 : !> A = intensity * sum_k kvec_k r_k (MO basis), Hermitizes it, and propagates ρ by exp(-iA),
2701 : !> so Δρ^+_nm = i(f_n - f_m) A_nm excites only the OV/VO blocks.
2702 : !> \param rtbse_env RT-BSE environment
2703 : !> \author Stepan Marek (09.24)
2704 : !> \author Maximilian Graml - trafo to MO and linearized version following 10.1021/acs.jctc.2c00644 (03.26)
2705 : ! **************************************************************************************************
2706 44 : SUBROUTINE apply_delta_pulse_MO(rtbse_env)
2707 : TYPE(rtbse_env_type), POINTER :: rtbse_env
2708 :
2709 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_MO'
2710 :
2711 : INTEGER :: handle, i, k
2712 : REAL(kind=dp) :: intensity, metric
2713 : REAL(kind=dp), DIMENSION(3) :: kvec
2714 :
2715 44 : CALL timeset(routineN, handle)
2716 :
2717 : ! Report application
2718 44 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A28)') ' RTBSE| Applying delta pulse'
2719 : ! Extra minus for the propagation of density
2720 44 : intensity = -rtbse_env%dft_control%rtp_control%delta_pulse_scale
2721 : metric = 0.0_dp
2722 176 : kvec(:) = rtbse_env%dft_control%rtp_control%delta_pulse_direction(:)
2723 44 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A38,E14.4E3,E14.4E3,E14.4E3)') &
2724 88 : " RTBSE| Delta pulse elements (a.u.) : ", intensity*kvec(:)
2725 : ! Per-spin kick: each spin uses its own MO-active dipole operator (C_active(i_spin) basis)
2726 94 : DO i = 1, rtbse_env%n_spin
2727 50 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(1), 0.0_dp)
2728 200 : DO k = 1, 3
2729 : CALL cp_fm_scale_and_add(1.0_dp, rtbse_env%real_workspace_mo(1), &
2730 200 : kvec(k), rtbse_env%moments_field(k, i))
2731 : END DO
2732 : ! enforce hermiticity of the effective Hamiltonian
2733 50 : CALL cp_fm_transpose(rtbse_env%real_workspace_mo(1), rtbse_env%real_workspace_mo(2))
2734 : CALL cp_fm_scale_and_add(0.5_dp, rtbse_env%real_workspace_mo(1), &
2735 50 : 0.5_dp, rtbse_env%real_workspace_mo(2))
2736 : ! multiply by intensity, set as the imaginary exponent for this spin
2737 50 : CALL cp_fm_scale(intensity, rtbse_env%real_workspace_mo(1))
2738 94 : CALL cp_fm_to_cfm(msourcei=rtbse_env%real_workspace_mo(1), mtarget=rtbse_env%ham_workspace(i))
2739 : END DO
2740 : ! Propagate the density by the effect of the delta pulse
2741 44 : CALL propagate_density(rtbse_env, rtbse_env%ham_workspace, rtbse_env%rho, rtbse_env%rho_new)
2742 44 : metric = rho_metric(rtbse_env%rho_new, rtbse_env%rho, rtbse_env%n_spin)
2743 44 : IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, ('(A42,E38.8E3)')) " RTBSE| Metric difference after delta kick", metric
2744 : ! Copy the new density to the old density
2745 94 : DO i = 1, rtbse_env%n_spin
2746 94 : CALL cp_cfm_to_cfm(rtbse_env%rho_new(i), rtbse_env%rho(i))
2747 : END DO
2748 :
2749 44 : CALL timestop(handle)
2750 44 : END SUBROUTINE apply_delta_pulse_MO
2751 :
2752 : ! **************************************************************************************************
2753 : !> \brief Zero the OO and VV blocks of an MO-basis derivative cfm in place. Used to enforce
2754 : !> strict linear response on the RK4 derivatives so that the linearized propagator
2755 : !> only carries the OV/VO branches and the OO/VV orbital-energy spreads do not enter
2756 : !> the RK4 stability bound.
2757 : !> \param rtbse_env Entry point - rtbse environment
2758 : !> \param cfm Derivative-like cfm in MO basis for one spin channel (modified in place)
2759 : !> \param i_spin Spin index
2760 : ! **************************************************************************************************
2761 2680 : SUBROUTINE project_drho_to_ov(rtbse_env, cfm, i_spin)
2762 : TYPE(rtbse_env_type), INTENT(IN) :: rtbse_env
2763 : TYPE(cp_cfm_type), INTENT(INOUT) :: cfm
2764 : INTEGER, INTENT(IN) :: i_spin
2765 :
2766 2680 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: local_data
2767 : INTEGER :: i_global, i_global_mo, i_local, &
2768 : j_global, j_global_mo, j_local, n_occ, &
2769 : ncol_local, nrow_local, shift
2770 2680 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2771 :
2772 2680 : n_occ = rtbse_env%n_occ(i_spin)
2773 2680 : shift = rtbse_env%first_active_mo - 1
2774 :
2775 : CALL cp_cfm_get_info(matrix=cfm, &
2776 : nrow_local=nrow_local, ncol_local=ncol_local, &
2777 2680 : row_indices=row_indices, col_indices=col_indices)
2778 2680 : local_data => cfm%local_data
2779 :
2780 : ! keep OV/VO, zero OO/VV: project Δρ onto the δ-kick sectors
2781 34840 : DO j_local = 1, ncol_local
2782 32160 : j_global = col_indices(j_local)
2783 32160 : j_global_mo = j_global + shift
2784 227800 : DO i_local = 1, nrow_local
2785 192960 : i_global = row_indices(i_local)
2786 192960 : i_global_mo = i_global + shift
2787 192960 : IF ((i_global_mo <= n_occ .AND. j_global_mo <= n_occ) .OR. &
2788 32160 : (i_global_mo > n_occ .AND. j_global_mo > n_occ)) THEN
2789 163480 : local_data(i_local, j_local) = CMPLX(0.0_dp, 0.0_dp, kind=dp)
2790 : END IF
2791 : END DO
2792 : END DO
2793 2680 : END SUBROUTINE project_drho_to_ov
2794 :
2795 : ! **************************************************************************************************
2796 : !> \brief Per-spin electron numbers from the MO density: N_e^σ = spin_degeneracy * Re Tr[ρ^σ],
2797 : !> returned as one entry per spin channel (alpha/beta). The node-local diagonal partial
2798 : !> sums are reduced over the BLACS grid before scaling; the imaginary trace is a
2799 : !> non-Hermiticity diagnostic.
2800 : !> \param rtbse_env Entry point - rtbse environment
2801 : !> \param rho Density matrix in MO basis (per spin)
2802 : !> \param electron_n_re Real electron number per spin channel (size n_spin)
2803 : !> \param electron_n_im Imaginary electron number per spin channel (numerical non-hermiticity)
2804 : ! **************************************************************************************************
2805 830 : SUBROUTINE get_electron_number_MO(rtbse_env, rho, electron_n_re, electron_n_im)
2806 : TYPE(rtbse_env_type) :: rtbse_env
2807 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho
2808 : REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: electron_n_re, electron_n_im
2809 :
2810 : CHARACTER(len=*), PARAMETER :: routineN = 'get_electron_number_MO'
2811 :
2812 830 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: local_data
2813 : INTEGER :: handle, i_global, i_local, j, j_global, &
2814 : j_local, ncol_local, nrow_local
2815 830 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
2816 :
2817 830 : CALL timeset(routineN, handle)
2818 1780 : electron_n_re(:) = 0.0_dp
2819 1780 : electron_n_im(:) = 0.0_dp
2820 1780 : DO j = 1, rtbse_env%n_spin
2821 : CALL cp_cfm_get_info(matrix=rho(j), &
2822 : nrow_local=nrow_local, &
2823 : ncol_local=ncol_local, &
2824 : row_indices=row_indices, &
2825 950 : col_indices=col_indices)
2826 950 : local_data => rho(j)%local_data
2827 : ! accumulate Tr[ρ^σ] = sum_m ρ^σ_mm (real and imaginary parts separately)
2828 6530 : DO i_local = 1, nrow_local
2829 5580 : i_global = row_indices(i_local)
2830 : ! Search column indices for the diagonal position
2831 36860 : DO j_local = 1, ncol_local
2832 35910 : j_global = col_indices(j_local)
2833 35910 : IF (j_global == i_global) THEN
2834 : ! Found diagonal element
2835 5580 : electron_n_re(j) = electron_n_re(j) + REAL(local_data(i_local, j_local), kind=dp)
2836 5580 : electron_n_im(j) = electron_n_im(j) + AIMAG(local_data(i_local, j_local))
2837 5580 : EXIT
2838 : END IF
2839 : END DO
2840 : END DO
2841 : ! reduce the per-rank partial traces over the process grid (the MO diagonal is distributed)
2842 950 : CALL rho(j)%matrix_struct%para_env%sum(electron_n_re(j))
2843 950 : CALL rho(j)%matrix_struct%para_env%sum(electron_n_im(j))
2844 : ! N_e^σ = spin_degeneracy * Tr[ρ^σ] (g=2 closed shell; 1 per channel open shell)
2845 950 : electron_n_re(j) = electron_n_re(j)*rtbse_env%spin_degeneracy
2846 2730 : electron_n_im(j) = electron_n_im(j)*rtbse_env%spin_degeneracy
2847 : END DO
2848 :
2849 830 : CALL timestop(handle)
2850 830 : END SUBROUTINE get_electron_number_MO
2851 :
2852 : END MODULE rt_bse_linearized
|
