Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Data storage and other types for propagation via RT-BSE method.
10 : !> \author Stepan Marek (01.24)
11 : ! **************************************************************************************************
12 :
13 : MODULE rt_bse_types
14 :
15 : USE kinds, ONLY: dp
16 : USE cp_fm_types, ONLY: cp_fm_type, &
17 : cp_fm_release, &
18 : cp_fm_create, &
19 : cp_fm_set_all
20 : USE cp_fm_struct, ONLY: cp_fm_struct_type, &
21 : cp_fm_struct_create, &
22 : cp_fm_struct_release
23 : USE cp_cfm_types, ONLY: cp_cfm_type, &
24 : cp_cfm_set_all, &
25 : cp_cfm_create, &
26 : cp_fm_to_cfm, &
27 : cp_cfm_to_fm, &
28 : cp_cfm_release
29 : USE cp_dbcsr_api, ONLY: dbcsr_type, &
30 : dbcsr_p_type, &
31 : dbcsr_create, &
32 : dbcsr_release, &
33 : dbcsr_get_info
34 : USE parallel_gemm_api, ONLY: parallel_gemm
35 : USE dbt_api, ONLY: dbt_type, &
36 : dbt_create, &
37 : dbt_destroy
38 : USE qs_mo_types, ONLY: mo_set_type
39 : USE basis_set_types, ONLY: gto_basis_set_p_type
40 : USE cp_control_types, ONLY: dft_control_type
41 : USE qs_environment_types, ONLY: qs_environment_type, &
42 : get_qs_env
43 : USE force_env_types, ONLY: force_env_type
44 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type, &
45 : eps_qp_gap, &
46 : max_qp_gap
47 : USE rt_propagation_types, ONLY: rt_prop_type
48 : USE rt_propagation_utils, ONLY: warn_section_unused
49 : USE gw_integrals, ONLY: build_3c_integral_block
50 : USE gw_large_cell_Gamma, ONLY: compute_3c_integrals
51 : USE gw_utils, ONLY: rtbse_resolve_rirs_flag
52 : USE qs_tensors, ONLY: neighbor_list_3c_destroy
53 : USE libint_2c_3c, ONLY: libint_potential_type
54 : USE input_constants, ONLY: use_mom_ref_coac, &
55 : do_bch, &
56 : do_exact, &
57 : rtp_bse_ham_g0w0, &
58 : rtp_method_bse_linearized
59 : USE bse_util, ONLY: determine_cutoff_indices
60 : USE cp_log_handling, ONLY: cp_to_string
61 : USE physcon, ONLY: evolt
62 : USE mathconstants, ONLY: z_zero
63 : USE input_section_types, ONLY: section_vals_type, &
64 : section_vals_val_get, &
65 : section_vals_get_subs_vals
66 :
67 : #include "../base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = "rt_bse"
74 :
75 : #:include "rt_bse_macros.fypp"
76 :
77 : PUBLIC :: rtbse_env_type, &
78 : create_rtbse_env, &
79 : release_rtbse_env, &
80 : multiply_cfm_fm, &
81 : multiply_fm_cfm, &
82 : create_hartree_ri_3c, &
83 : create_sigma_workspace_qs_only
84 :
85 : ! ! Created so that we can have an array of pointers to arrays
86 : ! TYPE series_real_type
87 : ! REAL(kind=dp), DIMENSION(:), POINTER :: series => NULL()
88 : ! END TYPE series_real_type
89 : ! TYPE series_complex_type
90 : ! COMPLEX(kind=dp), DIMENSION(:), POINTER :: series => NULL()
91 : ! END TYPE series_complex_type
92 :
93 : ! **************************************************************************************************
94 : !> \param n_spin Number of spin channels that are present
95 : !> \param n_ao Number of atomic orbitals
96 : !> \param n_RI Number of RI orbitals
97 : !> \param n_occ Number of occupied orbitals, spin dependent
98 : !> \param spin_degeneracy Number of electrons per orbital
99 : !> \param field Electric field calculated at the given timestep
100 : !> \param moments Moment operators (2nd index = spin) along cartesian directions - centered at zero charge - used for plotting
101 : !> \param moments_field Moment operators (2nd index = spin) along cartesian directions - used to coupling to the field -
102 : !> origin bound to unit cell
103 : !> \param sim_step Current step of the simulation
104 : !> \param sim_start Starting step of the simulation
105 : !> \param sim_nsteps Number of steps of the simulation
106 : !> \param sim_time Current time of the simulation
107 : !> \param sim_dt Timestep of the simulation
108 : !> \param sim_dt_restart Original-run timestep read from the trace on restart (< 0 when not a restart)
109 : !> \param etrs_threshold Self-consistency threshold for enforced time reversal symmetry propagation
110 : !> \param exp_accuracy Threshold for matrix exponential calculation
111 : !> \param dft_control DFT control parameters
112 : !> \param ham_effective Real and imaginary part of the effective Hamiltonian used to propagate
113 : !> the density matrix
114 : !> \param ham_reference Reference Hamiltonian, which does not change in the
115 : !> propagation = DFT+G0W0 - initial Hartree - initial COHSEX
116 : !> \param ham_workspace Workspace matrices for use with the Hamiltonian propagation - storage of
117 : !> exponential propagators etc.
118 : !> \param rho Density matrix at the current time step
119 : !> \param rho_new Density matrix - workspace in ETRS
120 : !> \param rho_last Density matrix - workspace in ETRS
121 : !> \param rho_new_last Density matrix - workspace in ETRS
122 : !> \param rho_M Density matrix - workspace in ETRS
123 : !> \param S_inv_fm Inverse overlap matrix, full matrix
124 : !> \param S_fm Overlap matrix, full matrix
125 : !> \param S_inv Inverse overlap matrix, sparse matrix
126 : !> \param rho_dbcsr Density matrix, sparse matrix
127 : !> \param rho_workspace Matrices for storage of density matrix at different timesteps for
128 : !> interpolation and self-consistency checks etc.
129 : !> \param complex_workspace Workspace for complex density (exact diagonalisation)
130 : !> \param complex_s Complex overlap matrix (exact diagonalisation)
131 : !> \param real_eigvals Eigenvalues of hermitian matrix (exact diagonalisation)
132 : !> \param exp_eigvals Exponentiated eigenvalues (exact diagonalisation)
133 : !> \param v_dbcsr Sparse matrix with bare Coulomb in RI basis
134 : !> \param w_dbcsr Sparse matrix with correlation part of dressed Coulomb in RI basis (without bare Coulomb)
135 : !> \param screened_dbt Tensor for screened Coulomb interaction
136 : !> \param greens_dbt Tensor for greens function/density matrix
137 : !> \param t_3c_w Tensor containing 3c integrals
138 : !> \param t_3c_work_RI_AO__AO Tensor sigma contraction
139 : !> \param t_3c_work2_RI_AO__AO Tensor sigma contraction
140 : !> \param sigma_SEX Screened exchange self-energy
141 : !> \param sigma_COH Coulomb hole self-energy
142 : !> \param hartree_curr Current Hartree matrix
143 : !> \param etrs_max_iter Maximum number of ETRS iterations
144 : !> \param ham_reference_type Which Hamiltonian to use as single particle basis
145 : !> \param mat_exp_method Which method to use for matrix exponentiation
146 : !> \param unit_nr Number of output unit
147 : !> \param int_3c_array Array containing the local 3c integrals
148 : !> \author Stepan Marek (01.24)
149 : ! **************************************************************************************************
150 : TYPE rtbse_env_type
151 : INTEGER :: n_spin = 1, &
152 : n_ao = -1, &
153 : n_RI = -1
154 : INTEGER, DIMENSION(2) :: n_occ = -1
155 : ! Active MO window for linearized RT-BSE truncation. When no truncation is requested,
156 : ! first_active_mo=1, last_active_mo=n_ao, and mo_active=n_ao. The window is the
157 : ! combined inclusive bound that covers both spin channels.
158 : INTEGER :: first_active_mo = 1, &
159 : last_active_mo = -1, &
160 : mo_active = -1
161 : REAL(KIND=dp) :: rtbse_energy_cutoff_occ = -1.0_dp, &
162 : rtbse_energy_cutoff_empty = -1.0_dp
163 : LOGICAL :: active_mo_truncation = .FALSE.
164 : LOGICAL :: linearized = .FALSE.
165 : ! Tamm-Dancoff approximation switch (linearized RT-BSE only).
166 : LOGICAL :: tda_active = .FALSE.
167 : ! First-peak shift for the TDA path (linearized RT-BSE only).
168 : ! Shifts active-MO single-particle diagonals by +Omega_0/2 (occ) / -Omega_0/2 (virt)
169 : ! with Omega_0 = eps_min_ai so the lowest active OV mode oscillates at zero in the
170 : ! rotating frame (RK4-exact for peak 1). omega_max becomes the full active OV width
171 : ! Delta = eps_max_ai - eps_min_ai. The resulting rotating-frame density is undone
172 : ! at I/O so observables stay lab-frame.
173 : LOGICAL :: tda_shift_to_first_peak = .FALSE.
174 : REAL(kind=dp) :: omega_shift = 0.0_dp
175 : ! Debug-only kernel switches shared by initialization and propagation.
176 : LOGICAL :: debug_disable_hartree = .FALSE., &
177 : debug_disable_sex = .FALSE.
178 : ! RI framework for the linRTBSE Hartree + screened-exchange kernels, set by the KERNEL_RI
179 : ! input keyword (DEFAULT inherits bs_env%do_gw_ri_rs; RS/AO force; full RT-BSE forced AO).
180 : ! .TRUE. = RI-RS grid kernels, .FALSE. = AO-RI. The required grid and V_grid/W0_grid
181 : ! kernels are built on demand and reused across steps.
182 : LOGICAL :: rirs_kernel = .FALSE.
183 : ! Liouvillian eigenvalue diagnostic (TDA + n_spin=1 only). When .TRUE., at job
184 : ! init the linearized RT-BSE assembles the OV-subspace Liouvillian by probing
185 : ! apply_liouvillian_to_drho with canonical OV basis vectors and diagonalizes
186 : ! via cp_cfm_heevd. In TDA this equals the Casida-A eigenvalue problem.
187 : ! Run once, no propagation impact.
188 : LOGICAL :: diagnose_liouvillian_eig = .FALSE.
189 : ! Whether to enforce max_dt within stability region of rk4
190 : LOGICAL :: enforce_max_dt = .FALSE.
191 : ! Owned fm structures sized to the active MO window. Equal to the full n_ao x n_ao
192 : ! when no truncation is active.
193 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_mo_active => NULL()
194 : ! n_ao x mo_active fm struct used for C_active and AO<->MO rectangular intermediates.
195 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_ao_mo_active => NULL()
196 : ! Liouvillian-diagnostic struct: (N_OV_joint x N_OV_joint), spin blocks stacked, on the
197 : ! same BLACS context as fm_struct_mo_active. Allocated only when diagnose_liouvillian_eig=.TRUE..
198 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_ov_pairs => NULL()
199 : ! Truncated MO coefficient slabs C_active(:,:) of size n_ao x mo_active for each spin
200 : ! (only allocated for the linearized path).
201 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: C_active => NULL()
202 : ! Rectangular n_ao x mo_active scratch used by linearized AO<->MO transforms.
203 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: ao_mo_workspace => NULL()
204 : REAL(kind=dp) :: spin_degeneracy = 2
205 : REAL(kind=dp), DIMENSION(3) :: field = 0.0_dp
206 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: moments => NULL(), &
207 : moments_field => NULL()
208 : INTEGER :: sim_step = 0, &
209 : sim_start = 0, &
210 : ! Needed to continue runs by loading previous moments trace
211 : sim_start_orig = 0, &
212 : sim_nsteps = -1, &
213 : ! Default reference point type for output moments
214 : ! Field moments always use zero reference
215 : moment_ref_type = use_mom_ref_coac
216 : ! Restart output bookkeeping: RESTART.trace header+prefix and (linearized) C_active are
217 : ! (re)written on the first output_restart call of a run, then .trace is appended per step.
218 : LOGICAL :: restart_trace_written = .FALSE., &
219 : restart_C_written = .FALSE.
220 : REAL(kind=dp), DIMENSION(:), POINTER :: user_moment_ref_point => NULL()
221 : REAL(kind=dp) :: sim_time = 0.0_dp, &
222 : sim_dt = 0.1_dp, &
223 : ! Original-run dt from the trace header on restart (< 0 when not
224 : ! a restart); ENFORCE_MAX_DT reuses it instead of recomputing dt
225 : sim_dt_restart = -1.0_dp, &
226 : maximum_timestep = -1.0_dp, &
227 : omega_max = -1.0_dp, &
228 : etrs_threshold = 1.0e-7_dp, &
229 : exp_accuracy = 1.0e-10_dp, &
230 : ft_damping = 0.0_dp, &
231 : ft_start = 0.0_dp
232 : ! Which element of polarizability to print out
233 : INTEGER, DIMENSION(:, :), POINTER :: pol_elements => NULL()
234 : TYPE(dft_control_type), POINTER :: dft_control => NULL()
235 : ! DEBUG : Trying keeping the reference to previous environments inside this one
236 : TYPE(qs_environment_type), POINTER :: qs_env => NULL()
237 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env => NULL()
238 : ! Stores data needed for reading/writing to the restart files
239 : TYPE(section_vals_type), POINTER :: restart_section => NULL(), &
240 : field_section => NULL(), &
241 : rho_section => NULL(), &
242 : ft_section => NULL(), &
243 : pol_section => NULL(), &
244 : eig_section => NULL(), &
245 : moments_section => NULL(), &
246 : rtp_section => NULL()
247 : LOGICAL :: restart_extracted = .FALSE.
248 :
249 : ! Different indices signify different spins
250 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_effective => NULL()
251 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_reference => NULL()
252 : !Only for linearised RTBSE
253 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_reference_singleparticle => NULL()
254 : ! Active single-particle energies = diag(ham_reference_singleparticle), replicated over all
255 : ! ranks. Drives the [H^0,rho]_mn = (eps_m - eps_n) rho_mn element-wise commutator (mo_active x n_spin).
256 : REAL(kind=dp), DIMENSION(:, :), POINTER :: eps_active => NULL()
257 : ! Original run's active eigenvalues stashed from RESTART.trace at read time; compared against
258 : ! the recomputed eps_active once the Hamiltonian is built (consistency heads-up), then freed.
259 : REAL(kind=dp), DIMENSION(:, :), POINTER :: eps_active_restart => NULL()
260 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: ham_workspace => NULL()
261 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: sigma_SEX => NULL()
262 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: sigma_COH => NULL(), &
263 : hartree_curr => NULL()
264 : ! AO-sized scratch buffers used in the linearized RT-BSE path so that the MO-sized
265 : ! sigma_COH/sigma_SEX/hartree_curr matrices above can be allocated on
266 : ! fm_struct_mo_active. Only allocated when linearized=.TRUE..
267 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: sigma_SEX_ao => NULL()
268 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: hartree_curr_ao => NULL()
269 : ! Open-shell cross-spin Hartree (TDA AO-RI), single shared AO buffers built once per
270 : ! RK4 stage: rho_total_ao_scratch (fm_s struct, matches rho_ao_scratch for the spin sum);
271 : ! hartree_total_ao = V_H[sum] (fm_ks struct, matches sigma_SEX_ao as the Hartree output).
272 : TYPE(cp_cfm_type) :: rho_total_ao_scratch = cp_cfm_type(), &
273 : hartree_total_ao = cp_cfm_type()
274 : ! RI-RS Hartree diagonal reuse: spin-summed grid density n_l = diag(phi.rho.phi^T), harvested
275 : ! from the SEX rho_grid (before its Hadamard, scaled by spin_degeneracy) and consumed by
276 : ! compute_hartree_ri_rs_from_diag so Hartree never rebuilds the grid product. Re/Im; sized
277 : ! n_grid; allocated in initialize_hartree_potential when rirs_kernel. rtbse_env-owned (not
278 : ! bs_env%ri_rs) to avoid aliasing bs_env when passed through get_sigma_complex.
279 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: hartree_diag_re, hartree_diag_im
280 : ! mo_active x mo_active real workspace pair used in MO-side transforms (linearized only).
281 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: real_workspace_mo => NULL()
282 : ! Liouvillian-diagnostic scratch (allocated only when diagnose_liouvillian_eig=.TRUE.).
283 : ! drho_probe(:) / L_drho(:) are per-spin arrays on fm_struct_mo_active: the joint TDA
284 : ! diagnostic probes one spin and reads the Liouvillian response on every spin block.
285 : ! L_pairs / eigvecs_pairs live on fm_struct_ov_pairs (N_OV_joint x N_OV_joint, the spin
286 : ! blocks stacked). A_mat / B_mat / AmB_scratch / ApB_scratch are ABBA-only blocks holding
287 : ! A, B, (A-B) -> (A-B)^{1/2}, and (A+B) for the Furche reduction; allocated only
288 : ! when .NOT. tda_active (n_spin=1; TDA path uses L_pairs alone).
289 : ! eigenvalues_liouvillian holds the N_OV_joint real eigenvalues from cp_cfm_heevd.
290 : TYPE(cp_cfm_type) :: L_pairs = cp_cfm_type(), &
291 : eigvecs_pairs = cp_cfm_type(), &
292 : A_mat = cp_cfm_type(), &
293 : B_mat = cp_cfm_type(), &
294 : AmB_scratch = cp_cfm_type(), &
295 : ApB_scratch = cp_cfm_type()
296 : REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues_liouvillian => NULL()
297 :
298 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho => NULL(), &
299 : rho_new => NULL(), &
300 : rho_new_last => NULL(), &
301 : rho_M => NULL(), &
302 : rho_orig => NULL(), &
303 : rho_ao_scratch => NULL(), &
304 : rho_delta_mo => NULL(), &
305 : drho_probe => NULL(), &
306 : L_drho => NULL()
307 : ! Workspace for rk4 in linearized RTBSE
308 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rk4_coefficients => NULL()
309 :
310 : TYPE(cp_fm_type) :: S_inv_fm = cp_fm_type(), &
311 : S_fm = cp_fm_type()
312 : ! Many routines require overlap in the complex format
313 : TYPE(cp_cfm_type) :: S_cfm = cp_cfm_type()
314 : TYPE(dbcsr_type) :: rho_dbcsr = dbcsr_type(), &
315 : v_ao_dbcsr = dbcsr_type()
316 : ! Indices only correspond to different workspaces
317 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: rho_workspace => NULL()
318 : ! Many methods use real and imaginary parts separately - prevent unnecessary reallocation
319 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: real_workspace => NULL()
320 : ! AO-sized complex scratch used to stage the real part in get_sigma_complex.
321 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: sigma_complex_workspace => NULL()
322 : ! Workspace required for exact matrix exponentiation
323 : REAL(kind=dp), DIMENSION(:), POINTER :: real_eigvals => NULL()
324 : COMPLEX(kind=dp), DIMENSION(:), POINTER :: exp_eigvals => NULL()
325 : ! Workspace for saving the values for FT
326 : ! TODO : Change back to multi-dimensional arrays
327 : ! Index 1 : spin, Index 2 : direction, Index 3 : time point
328 : COMPLEX(kind=dp), DIMENSION(:, :, :), POINTER :: moments_trace => NULL()
329 : REAL(kind=dp), DIMENSION(:), POINTER :: time_trace => NULL()
330 : ! Index 1 : direction, Index 2 : time point
331 : COMPLEX(kind=dp), DIMENSION(:, :), POINTER :: field_trace => NULL()
332 : ! Workspace required for hartree_pw
333 : TYPE(dbcsr_type) :: v_dbcsr = dbcsr_type(), &
334 : w_dbcsr = dbcsr_type()
335 : #if defined(FTN_NO_DEFAULT_INIT)
336 : TYPE(dbt_type) :: screened_dbt, &
337 : greens_dbt, &
338 : t_3c_w, &
339 : t_3c_work_RI_AO__AO, &
340 : t_3c_work2_RI_AO__AO
341 : #else
342 : TYPE(dbt_type) :: screened_dbt = dbt_type(), &
343 : greens_dbt = dbt_type(), &
344 : t_3c_w = dbt_type(), &
345 : t_3c_work_RI_AO__AO = dbt_type(), &
346 : t_3c_work2_RI_AO__AO = dbt_type()
347 : #endif
348 : ! These matrices are always real
349 : INTEGER :: etrs_max_iter = 10
350 : INTEGER :: ham_reference_type = 2
351 : INTEGER :: mat_exp_method = 4
352 : INTEGER :: unit_nr = -1
353 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c_array => NULL()
354 : ! Parameters for Padé refinement
355 : REAL(kind=dp) :: pade_e_min = 0.0_dp, &
356 : pade_e_max = 100.0_dp, &
357 : pade_e_step = 0.05_dp, &
358 : pade_fit_e_min = 0.0_dp, &
359 : pade_fit_e_max = -1.0_dp
360 : INTEGER :: pade_npoints = 0
361 : LOGICAL :: pade_requested = .FALSE.
362 : COMPLEX(kind=dp), DIMENSION(:), POINTER :: pade_x_eval => NULL()
363 :
364 : END TYPE rtbse_env_type
365 :
366 : CONTAINS
367 :
368 : ! **************************************************************************************************
369 : !> \brief Allocates structures and prepares rtbse_env for run
370 : !> \param rtbse_env rtbse_env_type that is initialised
371 : !> \param force_env Force environment - entry point of the calculation
372 : !> \param linearized Optional; when present and .TRUE., configure the environment for the linearized RT-BSE path
373 : !> \author Stepan Marek
374 : !> \date 02.2024
375 : ! **************************************************************************************************
376 66 : SUBROUTINE create_rtbse_env(rtbse_env, force_env, linearized)
377 : TYPE(rtbse_env_type), POINTER :: rtbse_env
378 : TYPE(force_env_type), POINTER :: force_env
379 : LOGICAL, OPTIONAL :: linearized
380 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
381 : TYPE(rt_prop_type), POINTER :: rtp
382 66 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
383 66 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
384 : INTEGER :: i, k, n_ov, i_spin
385 : TYPE(section_vals_type), POINTER :: input, bs_sec, md_sec
386 : TYPE(cp_fm_struct_type), POINTER :: mo_struct
387 :
388 : ! Allocate the storage for the gwbse environment
389 66 : NULLIFY (rtbse_env, mo_struct)
390 3102 : ALLOCATE (rtbse_env)
391 66 : IF (PRESENT(linearized)) rtbse_env%linearized = linearized
392 : ! Extract the other types first
393 : CALL get_qs_env(force_env%qs_env, &
394 : bs_env=bs_env, &
395 : rtp=rtp, &
396 : matrix_s=matrix_s, &
397 : mos=mos, &
398 : dft_control=rtbse_env%dft_control, &
399 66 : input=input)
400 66 : bs_sec => section_vals_get_subs_vals(input, "PROPERTIES%BANDSTRUCTURE")
401 66 : IF (.NOT. ASSOCIATED(bs_env)) THEN
402 0 : CPABORT("Cannot run RT-BSE without running GW calculation (PROPERTIES) before")
403 : END IF
404 : ! Number of spins
405 66 : rtbse_env%n_spin = bs_env%n_spin
406 : ! Open shell (n_spin>1) is only implemented and tested for the linearized
407 : ! propagation; the full RT-BSE open-shell path is untested.
408 66 : IF (rtbse_env%n_spin > 1 .AND. .NOT. rtbse_env%linearized) THEN
409 : CALL cp_abort(__LOCATION__, &
410 : "Open-shell (n_spin>1) RT-BSE is only implemented and tested for the "// &
411 : "linearized propagation. Set DFT%REAL_TIME_PROPAGATION%RTBSE%LRRTBSE "// &
412 0 : ".TRUE.; the full (non-linearized) open-shell RT-BSE path is untested.")
413 : END IF
414 : ! Number of atomic orbitals
415 66 : rtbse_env%n_ao = bs_env%n_ao
416 : ! Number of auxiliary basis orbitals
417 66 : rtbse_env%n_RI = bs_env%n_RI
418 : ! Number of occupied orbitals - for closed shell equals to half the number of electrons
419 396 : rtbse_env%n_occ(:) = bs_env%n_occ(:)
420 : ! Spin degeneracy - number of spins per orbital
421 66 : rtbse_env%spin_degeneracy = bs_env%spin_degeneracy
422 : ! Default field is zero
423 264 : rtbse_env%field(:) = 0.0_dp
424 : ! Default time is zero
425 66 : rtbse_env%sim_step = 0
426 66 : rtbse_env%sim_time = 0
427 : ! Time step is taken from rtp
428 66 : md_sec => section_vals_get_subs_vals(force_env%root_section, "MOTION%MD")
429 66 : CALL section_vals_val_get(md_sec, "TIMESTEP", r_val=rtbse_env%sim_dt)
430 : ! rtbse_env%sim_dt = rtp%dt
431 : ! Threshold for etrs is taken from the eps_energy from RT propagation
432 66 : rtbse_env%etrs_threshold = rtbse_env%dft_control%rtp_control%eps_ener
433 66 : rtbse_env%exp_accuracy = rtbse_env%dft_control%rtp_control%eps_exp
434 : ! Recover custom options
435 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%RTBSE_HAMILTONIAN", &
436 66 : i_val=rtbse_env%ham_reference_type)
437 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAX_ITER", &
438 66 : i_val=rtbse_env%etrs_max_iter)
439 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAT_EXP", &
440 66 : i_val=rtbse_env%mat_exp_method)
441 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%ENERGY_CUTOFF_OCC", &
442 66 : r_val=rtbse_env%rtbse_energy_cutoff_occ)
443 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%ENERGY_CUTOFF_EMPTY", &
444 66 : r_val=rtbse_env%rtbse_energy_cutoff_empty)
445 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%TDA", &
446 66 : l_val=rtbse_env%tda_active)
447 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%TDA_SHIFT_TO_FIRST_PEAK", &
448 66 : l_val=rtbse_env%tda_shift_to_first_peak)
449 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%ENFORCE_MAX_DT", &
450 66 : l_val=rtbse_env%enforce_max_dt)
451 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%DEBUG_DISABLE_HARTREE", &
452 66 : l_val=rtbse_env%debug_disable_hartree)
453 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%DEBUG_DISABLE_SEX", &
454 66 : l_val=rtbse_env%debug_disable_sex)
455 : ! RI-RS kernel switch (KERNEL_RI: DEFAULT=do_gw_ri_rs, RS/AO force, full-RTBSE force-off)
456 : ! is resolved by the shared helper so de_init_bs_env reaches the same verdict when
457 : ! deciding whether to retain nl_3c.
458 66 : CALL rtbse_resolve_rirs_flag(force_env%qs_env, bs_env, rirs_kernel=rtbse_env%rirs_kernel)
459 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%DIAGNOSE_LIOUVILLIAN_EIG", &
460 66 : l_val=rtbse_env%diagnose_liouvillian_eig)
461 :
462 66 : IF (.NOT. rtbse_env%dft_control%rtp_control%rtp_method == rtp_method_bse_linearized) THEN
463 14 : rtbse_env%rtbse_energy_cutoff_occ = -1.0_dp
464 14 : rtbse_env%rtbse_energy_cutoff_empty = -1.0_dp
465 14 : rtbse_env%enforce_max_dt = .FALSE.
466 14 : rtbse_env%debug_disable_hartree = .FALSE.
467 14 : rtbse_env%debug_disable_sex = .FALSE.
468 14 : rtbse_env%tda_shift_to_first_peak = .FALSE.
469 : ! rirs_kernel is already forced .FALSE. here by rtbse_resolve_rirs_flag.
470 14 : rtbse_env%diagnose_liouvillian_eig = .FALSE.
471 : END IF
472 : ! First-peak shift only makes sense within TDA; force-disable otherwise.
473 66 : IF (.NOT. rtbse_env%tda_active) rtbse_env%tda_shift_to_first_peak = .FALSE.
474 66 : rtbse_env%omega_shift = 0.0_dp
475 :
476 66 : IF (rtbse_env%tda_active .AND. .NOT. rtbse_env%linearized) THEN
477 0 : CPABORT("RTBSE TDA keyword requires LINEARIZED_BSE_PROPAGATION=.TRUE.")
478 : END IF
479 : ! Open shell: omega_shift would be referenced to a non-physical cross-spin
480 : ! pseudo-gap (global MIN/MAX over both spins); abort until made per-spin.
481 66 : IF (rtbse_env%tda_shift_to_first_peak .AND. rtbse_env%n_spin > 1) THEN
482 : CALL cp_abort(__LOCATION__, &
483 : "TDA_SHIFT_TO_FIRST_PEAK is not implemented for open-shell (n_spin>1) "// &
484 : "systems - the first-peak gap estimate would mix spin channels. "// &
485 0 : "Set TDA_SHIFT_TO_FIRST_PEAK=.FALSE. for open-shell runs.")
486 : END IF
487 66 : CALL check_qp_gap_sanity(rtbse_env, bs_env)
488 66 : CALL determine_active_mo_window(rtbse_env, bs_env)
489 : ! Owned active-MO matrix structure (currently identical to the full n_ao x n_ao struct
490 : ! when no truncation is active; will be used by the linearized RT-BSE allocation path).
491 66 : NULLIFY (rtbse_env%fm_struct_mo_active)
492 : CALL cp_fm_struct_create(rtbse_env%fm_struct_mo_active, &
493 : bs_env%fm_ks_Gamma(1)%matrix_struct%para_env, &
494 : bs_env%fm_ks_Gamma(1)%matrix_struct%context, &
495 66 : rtbse_env%mo_active, rtbse_env%mo_active)
496 : ! Rectangular n_ao x mo_active struct used for C_active and AO<->MO intermediates
497 66 : NULLIFY (rtbse_env%fm_struct_ao_mo_active)
498 : CALL cp_fm_struct_create(rtbse_env%fm_struct_ao_mo_active, &
499 : bs_env%fm_ks_Gamma(1)%matrix_struct%para_env, &
500 : bs_env%fm_ks_Gamma(1)%matrix_struct%context, &
501 66 : rtbse_env%n_ao, rtbse_env%mo_active)
502 : ! Choose the matrix struct used for MO-side persistent matrices.
503 : ! Linearized RT-BSE: mo_active x mo_active. Full RT-BSE: full AO struct (unchanged).
504 66 : IF (rtbse_env%linearized) THEN
505 52 : mo_struct => rtbse_env%fm_struct_mo_active
506 : ELSE
507 14 : mo_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
508 : END IF
509 : ! Output unit number, recovered from the post_scf_bandstructure_type
510 66 : rtbse_env%unit_nr = bs_env%unit_nr
511 : ! Sim start index and total number of steps as well
512 66 : CALL section_vals_val_get(md_sec, "STEP_START_VAL", i_val=rtbse_env%sim_start)
513 : ! Copy this value to sim_start_orig for continuation runs
514 66 : rtbse_env%sim_start_orig = rtbse_env%sim_start
515 66 : CALL section_vals_val_get(md_sec, "STEPS", i_val=rtbse_env%sim_nsteps)
516 : ! Get the values for the FT
517 66 : rtbse_env%ft_damping = rtbse_env%dft_control%rtp_control%ft_damping
518 66 : rtbse_env%ft_damping = rtbse_env%dft_control%rtp_control%ft_t0
519 66 : rtbse_env%pol_elements => rtbse_env%dft_control%rtp_control%print_pol_elements
520 :
521 66 : rtbse_env%rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
522 : ! Get the restart section
523 66 : rtbse_env%restart_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%RESTART")
524 66 : rtbse_env%restart_extracted = .FALSE.
525 66 : rtbse_env%field_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%FIELD")
526 66 : rtbse_env%moments_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%MOMENTS")
527 : ! Moment specification
528 : CALL section_vals_val_get(rtbse_env%rtp_section, "PRINT%MOMENTS%REFERENCE", &
529 66 : i_val=rtbse_env%moment_ref_type)
530 : CALL section_vals_val_get(rtbse_env%rtp_section, "PRINT%MOMENTS%REFERENCE_POINT", &
531 66 : r_vals=rtbse_env%user_moment_ref_point)
532 66 : rtbse_env%rho_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%DENSITY_MATRIX")
533 66 : rtbse_env%ft_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%MOMENTS_FT")
534 66 : rtbse_env%pol_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%POLARIZABILITY")
535 66 : rtbse_env%eig_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%LIOUVILLIAN_EIG")
536 : ! Warn the user about print sections which are not yet implemented in the RTBSE run
537 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%CURRENT", &
538 66 : "CURRENT print section not yet implemented for RTBSE.")
539 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%E_CONSTITUENTS", &
540 66 : "E_CONSTITUENTS print section not yet implemented for RTBSE.")
541 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%PROGRAM_RUN_INFO", &
542 66 : "PROGRAM_RUN_INFO print section not yet implemented for RTBSE.")
543 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%PROJECTION_MO", &
544 66 : "PROJECTION_MO print section not yet implemented for RTBSE.")
545 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%RESTART_HISTORY", &
546 66 : "RESTART_HISTORY print section not yet implemented for RTBSE.")
547 : ! References to the parent qs_env / bs_env
548 66 : rtbse_env%qs_env => force_env%qs_env
549 66 : rtbse_env%bs_env => bs_env
550 : ! Padé refinement
551 66 : rtbse_env%pade_requested = rtbse_env%dft_control%rtp_control%pade_requested
552 66 : rtbse_env%pade_e_min = rtbse_env%dft_control%rtp_control%pade_e_min
553 66 : rtbse_env%pade_e_step = rtbse_env%dft_control%rtp_control%pade_e_step
554 66 : rtbse_env%pade_e_max = rtbse_env%dft_control%rtp_control%pade_e_max
555 66 : rtbse_env%pade_fit_e_min = rtbse_env%dft_control%rtp_control%pade_fit_e_min
556 66 : rtbse_env%pade_fit_e_max = rtbse_env%dft_control%rtp_control%pade_fit_e_max
557 66 : rtbse_env%pade_npoints = INT((rtbse_env%pade_e_max - rtbse_env%pade_e_min)/rtbse_env%pade_e_step)
558 : ! Evaluate the evaluation grid
559 66 : IF (rtbse_env%pade_requested) THEN
560 2 : NULLIFY (rtbse_env%pade_x_eval)
561 6 : ALLOCATE (rtbse_env%pade_x_eval(rtbse_env%pade_npoints))
562 2000 : DO i = 1, rtbse_env%pade_npoints
563 2000 : rtbse_env%pade_x_eval(i) = CMPLX(rtbse_env%pade_e_step*REAL(i - 1, kind=dp), 0.0, kind=dp)
564 : END DO
565 : END IF
566 :
567 : ! Allocate moments matrices.
568 : ! In linearized RT-BSE these store the MO-active transformed dipole moments;
569 : ! in full RT-BSE they remain AO-sized (initialized from overlap template).
570 66 : NULLIFY (rtbse_env%moments)
571 494 : ALLOCATE (rtbse_env%moments(3, rtbse_env%n_spin))
572 66 : NULLIFY (rtbse_env%moments_field)
573 494 : ALLOCATE (rtbse_env%moments_field(3, rtbse_env%n_spin))
574 140 : DO i_spin = 1, rtbse_env%n_spin
575 362 : DO k = 1, 3
576 222 : CALL cp_fm_create(rtbse_env%moments(k, i_spin), mo_struct)
577 296 : CALL cp_fm_create(rtbse_env%moments_field(k, i_spin), mo_struct)
578 : END DO
579 : END DO
580 :
581 : ! Allocate space for density propagation and other operations.
582 : ! In linearized RT-BSE these workspaces are MO-active sized; in full RT-BSE
583 : ! they remain at the full AO size.
584 66 : NULLIFY (rtbse_env%rho_workspace)
585 330 : ALLOCATE (rtbse_env%rho_workspace(4))
586 330 : DO i = 1, SIZE(rtbse_env%rho_workspace)
587 264 : CALL cp_cfm_create(rtbse_env%rho_workspace(i), mo_struct)
588 330 : CALL cp_cfm_set_all(rtbse_env%rho_workspace(i), CMPLX(0.0, 0.0, kind=dp))
589 : END DO
590 :
591 : ! TODO: gate workspace allocation so methods skip workspaces they don't need
592 :
593 : ! Allocate real workspace
594 66 : NULLIFY (rtbse_env%real_workspace)
595 66 : IF (rtbse_env%linearized) THEN
596 156 : ALLOCATE (rtbse_env%real_workspace(2))
597 : ELSE
598 14 : SELECT CASE (rtbse_env%mat_exp_method)
599 : CASE (do_exact)
600 0 : ALLOCATE (rtbse_env%real_workspace(4))
601 : CASE (do_bch)
602 42 : ALLOCATE (rtbse_env%real_workspace(2))
603 : CASE DEFAULT
604 14 : CPABORT("Only exact and BCH matrix propagation implemented in RT-BSE")
605 : END SELECT
606 : END IF
607 198 : DO i = 1, SIZE(rtbse_env%real_workspace)
608 132 : CALL cp_fm_create(rtbse_env%real_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
609 198 : CALL cp_fm_set_all(rtbse_env%real_workspace(i), 0.0_dp)
610 : END DO
611 66 : NULLIFY (rtbse_env%sigma_complex_workspace)
612 132 : ALLOCATE (rtbse_env%sigma_complex_workspace(1))
613 66 : CALL cp_cfm_create(rtbse_env%sigma_complex_workspace(1), bs_env%fm_ks_Gamma(1)%matrix_struct)
614 66 : CALL cp_cfm_set_all(rtbse_env%sigma_complex_workspace(1), CMPLX(0.0_dp, 0.0_dp, kind=dp))
615 : ! Allocate density matrix (MO-active sized when linearized; AO-sized otherwise)
616 66 : NULLIFY (rtbse_env%rho)
617 272 : ALLOCATE (rtbse_env%rho(rtbse_env%n_spin))
618 140 : DO i = 1, rtbse_env%n_spin
619 140 : CALL cp_cfm_create(rtbse_env%rho(i), matrix_struct=mo_struct)
620 : END DO
621 : ! Allocate additional space for AO density matrix
622 : ! in linearised RTBSE, where default is MO
623 66 : IF (rtbse_env%linearized) THEN
624 52 : NULLIFY (rtbse_env%rho_ao_scratch)
625 216 : ALLOCATE (rtbse_env%rho_ao_scratch(rtbse_env%n_spin))
626 112 : DO i = 1, rtbse_env%n_spin
627 112 : CALL cp_cfm_create(rtbse_env%rho_ao_scratch(i), matrix_struct=bs_env%fm_s_Gamma%matrix_struct)
628 : END DO
629 : END IF
630 : ! Create the inverse overlap matrix, for use in density propagation
631 : ! Start by creating the actual overlap matrix
632 66 : CALL cp_fm_create(rtbse_env%S_fm, bs_env%fm_s_Gamma%matrix_struct)
633 66 : CALL cp_fm_create(rtbse_env%S_inv_fm, bs_env%fm_s_Gamma%matrix_struct)
634 66 : CALL cp_cfm_create(rtbse_env%S_cfm, bs_env%fm_s_Gamma%matrix_struct)
635 :
636 : ! Create the single particle hamiltonian
637 : ! Allocate workspace (MO-active sized in linearized RT-BSE; AO sized otherwise)
638 66 : NULLIFY (rtbse_env%ham_workspace)
639 272 : ALLOCATE (rtbse_env%ham_workspace(rtbse_env%n_spin))
640 140 : DO i = 1, rtbse_env%n_spin
641 74 : CALL cp_cfm_create(rtbse_env%ham_workspace(i), mo_struct)
642 140 : CALL cp_cfm_set_all(rtbse_env%ham_workspace(i), CMPLX(0.0, 0.0, kind=dp))
643 : END DO
644 : ! Now onto the Hamiltonian itself
645 : ! full RTBSE: Contains energy differences and Hartree/COHSEX ρ_0 parts
646 : ! linearised RTBSE: Contains only the Hartree/SEX ρ_0 parts as Δε * Δρ(t) need to be updated
647 66 : NULLIFY (rtbse_env%ham_reference)
648 272 : ALLOCATE (rtbse_env%ham_reference(rtbse_env%n_spin))
649 140 : DO i = 1, rtbse_env%n_spin
650 140 : CALL cp_cfm_create(rtbse_env%ham_reference(i), mo_struct)
651 : END DO
652 : ! Single particle Hamiltonian (Δε * Δρ(t)) for updates during timesteps in LR-RTBSE
653 66 : IF (rtbse_env%linearized) THEN
654 52 : NULLIFY (rtbse_env%ham_reference_singleparticle)
655 216 : ALLOCATE (rtbse_env%ham_reference_singleparticle(rtbse_env%n_spin))
656 112 : DO i = 1, rtbse_env%n_spin
657 112 : CALL cp_cfm_create(rtbse_env%ham_reference_singleparticle(i), mo_struct)
658 : END DO
659 52 : NULLIFY (rtbse_env%eps_active)
660 208 : ALLOCATE (rtbse_env%eps_active(rtbse_env%mo_active, rtbse_env%n_spin))
661 820 : rtbse_env%eps_active(:, :) = 0.0_dp
662 : END IF
663 :
664 : ! Create the matrices and workspaces for ETRS propagation
665 66 : NULLIFY (rtbse_env%ham_effective)
666 66 : NULLIFY (rtbse_env%rho_new)
667 66 : NULLIFY (rtbse_env%rho_new_last)
668 66 : NULLIFY (rtbse_env%rho_M)
669 66 : NULLIFY (rtbse_env%rho_orig)
670 272 : ALLOCATE (rtbse_env%ham_effective(rtbse_env%n_spin))
671 272 : ALLOCATE (rtbse_env%rho_new(rtbse_env%n_spin))
672 272 : ALLOCATE (rtbse_env%rho_new_last(rtbse_env%n_spin))
673 272 : ALLOCATE (rtbse_env%rho_M(rtbse_env%n_spin))
674 272 : ALLOCATE (rtbse_env%rho_orig(rtbse_env%n_spin))
675 140 : DO i = 1, rtbse_env%n_spin
676 74 : CALL cp_cfm_create(rtbse_env%ham_effective(i), mo_struct)
677 74 : CALL cp_cfm_set_all(rtbse_env%ham_effective(i), CMPLX(0.0, 0.0, kind=dp))
678 74 : CALL cp_cfm_create(rtbse_env%rho_new(i), mo_struct)
679 74 : CALL cp_cfm_set_all(rtbse_env%rho_new(i), CMPLX(0.0, 0.0, kind=dp))
680 74 : CALL cp_cfm_create(rtbse_env%rho_new_last(i), mo_struct)
681 74 : CALL cp_cfm_set_all(rtbse_env%rho_new_last(i), CMPLX(0.0, 0.0, kind=dp))
682 74 : CALL cp_cfm_create(rtbse_env%rho_M(i), mo_struct)
683 74 : CALL cp_cfm_set_all(rtbse_env%rho_M(i), CMPLX(0.0, 0.0, kind=dp))
684 140 : CALL cp_cfm_create(rtbse_env%rho_orig(i), mo_struct)
685 : END DO
686 :
687 : !For LR-RTBSE we need RK4 coefficients - create new workspace
688 66 : IF (rtbse_env%linearized) THEN
689 : ! Indexed by SPIN, not by RK4 stage: the spin loop is inner to each stage (do_rk4_stage), so every
690 : ! spin's current-stage k must be live at once, but only one stage's k per spin - each is folded
691 : ! into rho_end and the next stage density before the next stage overwrites it. Hence size n_spin.
692 52 : NULLIFY (rtbse_env%rk4_coefficients)
693 216 : ALLOCATE (rtbse_env%rk4_coefficients(rtbse_env%n_spin))
694 112 : DO i = 1, rtbse_env%n_spin
695 60 : CALL cp_cfm_create(rtbse_env%rk4_coefficients(i), mo_struct)
696 112 : CALL cp_cfm_set_all(rtbse_env%rk4_coefficients(i), CMPLX(0.0, 0.0, kind=dp))
697 : END DO
698 : END IF
699 :
700 : ! Fields for exact diagonalisation
701 66 : NULLIFY (rtbse_env%real_eigvals)
702 198 : ALLOCATE (rtbse_env%real_eigvals(rtbse_env%n_ao))
703 718 : rtbse_env%real_eigvals(:) = 0.0_dp
704 66 : NULLIFY (rtbse_env%exp_eigvals)
705 198 : ALLOCATE (rtbse_env%exp_eigvals(rtbse_env%n_ao))
706 718 : rtbse_env%exp_eigvals(:) = CMPLX(0.0, 0.0, kind=dp)
707 :
708 : ! Workspace for FT - includes (in principle) the zeroth step and the extra last step
709 66 : NULLIFY (rtbse_env%moments_trace)
710 : ! TODO : Unite the number of steps with TD-DFT
711 18950 : ALLOCATE (rtbse_env%moments_trace(rtbse_env%n_spin, 3, rtbse_env%sim_nsteps + 1), source=z_zero)
712 66 : NULLIFY (rtbse_env%field_trace)
713 10622 : ALLOCATE (rtbse_env%field_trace(3, rtbse_env%sim_nsteps + 1), source=z_zero)
714 66 : NULLIFY (rtbse_env%time_trace)
715 2804 : ALLOCATE (rtbse_env%time_trace(rtbse_env%sim_nsteps + 1), source=0.0_dp)
716 :
717 : ! Allocate self-energy parts and dynamic Hartree potential.
718 : ! In linearized RT-BSE these matrices hold the MO-active-sized result of the
719 : ! AO->MO transform; the AO-sized buffer is allocated as sigma_*_ao below.
720 66 : NULLIFY (rtbse_env%hartree_curr)
721 66 : NULLIFY (rtbse_env%sigma_SEX)
722 66 : NULLIFY (rtbse_env%sigma_COH)
723 272 : ALLOCATE (rtbse_env%hartree_curr(rtbse_env%n_spin))
724 272 : ALLOCATE (rtbse_env%sigma_SEX(rtbse_env%n_spin))
725 272 : ALLOCATE (rtbse_env%sigma_COH(rtbse_env%n_spin))
726 140 : DO i = 1, rtbse_env%n_spin
727 74 : CALL cp_fm_create(rtbse_env%sigma_COH(i), mo_struct)
728 74 : CALL cp_cfm_create(rtbse_env%sigma_SEX(i), mo_struct)
729 74 : CALL cp_fm_create(rtbse_env%hartree_curr(i), mo_struct)
730 74 : CALL cp_fm_set_all(rtbse_env%sigma_COH(i), 0.0_dp)
731 74 : CALL cp_cfm_set_all(rtbse_env%sigma_SEX(i), CMPLX(0.0, 0.0, kind=dp))
732 140 : CALL cp_fm_set_all(rtbse_env%hartree_curr(i), 0.0_dp)
733 : END DO
734 : ! AO-sized scratch buffers used by the linearized RT-BSE path
735 66 : IF (rtbse_env%linearized) THEN
736 52 : NULLIFY (rtbse_env%hartree_curr_ao)
737 52 : NULLIFY (rtbse_env%sigma_SEX_ao)
738 216 : ALLOCATE (rtbse_env%hartree_curr_ao(rtbse_env%n_spin))
739 216 : ALLOCATE (rtbse_env%sigma_SEX_ao(rtbse_env%n_spin))
740 112 : DO i = 1, rtbse_env%n_spin
741 60 : CALL cp_cfm_create(rtbse_env%sigma_SEX_ao(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
742 60 : CALL cp_fm_create(rtbse_env%hartree_curr_ao(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
743 60 : CALL cp_cfm_set_all(rtbse_env%sigma_SEX_ao(i), CMPLX(0.0, 0.0, kind=dp))
744 112 : CALL cp_fm_set_all(rtbse_env%hartree_curr_ao(i), 0.0_dp)
745 : END DO
746 : ! mo_active x mo_active real workspace pair for MO-side intermediates
747 52 : NULLIFY (rtbse_env%real_workspace_mo)
748 156 : ALLOCATE (rtbse_env%real_workspace_mo(2))
749 156 : DO i = 1, SIZE(rtbse_env%real_workspace_mo)
750 104 : CALL cp_fm_create(rtbse_env%real_workspace_mo(i), rtbse_env%fm_struct_mo_active)
751 156 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(i), 0.0_dp)
752 : END DO
753 52 : NULLIFY (rtbse_env%ao_mo_workspace)
754 104 : ALLOCATE (rtbse_env%ao_mo_workspace(1))
755 52 : CALL cp_fm_create(rtbse_env%ao_mo_workspace(1), rtbse_env%fm_struct_ao_mo_active)
756 52 : CALL cp_fm_set_all(rtbse_env%ao_mo_workspace(1), 0.0_dp)
757 : ! Truncated MO coefficient slabs C_active (n_ao x mo_active) for each spin.
758 : ! Filled in initialize_rtbse_env from bs_env%fm_mo_coeff_Gamma via submatrix copy.
759 52 : NULLIFY (rtbse_env%C_active)
760 216 : ALLOCATE (rtbse_env%C_active(rtbse_env%n_spin))
761 112 : DO i = 1, rtbse_env%n_spin
762 60 : CALL cp_fm_create(rtbse_env%C_active(i), rtbse_env%fm_struct_ao_mo_active)
763 112 : CALL cp_fm_set_all(rtbse_env%C_active(i), 0.0_dp)
764 : END DO
765 : ! Masked-copy staging scratch for the builder (all propagation paths including closed-shell ABBA).
766 : ! Also used as conjugate-transpose scratch in the TDA consumer.
767 52 : NULLIFY (rtbse_env%rho_delta_mo)
768 216 : ALLOCATE (rtbse_env%rho_delta_mo(rtbse_env%n_spin))
769 112 : DO i = 1, rtbse_env%n_spin
770 60 : CALL cp_cfm_create(rtbse_env%rho_delta_mo(i), rtbse_env%fm_struct_mo_active)
771 112 : CALL cp_cfm_set_all(rtbse_env%rho_delta_mo(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
772 : END DO
773 : ! Shared AO Hartree buffers: widened from (tda_active .OR. n_spin>1) so n_spin=1 ABBA
774 : ! (diagnostic + propagator) gets a dedicated buffer instead of aliasing sigma_SEX_ao.
775 52 : IF (.NOT. rtbse_env%debug_disable_hartree) THEN
776 48 : CALL cp_cfm_create(rtbse_env%rho_total_ao_scratch, bs_env%fm_s_Gamma%matrix_struct)
777 48 : CALL cp_cfm_create(rtbse_env%hartree_total_ao, bs_env%fm_ks_Gamma(1)%matrix_struct)
778 48 : CALL cp_cfm_set_all(rtbse_env%rho_total_ao_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
779 48 : CALL cp_cfm_set_all(rtbse_env%hartree_total_ao, CMPLX(0.0_dp, 0.0_dp, kind=dp))
780 : END IF
781 : ! Liouvillian-eigenvalue diagnostic state. n_spin=1 enforced upstream.
782 : ! Shared scratch (TDA + ABBA) allocated on diagnose_liouvillian_eig=T;
783 : ! ABBA-only A/B/A±B blocks added below under .NOT. tda_active. Mirror of the
784 : ! existing shared-scratch pattern, with the extra tda_active gate as the deviation.
785 52 : IF (rtbse_env%diagnose_liouvillian_eig) THEN
786 : ! Joint OV dimension: the spin blocks are stacked (n_spin=1 -> the old single-spin
787 : ! size). drho_probe/L_drho stay mo_active-sized per spin; L_pairs is N_OV_joint.
788 52 : n_ov = 0
789 112 : DO i = 1, rtbse_env%n_spin
790 : n_ov = n_ov + (rtbse_env%n_occ(i) - rtbse_env%first_active_mo + 1)* &
791 112 : (rtbse_env%last_active_mo - rtbse_env%n_occ(i))
792 : END DO
793 52 : NULLIFY (rtbse_env%fm_struct_ov_pairs)
794 : CALL cp_fm_struct_create(rtbse_env%fm_struct_ov_pairs, &
795 : bs_env%fm_ks_Gamma(1)%matrix_struct%para_env, &
796 : bs_env%fm_ks_Gamma(1)%matrix_struct%context, &
797 52 : n_ov, n_ov)
798 216 : ALLOCATE (rtbse_env%drho_probe(rtbse_env%n_spin))
799 216 : ALLOCATE (rtbse_env%L_drho(rtbse_env%n_spin))
800 112 : DO i = 1, rtbse_env%n_spin
801 60 : CALL cp_cfm_create(rtbse_env%drho_probe(i), rtbse_env%fm_struct_mo_active)
802 60 : CALL cp_cfm_create(rtbse_env%L_drho(i), rtbse_env%fm_struct_mo_active)
803 60 : CALL cp_cfm_set_all(rtbse_env%drho_probe(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
804 112 : CALL cp_cfm_set_all(rtbse_env%L_drho(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
805 : END DO
806 52 : CALL cp_cfm_create(rtbse_env%L_pairs, rtbse_env%fm_struct_ov_pairs)
807 52 : CALL cp_cfm_create(rtbse_env%eigvecs_pairs, rtbse_env%fm_struct_ov_pairs)
808 52 : CALL cp_cfm_set_all(rtbse_env%L_pairs, CMPLX(0.0_dp, 0.0_dp, kind=dp))
809 52 : CALL cp_cfm_set_all(rtbse_env%eigvecs_pairs, CMPLX(0.0_dp, 0.0_dp, kind=dp))
810 52 : NULLIFY (rtbse_env%eigenvalues_liouvillian)
811 156 : ALLOCATE (rtbse_env%eigenvalues_liouvillian(n_ov))
812 700 : rtbse_env%eigenvalues_liouvillian = 0.0_dp
813 : ! ABBA-only Furche-reduction scratch (A, B, A-B->sqrt, A+B).
814 52 : IF (.NOT. rtbse_env%tda_active) THEN
815 16 : CALL cp_cfm_create(rtbse_env%A_mat, rtbse_env%fm_struct_ov_pairs)
816 16 : CALL cp_cfm_create(rtbse_env%B_mat, rtbse_env%fm_struct_ov_pairs)
817 16 : CALL cp_cfm_create(rtbse_env%AmB_scratch, rtbse_env%fm_struct_ov_pairs)
818 16 : CALL cp_cfm_create(rtbse_env%ApB_scratch, rtbse_env%fm_struct_ov_pairs)
819 16 : CALL cp_cfm_set_all(rtbse_env%A_mat, CMPLX(0.0_dp, 0.0_dp, kind=dp))
820 16 : CALL cp_cfm_set_all(rtbse_env%B_mat, CMPLX(0.0_dp, 0.0_dp, kind=dp))
821 16 : CALL cp_cfm_set_all(rtbse_env%AmB_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
822 16 : CALL cp_cfm_set_all(rtbse_env%ApB_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
823 : END IF
824 : END IF
825 : END IF
826 :
827 : ! Allocate workspaces for get_sigma
828 66 : CALL create_sigma_workspace(rtbse_env)
829 :
830 : ! Depending on the chosen methods, allocate extra workspace
831 66 : CALL create_hartree_ri_workspace(rtbse_env)
832 :
833 66 : END SUBROUTINE create_rtbse_env
834 :
835 : ! **************************************************************************************************
836 : !> \brief Simple reimplementation of cp_fm_release_pp1 for complex matrices
837 : !> \param matrices cp_cfm_type(:)
838 : !> \author Stepan Marek
839 : !> \date 02.2024
840 : ! **************************************************************************************************
841 1090 : SUBROUTINE cp_cfm_release_pa1(matrices)
842 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: matrices
843 : INTEGER :: i
844 :
845 2506 : DO i = 1, SIZE(matrices)
846 2506 : CALL cp_cfm_release(matrices(i))
847 : END DO
848 1090 : DEALLOCATE (matrices)
849 : NULLIFY (matrices)
850 1090 : END SUBROUTINE cp_cfm_release_pa1
851 :
852 : ! **************************************************************************************************
853 : !> \brief Releases the environment allocated structures
854 : !> \param rtbse_env
855 : !> \author Stepan Marek
856 : !> \date 02.2024
857 : ! **************************************************************************************************
858 66 : SUBROUTINE release_rtbse_env(rtbse_env)
859 : TYPE(rtbse_env_type), POINTER :: rtbse_env
860 :
861 66 : CALL cp_cfm_release_pa1(rtbse_env%ham_effective)
862 66 : CALL cp_cfm_release_pa1(rtbse_env%ham_workspace)
863 66 : CALL cp_fm_release(rtbse_env%sigma_COH)
864 66 : CALL cp_cfm_release_pa1(rtbse_env%sigma_SEX)
865 66 : CALL cp_fm_release(rtbse_env%hartree_curr)
866 66 : CALL cp_cfm_release_pa1(rtbse_env%ham_reference)
867 66 : IF (ASSOCIATED(rtbse_env%ham_reference_singleparticle)) THEN
868 52 : CALL cp_cfm_release_pa1(rtbse_env%ham_reference_singleparticle)
869 : END IF
870 66 : IF (ASSOCIATED(rtbse_env%eps_active)) DEALLOCATE (rtbse_env%eps_active)
871 66 : IF (ASSOCIATED(rtbse_env%eps_active_restart)) DEALLOCATE (rtbse_env%eps_active_restart)
872 66 : CALL cp_cfm_release_pa1(rtbse_env%rho)
873 66 : CALL cp_cfm_release_pa1(rtbse_env%rho_workspace)
874 66 : CALL cp_cfm_release_pa1(rtbse_env%rho_new)
875 66 : CALL cp_cfm_release_pa1(rtbse_env%rho_new_last)
876 66 : CALL cp_cfm_release_pa1(rtbse_env%rho_M)
877 66 : CALL cp_cfm_release_pa1(rtbse_env%rho_orig)
878 66 : IF (ASSOCIATED(rtbse_env%rk4_coefficients)) THEN
879 52 : CALL cp_cfm_release_pa1(rtbse_env%rk4_coefficients)
880 : END IF
881 66 : CALL cp_fm_release(rtbse_env%real_workspace)
882 66 : IF (ASSOCIATED(rtbse_env%sigma_complex_workspace)) CALL cp_cfm_release_pa1(rtbse_env%sigma_complex_workspace)
883 66 : CALL cp_fm_release(rtbse_env%S_inv_fm)
884 66 : CALL cp_fm_release(rtbse_env%S_fm)
885 66 : CALL cp_cfm_release(rtbse_env%S_cfm)
886 :
887 66 : CALL cp_fm_release(rtbse_env%moments)
888 66 : CALL cp_fm_release(rtbse_env%moments_field)
889 :
890 66 : CALL release_sigma_workspace(rtbse_env)
891 :
892 66 : CALL release_hartree_ri_workspace(rtbse_env)
893 :
894 66 : DEALLOCATE (rtbse_env%real_eigvals)
895 66 : DEALLOCATE (rtbse_env%exp_eigvals)
896 66 : DEALLOCATE (rtbse_env%moments_trace)
897 66 : DEALLOCATE (rtbse_env%field_trace)
898 66 : DEALLOCATE (rtbse_env%time_trace)
899 :
900 66 : IF (ASSOCIATED(rtbse_env%pol_elements)) DEALLOCATE (rtbse_env%pol_elements)
901 66 : IF (ASSOCIATED(rtbse_env%pade_x_eval)) DEALLOCATE (rtbse_env%pade_x_eval)
902 :
903 : ! Deallocate the neighbour list that is not deallocated in gw anymore
904 66 : IF (ASSOCIATED(rtbse_env%bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(rtbse_env%bs_env%nl_3c)
905 : ! Release linearized-only AO scratches and MO-side workspaces
906 66 : IF (ASSOCIATED(rtbse_env%rho_ao_scratch)) CALL cp_cfm_release_pa1(rtbse_env%rho_ao_scratch)
907 66 : IF (ASSOCIATED(rtbse_env%sigma_SEX_ao)) CALL cp_cfm_release_pa1(rtbse_env%sigma_SEX_ao)
908 66 : IF (ASSOCIATED(rtbse_env%hartree_curr_ao)) CALL cp_fm_release(rtbse_env%hartree_curr_ao)
909 66 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) CALL cp_fm_release(rtbse_env%real_workspace_mo)
910 66 : IF (ASSOCIATED(rtbse_env%ao_mo_workspace)) CALL cp_fm_release(rtbse_env%ao_mo_workspace)
911 66 : IF (ASSOCIATED(rtbse_env%C_active)) CALL cp_fm_release(rtbse_env%C_active)
912 66 : IF (ASSOCIATED(rtbse_env%rho_delta_mo)) CALL cp_cfm_release_pa1(rtbse_env%rho_delta_mo)
913 : ! Release shared bare-Hartree scratch. Mirror the alloc gate exactly (linearized .AND.
914 : ! .NOT. debug_disable_hartree, every shell incl closed-shell ABBA) — the old
915 : ! (tda_active .OR. n_spin>1) gate leaked both buffers on the closed-shell ABBA path.
916 66 : IF (rtbse_env%linearized .AND. .NOT. rtbse_env%debug_disable_hartree) THEN
917 48 : CALL cp_cfm_release(rtbse_env%rho_total_ao_scratch)
918 48 : CALL cp_cfm_release(rtbse_env%hartree_total_ao)
919 : END IF
920 : ! Release the RI-RS Hartree diagonal-reuse accumulators (allocated in initialize_hartree_potential).
921 66 : IF (ALLOCATED(rtbse_env%hartree_diag_re)) DEALLOCATE (rtbse_env%hartree_diag_re)
922 66 : IF (ALLOCATED(rtbse_env%hartree_diag_im)) DEALLOCATE (rtbse_env%hartree_diag_im)
923 : ! Release Liouvillian-diagnostic scratch (only when the diagnostic was requested).
924 66 : IF (rtbse_env%diagnose_liouvillian_eig) THEN
925 52 : IF (ASSOCIATED(rtbse_env%drho_probe)) CALL cp_cfm_release_pa1(rtbse_env%drho_probe)
926 52 : IF (ASSOCIATED(rtbse_env%L_drho)) CALL cp_cfm_release_pa1(rtbse_env%L_drho)
927 52 : CALL cp_cfm_release(rtbse_env%L_pairs)
928 52 : CALL cp_cfm_release(rtbse_env%eigvecs_pairs)
929 52 : IF (ASSOCIATED(rtbse_env%eigenvalues_liouvillian)) DEALLOCATE (rtbse_env%eigenvalues_liouvillian)
930 52 : IF (.NOT. rtbse_env%tda_active) THEN
931 16 : CALL cp_cfm_release(rtbse_env%A_mat)
932 16 : CALL cp_cfm_release(rtbse_env%B_mat)
933 16 : CALL cp_cfm_release(rtbse_env%AmB_scratch)
934 16 : CALL cp_cfm_release(rtbse_env%ApB_scratch)
935 : END IF
936 52 : IF (ASSOCIATED(rtbse_env%fm_struct_ov_pairs)) THEN
937 52 : CALL cp_fm_struct_release(rtbse_env%fm_struct_ov_pairs)
938 : END IF
939 : END IF
940 : ! Release owned active-MO matrix structures
941 66 : IF (ASSOCIATED(rtbse_env%fm_struct_mo_active)) THEN
942 66 : CALL cp_fm_struct_release(rtbse_env%fm_struct_mo_active)
943 : END IF
944 66 : IF (ASSOCIATED(rtbse_env%fm_struct_ao_mo_active)) THEN
945 66 : CALL cp_fm_struct_release(rtbse_env%fm_struct_ao_mo_active)
946 : END IF
947 : ! Deallocate the storage for the environment itself
948 66 : DEALLOCATE (rtbse_env)
949 : ! Nullify to make sure it is not used again
950 : NULLIFY (rtbse_env)
951 :
952 66 : END SUBROUTINE release_rtbse_env
953 :
954 : ! **************************************************************************************************
955 : !> \brief Abort if the quasiparticle spectrum handed to the propagator is inverted or has diverged.
956 : !>
957 : !> Tests the fundamental gap per spin channel - not E(HOMO+1) - E(HOMO), since G0W0 reorders levels -
958 : !> on the very array the propagator consumes. Under RTBSE_HAMILTONIAN KS the quasiparticle energies
959 : !> never enter the propagator, so a broken G0W0 spectrum is irrelevant there and does not abort.
960 : !> \param rtbse_env RT-BSE environment with n_ao, n_occ, n_spin, ham_reference_type populated.
961 : !> \param bs_env Bandstructure environment providing the eigenvalues.
962 : ! **************************************************************************************************
963 66 : SUBROUTINE check_qp_gap_sanity(rtbse_env, bs_env)
964 : TYPE(rtbse_env_type), POINTER :: rtbse_env
965 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
966 :
967 : INTEGER :: homo, ispin
968 : REAL(KIND=dp) :: gap, gap_scf
969 :
970 66 : IF (rtbse_env%ham_reference_type /= rtp_bse_ham_g0w0) RETURN
971 :
972 140 : DO ispin = 1, rtbse_env%n_spin
973 74 : homo = rtbse_env%n_occ(ispin)
974 74 : IF (homo < 1 .OR. homo >= rtbse_env%n_ao) CYCLE
975 :
976 : gap = MINVAL(bs_env%eigenval_G0W0(homo + 1:rtbse_env%n_ao, 1, ispin)) - &
977 822 : MAXVAL(bs_env%eigenval_G0W0(1:homo, 1, ispin))
978 : gap_scf = MINVAL(bs_env%eigenval_scf_Gamma(homo + 1:rtbse_env%n_ao, ispin)) - &
979 822 : MAXVAL(bs_env%eigenval_scf_Gamma(1:homo, ispin))
980 :
981 : ! requiring a healthy SCF gap keeps the inversion test from firing on a genuine metal
982 140 : IF (gap < -eps_qp_gap .AND. gap_scf > eps_qp_gap) THEN
983 : CALL cp_abort(__LOCATION__, &
984 : "RTBSE: G0W0 gap of spin "//TRIM(ADJUSTL(cp_to_string(ispin)))// &
985 : " is negative ("//TRIM(ADJUSTL(cp_to_string(gap*evolt, '(F12.3)')))// &
986 : " eV): propagating an inverted spectrum is meaningless. Check the GW "// &
987 0 : "numerical parameters, or use RTBSE_HAMILTONIAN KS.")
988 74 : ELSE IF (ABS(gap) > max_qp_gap) THEN
989 : CALL cp_abort(__LOCATION__, &
990 : "RTBSE: G0W0 gap of spin "//TRIM(ADJUSTL(cp_to_string(ispin)))// &
991 : " is implausibly large ("// &
992 : TRIM(ADJUSTL(cp_to_string(gap*evolt, '(F12.3)')))//" eV): the GW step "// &
993 : "has likely diverged. Check the GW numerical parameters, or use "// &
994 0 : "RTBSE_HAMILTONIAN KS.")
995 : END IF
996 : END DO
997 :
998 : END SUBROUTINE check_qp_gap_sanity
999 :
1000 : ! **************************************************************************************************
1001 : !> \brief Determine the combined active MO window for linearized RT-BSE truncation.
1002 : !>
1003 : !> Evaluates BSE-like cutoff indices per spin from the requested single-particle spectrum
1004 : !> (G0W0 or KS Gamma-point eigenvalues) and collapses them into a single combined window
1005 : !> covering both spin channels by choosing the most inclusive bounds. Issues a CPWARN if the
1006 : !> spin-resolved cutoff candidates differ. When cutoffs are disabled (or the run is not
1007 : !> linearized RT-BSE), the window is set to the full MO range.
1008 : !> \param rtbse_env RT-BSE environment with cutoff values, n_ao, n_occ, n_spin, ham_reference_type
1009 : !> already populated.
1010 : !> \param bs_env Bandstructure environment providing the eigenvalues.
1011 : ! **************************************************************************************************
1012 66 : SUBROUTINE determine_active_mo_window(rtbse_env, bs_env)
1013 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1014 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1015 : CHARACTER(LEN=*), PARAMETER :: routineN = "determine_active_mo_window"
1016 :
1017 : INTEGER :: handle, ispin, n_ao_full, n_virt
1018 : INTEGER :: homo_red, virt_red, homo_incl, virt_incl
1019 : INTEGER :: combined_first_occ, combined_last_virt
1020 : INTEGER :: first_occ_prev, last_virt_prev
1021 : LOGICAL :: spins_differ, do_truncation
1022 : REAL(KIND=dp) :: cutoff_occ, cutoff_empty
1023 :
1024 66 : CALL timeset(routineN, handle)
1025 :
1026 66 : n_ao_full = rtbse_env%n_ao
1027 66 : cutoff_occ = rtbse_env%rtbse_energy_cutoff_occ
1028 66 : cutoff_empty = rtbse_env%rtbse_energy_cutoff_empty
1029 66 : do_truncation = rtbse_env%linearized .AND. (cutoff_occ > 0.0_dp .OR. cutoff_empty > 0.0_dp)
1030 :
1031 : ! Default: full MO window
1032 66 : rtbse_env%first_active_mo = 1
1033 66 : rtbse_env%last_active_mo = n_ao_full
1034 66 : rtbse_env%mo_active = n_ao_full
1035 66 : rtbse_env%active_mo_truncation = .FALSE.
1036 :
1037 66 : IF (.NOT. do_truncation) THEN
1038 64 : CALL timestop(handle)
1039 64 : RETURN
1040 : END IF
1041 :
1042 2 : combined_first_occ = n_ao_full
1043 2 : combined_last_virt = 1
1044 2 : first_occ_prev = -1
1045 2 : last_virt_prev = -1
1046 2 : spins_differ = .FALSE.
1047 :
1048 4 : DO ispin = 1, rtbse_env%n_spin
1049 2 : n_virt = n_ao_full - rtbse_env%n_occ(ispin)
1050 : ! Cut on the DFT axis, as LRBSE does: it is ascending by construction, so the window is a
1051 : ! well-defined contiguous MO range, which is all C_active can extract. The G0W0 axis is
1052 : ! not ordered.
1053 : CALL determine_cutoff_indices(bs_env%eigenval_scf_Gamma(:, ispin), &
1054 : rtbse_env%n_occ(ispin), n_virt, &
1055 : homo_red, virt_red, homo_incl, virt_incl, &
1056 2 : cutoff_occ, cutoff_empty)
1057 : ! Translate the per-spin candidate to global MO indices [homo_incl, homo + virt_incl]
1058 2 : IF (ispin > 1) THEN
1059 0 : IF (homo_incl /= first_occ_prev .OR. (rtbse_env%n_occ(ispin) + virt_incl) /= last_virt_prev) THEN
1060 0 : spins_differ = .TRUE.
1061 : END IF
1062 : END IF
1063 2 : first_occ_prev = homo_incl
1064 2 : last_virt_prev = rtbse_env%n_occ(ispin) + virt_incl
1065 2 : combined_first_occ = MIN(combined_first_occ, homo_incl)
1066 4 : combined_last_virt = MAX(combined_last_virt, rtbse_env%n_occ(ispin) + virt_incl)
1067 : END DO
1068 :
1069 2 : IF (spins_differ) THEN
1070 0 : CPWARN("RTBSE: spin-resolved active MO cutoff candidates differ; using combined window.")
1071 : END IF
1072 :
1073 2 : rtbse_env%first_active_mo = combined_first_occ
1074 2 : rtbse_env%last_active_mo = combined_last_virt
1075 2 : rtbse_env%mo_active = combined_last_virt - combined_first_occ + 1
1076 2 : rtbse_env%active_mo_truncation = (rtbse_env%mo_active < n_ao_full)
1077 :
1078 2 : CALL timestop(handle)
1079 : END SUBROUTINE determine_active_mo_window
1080 :
1081 : ! **************************************************************************************************
1082 : !> \brief Allocates the workspaces for Hartree RI method
1083 : !> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors
1084 : !> \param rtbse_env
1085 : !> \author Stepan Marek
1086 : !> \date 05.2024
1087 : ! **************************************************************************************************
1088 66 : SUBROUTINE create_hartree_ri_workspace(rtbse_env)
1089 : TYPE(rtbse_env_type) :: rtbse_env
1090 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1091 :
1092 : ! Skip the AO-RI Hartree scratch when the RT-BSE Hartree path is fully RI-RS.
1093 : ! In that case rho_dbcsr / v_ao_dbcsr / int_3c_array are never read.
1094 : ! get_sigma_real (AO-RI SX) used to borrow rho_dbcsr as a workspace; that
1095 : ! cross-dependency was removed by giving get_sigma_real its own local
1096 : ! dbcsr scratch (see rt_bse.F::get_sigma_real). rho_dbcsr is now AO-RI
1097 : ! Hartree only, as its name suggests.
1098 66 : IF (rtbse_env%rirs_kernel) RETURN
1099 :
1100 50 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
1101 :
1102 50 : CALL dbcsr_create(rtbse_env%rho_dbcsr, name="Sparse density", template=bs_env%mat_ao_ao%matrix)
1103 50 : CALL dbcsr_create(rtbse_env%v_ao_dbcsr, name="Sparse Hartree", template=bs_env%mat_ao_ao%matrix)
1104 :
1105 : CALL create_hartree_ri_3c(rtbse_env%rho_dbcsr, rtbse_env%int_3c_array, rtbse_env%n_ao, rtbse_env%n_RI, &
1106 : bs_env%basis_set_AO, bs_env%basis_set_RI, bs_env%i_RI_start_from_atom, &
1107 50 : bs_env%ri_metric, rtbse_env%qs_env, rtbse_env%unit_nr)
1108 : END SUBROUTINE create_hartree_ri_workspace
1109 : ! **************************************************************************************************
1110 : !> \brief Separated method for allocating the 3c integrals for RI Hartree
1111 : !> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors
1112 : !> \param rho_dbcsr matrix used for the description of shape of 3c array
1113 : !> \param int_3c 3-center integral array to be allocated and filled
1114 : !> \param n_ao Number of atomic orbitals
1115 : !> \param n_RI Number of auxiliary RI orbitals
1116 : !> \param basis_set_AO AO basis set
1117 : !> \param basis_set_RI RI auxiliary basis set
1118 : !> \param i_RI_start_from_atom Array of indices where functions of a given atom in RI basis start
1119 : !> \param unit_nr Unit number used for printing information about the size of int_3c
1120 : !> \author Stepan Marek
1121 : !> \date 01.2025
1122 : ! **************************************************************************************************
1123 50 : SUBROUTINE create_hartree_ri_3c(rho_dbcsr, int_3c, n_ao, n_RI, basis_set_AO, basis_set_RI, &
1124 50 : i_RI_start_from_atom, ri_metric, qs_env, unit_nr)
1125 : TYPE(dbcsr_type) :: rho_dbcsr
1126 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c
1127 : INTEGER :: n_ao, n_RI
1128 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_set_AO, &
1129 : basis_set_RI
1130 : INTEGER, DIMENSION(:) :: i_RI_start_from_atom
1131 : TYPE(libint_potential_type) :: ri_metric
1132 : TYPE(qs_environment_type), POINTER :: qs_env
1133 : INTEGER :: unit_nr
1134 : REAL(kind=dp) :: size_mb
1135 : INTEGER :: nblkrows_local, &
1136 : nblkcols_local, &
1137 : i_blk_local, &
1138 : j_blk_local, &
1139 : nrows_local, &
1140 : ncols_local, &
1141 : col_local_offset, &
1142 : row_local_offset, &
1143 : start_col_index, &
1144 : end_col_index, &
1145 : start_row_index, &
1146 : end_row_index
1147 50 : INTEGER, DIMENSION(:), POINTER :: local_blk_rows, &
1148 50 : local_blk_cols, &
1149 50 : row_blk_size, &
1150 50 : col_blk_size
1151 : ! TODO : Implement option/decision to not precompute all the 3c integrals
1152 : size_mb = REAL(n_ao, kind=dp)*REAL(n_ao, kind=dp)*REAL(n_RI, kind=dp)* &
1153 50 : REAL(STORAGE_SIZE(size_mb), kind=dp)/8.0_dp/1024.0_dp/1024.0_dp
1154 50 : IF (unit_nr > 0) WRITE (unit_nr, '(A44,E32.2E3,A4)') &
1155 25 : " RTBSE| Approximate size of the 3c integrals", size_mb, " MiB"
1156 :
1157 : ! Get the number of block rows and columns
1158 50 : CALL dbcsr_get_info(rho_dbcsr, nblkrows_local=nblkrows_local, nblkcols_local=nblkcols_local)
1159 : ! Get the global indices of local rows and columns
1160 50 : CALL dbcsr_get_info(rho_dbcsr, local_rows=local_blk_rows, local_cols=local_blk_cols)
1161 : ! Get the sizes of all blocks
1162 50 : CALL dbcsr_get_info(rho_dbcsr, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1163 :
1164 : ! Get the total required local rows and cols
1165 50 : nrows_local = 0
1166 100 : DO i_blk_local = 1, nblkrows_local
1167 100 : nrows_local = nrows_local + row_blk_size(local_blk_rows(i_blk_local))
1168 : END DO
1169 50 : ncols_local = 0
1170 150 : DO j_blk_local = 1, nblkcols_local
1171 150 : ncols_local = ncols_local + col_blk_size(local_blk_cols(j_blk_local))
1172 : END DO
1173 :
1174 : ! Allocate the appropriate storage
1175 250 : ALLOCATE (int_3c(nrows_local, ncols_local, n_RI))
1176 :
1177 : ! Fill the storage with appropriate values, block by block
1178 50 : row_local_offset = 1
1179 100 : DO i_blk_local = 1, nblkrows_local
1180 : col_local_offset = 1
1181 150 : DO j_blk_local = 1, nblkcols_local
1182 100 : start_row_index = row_local_offset
1183 100 : end_row_index = start_row_index + row_blk_size(local_blk_rows(i_blk_local)) - 1
1184 100 : start_col_index = col_local_offset
1185 100 : end_col_index = start_col_index + col_blk_size(local_blk_cols(j_blk_local)) - 1
1186 : CALL build_3c_integral_block(int_3c(start_row_index:end_row_index, &
1187 : start_col_index:end_col_index, &
1188 : 1:n_RI), &
1189 : qs_env, potential_parameter=ri_metric, &
1190 : basis_j=basis_set_AO, basis_k=basis_set_AO, &
1191 : basis_i=basis_set_RI, &
1192 : atom_j=local_blk_rows(i_blk_local), &
1193 : atom_k=local_blk_cols(j_blk_local), &
1194 100 : i_bf_start_from_atom=i_RI_start_from_atom)
1195 150 : col_local_offset = col_local_offset + col_blk_size(local_blk_cols(j_blk_local))
1196 : END DO
1197 100 : row_local_offset = row_local_offset + row_blk_size(local_blk_rows(i_blk_local))
1198 : END DO
1199 75 : END SUBROUTINE create_hartree_ri_3c
1200 : ! **************************************************************************************************
1201 : !> \brief Releases the workspace for the Hartree RI method
1202 : !> \param rtbse_env RT-BSE Environment, containing specific RI Hartree storage
1203 : !> \author Stepan Marek
1204 : !> \date 09.2024
1205 : ! **************************************************************************************************
1206 66 : SUBROUTINE release_hartree_ri_workspace(rtbse_env)
1207 : TYPE(rtbse_env_type) :: rtbse_env
1208 :
1209 : ! Mirror the gate in create_hartree_ri_workspace and the v_dbcsr gate in
1210 : ! initialize_hartree_potential. With one KERNEL_RI switch the AO-RI Hartree
1211 : ! scratch (3c integrals + dbcsr work + v_dbcsr) is created iff `.NOT. rirs_kernel`.
1212 66 : IF (.NOT. rtbse_env%rirs_kernel) THEN
1213 50 : DEALLOCATE (rtbse_env%int_3c_array)
1214 50 : CALL dbcsr_release(rtbse_env%rho_dbcsr)
1215 50 : CALL dbcsr_release(rtbse_env%v_ao_dbcsr)
1216 50 : CALL dbcsr_release(rtbse_env%v_dbcsr)
1217 : END IF
1218 66 : END SUBROUTINE release_hartree_ri_workspace
1219 : ! **************************************************************************************************
1220 : !> \brief Allocates the workspaces for self-energy determination routine
1221 : !> \param rtbse_env Structure for holding information and workspace structures
1222 : !> \author Stepan Marek
1223 : !> \date 02.2024
1224 : ! **************************************************************************************************
1225 66 : SUBROUTINE create_sigma_workspace(rtbse_env)
1226 : TYPE(rtbse_env_type) :: rtbse_env
1227 :
1228 : ! Skip the AO-RI sigma scratch (W matrix + 3c integrals + work tensors)
1229 : ! when the RT-BSE SEX path is fully RI-RS — these workspaces are read
1230 : ! only by get_sigma in the `.NOT. rirs_kernel` branches.
1231 66 : IF (rtbse_env%rirs_kernel) RETURN
1232 :
1233 : CALL create_sigma_workspace_qs_only(rtbse_env%qs_env, rtbse_env%screened_dbt, rtbse_env%w_dbcsr, &
1234 : rtbse_env%t_3c_w, rtbse_env%t_3c_work_RI_AO__AO, &
1235 50 : rtbse_env%t_3c_work2_RI_AO__AO, rtbse_env%greens_dbt)
1236 : END SUBROUTINE create_sigma_workspace
1237 : ! **************************************************************************************************
1238 : !> \brief Allocates the workspaces for self-energy determination routine
1239 : !> \note Does so without referencing the rtbse_env
1240 : !> \note References bs_env
1241 : !> \param rtbse_env Structure for holding information and workspace structures
1242 : !> \param qs_env Quickstep environment - entry point of calculation
1243 : !> \author Stepan Marek
1244 : !> \date 02.2024
1245 : ! **************************************************************************************************
1246 50 : SUBROUTINE create_sigma_workspace_qs_only(qs_env, screened_dbt, screened_dbcsr, int_3c_dbt, &
1247 : work_dbt_3c_1, work_dbt_3c_2, work_dbt_2c)
1248 : TYPE(qs_environment_type), POINTER :: qs_env
1249 : TYPE(dbcsr_type) :: screened_dbcsr
1250 : TYPE(dbt_type) :: screened_dbt, &
1251 : int_3c_dbt, &
1252 : work_dbt_3c_1, &
1253 : work_dbt_3c_2, &
1254 : work_dbt_2c
1255 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1256 :
1257 50 : CALL get_qs_env(qs_env, bs_env=bs_env)
1258 :
1259 : ! t_3c_w
1260 50 : CALL dbt_create(bs_env%t_RI__AO_AO, int_3c_dbt)
1261 : ! TODO : Provide option/decision whether to store the 3c integrals precomputed
1262 50 : CALL compute_3c_integrals(qs_env, bs_env, int_3c_dbt)
1263 : ! t_3c_work_RI_AO__AO
1264 50 : CALL dbt_create(bs_env%t_RI_AO__AO, work_dbt_3c_1)
1265 : ! t_3c_work2_RI_AO__AO
1266 50 : CALL dbt_create(bs_env%t_RI_AO__AO, work_dbt_3c_2)
1267 : ! t_W
1268 : ! Populate screened_dbt from gw run
1269 50 : CALL dbcsr_create(screened_dbcsr, name="W", template=bs_env%mat_RI_RI%matrix)
1270 50 : CALL dbt_create(screened_dbcsr, screened_dbt)
1271 : ! greens_dbt
1272 50 : CALL dbt_create(bs_env%mat_ao_ao%matrix, work_dbt_2c)
1273 50 : END SUBROUTINE create_sigma_workspace_qs_only
1274 : ! **************************************************************************************************
1275 : !> \brief Releases the workspaces for self-energy determination
1276 : !> \param rtbse_env
1277 : !> \author Stepan Marek
1278 : !> \date 02.2024
1279 : ! **************************************************************************************************
1280 66 : SUBROUTINE release_sigma_workspace(rtbse_env)
1281 : TYPE(rtbse_env_type) :: rtbse_env
1282 :
1283 : ! Mirror the gate in create_sigma_workspace.
1284 66 : IF (rtbse_env%rirs_kernel) RETURN
1285 :
1286 50 : CALL dbt_destroy(rtbse_env%t_3c_w)
1287 50 : CALL dbt_destroy(rtbse_env%t_3c_work_RI_AO__AO)
1288 50 : CALL dbt_destroy(rtbse_env%t_3c_work2_RI_AO__AO)
1289 50 : CALL dbt_destroy(rtbse_env%screened_dbt)
1290 50 : CALL dbt_destroy(rtbse_env%greens_dbt)
1291 50 : CALL dbcsr_release(rtbse_env%w_dbcsr)
1292 : END SUBROUTINE release_sigma_workspace
1293 : ! **************************************************************************************************
1294 : !> \brief Multiplies real matrix by a complex matrix from the right
1295 : !> \note So far only converts the real matrix to complex one, potentially doubling the work
1296 : !> \param rtbse_env
1297 : !> \author Stepan Marek
1298 : !> \date 09.2024
1299 : ! **************************************************************************************************
1300 15064 : SUBROUTINE multiply_fm_cfm(trans_r, trans_c, na, nb, nc, &
1301 : alpha, matrix_r, matrix_c, beta, res)
1302 : ! Transposition
1303 : CHARACTER(len=1) :: trans_r, trans_c
1304 : INTEGER :: na, nb, nc
1305 : ! accept real numbers
1306 : ! TODO : Just use complex numbers and import z_one, z_zero etc.
1307 : REAL(kind=dp) :: alpha, beta
1308 : TYPE(cp_fm_type) :: matrix_r
1309 : TYPE(cp_cfm_type) :: matrix_c, res
1310 : TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im
1311 : REAL(kind=dp) :: i_unit
1312 : CHARACTER(len=1) :: trans_cr
1313 :
1314 3766 : CALL cp_fm_create(work_re, matrix_c%matrix_struct)
1315 3766 : CALL cp_fm_create(work_im, matrix_c%matrix_struct)
1316 3766 : CALL cp_fm_create(res_re, res%matrix_struct)
1317 3766 : CALL cp_fm_create(res_im, res%matrix_struct)
1318 3766 : CALL cp_cfm_to_fm(matrix_c, work_re, work_im)
1319 0 : SELECT CASE (trans_c)
1320 : CASE ("C")
1321 0 : i_unit = -1.0_dp
1322 0 : trans_cr = "T"
1323 : CASE ("T")
1324 0 : i_unit = 1.0_dp
1325 0 : trans_cr = "T"
1326 : CASE default
1327 3766 : i_unit = 1.0_dp
1328 3766 : trans_cr = "N"
1329 : END SELECT
1330 : ! Actual multiplication
1331 : CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, &
1332 3766 : alpha, matrix_r, work_re, beta, res_re)
1333 : CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, &
1334 3766 : i_unit*alpha, matrix_r, work_im, beta, res_im)
1335 3766 : CALL cp_fm_to_cfm(res_re, res_im, res)
1336 3766 : CALL cp_fm_release(work_re)
1337 3766 : CALL cp_fm_release(work_im)
1338 3766 : CALL cp_fm_release(res_re)
1339 3766 : CALL cp_fm_release(res_im)
1340 :
1341 3766 : END SUBROUTINE multiply_fm_cfm
1342 : ! **************************************************************************************************
1343 : !> \brief Multiplies complex matrix by a real matrix from the right
1344 : !> \note So far only converts the real matrix to complex one, potentially doubling the work
1345 : !> \param rtbse_env
1346 : !> \author Stepan Marek
1347 : !> \date 09.2024
1348 : ! **************************************************************************************************
1349 5424 : SUBROUTINE multiply_cfm_fm(trans_c, trans_r, na, nb, nc, &
1350 : alpha, matrix_c, matrix_r, beta, res)
1351 : ! Transposition
1352 : CHARACTER(len=1) :: trans_c, trans_r
1353 : INTEGER :: na, nb, nc
1354 : ! accept real numbers
1355 : ! TODO : complex number support via interface?
1356 : REAL(kind=dp) :: alpha, beta
1357 : TYPE(cp_cfm_type) :: matrix_c, res
1358 : TYPE(cp_fm_type) :: matrix_r
1359 : TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im
1360 : REAL(kind=dp) :: i_unit
1361 : CHARACTER(len=1) :: trans_cr
1362 :
1363 1356 : CALL cp_fm_create(work_re, matrix_c%matrix_struct)
1364 1356 : CALL cp_fm_create(work_im, matrix_c%matrix_struct)
1365 1356 : CALL cp_fm_create(res_re, res%matrix_struct)
1366 1356 : CALL cp_fm_create(res_im, res%matrix_struct)
1367 1356 : CALL cp_cfm_to_fm(matrix_c, work_re, work_im)
1368 0 : SELECT CASE (trans_c)
1369 : CASE ("C")
1370 0 : i_unit = -1.0_dp
1371 0 : trans_cr = "T"
1372 : CASE ("T")
1373 0 : i_unit = 1.0_dp
1374 0 : trans_cr = "T"
1375 : CASE default
1376 1356 : i_unit = 1.0_dp
1377 1356 : trans_cr = "N"
1378 : END SELECT
1379 : ! Actual multiplication
1380 : CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, &
1381 1356 : alpha, work_re, matrix_r, beta, res_re)
1382 : CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, &
1383 1356 : i_unit*alpha, work_im, matrix_r, beta, res_im)
1384 1356 : CALL cp_fm_to_cfm(res_re, res_im, res)
1385 1356 : CALL cp_fm_release(work_re)
1386 1356 : CALL cp_fm_release(work_im)
1387 1356 : CALL cp_fm_release(res_re)
1388 1356 : CALL cp_fm_release(res_im)
1389 :
1390 1356 : END SUBROUTINE multiply_cfm_fm
1391 0 : END MODULE rt_bse_types
|
