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_gw, &
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 72 : 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 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
383 72 : 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 72 : NULLIFY (rtbse_env, mo_struct)
390 3384 : ALLOCATE (rtbse_env)
391 72 : 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 72 : input=input)
400 72 : bs_sec => section_vals_get_subs_vals(input, "PROPERTIES%BANDSTRUCTURE")
401 72 : 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 72 : 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 72 : 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 72 : rtbse_env%n_ao = bs_env%n_ao
416 : ! Number of auxiliary basis orbitals
417 72 : rtbse_env%n_RI = bs_env%n_RI
418 : ! Number of occupied orbitals - for closed shell equals to half the number of electrons
419 432 : rtbse_env%n_occ(:) = bs_env%n_occ(:)
420 : ! Spin degeneracy - number of spins per orbital
421 72 : rtbse_env%spin_degeneracy = bs_env%spin_degeneracy
422 : ! Default field is zero
423 288 : rtbse_env%field(:) = 0.0_dp
424 : ! Default time is zero
425 72 : rtbse_env%sim_step = 0
426 72 : rtbse_env%sim_time = 0
427 : ! Time step is taken from rtp
428 72 : md_sec => section_vals_get_subs_vals(force_env%root_section, "MOTION%MD")
429 72 : 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 72 : rtbse_env%etrs_threshold = rtbse_env%dft_control%rtp_control%eps_ener
433 72 : 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 72 : i_val=rtbse_env%ham_reference_type)
437 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAX_ITER", &
438 72 : i_val=rtbse_env%etrs_max_iter)
439 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%MAT_EXP", &
440 72 : i_val=rtbse_env%mat_exp_method)
441 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%ENERGY_CUTOFF_OCC", &
442 72 : r_val=rtbse_env%rtbse_energy_cutoff_occ)
443 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%ENERGY_CUTOFF_EMPTY", &
444 72 : r_val=rtbse_env%rtbse_energy_cutoff_empty)
445 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%TDA", &
446 72 : l_val=rtbse_env%tda_active)
447 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%TDA_SHIFT_TO_FIRST_PEAK", &
448 72 : 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 72 : l_val=rtbse_env%enforce_max_dt)
451 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%DEBUG_DISABLE_HARTREE", &
452 72 : l_val=rtbse_env%debug_disable_hartree)
453 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%DEBUG_DISABLE_SEX", &
454 72 : 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 72 : 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 72 : l_val=rtbse_env%diagnose_liouvillian_eig)
461 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%EPS_FILTER_RHO", &
462 72 : r_val=bs_env%ri_rs%eps_filter_rho)
463 : CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%RTBSE%CUTOFF_RADIUS_W0", &
464 72 : r_val=bs_env%ri_rs%cutoff_radius_w0)
465 :
466 72 : IF (.NOT. rtbse_env%dft_control%rtp_control%rtp_method == rtp_method_bse_linearized) THEN
467 14 : rtbse_env%rtbse_energy_cutoff_occ = -1.0_dp
468 14 : rtbse_env%rtbse_energy_cutoff_empty = -1.0_dp
469 14 : rtbse_env%enforce_max_dt = .FALSE.
470 14 : rtbse_env%debug_disable_hartree = .FALSE.
471 14 : rtbse_env%debug_disable_sex = .FALSE.
472 14 : rtbse_env%tda_shift_to_first_peak = .FALSE.
473 : ! rirs_kernel is already forced .FALSE. here by rtbse_resolve_rirs_flag.
474 14 : rtbse_env%diagnose_liouvillian_eig = .FALSE.
475 14 : bs_env%ri_rs%eps_filter_rho = -1.0_dp
476 14 : bs_env%ri_rs%cutoff_radius_w0 = -1.0_dp
477 : END IF
478 : ! First-peak shift only makes sense within TDA; force-disable otherwise.
479 72 : IF (.NOT. rtbse_env%tda_active) rtbse_env%tda_shift_to_first_peak = .FALSE.
480 72 : rtbse_env%omega_shift = 0.0_dp
481 :
482 72 : IF (rtbse_env%tda_active .AND. .NOT. rtbse_env%linearized) THEN
483 0 : CPABORT("RTBSE TDA keyword requires LINEARIZED_BSE_PROPAGATION=.TRUE.")
484 : END IF
485 : ! Open shell: omega_shift would be referenced to a non-physical cross-spin
486 : ! pseudo-gap (global MIN/MAX over both spins); abort until made per-spin.
487 72 : IF (rtbse_env%tda_shift_to_first_peak .AND. rtbse_env%n_spin > 1) THEN
488 : CALL cp_abort(__LOCATION__, &
489 : "TDA_SHIFT_TO_FIRST_PEAK is not implemented for open-shell (n_spin>1) "// &
490 : "systems - the first-peak gap estimate would mix spin channels. "// &
491 0 : "Set TDA_SHIFT_TO_FIRST_PEAK=.FALSE. for open-shell runs.")
492 : END IF
493 72 : CALL check_qp_gap_sanity(rtbse_env, bs_env)
494 72 : CALL determine_active_mo_window(rtbse_env, bs_env)
495 : ! Owned active-MO matrix structure (currently identical to the full n_ao x n_ao struct
496 : ! when no truncation is active; will be used by the linearized RT-BSE allocation path).
497 72 : NULLIFY (rtbse_env%fm_struct_mo_active)
498 : CALL cp_fm_struct_create(rtbse_env%fm_struct_mo_active, &
499 : bs_env%fm_ks_Gamma(1)%matrix_struct%para_env, &
500 : bs_env%fm_ks_Gamma(1)%matrix_struct%context, &
501 72 : rtbse_env%mo_active, rtbse_env%mo_active)
502 : ! Rectangular n_ao x mo_active struct used for C_active and AO<->MO intermediates
503 72 : NULLIFY (rtbse_env%fm_struct_ao_mo_active)
504 : CALL cp_fm_struct_create(rtbse_env%fm_struct_ao_mo_active, &
505 : bs_env%fm_ks_Gamma(1)%matrix_struct%para_env, &
506 : bs_env%fm_ks_Gamma(1)%matrix_struct%context, &
507 72 : rtbse_env%n_ao, rtbse_env%mo_active)
508 : ! Choose the matrix struct used for MO-side persistent matrices.
509 : ! Linearized RT-BSE: mo_active x mo_active. Full RT-BSE: full AO struct (unchanged).
510 72 : IF (rtbse_env%linearized) THEN
511 58 : mo_struct => rtbse_env%fm_struct_mo_active
512 : ELSE
513 14 : mo_struct => bs_env%fm_ks_Gamma(1)%matrix_struct
514 : END IF
515 : ! Output unit number, recovered from the post_scf_bandstructure_type
516 72 : rtbse_env%unit_nr = bs_env%unit_nr
517 : ! Sim start index and total number of steps as well
518 72 : CALL section_vals_val_get(md_sec, "STEP_START_VAL", i_val=rtbse_env%sim_start)
519 : ! Copy this value to sim_start_orig for continuation runs
520 72 : rtbse_env%sim_start_orig = rtbse_env%sim_start
521 72 : CALL section_vals_val_get(md_sec, "STEPS", i_val=rtbse_env%sim_nsteps)
522 : ! Get the values for the FT
523 72 : rtbse_env%ft_damping = rtbse_env%dft_control%rtp_control%ft_damping
524 72 : rtbse_env%ft_damping = rtbse_env%dft_control%rtp_control%ft_t0
525 72 : rtbse_env%pol_elements => rtbse_env%dft_control%rtp_control%print_pol_elements
526 :
527 72 : rtbse_env%rtp_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION")
528 : ! Get the restart section
529 72 : rtbse_env%restart_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%RESTART")
530 72 : rtbse_env%restart_extracted = .FALSE.
531 72 : rtbse_env%field_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%FIELD")
532 72 : rtbse_env%moments_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%MOMENTS")
533 : ! Moment specification
534 : CALL section_vals_val_get(rtbse_env%rtp_section, "PRINT%MOMENTS%REFERENCE", &
535 72 : i_val=rtbse_env%moment_ref_type)
536 : CALL section_vals_val_get(rtbse_env%rtp_section, "PRINT%MOMENTS%REFERENCE_POINT", &
537 72 : r_vals=rtbse_env%user_moment_ref_point)
538 72 : rtbse_env%rho_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%DENSITY_MATRIX")
539 72 : rtbse_env%ft_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%MOMENTS_FT")
540 72 : rtbse_env%pol_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%POLARIZABILITY")
541 72 : rtbse_env%eig_section => section_vals_get_subs_vals(rtbse_env%rtp_section, "PRINT%LIOUVILLIAN_EIG")
542 : ! Warn the user about print sections which are not yet implemented in the RTBSE run
543 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%CURRENT", &
544 72 : "CURRENT print section not yet implemented for RTBSE.")
545 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%E_CONSTITUENTS", &
546 72 : "E_CONSTITUENTS print section not yet implemented for RTBSE.")
547 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%PROGRAM_RUN_INFO", &
548 72 : "PROGRAM_RUN_INFO print section not yet implemented for RTBSE.")
549 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%PROJECTION_MO", &
550 72 : "PROJECTION_MO print section not yet implemented for RTBSE.")
551 : CALL warn_section_unused(rtbse_env%rtp_section, "PRINT%RESTART_HISTORY", &
552 72 : "RESTART_HISTORY print section not yet implemented for RTBSE.")
553 : ! References to the parent qs_env / bs_env
554 72 : rtbse_env%qs_env => force_env%qs_env
555 72 : rtbse_env%bs_env => bs_env
556 : ! Padé refinement
557 72 : rtbse_env%pade_requested = rtbse_env%dft_control%rtp_control%pade_requested
558 72 : rtbse_env%pade_e_min = rtbse_env%dft_control%rtp_control%pade_e_min
559 72 : rtbse_env%pade_e_step = rtbse_env%dft_control%rtp_control%pade_e_step
560 72 : rtbse_env%pade_e_max = rtbse_env%dft_control%rtp_control%pade_e_max
561 72 : rtbse_env%pade_fit_e_min = rtbse_env%dft_control%rtp_control%pade_fit_e_min
562 72 : rtbse_env%pade_fit_e_max = rtbse_env%dft_control%rtp_control%pade_fit_e_max
563 72 : rtbse_env%pade_npoints = INT((rtbse_env%pade_e_max - rtbse_env%pade_e_min)/rtbse_env%pade_e_step)
564 : ! Evaluate the evaluation grid
565 72 : IF (rtbse_env%pade_requested) THEN
566 2 : NULLIFY (rtbse_env%pade_x_eval)
567 6 : ALLOCATE (rtbse_env%pade_x_eval(rtbse_env%pade_npoints))
568 2000 : DO i = 1, rtbse_env%pade_npoints
569 2000 : rtbse_env%pade_x_eval(i) = CMPLX(rtbse_env%pade_e_step*REAL(i - 1, kind=dp), 0.0, kind=dp)
570 : END DO
571 : END IF
572 :
573 : ! Allocate moments matrices.
574 : ! In linearized RT-BSE these store the MO-active transformed dipole moments;
575 : ! in full RT-BSE they remain AO-sized (initialized from overlap template).
576 72 : NULLIFY (rtbse_env%moments)
577 536 : ALLOCATE (rtbse_env%moments(3, rtbse_env%n_spin))
578 72 : NULLIFY (rtbse_env%moments_field)
579 536 : ALLOCATE (rtbse_env%moments_field(3, rtbse_env%n_spin))
580 152 : DO i_spin = 1, rtbse_env%n_spin
581 392 : DO k = 1, 3
582 240 : CALL cp_fm_create(rtbse_env%moments(k, i_spin), mo_struct)
583 320 : CALL cp_fm_create(rtbse_env%moments_field(k, i_spin), mo_struct)
584 : END DO
585 : END DO
586 :
587 : ! Allocate space for density propagation and other operations.
588 : ! In linearized RT-BSE these workspaces are MO-active sized; in full RT-BSE
589 : ! they remain at the full AO size.
590 72 : NULLIFY (rtbse_env%rho_workspace)
591 360 : ALLOCATE (rtbse_env%rho_workspace(4))
592 360 : DO i = 1, SIZE(rtbse_env%rho_workspace)
593 288 : CALL cp_cfm_create(rtbse_env%rho_workspace(i), mo_struct)
594 360 : CALL cp_cfm_set_all(rtbse_env%rho_workspace(i), CMPLX(0.0, 0.0, kind=dp))
595 : END DO
596 :
597 : ! TODO: gate workspace allocation so methods skip workspaces they don't need
598 :
599 : ! Allocate real workspace
600 72 : NULLIFY (rtbse_env%real_workspace)
601 72 : IF (rtbse_env%linearized) THEN
602 174 : ALLOCATE (rtbse_env%real_workspace(2))
603 : ELSE
604 14 : SELECT CASE (rtbse_env%mat_exp_method)
605 : CASE (do_exact)
606 0 : ALLOCATE (rtbse_env%real_workspace(4))
607 : CASE (do_bch)
608 42 : ALLOCATE (rtbse_env%real_workspace(2))
609 : CASE DEFAULT
610 14 : CPABORT("Only exact and BCH matrix propagation implemented in RT-BSE")
611 : END SELECT
612 : END IF
613 216 : DO i = 1, SIZE(rtbse_env%real_workspace)
614 144 : CALL cp_fm_create(rtbse_env%real_workspace(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
615 216 : CALL cp_fm_set_all(rtbse_env%real_workspace(i), 0.0_dp)
616 : END DO
617 72 : NULLIFY (rtbse_env%sigma_complex_workspace)
618 144 : ALLOCATE (rtbse_env%sigma_complex_workspace(1))
619 72 : CALL cp_cfm_create(rtbse_env%sigma_complex_workspace(1), bs_env%fm_ks_Gamma(1)%matrix_struct)
620 72 : CALL cp_cfm_set_all(rtbse_env%sigma_complex_workspace(1), CMPLX(0.0_dp, 0.0_dp, kind=dp))
621 : ! Allocate density matrix (MO-active sized when linearized; AO-sized otherwise)
622 72 : NULLIFY (rtbse_env%rho)
623 296 : ALLOCATE (rtbse_env%rho(rtbse_env%n_spin))
624 152 : DO i = 1, rtbse_env%n_spin
625 152 : CALL cp_cfm_create(rtbse_env%rho(i), matrix_struct=mo_struct)
626 : END DO
627 : ! Allocate additional space for AO density matrix
628 : ! in linearised RTBSE, where default is MO
629 72 : IF (rtbse_env%linearized) THEN
630 58 : NULLIFY (rtbse_env%rho_ao_scratch)
631 240 : ALLOCATE (rtbse_env%rho_ao_scratch(rtbse_env%n_spin))
632 124 : DO i = 1, rtbse_env%n_spin
633 124 : CALL cp_cfm_create(rtbse_env%rho_ao_scratch(i), matrix_struct=bs_env%fm_s_Gamma%matrix_struct)
634 : END DO
635 : END IF
636 : ! Create the inverse overlap matrix, for use in density propagation
637 : ! Start by creating the actual overlap matrix
638 72 : CALL cp_fm_create(rtbse_env%S_fm, bs_env%fm_s_Gamma%matrix_struct)
639 72 : CALL cp_fm_create(rtbse_env%S_inv_fm, bs_env%fm_s_Gamma%matrix_struct)
640 72 : CALL cp_cfm_create(rtbse_env%S_cfm, bs_env%fm_s_Gamma%matrix_struct)
641 :
642 : ! Create the single particle hamiltonian
643 : ! Allocate workspace (MO-active sized in linearized RT-BSE; AO sized otherwise)
644 72 : NULLIFY (rtbse_env%ham_workspace)
645 296 : ALLOCATE (rtbse_env%ham_workspace(rtbse_env%n_spin))
646 152 : DO i = 1, rtbse_env%n_spin
647 80 : CALL cp_cfm_create(rtbse_env%ham_workspace(i), mo_struct)
648 152 : CALL cp_cfm_set_all(rtbse_env%ham_workspace(i), CMPLX(0.0, 0.0, kind=dp))
649 : END DO
650 : ! Now onto the Hamiltonian itself
651 : ! full RTBSE: Contains energy differences and Hartree/COHSEX ρ_0 parts
652 : ! linearised RTBSE: Contains only the Hartree/SEX ρ_0 parts as Δε * Δρ(t) need to be updated
653 72 : NULLIFY (rtbse_env%ham_reference)
654 296 : ALLOCATE (rtbse_env%ham_reference(rtbse_env%n_spin))
655 152 : DO i = 1, rtbse_env%n_spin
656 152 : CALL cp_cfm_create(rtbse_env%ham_reference(i), mo_struct)
657 : END DO
658 : ! Single particle Hamiltonian (Δε * Δρ(t)) for updates during timesteps in LR-RTBSE
659 72 : IF (rtbse_env%linearized) THEN
660 58 : NULLIFY (rtbse_env%ham_reference_singleparticle)
661 240 : ALLOCATE (rtbse_env%ham_reference_singleparticle(rtbse_env%n_spin))
662 124 : DO i = 1, rtbse_env%n_spin
663 124 : CALL cp_cfm_create(rtbse_env%ham_reference_singleparticle(i), mo_struct)
664 : END DO
665 58 : NULLIFY (rtbse_env%eps_active)
666 232 : ALLOCATE (rtbse_env%eps_active(rtbse_env%mo_active, rtbse_env%n_spin))
667 904 : rtbse_env%eps_active(:, :) = 0.0_dp
668 : END IF
669 :
670 : ! Create the matrices and workspaces for ETRS propagation
671 72 : NULLIFY (rtbse_env%ham_effective)
672 72 : NULLIFY (rtbse_env%rho_new)
673 72 : NULLIFY (rtbse_env%rho_new_last)
674 72 : NULLIFY (rtbse_env%rho_M)
675 72 : NULLIFY (rtbse_env%rho_orig)
676 296 : ALLOCATE (rtbse_env%ham_effective(rtbse_env%n_spin))
677 296 : ALLOCATE (rtbse_env%rho_new(rtbse_env%n_spin))
678 296 : ALLOCATE (rtbse_env%rho_new_last(rtbse_env%n_spin))
679 296 : ALLOCATE (rtbse_env%rho_M(rtbse_env%n_spin))
680 296 : ALLOCATE (rtbse_env%rho_orig(rtbse_env%n_spin))
681 152 : DO i = 1, rtbse_env%n_spin
682 80 : CALL cp_cfm_create(rtbse_env%ham_effective(i), mo_struct)
683 80 : CALL cp_cfm_set_all(rtbse_env%ham_effective(i), CMPLX(0.0, 0.0, kind=dp))
684 80 : CALL cp_cfm_create(rtbse_env%rho_new(i), mo_struct)
685 80 : CALL cp_cfm_set_all(rtbse_env%rho_new(i), CMPLX(0.0, 0.0, kind=dp))
686 80 : CALL cp_cfm_create(rtbse_env%rho_new_last(i), mo_struct)
687 80 : CALL cp_cfm_set_all(rtbse_env%rho_new_last(i), CMPLX(0.0, 0.0, kind=dp))
688 80 : CALL cp_cfm_create(rtbse_env%rho_M(i), mo_struct)
689 80 : CALL cp_cfm_set_all(rtbse_env%rho_M(i), CMPLX(0.0, 0.0, kind=dp))
690 152 : CALL cp_cfm_create(rtbse_env%rho_orig(i), mo_struct)
691 : END DO
692 :
693 : !For LR-RTBSE we need RK4 coefficients - create new workspace
694 72 : IF (rtbse_env%linearized) THEN
695 : ! Indexed by SPIN, not by RK4 stage: the spin loop is inner to each stage (do_rk4_stage), so every
696 : ! spin's current-stage k must be live at once, but only one stage's k per spin - each is folded
697 : ! into rho_end and the next stage density before the next stage overwrites it. Hence size n_spin.
698 58 : NULLIFY (rtbse_env%rk4_coefficients)
699 240 : ALLOCATE (rtbse_env%rk4_coefficients(rtbse_env%n_spin))
700 124 : DO i = 1, rtbse_env%n_spin
701 66 : CALL cp_cfm_create(rtbse_env%rk4_coefficients(i), mo_struct)
702 124 : CALL cp_cfm_set_all(rtbse_env%rk4_coefficients(i), CMPLX(0.0, 0.0, kind=dp))
703 : END DO
704 : END IF
705 :
706 : ! Fields for exact diagonalisation
707 72 : NULLIFY (rtbse_env%real_eigvals)
708 216 : ALLOCATE (rtbse_env%real_eigvals(rtbse_env%n_ao))
709 796 : rtbse_env%real_eigvals(:) = 0.0_dp
710 72 : NULLIFY (rtbse_env%exp_eigvals)
711 216 : ALLOCATE (rtbse_env%exp_eigvals(rtbse_env%n_ao))
712 796 : rtbse_env%exp_eigvals(:) = CMPLX(0.0, 0.0, kind=dp)
713 :
714 : ! Workspace for FT - includes (in principle) the zeroth step and the extra last step
715 72 : NULLIFY (rtbse_env%moments_trace)
716 : ! TODO : Unite the number of steps with TD-DFT
717 19856 : ALLOCATE (rtbse_env%moments_trace(rtbse_env%n_spin, 3, rtbse_env%sim_nsteps + 1), source=z_zero)
718 72 : NULLIFY (rtbse_env%field_trace)
719 11144 : ALLOCATE (rtbse_env%field_trace(3, rtbse_env%sim_nsteps + 1), source=z_zero)
720 72 : NULLIFY (rtbse_env%time_trace)
721 2948 : ALLOCATE (rtbse_env%time_trace(rtbse_env%sim_nsteps + 1), source=0.0_dp)
722 :
723 : ! Allocate self-energy parts and dynamic Hartree potential.
724 : ! In linearized RT-BSE these matrices hold the MO-active-sized result of the
725 : ! AO->MO transform; the AO-sized buffer is allocated as sigma_*_ao below.
726 72 : NULLIFY (rtbse_env%hartree_curr)
727 72 : NULLIFY (rtbse_env%sigma_SEX)
728 72 : NULLIFY (rtbse_env%sigma_COH)
729 296 : ALLOCATE (rtbse_env%hartree_curr(rtbse_env%n_spin))
730 296 : ALLOCATE (rtbse_env%sigma_SEX(rtbse_env%n_spin))
731 296 : ALLOCATE (rtbse_env%sigma_COH(rtbse_env%n_spin))
732 152 : DO i = 1, rtbse_env%n_spin
733 80 : CALL cp_fm_create(rtbse_env%sigma_COH(i), mo_struct)
734 80 : CALL cp_cfm_create(rtbse_env%sigma_SEX(i), mo_struct)
735 80 : CALL cp_fm_create(rtbse_env%hartree_curr(i), mo_struct)
736 80 : CALL cp_fm_set_all(rtbse_env%sigma_COH(i), 0.0_dp)
737 80 : CALL cp_cfm_set_all(rtbse_env%sigma_SEX(i), CMPLX(0.0, 0.0, kind=dp))
738 152 : CALL cp_fm_set_all(rtbse_env%hartree_curr(i), 0.0_dp)
739 : END DO
740 : ! AO-sized scratch buffers used by the linearized RT-BSE path
741 72 : IF (rtbse_env%linearized) THEN
742 58 : NULLIFY (rtbse_env%hartree_curr_ao)
743 58 : NULLIFY (rtbse_env%sigma_SEX_ao)
744 240 : ALLOCATE (rtbse_env%hartree_curr_ao(rtbse_env%n_spin))
745 240 : ALLOCATE (rtbse_env%sigma_SEX_ao(rtbse_env%n_spin))
746 124 : DO i = 1, rtbse_env%n_spin
747 66 : CALL cp_cfm_create(rtbse_env%sigma_SEX_ao(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
748 66 : CALL cp_fm_create(rtbse_env%hartree_curr_ao(i), bs_env%fm_ks_Gamma(1)%matrix_struct)
749 66 : CALL cp_cfm_set_all(rtbse_env%sigma_SEX_ao(i), CMPLX(0.0, 0.0, kind=dp))
750 124 : CALL cp_fm_set_all(rtbse_env%hartree_curr_ao(i), 0.0_dp)
751 : END DO
752 : ! mo_active x mo_active real workspace pair for MO-side intermediates
753 58 : NULLIFY (rtbse_env%real_workspace_mo)
754 174 : ALLOCATE (rtbse_env%real_workspace_mo(2))
755 174 : DO i = 1, SIZE(rtbse_env%real_workspace_mo)
756 116 : CALL cp_fm_create(rtbse_env%real_workspace_mo(i), rtbse_env%fm_struct_mo_active)
757 174 : CALL cp_fm_set_all(rtbse_env%real_workspace_mo(i), 0.0_dp)
758 : END DO
759 58 : NULLIFY (rtbse_env%ao_mo_workspace)
760 116 : ALLOCATE (rtbse_env%ao_mo_workspace(1))
761 58 : CALL cp_fm_create(rtbse_env%ao_mo_workspace(1), rtbse_env%fm_struct_ao_mo_active)
762 58 : CALL cp_fm_set_all(rtbse_env%ao_mo_workspace(1), 0.0_dp)
763 : ! Truncated MO coefficient slabs C_active (n_ao x mo_active) for each spin.
764 : ! Filled in initialize_rtbse_env from bs_env%fm_mo_coeff_Gamma via submatrix copy.
765 58 : NULLIFY (rtbse_env%C_active)
766 240 : ALLOCATE (rtbse_env%C_active(rtbse_env%n_spin))
767 124 : DO i = 1, rtbse_env%n_spin
768 66 : CALL cp_fm_create(rtbse_env%C_active(i), rtbse_env%fm_struct_ao_mo_active)
769 124 : CALL cp_fm_set_all(rtbse_env%C_active(i), 0.0_dp)
770 : END DO
771 : ! Masked-copy staging scratch for the builder (all propagation paths including closed-shell ABBA).
772 : ! Also used as conjugate-transpose scratch in the TDA consumer.
773 58 : NULLIFY (rtbse_env%rho_delta_mo)
774 240 : ALLOCATE (rtbse_env%rho_delta_mo(rtbse_env%n_spin))
775 124 : DO i = 1, rtbse_env%n_spin
776 66 : CALL cp_cfm_create(rtbse_env%rho_delta_mo(i), rtbse_env%fm_struct_mo_active)
777 124 : CALL cp_cfm_set_all(rtbse_env%rho_delta_mo(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
778 : END DO
779 : ! Shared AO Hartree buffers: widened from (tda_active .OR. n_spin>1) so n_spin=1 ABBA
780 : ! (diagnostic + propagator) gets a dedicated buffer instead of aliasing sigma_SEX_ao.
781 58 : IF (.NOT. rtbse_env%debug_disable_hartree) THEN
782 54 : CALL cp_cfm_create(rtbse_env%rho_total_ao_scratch, bs_env%fm_s_Gamma%matrix_struct)
783 54 : CALL cp_cfm_create(rtbse_env%hartree_total_ao, bs_env%fm_ks_Gamma(1)%matrix_struct)
784 54 : CALL cp_cfm_set_all(rtbse_env%rho_total_ao_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
785 54 : CALL cp_cfm_set_all(rtbse_env%hartree_total_ao, CMPLX(0.0_dp, 0.0_dp, kind=dp))
786 : END IF
787 : ! Liouvillian-eigenvalue diagnostic state. n_spin=1 enforced upstream.
788 : ! Shared scratch (TDA + ABBA) allocated on diagnose_liouvillian_eig=T;
789 : ! ABBA-only A/B/A±B blocks added below under .NOT. tda_active. Mirror of the
790 : ! existing shared-scratch pattern, with the extra tda_active gate as the deviation.
791 58 : IF (rtbse_env%diagnose_liouvillian_eig) THEN
792 : ! Joint OV dimension: the spin blocks are stacked (n_spin=1 -> the old single-spin
793 : ! size). drho_probe/L_drho stay mo_active-sized per spin; L_pairs is N_OV_joint.
794 58 : n_ov = 0
795 124 : DO i = 1, rtbse_env%n_spin
796 : n_ov = n_ov + (rtbse_env%n_occ(i) - rtbse_env%first_active_mo + 1)* &
797 124 : (rtbse_env%last_active_mo - rtbse_env%n_occ(i))
798 : END DO
799 58 : NULLIFY (rtbse_env%fm_struct_ov_pairs)
800 : CALL cp_fm_struct_create(rtbse_env%fm_struct_ov_pairs, &
801 : bs_env%fm_ks_Gamma(1)%matrix_struct%para_env, &
802 : bs_env%fm_ks_Gamma(1)%matrix_struct%context, &
803 58 : n_ov, n_ov)
804 240 : ALLOCATE (rtbse_env%drho_probe(rtbse_env%n_spin))
805 240 : ALLOCATE (rtbse_env%L_drho(rtbse_env%n_spin))
806 124 : DO i = 1, rtbse_env%n_spin
807 66 : CALL cp_cfm_create(rtbse_env%drho_probe(i), rtbse_env%fm_struct_mo_active)
808 66 : CALL cp_cfm_create(rtbse_env%L_drho(i), rtbse_env%fm_struct_mo_active)
809 66 : CALL cp_cfm_set_all(rtbse_env%drho_probe(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
810 124 : CALL cp_cfm_set_all(rtbse_env%L_drho(i), CMPLX(0.0_dp, 0.0_dp, kind=dp))
811 : END DO
812 58 : CALL cp_cfm_create(rtbse_env%L_pairs, rtbse_env%fm_struct_ov_pairs)
813 58 : CALL cp_cfm_create(rtbse_env%eigvecs_pairs, rtbse_env%fm_struct_ov_pairs)
814 58 : CALL cp_cfm_set_all(rtbse_env%L_pairs, CMPLX(0.0_dp, 0.0_dp, kind=dp))
815 58 : CALL cp_cfm_set_all(rtbse_env%eigvecs_pairs, CMPLX(0.0_dp, 0.0_dp, kind=dp))
816 58 : NULLIFY (rtbse_env%eigenvalues_liouvillian)
817 174 : ALLOCATE (rtbse_env%eigenvalues_liouvillian(n_ov))
818 772 : rtbse_env%eigenvalues_liouvillian = 0.0_dp
819 : ! ABBA-only Furche-reduction scratch (A, B, A-B->sqrt, A+B).
820 58 : IF (.NOT. rtbse_env%tda_active) THEN
821 18 : CALL cp_cfm_create(rtbse_env%A_mat, rtbse_env%fm_struct_ov_pairs)
822 18 : CALL cp_cfm_create(rtbse_env%B_mat, rtbse_env%fm_struct_ov_pairs)
823 18 : CALL cp_cfm_create(rtbse_env%AmB_scratch, rtbse_env%fm_struct_ov_pairs)
824 18 : CALL cp_cfm_create(rtbse_env%ApB_scratch, rtbse_env%fm_struct_ov_pairs)
825 18 : CALL cp_cfm_set_all(rtbse_env%A_mat, CMPLX(0.0_dp, 0.0_dp, kind=dp))
826 18 : CALL cp_cfm_set_all(rtbse_env%B_mat, CMPLX(0.0_dp, 0.0_dp, kind=dp))
827 18 : CALL cp_cfm_set_all(rtbse_env%AmB_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
828 18 : CALL cp_cfm_set_all(rtbse_env%ApB_scratch, CMPLX(0.0_dp, 0.0_dp, kind=dp))
829 : END IF
830 : END IF
831 : END IF
832 :
833 : ! Allocate workspaces for get_sigma
834 72 : CALL create_sigma_workspace(rtbse_env)
835 :
836 : ! Depending on the chosen methods, allocate extra workspace
837 72 : CALL create_hartree_ri_workspace(rtbse_env)
838 :
839 72 : END SUBROUTINE create_rtbse_env
840 :
841 : ! **************************************************************************************************
842 : !> \brief Simple reimplementation of cp_fm_release_pp1 for complex matrices
843 : !> \param matrices cp_cfm_type(:)
844 : !> \author Stepan Marek
845 : !> \date 02.2024
846 : ! **************************************************************************************************
847 1198 : SUBROUTINE cp_cfm_release_pa1(matrices)
848 : TYPE(cp_cfm_type), DIMENSION(:), POINTER :: matrices
849 : INTEGER :: i
850 :
851 2740 : DO i = 1, SIZE(matrices)
852 2740 : CALL cp_cfm_release(matrices(i))
853 : END DO
854 1198 : DEALLOCATE (matrices)
855 : NULLIFY (matrices)
856 1198 : END SUBROUTINE cp_cfm_release_pa1
857 :
858 : ! **************************************************************************************************
859 : !> \brief Releases the environment allocated structures
860 : !> \param rtbse_env
861 : !> \author Stepan Marek
862 : !> \date 02.2024
863 : ! **************************************************************************************************
864 72 : SUBROUTINE release_rtbse_env(rtbse_env)
865 : TYPE(rtbse_env_type), POINTER :: rtbse_env
866 :
867 72 : CALL cp_cfm_release_pa1(rtbse_env%ham_effective)
868 72 : CALL cp_cfm_release_pa1(rtbse_env%ham_workspace)
869 72 : CALL cp_fm_release(rtbse_env%sigma_COH)
870 72 : CALL cp_cfm_release_pa1(rtbse_env%sigma_SEX)
871 72 : CALL cp_fm_release(rtbse_env%hartree_curr)
872 72 : CALL cp_cfm_release_pa1(rtbse_env%ham_reference)
873 72 : IF (ASSOCIATED(rtbse_env%ham_reference_singleparticle)) THEN
874 58 : CALL cp_cfm_release_pa1(rtbse_env%ham_reference_singleparticle)
875 : END IF
876 72 : IF (ASSOCIATED(rtbse_env%eps_active)) DEALLOCATE (rtbse_env%eps_active)
877 72 : IF (ASSOCIATED(rtbse_env%eps_active_restart)) DEALLOCATE (rtbse_env%eps_active_restart)
878 72 : CALL cp_cfm_release_pa1(rtbse_env%rho)
879 72 : CALL cp_cfm_release_pa1(rtbse_env%rho_workspace)
880 72 : CALL cp_cfm_release_pa1(rtbse_env%rho_new)
881 72 : CALL cp_cfm_release_pa1(rtbse_env%rho_new_last)
882 72 : CALL cp_cfm_release_pa1(rtbse_env%rho_M)
883 72 : CALL cp_cfm_release_pa1(rtbse_env%rho_orig)
884 72 : IF (ASSOCIATED(rtbse_env%rk4_coefficients)) THEN
885 58 : CALL cp_cfm_release_pa1(rtbse_env%rk4_coefficients)
886 : END IF
887 72 : CALL cp_fm_release(rtbse_env%real_workspace)
888 72 : IF (ASSOCIATED(rtbse_env%sigma_complex_workspace)) CALL cp_cfm_release_pa1(rtbse_env%sigma_complex_workspace)
889 72 : CALL cp_fm_release(rtbse_env%S_inv_fm)
890 72 : CALL cp_fm_release(rtbse_env%S_fm)
891 72 : CALL cp_cfm_release(rtbse_env%S_cfm)
892 :
893 72 : CALL cp_fm_release(rtbse_env%moments)
894 72 : CALL cp_fm_release(rtbse_env%moments_field)
895 :
896 72 : CALL release_sigma_workspace(rtbse_env)
897 :
898 72 : CALL release_hartree_ri_workspace(rtbse_env)
899 :
900 72 : DEALLOCATE (rtbse_env%real_eigvals)
901 72 : DEALLOCATE (rtbse_env%exp_eigvals)
902 72 : DEALLOCATE (rtbse_env%moments_trace)
903 72 : DEALLOCATE (rtbse_env%field_trace)
904 72 : DEALLOCATE (rtbse_env%time_trace)
905 :
906 72 : IF (ASSOCIATED(rtbse_env%pol_elements)) DEALLOCATE (rtbse_env%pol_elements)
907 72 : IF (ASSOCIATED(rtbse_env%pade_x_eval)) DEALLOCATE (rtbse_env%pade_x_eval)
908 :
909 : ! Deallocate the neighbour list that is not deallocated in gw anymore
910 72 : IF (ASSOCIATED(rtbse_env%bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(rtbse_env%bs_env%nl_3c)
911 : ! Release linearized-only AO scratches and MO-side workspaces
912 72 : IF (ASSOCIATED(rtbse_env%rho_ao_scratch)) CALL cp_cfm_release_pa1(rtbse_env%rho_ao_scratch)
913 72 : IF (ASSOCIATED(rtbse_env%sigma_SEX_ao)) CALL cp_cfm_release_pa1(rtbse_env%sigma_SEX_ao)
914 72 : IF (ASSOCIATED(rtbse_env%hartree_curr_ao)) CALL cp_fm_release(rtbse_env%hartree_curr_ao)
915 72 : IF (ASSOCIATED(rtbse_env%real_workspace_mo)) CALL cp_fm_release(rtbse_env%real_workspace_mo)
916 72 : IF (ASSOCIATED(rtbse_env%ao_mo_workspace)) CALL cp_fm_release(rtbse_env%ao_mo_workspace)
917 72 : IF (ASSOCIATED(rtbse_env%C_active)) CALL cp_fm_release(rtbse_env%C_active)
918 72 : IF (ASSOCIATED(rtbse_env%rho_delta_mo)) CALL cp_cfm_release_pa1(rtbse_env%rho_delta_mo)
919 : ! Release shared bare-Hartree scratch. Mirror the alloc gate exactly (linearized .AND.
920 : ! .NOT. debug_disable_hartree, every shell incl closed-shell ABBA) — the old
921 : ! (tda_active .OR. n_spin>1) gate leaked both buffers on the closed-shell ABBA path.
922 72 : IF (rtbse_env%linearized .AND. .NOT. rtbse_env%debug_disable_hartree) THEN
923 54 : CALL cp_cfm_release(rtbse_env%rho_total_ao_scratch)
924 54 : CALL cp_cfm_release(rtbse_env%hartree_total_ao)
925 : END IF
926 : ! Release the RI-RS Hartree diagonal-reuse accumulators (allocated in initialize_hartree_potential).
927 72 : IF (ALLOCATED(rtbse_env%hartree_diag_re)) DEALLOCATE (rtbse_env%hartree_diag_re)
928 72 : IF (ALLOCATED(rtbse_env%hartree_diag_im)) DEALLOCATE (rtbse_env%hartree_diag_im)
929 : ! Release Liouvillian-diagnostic scratch (only when the diagnostic was requested).
930 72 : IF (rtbse_env%diagnose_liouvillian_eig) THEN
931 58 : IF (ASSOCIATED(rtbse_env%drho_probe)) CALL cp_cfm_release_pa1(rtbse_env%drho_probe)
932 58 : IF (ASSOCIATED(rtbse_env%L_drho)) CALL cp_cfm_release_pa1(rtbse_env%L_drho)
933 58 : CALL cp_cfm_release(rtbse_env%L_pairs)
934 58 : CALL cp_cfm_release(rtbse_env%eigvecs_pairs)
935 58 : IF (ASSOCIATED(rtbse_env%eigenvalues_liouvillian)) DEALLOCATE (rtbse_env%eigenvalues_liouvillian)
936 58 : IF (.NOT. rtbse_env%tda_active) THEN
937 18 : CALL cp_cfm_release(rtbse_env%A_mat)
938 18 : CALL cp_cfm_release(rtbse_env%B_mat)
939 18 : CALL cp_cfm_release(rtbse_env%AmB_scratch)
940 18 : CALL cp_cfm_release(rtbse_env%ApB_scratch)
941 : END IF
942 58 : IF (ASSOCIATED(rtbse_env%fm_struct_ov_pairs)) THEN
943 58 : CALL cp_fm_struct_release(rtbse_env%fm_struct_ov_pairs)
944 : END IF
945 : END IF
946 : ! Release owned active-MO matrix structures
947 72 : IF (ASSOCIATED(rtbse_env%fm_struct_mo_active)) THEN
948 72 : CALL cp_fm_struct_release(rtbse_env%fm_struct_mo_active)
949 : END IF
950 72 : IF (ASSOCIATED(rtbse_env%fm_struct_ao_mo_active)) THEN
951 72 : CALL cp_fm_struct_release(rtbse_env%fm_struct_ao_mo_active)
952 : END IF
953 : ! Deallocate the storage for the environment itself
954 72 : DEALLOCATE (rtbse_env)
955 : ! Nullify to make sure it is not used again
956 : NULLIFY (rtbse_env)
957 :
958 72 : END SUBROUTINE release_rtbse_env
959 :
960 : ! **************************************************************************************************
961 : !> \brief Abort if the quasiparticle spectrum handed to the propagator is inverted or has diverged.
962 : !>
963 : !> Tests the fundamental gap per spin channel - not E(HOMO+1) - E(HOMO), since G0W0 reorders levels -
964 : !> on the very array the propagator consumes. Under RTBSE_HAMILTONIAN KS the quasiparticle energies
965 : !> never enter the propagator, so a broken G0W0 spectrum is irrelevant there and does not abort.
966 : !> \param rtbse_env RT-BSE environment with n_ao, n_occ, n_spin, ham_reference_type populated.
967 : !> \param bs_env Bandstructure environment providing the eigenvalues.
968 : ! **************************************************************************************************
969 72 : SUBROUTINE check_qp_gap_sanity(rtbse_env, bs_env)
970 : TYPE(rtbse_env_type), POINTER :: rtbse_env
971 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
972 :
973 : INTEGER :: homo, ispin
974 : REAL(KIND=dp) :: gap, gap_scf
975 :
976 72 : IF (rtbse_env%ham_reference_type /= rtp_bse_ham_gw) RETURN
977 :
978 152 : DO ispin = 1, rtbse_env%n_spin
979 80 : homo = rtbse_env%n_occ(ispin)
980 80 : IF (homo < 1 .OR. homo >= rtbse_env%n_ao) CYCLE
981 :
982 : gap = MINVAL(bs_env%eigenval_GW(homo + 1:rtbse_env%n_ao, 1, ispin)) - &
983 900 : MAXVAL(bs_env%eigenval_GW(1:homo, 1, ispin))
984 : gap_scf = MINVAL(bs_env%eigenval_scf_Gamma(homo + 1:rtbse_env%n_ao, ispin)) - &
985 900 : MAXVAL(bs_env%eigenval_scf_Gamma(1:homo, ispin))
986 :
987 : ! requiring a healthy SCF gap keeps the inversion test from firing on a genuine metal
988 152 : IF (gap < -eps_qp_gap .AND. gap_scf > eps_qp_gap) THEN
989 : CALL cp_abort(__LOCATION__, &
990 : "RTBSE: G0W0 gap of spin "//TRIM(ADJUSTL(cp_to_string(ispin)))// &
991 : " is negative ("//TRIM(ADJUSTL(cp_to_string(gap*evolt, '(F12.3)')))// &
992 : " eV): propagating an inverted spectrum is meaningless. Check the GW "// &
993 0 : "numerical parameters, or use RTBSE_HAMILTONIAN KS.")
994 80 : ELSE IF (ABS(gap) > max_qp_gap) THEN
995 : CALL cp_abort(__LOCATION__, &
996 : "RTBSE: G0W0 gap of spin "//TRIM(ADJUSTL(cp_to_string(ispin)))// &
997 : " is implausibly large ("// &
998 : TRIM(ADJUSTL(cp_to_string(gap*evolt, '(F12.3)')))//" eV): the GW step "// &
999 : "has likely diverged. Check the GW numerical parameters, or use "// &
1000 0 : "RTBSE_HAMILTONIAN KS.")
1001 : END IF
1002 : END DO
1003 :
1004 : END SUBROUTINE check_qp_gap_sanity
1005 :
1006 : ! **************************************************************************************************
1007 : !> \brief Determine the combined active MO window for linearized RT-BSE truncation.
1008 : !>
1009 : !> Evaluates BSE-like cutoff indices per spin from the requested single-particle spectrum
1010 : !> (G0W0 or KS Gamma-point eigenvalues) and collapses them into a single combined window
1011 : !> covering both spin channels by choosing the most inclusive bounds. Issues a CPWARN if the
1012 : !> spin-resolved cutoff candidates differ. When cutoffs are disabled (or the run is not
1013 : !> linearized RT-BSE), the window is set to the full MO range.
1014 : !> \param rtbse_env RT-BSE environment with cutoff values, n_ao, n_occ, n_spin, ham_reference_type
1015 : !> already populated.
1016 : !> \param bs_env Bandstructure environment providing the eigenvalues.
1017 : ! **************************************************************************************************
1018 72 : SUBROUTINE determine_active_mo_window(rtbse_env, bs_env)
1019 : TYPE(rtbse_env_type), POINTER :: rtbse_env
1020 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1021 : CHARACTER(LEN=*), PARAMETER :: routineN = "determine_active_mo_window"
1022 :
1023 : INTEGER :: handle, ispin, n_ao_full, n_virt
1024 : INTEGER :: homo_red, virt_red, homo_incl, virt_incl
1025 : INTEGER :: combined_first_occ, combined_last_virt
1026 : INTEGER :: first_occ_prev, last_virt_prev
1027 : LOGICAL :: spins_differ, do_truncation
1028 : REAL(KIND=dp) :: cutoff_occ, cutoff_empty
1029 :
1030 72 : CALL timeset(routineN, handle)
1031 :
1032 72 : n_ao_full = rtbse_env%n_ao
1033 72 : cutoff_occ = rtbse_env%rtbse_energy_cutoff_occ
1034 72 : cutoff_empty = rtbse_env%rtbse_energy_cutoff_empty
1035 72 : do_truncation = rtbse_env%linearized .AND. (cutoff_occ > 0.0_dp .OR. cutoff_empty > 0.0_dp)
1036 :
1037 : ! Default: full MO window
1038 72 : rtbse_env%first_active_mo = 1
1039 72 : rtbse_env%last_active_mo = n_ao_full
1040 72 : rtbse_env%mo_active = n_ao_full
1041 72 : rtbse_env%active_mo_truncation = .FALSE.
1042 :
1043 72 : IF (.NOT. do_truncation) THEN
1044 70 : CALL timestop(handle)
1045 70 : RETURN
1046 : END IF
1047 :
1048 2 : combined_first_occ = n_ao_full
1049 2 : combined_last_virt = 1
1050 2 : first_occ_prev = -1
1051 2 : last_virt_prev = -1
1052 2 : spins_differ = .FALSE.
1053 :
1054 4 : DO ispin = 1, rtbse_env%n_spin
1055 2 : n_virt = n_ao_full - rtbse_env%n_occ(ispin)
1056 : ! Cut on the DFT axis, as LRBSE does: it is ascending by construction, so the window is a
1057 : ! well-defined contiguous MO range, which is all C_active can extract. The G0W0 axis is
1058 : ! not ordered.
1059 : CALL determine_cutoff_indices(bs_env%eigenval_scf_Gamma(:, ispin), &
1060 : rtbse_env%n_occ(ispin), n_virt, &
1061 : homo_red, virt_red, homo_incl, virt_incl, &
1062 2 : cutoff_occ, cutoff_empty)
1063 : ! Translate the per-spin candidate to global MO indices [homo_incl, homo + virt_incl]
1064 2 : IF (ispin > 1) THEN
1065 0 : IF (homo_incl /= first_occ_prev .OR. (rtbse_env%n_occ(ispin) + virt_incl) /= last_virt_prev) THEN
1066 0 : spins_differ = .TRUE.
1067 : END IF
1068 : END IF
1069 2 : first_occ_prev = homo_incl
1070 2 : last_virt_prev = rtbse_env%n_occ(ispin) + virt_incl
1071 2 : combined_first_occ = MIN(combined_first_occ, homo_incl)
1072 4 : combined_last_virt = MAX(combined_last_virt, rtbse_env%n_occ(ispin) + virt_incl)
1073 : END DO
1074 :
1075 2 : IF (spins_differ) THEN
1076 0 : CPWARN("RTBSE: spin-resolved active MO cutoff candidates differ; using combined window.")
1077 : END IF
1078 :
1079 2 : rtbse_env%first_active_mo = combined_first_occ
1080 2 : rtbse_env%last_active_mo = combined_last_virt
1081 2 : rtbse_env%mo_active = combined_last_virt - combined_first_occ + 1
1082 2 : rtbse_env%active_mo_truncation = (rtbse_env%mo_active < n_ao_full)
1083 :
1084 2 : CALL timestop(handle)
1085 : END SUBROUTINE determine_active_mo_window
1086 :
1087 : ! **************************************************************************************************
1088 : !> \brief Allocates the workspaces for Hartree RI method
1089 : !> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors
1090 : !> \param rtbse_env
1091 : !> \author Stepan Marek
1092 : !> \date 05.2024
1093 : ! **************************************************************************************************
1094 72 : SUBROUTINE create_hartree_ri_workspace(rtbse_env)
1095 : TYPE(rtbse_env_type) :: rtbse_env
1096 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1097 :
1098 : ! Skip the AO-RI Hartree scratch when the RT-BSE Hartree path is fully RI-RS.
1099 : ! In that case rho_dbcsr / v_ao_dbcsr / int_3c_array are never read.
1100 : ! get_sigma_real (AO-RI SX) used to borrow rho_dbcsr as a workspace; that
1101 : ! cross-dependency was removed by giving get_sigma_real its own local
1102 : ! dbcsr scratch (see rt_bse.F::get_sigma_real). rho_dbcsr is now AO-RI
1103 : ! Hartree only, as its name suggests.
1104 72 : IF (rtbse_env%rirs_kernel) RETURN
1105 :
1106 50 : CALL get_qs_env(rtbse_env%qs_env, bs_env=bs_env)
1107 :
1108 50 : CALL dbcsr_create(rtbse_env%rho_dbcsr, name="Sparse density", template=bs_env%mat_ao_ao%matrix)
1109 50 : CALL dbcsr_create(rtbse_env%v_ao_dbcsr, name="Sparse Hartree", template=bs_env%mat_ao_ao%matrix)
1110 :
1111 : CALL create_hartree_ri_3c(rtbse_env%rho_dbcsr, rtbse_env%int_3c_array, rtbse_env%n_ao, rtbse_env%n_RI, &
1112 : bs_env%basis_set_AO, bs_env%basis_set_RI, bs_env%i_RI_start_from_atom, &
1113 50 : bs_env%ri_metric, rtbse_env%qs_env, rtbse_env%unit_nr)
1114 : END SUBROUTINE create_hartree_ri_workspace
1115 : ! **************************************************************************************************
1116 : !> \brief Separated method for allocating the 3c integrals for RI Hartree
1117 : !> \note RI method calculates the Hartree contraction without the use of DBT, as it cannot emulate vectors
1118 : !> \param rho_dbcsr matrix used for the description of shape of 3c array
1119 : !> \param int_3c 3-center integral array to be allocated and filled
1120 : !> \param n_ao Number of atomic orbitals
1121 : !> \param n_RI Number of auxiliary RI orbitals
1122 : !> \param basis_set_AO AO basis set
1123 : !> \param basis_set_RI RI auxiliary basis set
1124 : !> \param i_RI_start_from_atom Array of indices where functions of a given atom in RI basis start
1125 : !> \param unit_nr Unit number used for printing information about the size of int_3c
1126 : !> \author Stepan Marek
1127 : !> \date 01.2025
1128 : ! **************************************************************************************************
1129 50 : SUBROUTINE create_hartree_ri_3c(rho_dbcsr, int_3c, n_ao, n_RI, basis_set_AO, basis_set_RI, &
1130 50 : i_RI_start_from_atom, ri_metric, qs_env, unit_nr)
1131 : TYPE(dbcsr_type) :: rho_dbcsr
1132 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: int_3c
1133 : INTEGER :: n_ao, n_RI
1134 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_set_AO, &
1135 : basis_set_RI
1136 : INTEGER, DIMENSION(:) :: i_RI_start_from_atom
1137 : TYPE(libint_potential_type) :: ri_metric
1138 : TYPE(qs_environment_type), POINTER :: qs_env
1139 : INTEGER :: unit_nr
1140 : REAL(kind=dp) :: size_mb
1141 : INTEGER :: nblkrows_local, &
1142 : nblkcols_local, &
1143 : i_blk_local, &
1144 : j_blk_local, &
1145 : nrows_local, &
1146 : ncols_local, &
1147 : col_local_offset, &
1148 : row_local_offset, &
1149 : start_col_index, &
1150 : end_col_index, &
1151 : start_row_index, &
1152 : end_row_index
1153 50 : INTEGER, DIMENSION(:), POINTER :: local_blk_rows, &
1154 50 : local_blk_cols, &
1155 50 : row_blk_size, &
1156 50 : col_blk_size
1157 : ! TODO : Implement option/decision to not precompute all the 3c integrals
1158 : size_mb = REAL(n_ao, kind=dp)*REAL(n_ao, kind=dp)*REAL(n_RI, kind=dp)* &
1159 50 : REAL(STORAGE_SIZE(size_mb), kind=dp)/8.0_dp/1024.0_dp/1024.0_dp
1160 50 : IF (unit_nr > 0) WRITE (unit_nr, '(A44,E32.2E3,A4)') &
1161 25 : " RTBSE| Approximate size of the 3c integrals", size_mb, " MiB"
1162 :
1163 : ! Get the number of block rows and columns
1164 50 : CALL dbcsr_get_info(rho_dbcsr, nblkrows_local=nblkrows_local, nblkcols_local=nblkcols_local)
1165 : ! Get the global indices of local rows and columns
1166 50 : CALL dbcsr_get_info(rho_dbcsr, local_rows=local_blk_rows, local_cols=local_blk_cols)
1167 : ! Get the sizes of all blocks
1168 50 : CALL dbcsr_get_info(rho_dbcsr, row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1169 :
1170 : ! Get the total required local rows and cols
1171 50 : nrows_local = 0
1172 100 : DO i_blk_local = 1, nblkrows_local
1173 100 : nrows_local = nrows_local + row_blk_size(local_blk_rows(i_blk_local))
1174 : END DO
1175 50 : ncols_local = 0
1176 150 : DO j_blk_local = 1, nblkcols_local
1177 150 : ncols_local = ncols_local + col_blk_size(local_blk_cols(j_blk_local))
1178 : END DO
1179 :
1180 : ! Allocate the appropriate storage
1181 250 : ALLOCATE (int_3c(nrows_local, ncols_local, n_RI))
1182 :
1183 : ! Fill the storage with appropriate values, block by block
1184 50 : row_local_offset = 1
1185 100 : DO i_blk_local = 1, nblkrows_local
1186 : col_local_offset = 1
1187 150 : DO j_blk_local = 1, nblkcols_local
1188 100 : start_row_index = row_local_offset
1189 100 : end_row_index = start_row_index + row_blk_size(local_blk_rows(i_blk_local)) - 1
1190 100 : start_col_index = col_local_offset
1191 100 : end_col_index = start_col_index + col_blk_size(local_blk_cols(j_blk_local)) - 1
1192 : CALL build_3c_integral_block(int_3c(start_row_index:end_row_index, &
1193 : start_col_index:end_col_index, &
1194 : 1:n_RI), &
1195 : qs_env, potential_parameter=ri_metric, &
1196 : basis_j=basis_set_AO, basis_k=basis_set_AO, &
1197 : basis_i=basis_set_RI, &
1198 : atom_j=local_blk_rows(i_blk_local), &
1199 : atom_k=local_blk_cols(j_blk_local), &
1200 100 : i_bf_start_from_atom=i_RI_start_from_atom)
1201 150 : col_local_offset = col_local_offset + col_blk_size(local_blk_cols(j_blk_local))
1202 : END DO
1203 100 : row_local_offset = row_local_offset + row_blk_size(local_blk_rows(i_blk_local))
1204 : END DO
1205 75 : END SUBROUTINE create_hartree_ri_3c
1206 : ! **************************************************************************************************
1207 : !> \brief Releases the workspace for the Hartree RI method
1208 : !> \param rtbse_env RT-BSE Environment, containing specific RI Hartree storage
1209 : !> \author Stepan Marek
1210 : !> \date 09.2024
1211 : ! **************************************************************************************************
1212 72 : SUBROUTINE release_hartree_ri_workspace(rtbse_env)
1213 : TYPE(rtbse_env_type) :: rtbse_env
1214 :
1215 : ! Mirror the gate in create_hartree_ri_workspace and the v_dbcsr gate in
1216 : ! initialize_hartree_potential. With one KERNEL_RI switch the AO-RI Hartree
1217 : ! scratch (3c integrals + dbcsr work + v_dbcsr) is created iff `.NOT. rirs_kernel`.
1218 72 : IF (.NOT. rtbse_env%rirs_kernel) THEN
1219 50 : DEALLOCATE (rtbse_env%int_3c_array)
1220 50 : CALL dbcsr_release(rtbse_env%rho_dbcsr)
1221 50 : CALL dbcsr_release(rtbse_env%v_ao_dbcsr)
1222 50 : CALL dbcsr_release(rtbse_env%v_dbcsr)
1223 : END IF
1224 72 : END SUBROUTINE release_hartree_ri_workspace
1225 : ! **************************************************************************************************
1226 : !> \brief Allocates the workspaces for self-energy determination routine
1227 : !> \param rtbse_env Structure for holding information and workspace structures
1228 : !> \author Stepan Marek
1229 : !> \date 02.2024
1230 : ! **************************************************************************************************
1231 72 : SUBROUTINE create_sigma_workspace(rtbse_env)
1232 : TYPE(rtbse_env_type) :: rtbse_env
1233 :
1234 : ! Skip the AO-RI sigma scratch (W matrix + 3c integrals + work tensors)
1235 : ! when the RT-BSE SEX path is fully RI-RS — these workspaces are read
1236 : ! only by get_sigma in the `.NOT. rirs_kernel` branches.
1237 72 : IF (rtbse_env%rirs_kernel) RETURN
1238 :
1239 : CALL create_sigma_workspace_qs_only(rtbse_env%qs_env, rtbse_env%screened_dbt, rtbse_env%w_dbcsr, &
1240 : rtbse_env%t_3c_w, rtbse_env%t_3c_work_RI_AO__AO, &
1241 50 : rtbse_env%t_3c_work2_RI_AO__AO, rtbse_env%greens_dbt)
1242 : END SUBROUTINE create_sigma_workspace
1243 : ! **************************************************************************************************
1244 : !> \brief Allocates the workspaces for self-energy determination routine
1245 : !> \note Does so without referencing the rtbse_env
1246 : !> \note References bs_env
1247 : !> \param rtbse_env Structure for holding information and workspace structures
1248 : !> \param qs_env Quickstep environment - entry point of calculation
1249 : !> \author Stepan Marek
1250 : !> \date 02.2024
1251 : ! **************************************************************************************************
1252 50 : SUBROUTINE create_sigma_workspace_qs_only(qs_env, screened_dbt, screened_dbcsr, int_3c_dbt, &
1253 : work_dbt_3c_1, work_dbt_3c_2, work_dbt_2c)
1254 : TYPE(qs_environment_type), POINTER :: qs_env
1255 : TYPE(dbcsr_type) :: screened_dbcsr
1256 : TYPE(dbt_type) :: screened_dbt, &
1257 : int_3c_dbt, &
1258 : work_dbt_3c_1, &
1259 : work_dbt_3c_2, &
1260 : work_dbt_2c
1261 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1262 :
1263 50 : CALL get_qs_env(qs_env, bs_env=bs_env)
1264 :
1265 : ! t_3c_w
1266 50 : CALL dbt_create(bs_env%t_RI__AO_AO, int_3c_dbt)
1267 : ! TODO : Provide option/decision whether to store the 3c integrals precomputed
1268 50 : CALL compute_3c_integrals(qs_env, bs_env, int_3c_dbt)
1269 : ! t_3c_work_RI_AO__AO
1270 50 : CALL dbt_create(bs_env%t_RI_AO__AO, work_dbt_3c_1)
1271 : ! t_3c_work2_RI_AO__AO
1272 50 : CALL dbt_create(bs_env%t_RI_AO__AO, work_dbt_3c_2)
1273 : ! t_W
1274 : ! Populate screened_dbt from gw run
1275 50 : CALL dbcsr_create(screened_dbcsr, name="W", template=bs_env%mat_RI_RI%matrix)
1276 50 : CALL dbt_create(screened_dbcsr, screened_dbt)
1277 : ! greens_dbt
1278 50 : CALL dbt_create(bs_env%mat_ao_ao%matrix, work_dbt_2c)
1279 50 : END SUBROUTINE create_sigma_workspace_qs_only
1280 : ! **************************************************************************************************
1281 : !> \brief Releases the workspaces for self-energy determination
1282 : !> \param rtbse_env
1283 : !> \author Stepan Marek
1284 : !> \date 02.2024
1285 : ! **************************************************************************************************
1286 72 : SUBROUTINE release_sigma_workspace(rtbse_env)
1287 : TYPE(rtbse_env_type) :: rtbse_env
1288 :
1289 : ! Mirror the gate in create_sigma_workspace.
1290 72 : IF (rtbse_env%rirs_kernel) RETURN
1291 :
1292 50 : CALL dbt_destroy(rtbse_env%t_3c_w)
1293 50 : CALL dbt_destroy(rtbse_env%t_3c_work_RI_AO__AO)
1294 50 : CALL dbt_destroy(rtbse_env%t_3c_work2_RI_AO__AO)
1295 50 : CALL dbt_destroy(rtbse_env%screened_dbt)
1296 50 : CALL dbt_destroy(rtbse_env%greens_dbt)
1297 50 : CALL dbcsr_release(rtbse_env%w_dbcsr)
1298 : END SUBROUTINE release_sigma_workspace
1299 : ! **************************************************************************************************
1300 : !> \brief Multiplies real matrix by a complex matrix from the right
1301 : !> \note So far only converts the real matrix to complex one, potentially doubling the work
1302 : !> \param rtbse_env
1303 : !> \author Stepan Marek
1304 : !> \date 09.2024
1305 : ! **************************************************************************************************
1306 15064 : SUBROUTINE multiply_fm_cfm(trans_r, trans_c, na, nb, nc, &
1307 : alpha, matrix_r, matrix_c, beta, res)
1308 : ! Transposition
1309 : CHARACTER(len=1) :: trans_r, trans_c
1310 : INTEGER :: na, nb, nc
1311 : ! accept real numbers
1312 : ! TODO : Just use complex numbers and import z_one, z_zero etc.
1313 : REAL(kind=dp) :: alpha, beta
1314 : TYPE(cp_fm_type) :: matrix_r
1315 : TYPE(cp_cfm_type) :: matrix_c, res
1316 : TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im
1317 : REAL(kind=dp) :: i_unit
1318 : CHARACTER(len=1) :: trans_cr
1319 :
1320 3766 : CALL cp_fm_create(work_re, matrix_c%matrix_struct)
1321 3766 : CALL cp_fm_create(work_im, matrix_c%matrix_struct)
1322 3766 : CALL cp_fm_create(res_re, res%matrix_struct)
1323 3766 : CALL cp_fm_create(res_im, res%matrix_struct)
1324 3766 : CALL cp_cfm_to_fm(matrix_c, work_re, work_im)
1325 0 : SELECT CASE (trans_c)
1326 : CASE ("C")
1327 0 : i_unit = -1.0_dp
1328 0 : trans_cr = "T"
1329 : CASE ("T")
1330 0 : i_unit = 1.0_dp
1331 0 : trans_cr = "T"
1332 : CASE default
1333 3766 : i_unit = 1.0_dp
1334 3766 : trans_cr = "N"
1335 : END SELECT
1336 : ! Actual multiplication
1337 : CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, &
1338 3766 : alpha, matrix_r, work_re, beta, res_re)
1339 : CALL parallel_gemm(trans_r, trans_cr, na, nb, nc, &
1340 3766 : i_unit*alpha, matrix_r, work_im, beta, res_im)
1341 3766 : CALL cp_fm_to_cfm(res_re, res_im, res)
1342 3766 : CALL cp_fm_release(work_re)
1343 3766 : CALL cp_fm_release(work_im)
1344 3766 : CALL cp_fm_release(res_re)
1345 3766 : CALL cp_fm_release(res_im)
1346 :
1347 3766 : END SUBROUTINE multiply_fm_cfm
1348 : ! **************************************************************************************************
1349 : !> \brief Multiplies complex matrix by a real matrix from the right
1350 : !> \note So far only converts the real matrix to complex one, potentially doubling the work
1351 : !> \param rtbse_env
1352 : !> \author Stepan Marek
1353 : !> \date 09.2024
1354 : ! **************************************************************************************************
1355 5424 : SUBROUTINE multiply_cfm_fm(trans_c, trans_r, na, nb, nc, &
1356 : alpha, matrix_c, matrix_r, beta, res)
1357 : ! Transposition
1358 : CHARACTER(len=1) :: trans_c, trans_r
1359 : INTEGER :: na, nb, nc
1360 : ! accept real numbers
1361 : ! TODO : complex number support via interface?
1362 : REAL(kind=dp) :: alpha, beta
1363 : TYPE(cp_cfm_type) :: matrix_c, res
1364 : TYPE(cp_fm_type) :: matrix_r
1365 : TYPE(cp_fm_type) :: work_re, work_im, res_re, res_im
1366 : REAL(kind=dp) :: i_unit
1367 : CHARACTER(len=1) :: trans_cr
1368 :
1369 1356 : CALL cp_fm_create(work_re, matrix_c%matrix_struct)
1370 1356 : CALL cp_fm_create(work_im, matrix_c%matrix_struct)
1371 1356 : CALL cp_fm_create(res_re, res%matrix_struct)
1372 1356 : CALL cp_fm_create(res_im, res%matrix_struct)
1373 1356 : CALL cp_cfm_to_fm(matrix_c, work_re, work_im)
1374 0 : SELECT CASE (trans_c)
1375 : CASE ("C")
1376 0 : i_unit = -1.0_dp
1377 0 : trans_cr = "T"
1378 : CASE ("T")
1379 0 : i_unit = 1.0_dp
1380 0 : trans_cr = "T"
1381 : CASE default
1382 1356 : i_unit = 1.0_dp
1383 1356 : trans_cr = "N"
1384 : END SELECT
1385 : ! Actual multiplication
1386 : CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, &
1387 1356 : alpha, work_re, matrix_r, beta, res_re)
1388 : CALL parallel_gemm(trans_cr, trans_r, na, nb, nc, &
1389 1356 : i_unit*alpha, work_im, matrix_r, beta, res_im)
1390 1356 : CALL cp_fm_to_cfm(res_re, res_im, res)
1391 1356 : CALL cp_fm_release(work_re)
1392 1356 : CALL cp_fm_release(work_im)
1393 1356 : CALL cp_fm_release(res_re)
1394 1356 : CALL cp_fm_release(res_im)
1395 :
1396 1356 : END SUBROUTINE multiply_cfm_fm
1397 0 : END MODULE rt_bse_types
|
