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