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