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