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