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