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 RT-BSE RI-RS kernels: SEX and Hartree evaluated by collocation on grid points r_l.
10 : !> Once-built grid objects:
11 : !> φ_µ(r_l) [grid×AO], Z_lP [grid×RI],
12 : !> V^aux_PQ = [M^-1 V^tr M^-1]_PQ [RI×RI], W^0_ll' = sum_PQ Z_lP (V + W^c(ω=0))_PQ Z_l'Q.
13 : !> Per-call kernels, collocation X = φ_µ(r_l) (AO domain, written φ_lµ below):
14 : !> SEX: ρ^grid_ll' = sum_µν φ_lµ Δρ_µν φ_l'ν ;
15 : !> Σ_µν = pref * sum_ll' φ_lµ [ρ^grid ∘ W^0]_ll' φ_l'ν
16 : !> Hartree: n_l = sum_µν φ_lµ Δρ_µν φ_lν ;
17 : !> v_l = sum_PQl' Z_lP V^aux_PQ Z_l'Q n_l' (applied factorized, stage 2) ;
18 : !> V^H_µν = sum_l φ_lµ v_l φ_lν (diagonal-only, no grid×grid)
19 : !> Independent of which GW variant produced bs_env%fm_W_MIC_freq_zero.
20 : !> \author Maximilian Graml (05.26)
21 : ! **************************************************************************************************
22 : MODULE rt_bse_ri_rs
23 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add
24 : USE cp_cfm_types, ONLY: cp_cfm_create,&
25 : cp_cfm_release,&
26 : cp_cfm_to_fm,&
27 : cp_cfm_type,&
28 : cp_fm_to_cfm
29 : USE cp_dbcsr_api, ONLY: &
30 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, &
31 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
32 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_release, dbcsr_scale, &
33 : dbcsr_type, dbcsr_type_no_symmetry
34 : USE cp_dbcsr_contrib, ONLY: dbcsr_get_diag
35 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
36 : copy_fm_to_dbcsr
37 : USE cp_fm_types, ONLY: cp_fm_create,&
38 : cp_fm_release,&
39 : cp_fm_type
40 : USE gw_large_cell_Gamma_ri_rs, ONLY: contract_A_B_A,&
41 : hadamard_product_inplace,&
42 : release_dbcsr_topology_and_matrices,&
43 : setup_square_topology
44 : USE gw_large_cell_gamma, ONLY: multiply_fm_W_MIC_time_with_Minv_Gamma
45 : USE gw_non_periodic_ri_rs, ONLY: atomic_basis_at_grid_point,&
46 : compute_coeff_Z_lP,&
47 : precompute_ri_rs_radii,&
48 : ri_rs_grid_assembler
49 : USE kinds, ONLY: dp
50 : USE machine, ONLY: m_walltime
51 : USE message_passing, ONLY: mp_para_env_type
52 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
53 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
54 : USE qs_environment_types, ONLY: qs_environment_type
55 : #include "../base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_bse_ri_rs'
62 :
63 : PUBLIC :: rt_bse_ri_rs_ensure_grid, &
64 : rt_bse_ri_rs_ensure_V_grid, &
65 : rt_bse_ri_rs_ensure_W0_grid, &
66 : compute_sigma_ri_rs, &
67 : compute_sigma_ri_rs_complex, &
68 : compute_hartree_ri_rs, &
69 : compute_hartree_ri_rs_complex, &
70 : compute_hartree_ri_rs_from_diag, &
71 : hartree_potential_from_diag_ri_rs
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief Make sure the AO collocation φ_µ(r_l) (mat_phi_mu_l) and the RI fit coefficients
77 : !> Z_lP (mat_Z_lP) are populated in memory.
78 : !> If GW was run with RI-RS the grid is already built; otherwise build it here so the
79 : !> AO-RI GW + RI-RS RT-BSE combination is possible.
80 : !> \param bs_env ...
81 : !> \param qs_env ...
82 : ! **************************************************************************************************
83 32 : SUBROUTINE rt_bse_ri_rs_ensure_grid(bs_env, qs_env)
84 :
85 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
86 : TYPE(qs_environment_type), POINTER :: qs_env
87 :
88 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_bse_ri_rs_ensure_grid'
89 :
90 : INTEGER :: handle
91 : REAL(KIND=dp) :: t1
92 :
93 32 : CALL timeset(routineN, handle)
94 :
95 32 : IF (bs_env%ri_rs%grid_built) THEN
96 24 : CALL timestop(handle)
97 24 : RETURN
98 : END IF
99 :
100 8 : t1 = m_walltime()
101 :
102 8 : CALL ri_rs_grid_assembler(qs_env, bs_env, bs_env%ri_rs%grid_points)
103 : ! Per-atom AO/RI screening radii required by the screened grid-fill and d_lP routines; the
104 : ! GW RI-RS driver populates these, but the standalone RT-BSE grid build must do so itself.
105 8 : IF (.NOT. ALLOCATED(bs_env%ri_rs%radius_ao_per_atom)) THEN
106 8 : CALL precompute_ri_rs_radii(qs_env, bs_env)
107 : END IF
108 : CALL atomic_basis_at_grid_point(qs_env, bs_env, bs_env%ri_rs%grid_points, &
109 8 : bs_env%ri_rs%mat_phi_mu_l)
110 : CALL compute_coeff_Z_lP(qs_env, bs_env, bs_env%ri_rs%grid_points, &
111 8 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
112 :
113 8 : bs_env%ri_rs%grid_built = .TRUE.
114 :
115 8 : IF (bs_env%unit_nr > 0) THEN
116 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
117 4 : 'Built RI-RS grid for RT-BSE (no GW_RI_RS used),', ' Execution time', &
118 8 : m_walltime() - t1, ' s'
119 4 : WRITE (bs_env%unit_nr, '(A)') ' '
120 : END IF
121 :
122 8 : CALL timestop(handle)
123 :
124 : END SUBROUTINE rt_bse_ri_rs_ensure_grid
125 :
126 : ! **************************************************************************************************
127 : !> \brief Build V^aux_PQ = [M^-1 V^tr M^-1]_PQ (truncated Coulomb in the RI basis, M^-1-sandwiched
128 : !> to match the W^MIC convention). The grid Coulomb V_ll' = sum_PQ Z_lP V^aux_PQ Z_l'Q is
129 : !> never materialized -- the Hartree stage applies it factorized:
130 : !> v_l = sum_PQl' Z_lP V^aux_PQ Z_l'Q n_l'
131 : !> Needed by RT-BSE RI-RS Hartree.
132 : !> \param bs_env ...
133 : !> \param qs_env ...
134 : ! **************************************************************************************************
135 16 : SUBROUTINE rt_bse_ri_rs_ensure_V_grid(bs_env, qs_env)
136 :
137 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
138 : TYPE(qs_environment_type), POINTER :: qs_env
139 :
140 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_bse_ri_rs_ensure_V_grid'
141 :
142 : INTEGER :: handle
143 16 : INTEGER, DIMENSION(:), POINTER :: blk_aux, dist_row_aux
144 : REAL(KIND=dp) :: t1
145 16 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
146 : TYPE(dbcsr_distribution_type) :: dist_aux_aux
147 :
148 16 : CALL timeset(routineN, handle)
149 :
150 16 : IF (bs_env%ri_rs%V_grid_built) THEN
151 0 : CALL timestop(handle)
152 0 : RETURN
153 : END IF
154 :
155 16 : CALL rt_bse_ri_rs_ensure_grid(bs_env, qs_env)
156 16 : t1 = m_walltime()
157 :
158 16 : CALL setup_square_topology(bs_env%ri_rs%mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
159 :
160 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
161 16 : bs_env%trunc_coulomb, do_kpoints=.FALSE.)
162 : ! Apply M^-1 sandwich to match the W^MIC convention; same scale as W^c when both used.
163 16 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_Vtr_Gamma(:, 1))
164 :
165 : ! Store the M^-1-sandwiched RI-basis Coulomb; the grid kernel Z V Z^T is applied factorized.
166 : CALL dbcsr_create(bs_env%ri_rs%mat_V_aux_rtbse, "V_aux_rtbse", dist_aux_aux, &
167 16 : dbcsr_type_no_symmetry, blk_aux, blk_aux)
168 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), bs_env%ri_rs%mat_V_aux_rtbse, &
169 16 : keep_sparsity=.FALSE.)
170 :
171 16 : bs_env%ri_rs%V_grid_built = .TRUE.
172 :
173 16 : IF (bs_env%unit_nr > 0) THEN
174 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
175 8 : 'Precomputed RT-BSE RI-RS V_aux kernel,', ' Execution time', &
176 16 : m_walltime() - t1, ' s'
177 8 : WRITE (bs_env%unit_nr, '(A)') ' '
178 : END IF
179 :
180 16 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
181 16 : CALL cp_fm_release(fm_Vtr_Gamma)
182 :
183 16 : CALL timestop(handle)
184 :
185 32 : END SUBROUTINE rt_bse_ri_rs_ensure_V_grid
186 :
187 : ! **************************************************************************************************
188 : !> \brief Build W^0_ll' = sum_PQ Z_lP (V + W^c(ω=0))_PQ Z_l'Q (statically screened W on the grid).
189 : !> Needed by RT-BSE RI-RS SEX/COH. Reuses bs_env%fm_W_MIC_freq_zero which must already
190 : !> contain M^-1 W^c(ω=0) M^-1 (built by either GW path under the BSE rtp_method gate).
191 : !> W^0 enters only through the Hadamard ρ^grid ∘ W^0 -- its elements are needed, so it is
192 : !> the one persistent grid×grid object of the kernel layer.
193 : !> \param bs_env ...
194 : !> \param qs_env ...
195 : ! **************************************************************************************************
196 16 : SUBROUTINE rt_bse_ri_rs_ensure_W0_grid(bs_env, qs_env)
197 :
198 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
199 : TYPE(qs_environment_type), POINTER :: qs_env
200 :
201 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_bse_ri_rs_ensure_W0_grid'
202 :
203 : INTEGER :: handle
204 16 : INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
205 16 : dist_row_aux
206 : REAL(KIND=dp) :: t1
207 16 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
208 : TYPE(dbcsr_distribution_type) :: dist_aux_aux, dist_grid_grid
209 : TYPE(dbcsr_type) :: matrix_V_aux, matrix_W_aux
210 :
211 16 : CALL timeset(routineN, handle)
212 :
213 16 : IF (bs_env%ri_rs%W0_grid_built) THEN
214 0 : CALL timestop(handle)
215 0 : RETURN
216 : END IF
217 :
218 : ! W(w=0) is built by the GW step only under its RTBSE rtp_method gate; reaching this consumer
219 : ! without it means gate and consumer disagree. Abort rather than read a never-created cp_fm.
220 16 : IF (.NOT. ASSOCIATED(bs_env%fm_W_MIC_freq_zero%matrix_struct)) THEN
221 : CALL cp_abort(__LOCATION__, &
222 : "RT-BSE RI-RS kernel needs the screened interaction W(w=0), which the GW "// &
223 : "step did not build. Select the RT-BSE propagator with '&RTBSE' or "// &
224 0 : "'&RTBSE RTBSE', not '&RTBSE TDDFT'.")
225 : END IF
226 :
227 16 : CALL rt_bse_ri_rs_ensure_grid(bs_env, qs_env)
228 16 : t1 = m_walltime()
229 :
230 : CALL setup_square_topology(bs_env%ri_rs%mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, &
231 16 : dist_col_grid)
232 16 : CALL setup_square_topology(bs_env%ri_rs%mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
233 :
234 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
235 16 : bs_env%trunc_coulomb, do_kpoints=.FALSE.)
236 16 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_Vtr_Gamma(:, 1))
237 :
238 : CALL dbcsr_create(matrix_V_aux, "V_aux_rtbse_W0", dist_aux_aux, dbcsr_type_no_symmetry, &
239 16 : blk_aux, blk_aux)
240 16 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.FALSE.)
241 :
242 : CALL dbcsr_create(matrix_W_aux, "W_aux_rtbse", dist_aux_aux, dbcsr_type_no_symmetry, &
243 16 : blk_aux, blk_aux)
244 16 : CALL copy_fm_to_dbcsr(bs_env%fm_W_MIC_freq_zero, matrix_W_aux, keep_sparsity=.FALSE.)
245 16 : CALL dbcsr_add(matrix_W_aux, matrix_V_aux, 1.0_dp, 1.0_dp)
246 :
247 : CALL dbcsr_create(bs_env%ri_rs%mat_W0_grid_rtbse, "W0_grid_rtbse", dist_grid_grid, &
248 16 : dbcsr_type_no_symmetry, blk_grid, blk_grid)
249 : CALL contract_A_B_A("N", "T", bs_env%ri_rs%mat_Z_lP, matrix_W_aux, &
250 16 : bs_env%ri_rs%mat_W0_grid_rtbse, bs_env%eps_filter)
251 :
252 16 : bs_env%ri_rs%W0_grid_built = .TRUE.
253 16 : bs_env%ri_rs%rtbse_kernels_ready = .TRUE.
254 :
255 16 : IF (bs_env%unit_nr > 0) THEN
256 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
257 8 : 'Precomputed RT-BSE RI-RS W0_grid kernel,', ' Execution time', &
258 16 : m_walltime() - t1, ' s'
259 8 : WRITE (bs_env%unit_nr, '(A)') ' '
260 : END IF
261 :
262 16 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid)
263 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux, &
264 16 : m1=matrix_V_aux, m2=matrix_W_aux)
265 16 : CALL cp_fm_release(fm_Vtr_Gamma)
266 :
267 16 : CALL timestop(handle)
268 :
269 32 : END SUBROUTINE rt_bse_ri_rs_ensure_W0_grid
270 :
271 : ! **************************************************************************************************
272 : !> \brief AO-domain SEX: Σ_µν = pref * sum_ll' φ_lµ [ρ^grid ∘ W^0]_ll' φ_l'ν,
273 : !> ρ^grid_ll' = sum_µν φ_lµ Δρ_µν φ_l'ν.
274 : !> The grid×grid ρ^grid is intrinsic to SEX -- the Hadamard needs W^0's elements, so no
275 : !> factorized application exists (unlike the Hartree V_ll'). Real input, real output;
276 : !> used for COH (input S^-1) and the init reference (ρ^0); dynamic Δρ goes through the
277 : !> complex variant. Mirrors the AO-RI get_sigma(rtbse_env, sigma_fm, prefactor, rho_fm) API.
278 : !> \param bs_env ...
279 : !> \param sigma_AO_fm result, AO x AO
280 : !> \param prefactor scaling applied to the final result
281 : !> \param rho_AO_fm input density-like matrix, AO x AO
282 : !> \param grid_diag_accum ...
283 : ! **************************************************************************************************
284 3504 : SUBROUTINE compute_sigma_ri_rs(bs_env, sigma_AO_fm, prefactor, rho_AO_fm, &
285 3504 : grid_diag_accum)
286 :
287 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
288 : TYPE(cp_fm_type), INTENT(INOUT) :: sigma_AO_fm
289 : REAL(KIND=dp), INTENT(IN) :: prefactor
290 : TYPE(cp_fm_type), INTENT(IN) :: rho_AO_fm
291 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: grid_diag_accum(:)
292 :
293 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_sigma_ri_rs'
294 :
295 : INTEGER :: handle, n_grid
296 3504 : INTEGER, DIMENSION(:), POINTER :: blk_ao, blk_grid, dist_col_grid, &
297 3504 : dist_row_ao
298 3504 : REAL(KIND=dp), ALLOCATABLE :: diag_local(:)
299 : TYPE(dbcsr_distribution_type) :: dist_ao_ao, dist_grid_grid
300 : TYPE(dbcsr_type) :: matrix_rho_AO, matrix_rho_grid, &
301 : matrix_Sigma_AO
302 :
303 3504 : CALL timeset(routineN, handle)
304 :
305 3504 : CPASSERT(bs_env%ri_rs%W0_grid_built)
306 :
307 3504 : CALL setup_square_topology(bs_env%ri_rs%mat_phi_mu_l, 'COL', dist_ao_ao, blk_ao, dist_row_ao)
308 : CALL setup_square_topology(bs_env%ri_rs%mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, &
309 3504 : dist_col_grid)
310 :
311 : CALL dbcsr_create(matrix_rho_AO, "rho_AO_ri_rs", dist_ao_ao, dbcsr_type_no_symmetry, &
312 3504 : blk_ao, blk_ao)
313 3504 : CALL copy_fm_to_dbcsr(rho_AO_fm, matrix_rho_AO, keep_sparsity=.FALSE.)
314 :
315 : ! ρ^grid_ll' = sum_µν φ_lµ Δρ_µν φ_l'ν (the SEX-intrinsic grid×grid transient)
316 : CALL dbcsr_create(matrix_rho_grid, "rho_grid_ri_rs", dist_grid_grid, &
317 3504 : dbcsr_type_no_symmetry, blk_grid, blk_grid)
318 : ! No eps_filter on the SEX projections: a Δρ-dependent block drop breaks kernel
319 : ! self-adjointness (L non-Hermitian). The W0 Hadamard below supplies the symmetric,
320 : ! Δρ-independent grid sparsity, so unfiltered projections stay sparse via W0.
321 : ! PERF(rirs-selfadjoint): a retain_sparsity-to-W0 build (dbcsr_multiply) gives identical
322 : ! numbers with W0-sparse memory; NOT YET APPLIED, needs a perf test.
323 : CALL contract_A_B_A("N", "T", bs_env%ri_rs%mat_phi_mu_l, matrix_rho_AO, &
324 3504 : matrix_rho_grid, 0.0_dp)
325 :
326 : ! Harvest n_l = diag(ρ^grid) BEFORE the in-place Hadamard destroys it; the Hartree reuses it
327 : ! (compute_hartree_ri_rs_from_diag) instead of rebuilding φρφ^T. Bare accumulate (caller
328 : ! pre-zeroes): the cross-spin Hartree density is the spin SUM of these, and spin_degeneracy is
329 : ! applied to the V_H OUTPUT (post-filter, bit-identical) -- never to the diagonal, which would
330 : ! shift the stage-3 eps_filter cut (coarse-filter sensitivity).
331 3504 : IF (PRESENT(grid_diag_accum)) THEN
332 2240 : n_grid = SIZE(grid_diag_accum)
333 6720 : ALLOCATE (diag_local(n_grid))
334 2240 : diag_local = 0.0_dp
335 2240 : CALL dbcsr_get_diag(matrix_rho_grid, diag_local)
336 2240 : CALL bs_env%para_env%sum(diag_local)
337 671040 : grid_diag_accum(:) = grid_diag_accum(:) + diag_local(:)
338 2240 : DEALLOCATE (diag_local)
339 : END IF
340 :
341 : ! ρ^grid_ll' <- ρ^grid_ll' * W^0_ll' (in place; blocks without a W^0 partner zeroed)
342 : CALL hadamard_product_inplace(matrix_rho_grid, bs_env%ri_rs%mat_W0_grid_rtbse, &
343 3504 : 1.0_dp)
344 :
345 : ! Σ_µν = sum_ll' φ_lµ [ρ^grid ∘ W^0]_ll' φ_l'ν
346 3504 : CALL dbcsr_create(matrix_Sigma_AO, template=matrix_rho_AO)
347 : ! Unfiltered for self-adjointness (see the forward projection above; PERF note there).
348 : CALL contract_A_B_A("T", "N", bs_env%ri_rs%mat_phi_mu_l, matrix_rho_grid, &
349 3504 : matrix_Sigma_AO, 0.0_dp)
350 :
351 3504 : CALL dbcsr_scale(matrix_Sigma_AO, prefactor)
352 3504 : CALL copy_dbcsr_to_fm(matrix_Sigma_AO, sigma_AO_fm)
353 :
354 : CALL release_dbcsr_topology_and_matrices(dist=dist_ao_ao, mapped_dist=dist_row_ao, &
355 3504 : m1=matrix_rho_AO, m2=matrix_Sigma_AO)
356 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid, &
357 3504 : m1=matrix_rho_grid)
358 :
359 3504 : CALL timestop(handle)
360 :
361 7008 : END SUBROUTINE compute_sigma_ri_rs
362 :
363 : ! **************************************************************************************************
364 : !> \brief Complex-input AO SEX via Re/Im split: the kernel is real, so complex linearity holds as
365 : !> Σ[Δρ] = Σ[Re Δρ] + i Σ[Im Δρ]. Required for non-Hermitian Δρ inputs
366 : !> (TDA OV-only / ABBA OV+VO).
367 : !> \param bs_env ...
368 : !> \param sigma_AO_cfm result, AO x AO (complex)
369 : !> \param prefactor scaling applied to the final result
370 : !> \param rho_AO_cfm input AO x AO complex matrix
371 : !> \param grid_diag_re_accum optional: accumulate diag(φ.Re(ρ).φ^T) (bare; for the Hartree reuse)
372 : !> \param grid_diag_im_accum optional: accumulate diag(φ.Im(ρ).φ^T) (bare; for the Hartree reuse)
373 : ! **************************************************************************************************
374 8760 : SUBROUTINE compute_sigma_ri_rs_complex(bs_env, sigma_AO_cfm, prefactor, rho_AO_cfm, &
375 1752 : grid_diag_re_accum, grid_diag_im_accum)
376 :
377 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
378 : TYPE(cp_cfm_type), INTENT(INOUT) :: sigma_AO_cfm
379 : REAL(KIND=dp), INTENT(IN) :: prefactor
380 : TYPE(cp_cfm_type), INTENT(IN) :: rho_AO_cfm
381 : REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: grid_diag_re_accum(:), &
382 : grid_diag_im_accum(:)
383 :
384 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_sigma_ri_rs_complex'
385 :
386 : INTEGER :: handle
387 : TYPE(cp_cfm_type) :: cfm_real_part
388 : TYPE(cp_fm_type) :: fm_rho, fm_sigma
389 :
390 1752 : CALL timeset(routineN, handle)
391 :
392 1752 : CALL cp_fm_create(fm_rho, bs_env%fm_s_Gamma%matrix_struct)
393 1752 : CALL cp_fm_create(fm_sigma, bs_env%fm_s_Gamma%matrix_struct)
394 1752 : CALL cp_cfm_create(cfm_real_part, bs_env%fm_s_Gamma%matrix_struct)
395 :
396 : ! Re/Im each harvest into their own accumulator; absent optionals propagate as absent.
397 1752 : CALL cp_cfm_to_fm(msource=rho_AO_cfm, mtargetr=fm_rho)
398 : CALL compute_sigma_ri_rs(bs_env, fm_sigma, prefactor, fm_rho, &
399 2064 : grid_diag_accum=grid_diag_re_accum)
400 1752 : CALL cp_fm_to_cfm(msourcer=fm_sigma, mtarget=cfm_real_part)
401 :
402 1752 : CALL cp_cfm_to_fm(msource=rho_AO_cfm, mtargeti=fm_rho)
403 : CALL compute_sigma_ri_rs(bs_env, fm_sigma, prefactor, fm_rho, &
404 2704 : grid_diag_accum=grid_diag_im_accum)
405 1752 : CALL cp_fm_to_cfm(msourcei=fm_sigma, mtarget=sigma_AO_cfm)
406 :
407 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), sigma_AO_cfm, &
408 1752 : CMPLX(1.0_dp, 0.0_dp, kind=dp), cfm_real_part)
409 :
410 1752 : CALL cp_fm_release(fm_rho)
411 1752 : CALL cp_fm_release(fm_sigma)
412 1752 : CALL cp_cfm_release(cfm_real_part)
413 :
414 1752 : CALL timestop(handle)
415 :
416 1752 : END SUBROUTINE compute_sigma_ri_rs_complex
417 :
418 : ! **************************************************************************************************
419 : !> \brief AO-domain Hartree via RI-RS:
420 : !> n_l = sum_µν φ_lµ Δρ_µν φ_lν = (φ ρ φ^T)_ll (diagonal of materialized grid×grid) ;
421 : !> v_l = sum_PQl' Z_lP V^aux_PQ Z_l'Q n_l' (factorized, stage 2) ;
422 : !> V^H_µν = sum_l φ_lµ v_l φ_lν (diagonal-only row-scale, stage 3)
423 : !> Real input, real output; complex inputs go through compute_hartree_ri_rs_complex.
424 : !> \param bs_env ...
425 : !> \param rho_AO_fm input AO x AO density matrix
426 : !> \param V_H_AO_fm output AO x AO Hartree potential
427 : ! **************************************************************************************************
428 444 : SUBROUTINE compute_hartree_ri_rs(bs_env, rho_AO_fm, V_H_AO_fm)
429 :
430 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
431 : TYPE(cp_fm_type), INTENT(IN) :: rho_AO_fm
432 : TYPE(cp_fm_type), INTENT(INOUT) :: V_H_AO_fm
433 :
434 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_hartree_ri_rs'
435 :
436 : INTEGER :: handle, n_grid
437 444 : INTEGER, DIMENSION(:), POINTER :: blk_ao, blk_grid, dist_col_grid, &
438 444 : dist_row_ao
439 : REAL(KIND=dp), ALLOCATABLE :: n_vec(:)
440 : TYPE(dbcsr_distribution_type) :: dist_ao_ao, dist_grid_grid
441 : TYPE(dbcsr_type) :: matrix_rho_AO, matrix_rho_grid
442 :
443 444 : CALL timeset(routineN, handle)
444 :
445 444 : CPASSERT(bs_env%ri_rs%V_grid_built)
446 :
447 444 : CALL setup_square_topology(bs_env%ri_rs%mat_phi_mu_l, 'COL', dist_ao_ao, blk_ao, dist_row_ao)
448 : CALL setup_square_topology(bs_env%ri_rs%mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, &
449 444 : dist_col_grid)
450 :
451 : CALL dbcsr_create(matrix_rho_AO, "rho_AO_hartree", dist_ao_ao, dbcsr_type_no_symmetry, &
452 444 : blk_ao, blk_ao)
453 444 : CALL copy_fm_to_dbcsr(rho_AO_fm, matrix_rho_AO, keep_sparsity=.FALSE.)
454 :
455 5236 : n_grid = SUM(blk_grid)
456 1332 : ALLOCATE (n_vec(n_grid))
457 :
458 : ! stage 1: n_l = (φ ρ φ^T)_ll. A grid×AO row-dot (iterate φρ, look up φ) silently drops
459 : ! off-rank pairs at ≥2 ranks — φ is sparse and the product is not co-located with φ — so
460 : ! materialize φρφ^T and take its diagonal, co-location-safe like the SEX grid kernel.
461 : CALL dbcsr_create(matrix_rho_grid, "rho_grid_hartree", dist_grid_grid, &
462 444 : dbcsr_type_no_symmetry, blk_grid, blk_grid)
463 : ! No eps_filter: a Δρ-dependent diagonal-block drop zeroes n_l inconsistently between OV
464 : ! pairs and breaks Hartree self-adjointness (L non-Hermitian).
465 : ! PERF(rirs-selfadjoint): this full grid×grid is built only for its diagonal n_l; a
466 : ! diagonal-only build would avoid it. NOT YET APPLIED, needs a perf test.
467 : CALL contract_A_B_A("N", "T", bs_env%ri_rs%mat_phi_mu_l, matrix_rho_AO, &
468 444 : matrix_rho_grid, 0.0_dp)
469 444 : n_vec = 0.0_dp
470 444 : CALL dbcsr_get_diag(matrix_rho_grid, n_vec)
471 444 : CALL bs_env%para_env%sum(n_vec)
472 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid, &
473 444 : m1=matrix_rho_grid)
474 : CALL release_dbcsr_topology_and_matrices(dist=dist_ao_ao, mapped_dist=dist_row_ao, &
475 444 : m1=matrix_rho_AO)
476 :
477 : ! stages 2-3: factorized Coulomb v = Z V^aux Z^T n, then V^H = φ^T diag(v) φ.
478 444 : CALL hartree_potential_from_diag_ri_rs(bs_env, n_vec, V_H_AO_fm)
479 :
480 444 : DEALLOCATE (n_vec)
481 :
482 444 : CALL timestop(handle)
483 :
484 888 : END SUBROUTINE compute_hartree_ri_rs
485 :
486 : ! **************************************************************************************************
487 : !> \brief Hartree stages 2-3 from a precomputed grid density n_l (skips the stage-1 φρφ^T build):
488 : !> v_l = sum_PQl' Z_lP V^aux_PQ Z_l'Q n_l' (factorized Coulomb) ;
489 : !> V^H_µν = sum_l φ_lµ v_l φ_lν (Φ_lν = v_l φ_lν rowscale, then V^H = φ^T Φ).
490 : !> n_l is harvested as diag(φρφ^T) inside compute_sigma_ri_rs (the SEX grid kernel), so the
491 : !> Hartree never rebuilds the grid×grid product. Real in/out.
492 : !> \param bs_env ...
493 : !> \param n_vec grid density n_l (length n_grid, replicated)
494 : !> \param V_H_AO_fm output AO x AO Hartree potential
495 : ! **************************************************************************************************
496 2204 : SUBROUTINE hartree_potential_from_diag_ri_rs(bs_env, n_vec, V_H_AO_fm)
497 :
498 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
499 : REAL(KIND=dp), INTENT(IN) :: n_vec(:)
500 : TYPE(cp_fm_type), INTENT(INOUT) :: V_H_AO_fm
501 :
502 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hartree_potential_from_diag_ri_rs'
503 :
504 : INTEGER :: handle
505 2204 : INTEGER, DIMENSION(:), POINTER :: blk_ao, dist_row_ao
506 : REAL(KIND=dp), ALLOCATABLE :: u_ri(:), v_vec(:), w_ri(:)
507 : TYPE(dbcsr_distribution_type) :: dist_ao_ao
508 : TYPE(dbcsr_type) :: matrix_Phi, matrix_V_H_AO
509 :
510 2204 : CALL timeset(routineN, handle)
511 :
512 2204 : CPASSERT(bs_env%ri_rs%V_grid_built)
513 :
514 2204 : CALL setup_square_topology(bs_env%ri_rs%mat_phi_mu_l, 'COL', dist_ao_ao, blk_ao, dist_row_ao)
515 6612 : ALLOCATE (v_vec(SIZE(n_vec)))
516 :
517 : ! stage 2: v_l = sum_PQl' Z_lP V^aux_PQ Z_l'Q n_l' (factorized; helper zeroes each output)
518 8816 : ALLOCATE (w_ri(bs_env%n_RI), u_ri(bs_env%n_RI))
519 : ! w_Q = sum_l Z_lQ n_l
520 : CALL dbcsr_matvec_replicated(bs_env%ri_rs%mat_Z_lP, n_vec, w_ri, bs_env%para_env, &
521 2204 : transposed=.TRUE.)
522 : ! u_P = sum_Q V^aux_PQ w_Q
523 2204 : CALL dbcsr_matvec_replicated(bs_env%ri_rs%mat_V_aux_rtbse, w_ri, u_ri, bs_env%para_env)
524 : ! v_l = sum_P Z_lP u_P
525 2204 : CALL dbcsr_matvec_replicated(bs_env%ri_rs%mat_Z_lP, u_ri, v_vec, bs_env%para_env)
526 2204 : DEALLOCATE (w_ri, u_ri)
527 :
528 : ! stage 3: V^H_µν = sum_l φ_lµ v_l φ_lν (Φ_lν = v_l φ_lν rowscale, then V^H = φ^T Φ)
529 2204 : CALL dbcsr_create(matrix_Phi, template=bs_env%ri_rs%mat_phi_mu_l)
530 2204 : CALL dbcsr_copy(matrix_Phi, bs_env%ri_rs%mat_phi_mu_l)
531 2204 : CALL dbcsr_scale_rows_replicated(matrix_Phi, v_vec)
532 :
533 : CALL dbcsr_create(matrix_V_H_AO, "V_H_AO_hartree", dist_ao_ao, dbcsr_type_no_symmetry, &
534 2204 : blk_ao, blk_ao)
535 : ! No eps_filter: a Δρ-dependent block drop breaks Hartree kernel self-adjointness.
536 : CALL dbcsr_multiply("T", "N", 1.0_dp, bs_env%ri_rs%mat_phi_mu_l, matrix_Phi, &
537 2204 : 0.0_dp, matrix_V_H_AO, filter_eps=0.0_dp)
538 2204 : CALL dbcsr_release(matrix_Phi)
539 :
540 2204 : CALL copy_dbcsr_to_fm(matrix_V_H_AO, V_H_AO_fm)
541 :
542 2204 : DEALLOCATE (v_vec)
543 :
544 : CALL release_dbcsr_topology_and_matrices(dist=dist_ao_ao, mapped_dist=dist_row_ao, &
545 2204 : m1=matrix_V_H_AO)
546 :
547 2204 : CALL timestop(handle)
548 :
549 4408 : END SUBROUTINE hartree_potential_from_diag_ri_rs
550 :
551 : ! **************************************************************************************************
552 : !> \brief Complex Hartree from precomputed grid diagonals: V^H = V^H[n_re] + i V^H[n_im], each via
553 : !> hartree_potential_from_diag_ri_rs (stages 2-3 only). n_re/n_im are the spin-summed grid
554 : !> densities harvested in the SEX kernel; this is the cross-spin / TDA complex consumer that
555 : !> replaces compute_hartree_ri_rs_complex when SEX already built the grid. n_im optional: when
556 : !> absent the result is purely real (matches the Re-only real-input Hartree).
557 : !> \param bs_env ...
558 : !> \param n_re grid density Re part (length n_grid, replicated)
559 : !> \param V_H_AO_cfm output AO x AO complex Hartree potential
560 : !> \param n_im optional grid density Im part (length n_grid, replicated)
561 : ! **************************************************************************************************
562 1120 : SUBROUTINE compute_hartree_ri_rs_from_diag(bs_env, n_re, V_H_AO_cfm, n_im)
563 :
564 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
565 : REAL(KIND=dp), INTENT(IN) :: n_re(:)
566 : TYPE(cp_cfm_type), INTENT(INOUT) :: V_H_AO_cfm
567 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: n_im(:)
568 :
569 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_hartree_ri_rs_from_diag'
570 :
571 : INTEGER :: handle
572 : TYPE(cp_cfm_type) :: cfm_real_part
573 : TYPE(cp_fm_type) :: fm_v
574 :
575 1120 : CALL timeset(routineN, handle)
576 :
577 1120 : CALL cp_fm_create(fm_v, bs_env%fm_s_Gamma%matrix_struct)
578 :
579 1120 : CALL hartree_potential_from_diag_ri_rs(bs_env, n_re, fm_v)
580 1120 : IF (PRESENT(n_im)) THEN
581 640 : CALL cp_cfm_create(cfm_real_part, bs_env%fm_s_Gamma%matrix_struct)
582 640 : CALL cp_fm_to_cfm(msourcer=fm_v, mtarget=cfm_real_part)
583 640 : CALL hartree_potential_from_diag_ri_rs(bs_env, n_im, fm_v)
584 640 : CALL cp_fm_to_cfm(msourcei=fm_v, mtarget=V_H_AO_cfm)
585 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), V_H_AO_cfm, &
586 640 : CMPLX(1.0_dp, 0.0_dp, kind=dp), cfm_real_part)
587 640 : CALL cp_cfm_release(cfm_real_part)
588 : ELSE
589 480 : CALL cp_fm_to_cfm(msourcer=fm_v, mtarget=V_H_AO_cfm)
590 : END IF
591 :
592 1120 : CALL cp_fm_release(fm_v)
593 :
594 1120 : CALL timestop(handle)
595 :
596 1120 : END SUBROUTINE compute_hartree_ri_rs_from_diag
597 :
598 : ! **************************************************************************************************
599 : !> \brief Complex-input Hartree potential via RI-RS. Re/Im split: feed each part to the real
600 : !> compute_hartree_ri_rs and reassemble. Real-input Hartree on a non-Hermitian input
601 : !> would silently drop Im and break Hermitian conjugacy of OV+VO contributions in TDA.
602 : !> \param bs_env ...
603 : !> \param rho_AO_cfm input AO x AO complex density-like matrix
604 : !> \param V_H_AO_cfm output AO x AO complex Hartree potential
605 : ! **************************************************************************************************
606 880 : SUBROUTINE compute_hartree_ri_rs_complex(bs_env, rho_AO_cfm, V_H_AO_cfm)
607 :
608 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
609 : TYPE(cp_cfm_type), INTENT(IN) :: rho_AO_cfm
610 : TYPE(cp_cfm_type), INTENT(INOUT) :: V_H_AO_cfm
611 :
612 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_hartree_ri_rs_complex'
613 :
614 : INTEGER :: handle
615 : TYPE(cp_cfm_type) :: cfm_real_part
616 : TYPE(cp_fm_type) :: fm_rho, fm_v
617 :
618 220 : CALL timeset(routineN, handle)
619 :
620 220 : CALL cp_fm_create(fm_rho, bs_env%fm_s_Gamma%matrix_struct)
621 220 : CALL cp_fm_create(fm_v, bs_env%fm_s_Gamma%matrix_struct)
622 220 : CALL cp_cfm_create(cfm_real_part, bs_env%fm_s_Gamma%matrix_struct)
623 :
624 220 : CALL cp_cfm_to_fm(msource=rho_AO_cfm, mtargetr=fm_rho)
625 220 : CALL compute_hartree_ri_rs(bs_env, fm_rho, fm_v)
626 220 : CALL cp_fm_to_cfm(msourcer=fm_v, mtarget=cfm_real_part)
627 :
628 220 : CALL cp_cfm_to_fm(msource=rho_AO_cfm, mtargeti=fm_rho)
629 220 : CALL compute_hartree_ri_rs(bs_env, fm_rho, fm_v)
630 220 : CALL cp_fm_to_cfm(msourcei=fm_v, mtarget=V_H_AO_cfm)
631 :
632 : CALL cp_cfm_scale_and_add(CMPLX(1.0_dp, 0.0_dp, kind=dp), V_H_AO_cfm, &
633 220 : CMPLX(1.0_dp, 0.0_dp, kind=dp), cfm_real_part)
634 :
635 220 : CALL cp_fm_release(fm_rho)
636 220 : CALL cp_fm_release(fm_v)
637 220 : CALL cp_cfm_release(cfm_real_part)
638 :
639 220 : CALL timestop(handle)
640 :
641 220 : END SUBROUTINE compute_hartree_ri_rs_complex
642 :
643 : ! **************************************************************************************************
644 : !> \brief Scale each row of a dbcsr matrix by a replicated full-length vector:
645 : !> block(ir,ic) <- vec(global_row(ir)) * block(ir,ic). Value mutation only.
646 : !> \param matrix ...
647 : !> \param vec full-length replicated row-scaling vector
648 : ! **************************************************************************************************
649 2204 : SUBROUTINE dbcsr_scale_rows_replicated(matrix, vec)
650 :
651 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
652 : REAL(KIND=dp), INTENT(IN) :: vec(:)
653 :
654 : INTEGER :: col_blk, ic, ir, nblkrows_total, &
655 : row_blk, row_off
656 2204 : INTEGER, ALLOCATABLE :: row_offset(:)
657 2204 : INTEGER, DIMENSION(:), POINTER :: row_blk_size
658 2204 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk
659 : TYPE(dbcsr_iterator_type) :: iter
660 :
661 2204 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_total, row_blk_size=row_blk_size)
662 6612 : ALLOCATE (row_offset(nblkrows_total + 1))
663 2204 : row_offset(1) = 0
664 25556 : DO ir = 1, nblkrows_total
665 25556 : row_offset(ir + 1) = row_offset(ir) + row_blk_size(ir)
666 : END DO
667 :
668 2204 : CALL dbcsr_iterator_start(iter, matrix)
669 25556 : DO WHILE (dbcsr_iterator_blocks_left(iter))
670 23352 : CALL dbcsr_iterator_next_block(iter, row_blk, col_blk, blk)
671 23352 : row_off = row_offset(row_blk)
672 165668 : DO ic = 1, SIZE(blk, 2)
673 4004424 : DO ir = 1, SIZE(blk, 1)
674 3981072 : blk(ir, ic) = vec(row_off + ir)*blk(ir, ic)
675 : END DO
676 : END DO
677 : END DO
678 2204 : CALL dbcsr_iterator_stop(iter)
679 :
680 2204 : DEALLOCATE (row_offset)
681 :
682 4408 : END SUBROUTINE dbcsr_scale_rows_replicated
683 :
684 : ! **************************************************************************************************
685 : !> \brief Replicated matvec: vec_out = matrix * vec_in (or matrix^T * vec_in if transposed) for a
686 : !> possibly rectangular distributed dbcsr matrix and replicated full-length vectors.
687 : !> Iterates over local blocks and reduces.
688 : !> \param matrix distributed dbcsr matrix (may be rectangular)
689 : !> \param vec_in full input vector, replicated on all ranks (column length, or row length if transposed)
690 : !> \param vec_out full output vector, replicated on all ranks (zeroed on entry; sum-reduced on exit)
691 : !> \param para_env ...
692 : !> \param transposed if .TRUE. compute vec_out = matrix^T * vec_in
693 : ! **************************************************************************************************
694 6612 : SUBROUTINE dbcsr_matvec_replicated(matrix, vec_in, vec_out, para_env, transposed)
695 :
696 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
697 : REAL(KIND=dp), INTENT(IN) :: vec_in(:)
698 : REAL(KIND=dp), INTENT(INOUT) :: vec_out(:)
699 : TYPE(mp_para_env_type), POINTER :: para_env
700 : LOGICAL, INTENT(IN), OPTIONAL :: transposed
701 :
702 : INTEGER :: col_blk, col_off, ic, ir, &
703 : nblkcols_total, nblkrows_total, &
704 : row_blk, row_off
705 6612 : INTEGER, ALLOCATABLE :: col_offset(:), row_offset(:)
706 6612 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
707 : LOGICAL :: my_trans
708 6612 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: block
709 : TYPE(dbcsr_iterator_type) :: iter
710 :
711 6612 : my_trans = .FALSE.
712 2204 : IF (PRESENT(transposed)) my_trans = transposed
713 :
714 : CALL dbcsr_get_info(matrix, nblkrows_total=nblkrows_total, &
715 : nblkcols_total=nblkcols_total, &
716 6612 : row_blk_size=row_blk_size, col_blk_size=col_blk_size)
717 :
718 33060 : ALLOCATE (row_offset(nblkrows_total + 1), col_offset(nblkcols_total + 1))
719 6612 : row_offset(1) = 0
720 57724 : DO ir = 1, nblkrows_total
721 57724 : row_offset(ir + 1) = row_offset(ir) + row_blk_size(ir)
722 : END DO
723 6612 : col_offset(1) = 0
724 19836 : DO ic = 1, nblkcols_total
725 19836 : col_offset(ic + 1) = col_offset(ic) + col_blk_size(ic)
726 : END DO
727 :
728 849540 : vec_out(:) = 0.0_dp
729 :
730 6612 : CALL dbcsr_iterator_start(iter, matrix)
731 57724 : DO WHILE (dbcsr_iterator_blocks_left(iter))
732 51112 : CALL dbcsr_iterator_next_block(iter, row_blk, col_blk, block)
733 51112 : row_off = row_offset(row_blk)
734 51112 : col_off = col_offset(col_blk)
735 57724 : IF (my_trans) THEN
736 560448 : DO ic = 1, SIZE(block, 2)
737 15284128 : DO ir = 1, SIZE(block, 1)
738 : vec_out(col_off + ic) = vec_out(col_off + ic) + &
739 15260776 : block(ir, ic)*vec_in(row_off + ir)
740 : END DO
741 : END DO
742 : ELSE
743 666240 : DO ic = 1, SIZE(block, 2)
744 17721752 : DO ir = 1, SIZE(block, 1)
745 : vec_out(row_off + ir) = vec_out(row_off + ir) + &
746 17693992 : block(ir, ic)*vec_in(col_off + ic)
747 : END DO
748 : END DO
749 : END IF
750 : END DO
751 6612 : CALL dbcsr_iterator_stop(iter)
752 :
753 6612 : DEALLOCATE (row_offset, col_offset)
754 :
755 1692468 : CALL para_env%sum(vec_out)
756 :
757 15428 : END SUBROUTINE dbcsr_matvec_replicated
758 :
759 : END MODULE rt_bse_ri_rs
|
