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 GW using RI-RS Approximation for molecules
10 : !> \par History
11 : !> 04.2026 created [Ritaj Tyagi]
12 : ! **************************************************************************************************
13 : MODULE gw_large_cell_Gamma_ri_rs
14 : USE atomic_kind_types, ONLY: atomic_kind_type
15 : USE basis_set_types, ONLY: gto_basis_set_type
16 : USE cell_types, ONLY: cell_type,&
17 : get_cell,&
18 : pbc
19 : USE constants_operator, ONLY: operator_coulomb
20 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
21 : cp_blacs_env_release,&
22 : cp_blacs_env_type
23 : USE cp_dbcsr_api, ONLY: &
24 : dbcsr_add, dbcsr_binary_read, dbcsr_binary_write, dbcsr_copy, dbcsr_create, &
25 : dbcsr_deallocate_matrix, dbcsr_distribution_get, dbcsr_distribution_new, &
26 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, &
27 : dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
28 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
29 : dbcsr_p_type, dbcsr_put_block, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
30 : dbcsr_type_no_symmetry
31 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_all_blocks
32 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
33 : copy_fm_to_dbcsr,&
34 : dbcsr_deallocate_matrix_set,&
35 : max_elements_per_block
36 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
37 : cp_fm_uplo_to_full
38 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
39 : cp_fm_cholesky_invert
40 : USE cp_fm_diag, ONLY: cp_fm_power
41 : USE cp_fm_struct, ONLY: cp_fm_struct_type
42 : USE cp_fm_types, ONLY: cp_fm_create,&
43 : cp_fm_get_info,&
44 : cp_fm_release,&
45 : cp_fm_set_all,&
46 : cp_fm_to_fm,&
47 : cp_fm_type
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_type
50 : USE cp_output_handling, ONLY: cp_p_file,&
51 : cp_print_key_should_output
52 : USE gw_integrals, ONLY: build_3c_integral_block_ctx,&
53 : gw_3c_ctx_create,&
54 : gw_3c_ctx_release,&
55 : gw_3c_ctx_type,&
56 : gw_3c_ws_create,&
57 : gw_3c_ws_release,&
58 : gw_3c_ws_type
59 : USE gw_large_cell_gamma, ONLY: &
60 : Fourier_transform_w_to_t, G_occ_vir, compute_QP_energies, compute_fm_chi_Gamma_freq, &
61 : create_fm_W_MIC_time, delete_unnecessary_files, fill_fm_Sigma_c_Gamma_time, fm_write, &
62 : get_W_MIC, multiply_fm_W_MIC_time_with_Minv_Gamma
63 : USE gw_non_periodic_ri_rs, ONLY: get_basis_offsets,&
64 : precompute_ri_rs_radii,&
65 : ri_rs_grid_assembler,&
66 : solve_D_lp_distributed
67 : USE gw_utils, ONLY: de_init_bs_env
68 : USE input_constants, ONLY: rtp_method_bse
69 : USE input_section_types, ONLY: section_vals_type
70 : USE kinds, ONLY: dp
71 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
72 : USE machine, ONLY: m_walltime
73 : USE message_passing, ONLY: mp_para_env_type
74 : USE mp2_ri_2c, ONLY: RI_2c_integral_mat
75 : USE orbital_pointers, ONLY: indco,&
76 : ncoset
77 : USE parallel_gemm_api, ONLY: parallel_gemm
78 : USE particle_types, ONLY: particle_type
79 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
80 : USE qs_environment_types, ONLY: get_qs_env,&
81 : qs_environment_type
82 : USE qs_kind_types, ONLY: get_qs_kind,&
83 : qs_kind_type
84 : #include "./base/base_uses.f90"
85 :
86 : IMPLICIT NONE
87 :
88 : PRIVATE
89 :
90 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_large_cell_Gamma_ri_rs'
91 :
92 : PUBLIC :: gw_calc_large_cell_Gamma_ri_rs, &
93 : contract_A_B_A, &
94 : hadamard_product_inplace, &
95 : release_dbcsr_topology_and_matrices, &
96 : setup_square_topology
97 :
98 : CONTAINS
99 :
100 : ! **************************************************************************************************
101 : !> \brief GW calculation using RI-RS formalism for molecules
102 : !> \param qs_env ...
103 : !> \param bs_env ...
104 : ! **************************************************************************************************
105 :
106 0 : SUBROUTINE gw_calc_large_cell_Gamma_ri_rs(qs_env, bs_env)
107 :
108 : TYPE(qs_environment_type), POINTER :: qs_env
109 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
110 :
111 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_large_cell_Gamma_ri_rs'
112 :
113 : INTEGER :: handle
114 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_W_time
115 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
116 :
117 0 : CALL timeset(routineN, handle)
118 :
119 : !!========================================================================
120 : !! 0. Precompute AO and RI Radii
121 : !! Per-atom cutoff radii from the most diffuse Gaussian primitives in
122 : !! the AO and RI auxiliary basis sets. Stored in bs_env%ri_rs%
123 : !! radius_ao_per_atom and radius_ri_per_atom, used for sphere-cutoff
124 : !! and phi_local screening.
125 : !!========================================================================
126 0 : CALL precompute_ri_rs_radii(qs_env, bs_env)
127 :
128 : !!========================================================================
129 : !! 1. Grid Generation for RI-RS
130 : !! (Modified Lebedev grids from Ivan Duchemin and Xavier Blase)
131 : !! Generate flattened 1D array of grid points for RI-RS.
132 : !! Equation: r_g(k) = R_A + r_g(A)
133 : !!========================================================================
134 0 : CALL ri_rs_grid_assembler(qs_env, bs_env, bs_env%ri_rs%grid_points)
135 :
136 : !!========================================================================
137 : !! 2. Atomic Basis Evaluation
138 : !! Compute values of spherical atomic basis functions at grid points.
139 : !! Expression: Φ_μl = Φ_μ(r_l) (mat_phi_mu_l)
140 : !!========================================================================
141 : CALL atomic_basis_at_grid_point(qs_env, bs_env, bs_env%ri_rs%grid_points, &
142 0 : bs_env%ri_rs%mat_phi_mu_l)
143 :
144 : !!========================================================================
145 : !! 3. Compute RI-RS Coefficients (Z_lp)
146 : !! Solve the regularized system for each atom P, where the grid domain
147 : !! is restricted to r_l within a cutoff distance of atom P:
148 : !! a. D_ll' = [ Σ_μ Φ_μ(r_l) Φ_μ(r_l') ]^2 (Equation 13)
149 : !! b. D_lP = Σ_{μν} Φ_μ(r_l) Φ_ν(r_l) (μν|P) (Equation 15)
150 : !! c. Conditioning:
151 : !! Dvec_l = 1 / sqrt(D_ll) (Diagonal scaling vector)
152 : !! D'_ll' = Dvec_l * D_ll' * Dvec_l' + λδ_ll'
153 : !! D'_lP = Dvec_l * D_lP
154 : !! d. Solve: Σ_l' D'_ll' * Z'_l'P = D'_lP (Equation 14)
155 : !! e. Rescale: Z_lP = Z'_lP * Dvec_l (Z_lP stored in mat_Z_lP)
156 : !!========================================================================
157 : CALL compute_coeff_Z_lP(qs_env, bs_env, bs_env%ri_rs%grid_points, &
158 0 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
159 0 : bs_env%ri_rs%grid_built = .TRUE.
160 :
161 : !!========================================================================
162 : !! 4. Compute Independent-Particle Polarizability (χ)
163 : !! G^occ_µλ(i|τ|) = sum_n^occ C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
164 : !! G^vir_µλ(i|τ|) = sum_n^vir C_µn e^(-|(ϵ_n-ϵ_F)τ|) C_λn
165 : !! G^occ_ll'(i|τ|) = sum_µν Φ_µ(r_l) G^occ_µν Φ_ν(r_l')
166 : !! G^vir_ll'(i|τ|) = sum_µν Φ_µ(r_l) G^vir_µν Φ_ν(r_l')
167 : !! χ_ll'(iτ) = G^occ_ll'(i|τ|) * G^vir_ll'(i|τ|)
168 : !! χ_PQ(iτ) = sum_ll' Z_lP χ_ll'(iτ) Z_l'Q
169 : !!========================================================================
170 : CALL get_mat_chi_Gamma_tau(bs_env, bs_env%mat_chi_Gamma_tau, &
171 0 : bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
172 :
173 : !!========================================================================
174 : !! 5. Compute Screened Interaction (W^MIC)
175 : !! χ_PQ(iτ) -> χ_PQ(iω) -> ε_PQ(iω) -> W_PQ(iω) -> W^MIC_PQ(iτ)
176 : !!========================================================================
177 0 : CALL get_W_MIC(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_time)
178 :
179 : !!========================================================================
180 : !! 6. Compute Exact Exchange Self-Energy (Σ^x)
181 : !! D_µν = sum_n^occ C_µn C_νn
182 : !! D_ll' = sum_µν Φ_µ(r_l) D_µν Φ_ν(r_l')
183 : !! V^trunc_ll' = sum_PQ Z_lP V^trunc_PQ Z_l'Q
184 : !! Σ^x_ll' = D_ll' * V^trunc_ll'
185 : !! Σ^x_λσ(k=0) = -sum_ll' Φ_λ(r_l) Σ^x_ll' Φ_σ(r_l')
186 : !!========================================================================
187 : CALL compute_Sigma_x(bs_env, qs_env, bs_env%ri_rs%mat_phi_mu_l, &
188 0 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_x_Gamma)
189 :
190 : !!========================================================================
191 : !! 7. Compute Correlation Self-Energy (Σ^c)
192 : !! W^MIC_ll'(iτ) = sum_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
193 : !! Σ^c_ll'(iτ) = -G^occ_ll'(i|τ|) * W^MIC_ll'(iτ), for τ < 0
194 : !! Σ^c_ll'(iτ) = G^vir_ll'(i|τ|) * W^MIC_ll'(iτ), for τ > 0
195 : !! Σ^c_λσ(iτ) = sum_ll' Φ_λ(r_l) Σ^c_ll'(iτ) Φ_σ(r_l')
196 : !!========================================================================
197 : CALL compute_Sigma_c(bs_env, fm_W_time, bs_env%ri_rs%mat_phi_mu_l, &
198 0 : bs_env%ri_rs%mat_Z_lP, fm_Sigma_c_Gamma_time)
199 :
200 : !!========================================================================
201 : !! 8. Compute Quasiparticle Energies
202 : !! Σ^c_λσ(iτ) -> Σ^c_nn(ϵ)
203 : !! ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ) + Σ^x_nn - v^xc_nn
204 : !!========================================================================
205 0 : CALL compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma_time)
206 :
207 0 : CALL de_init_bs_env(qs_env, bs_env)
208 :
209 0 : CALL timestop(handle)
210 :
211 0 : END SUBROUTINE gw_calc_large_cell_Gamma_ri_rs
212 :
213 : ! **************************************************************************************************
214 : !> \brief Evaluates atomic basis functions on a real-space grid and builds a sparse DBCSR matrix.
215 : !> \param qs_env ...
216 : !> \param bs_env ...
217 : !> \param ri_rs_grid_points ...
218 : !> \param mat_phi_mu_l ...
219 : ! **************************************************************************************************
220 :
221 0 : SUBROUTINE atomic_basis_at_grid_point(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l)
222 :
223 : TYPE(qs_environment_type), POINTER :: qs_env
224 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
225 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
226 : TYPE(dbcsr_type), INTENT(OUT) :: mat_phi_mu_l
227 :
228 : CHARACTER(LEN=*), PARAMETER :: routineN = 'atomic_basis_at_grid_point'
229 :
230 : INTEGER :: c_size, chunk_size, dimen_ORB, handle, i, i_blk, iatom, natom, npcol, nprow, &
231 : num_grid_chunks, r_end, r_start, total_grid_npts
232 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf
233 0 : INTEGER, DIMENSION(:), POINTER :: c_blk_sizes, col_dist, r_blk_sizes, &
234 0 : row_dist
235 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_col_buffer
236 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
237 : TYPE(cell_type), POINTER :: cell
238 : TYPE(dbcsr_distribution_type) :: dist
239 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist_ks
240 : TYPE(mp_para_env_type), POINTER :: para_env
241 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
242 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
243 :
244 0 : CALL timeset(routineN, handle)
245 :
246 : ! Setup Grid Blocking
247 0 : chunk_size = max_elements_per_block
248 :
249 : ! Extract environment variables
250 : CALL get_qs_env(qs_env, cell=cell, atomic_kind_set=atomic_kind_set, &
251 : qs_kind_set=qs_kind_set, particle_set=particle_set, &
252 0 : para_env=para_env)
253 :
254 0 : natom = SIZE(particle_set)
255 0 : total_grid_npts = SIZE(ri_rs_grid_points, 2)
256 :
257 : ! Map the starting indices of spherical gaussian functions (SGF) for each atom
258 0 : ALLOCATE (first_sgf(natom + 1))
259 0 : CALL get_basis_offsets(particle_set, qs_kind_set, first_sgf, dimen_ORB)
260 :
261 : ! =========================================================================
262 : ! 1. SETUP DBCSR MATRIX TOPOLOGY
263 : ! =========================================================================
264 :
265 : ! A. Define Column Block Sizes (1 Block = 1 Atom's full basis set)
266 0 : ALLOCATE (c_blk_sizes(natom))
267 0 : DO iatom = 1, natom
268 0 : c_blk_sizes(iatom) = first_sgf(iatom + 1) - first_sgf(iatom)
269 : END DO
270 :
271 : ! B. Define Row Block Sizes (Grid chunks of max size 256)
272 0 : num_grid_chunks = CEILING(REAL(total_grid_npts, KIND=dp)/REAL(chunk_size, KIND=dp))
273 0 : ALLOCATE (r_blk_sizes(num_grid_chunks))
274 0 : r_blk_sizes = chunk_size
275 0 : IF (MOD(total_grid_npts, chunk_size) /= 0) THEN
276 0 : r_blk_sizes(num_grid_chunks) = MOD(total_grid_npts, chunk_size)
277 : END IF
278 :
279 : ! C. Fetch CP2K's Default Process Grid Configuration
280 0 : CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist_ks)
281 0 : CALL dbcsr_distribution_get(dbcsr_dist_ks, nprows=nprow, npcols=npcol)
282 :
283 : ! D. Build Custom Mappings using Round-Robin across the 2D process grid
284 :
285 0 : ALLOCATE (row_dist(num_grid_chunks))
286 0 : DO i = 1, num_grid_chunks
287 0 : row_dist(i) = MOD(i - 1, nprow)
288 : END DO
289 :
290 0 : ALLOCATE (col_dist(natom))
291 0 : DO i = 1, natom
292 0 : col_dist(i) = MOD(i - 1, npcol)
293 : END DO
294 :
295 : ! E. Create the DBCSR Distribution and Initialize the Matrix
296 : CALL dbcsr_distribution_new(dist, template=dbcsr_dist_ks, &
297 0 : row_dist=row_dist, col_dist=col_dist)
298 :
299 : CALL dbcsr_create(mat_phi_mu_l, name="phi_val_sparse", dist=dist, &
300 : matrix_type=dbcsr_type_no_symmetry, &
301 0 : row_blk_size=r_blk_sizes, col_blk_size=c_blk_sizes)
302 :
303 : ! =========================================================================
304 : ! 2. STREAM DATA DIRECTLY INTO SPARSE MATRIX
305 : ! =========================================================================
306 : ! Iterate over the atoms assigned to this specific MPI rank
307 0 : DO iatom = para_env%mepos + 1, natom, para_env%num_pe
308 :
309 0 : c_size = c_blk_sizes(iatom)
310 :
311 : ! Allocate a temporary dense buffer just for this specific atom
312 0 : ALLOCATE (atom_col_buffer(total_grid_npts, c_size))
313 0 : atom_col_buffer = 0.0_dp
314 :
315 : ! Evaluate the basis functions on the grid. Skip grid points outside the spatial
316 : ! extent of the most diffuse AO Gaussian on iatom; beyond that radius the contribution
317 : ! is guaranteed below eps_filter.
318 : CALL fill_phi_for_atom(atom_col_buffer, ri_rs_grid_points, total_grid_npts, &
319 : iatom, particle_set, qs_kind_set, cell, &
320 0 : r2_threshold=bs_env%ri_rs%radius_ao_per_atom(iatom)**2)
321 :
322 : ! Slice the dense column into chunks and insert into DBCSR
323 0 : DO i_blk = 1, num_grid_chunks
324 0 : r_start = (i_blk - 1)*chunk_size + 1
325 0 : r_end = MIN(i_blk*chunk_size, total_grid_npts)
326 :
327 : ! Apply dynamic sparsity filtering: Only store blocks with physical significance
328 0 : IF (MAXVAL(ABS(atom_col_buffer(r_start:r_end, 1:c_size))) > bs_env%eps_filter) THEN
329 : CALL dbcsr_put_block(mat_phi_mu_l, row=i_blk, col=iatom, &
330 0 : block=atom_col_buffer(r_start:r_end, 1:c_size))
331 : END IF
332 : END DO
333 :
334 0 : DEALLOCATE (atom_col_buffer)
335 :
336 : END DO
337 :
338 : ! Finalize triggers internal MPI communication to route blocks to their correct 2D process owners
339 0 : CALL dbcsr_finalize(mat_phi_mu_l)
340 :
341 0 : IF (bs_env%unit_nr > 0) THEN
342 0 : WRITE (bs_env%unit_nr, *) "Done with evaluation of phi"
343 : END IF
344 :
345 : ! -------------------------------------------------------------------------
346 : ! CLEANUP
347 : ! -------------------------------------------------------------------------
348 0 : DEALLOCATE (first_sgf, r_blk_sizes, c_blk_sizes, row_dist, col_dist)
349 0 : CALL dbcsr_distribution_release(dist)
350 :
351 0 : CALL timestop(handle)
352 :
353 0 : END SUBROUTINE atomic_basis_at_grid_point
354 :
355 : ! **************************************************************************************************
356 : !> \brief Compute value of all basis functions for a single atom across all grid points.
357 : !> Sums contributions from periodic images of `iatom` (loop over (ix, iy, iz) cells gated
358 : !> by `cell%perd`). Each per-image squared distance is compared against `r2_threshold`
359 : !> (per-atom AO Gaussian extent²); images beyond that radius contribute below eps_filter
360 : !> and are skipped.
361 : !> \param phi_val ...
362 : !> \param ri_rs_grid ...
363 : !> \param npts ...
364 : !> \param iatom ...
365 : !> \param particle_set ...
366 : !> \param qs_kind_set ...
367 : !> \param cell ...
368 : !> \param r2_threshold per-image squared-distance threshold; CYCLE if r² > r2_threshold. Pass
369 : !> HUGE(1.0_dp) to disable.
370 : ! **************************************************************************************************
371 :
372 0 : SUBROUTINE fill_phi_for_atom(phi_val, ri_rs_grid, npts, iatom, &
373 : particle_set, qs_kind_set, cell, r2_threshold)
374 :
375 : REAL(KIND=dp), INTENT(INOUT) :: phi_val(:, :)
376 : INTEGER, INTENT(IN) :: npts
377 : REAL(KIND=dp), INTENT(IN) :: ri_rs_grid(3, npts)
378 : INTEGER, INTENT(IN) :: iatom
379 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
380 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
381 : TYPE(cell_type), POINTER :: cell
382 : REAL(KIND=dp), INTENT(IN) :: r2_threshold
383 :
384 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_phi_for_atom'
385 :
386 : INTEGER :: first_sgf, handle, i_pt, ico, iend_co, ikind, ipgf, iset, isgf, ishell, &
387 : istart_co, ix, ix_max, ix_min, iy, iy_max, iy_min, iz, iz_max, iz_min, l, last_sgf, lx, &
388 : ly, lz, n_cart_total, row_idx
389 : REAL(KIND=dp) :: alpha, cell_vector(3), dist_vec(3), &
390 : dist_vec_raw(3), exp_val, poly, r2, &
391 : r_atom(3), weight
392 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
393 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
394 :
395 0 : CALL timeset(routineN, handle)
396 :
397 : ! Get Atom Info
398 0 : ikind = particle_set(iatom)%atomic_kind%kind_number
399 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
400 0 : CALL get_cell(cell=cell, h=hmat)
401 :
402 0 : IF (.NOT. ASSOCIATED(orb_basis_set)) THEN
403 0 : CALL timestop(handle)
404 0 : RETURN
405 : END IF
406 :
407 0 : IF (cell%perd(1) == 1) THEN; ix_min = -1; ix_max = 1; ELSE; ix_min = 0; ix_max = 0
408 : END IF
409 0 : IF (cell%perd(2) == 1) THEN; iy_min = -1; iy_max = 1; ELSE; iy_min = 0; iy_max = 0
410 : END IF
411 0 : IF (cell%perd(3) == 1) THEN; iz_min = -1; iz_max = 1; ELSE; iz_min = 0; iz_max = 0
412 : END IF
413 :
414 0 : r_atom = particle_set(iatom)%r
415 :
416 : !$OMP PARALLEL DO DEFAULT(NONE) &
417 : !$OMP SHARED(phi_val, ri_rs_grid, npts, orb_basis_set, r_atom, hmat, &
418 : !$OMP ncoset, indco, cell, ix_min, ix_max, &
419 : !$OMP iy_min, iy_max, iz_min, iz_max, r2_threshold) &
420 : !$OMP PRIVATE(i_pt, dist_vec_raw, ix, iy, iz, cell_vector, dist_vec, r2, iset, &
421 : !$OMP n_cart_total, ishell, l, istart_co, iend_co, first_sgf, last_sgf, &
422 : !$OMP ipgf, alpha, exp_val, isgf, ico, row_idx, weight, lx, ly, lz, poly) &
423 0 : !$OMP SCHEDULE(DYNAMIC)
424 :
425 : DO i_pt = 1, npts
426 :
427 : dist_vec_raw = ri_rs_grid(:, i_pt) - r_atom
428 :
429 : DO ix = ix_min, ix_max
430 : DO iy = iy_min, iy_max
431 : DO iz = iz_min, iz_max
432 :
433 : cell_vector(1:3) = MATMUL(hmat, REAL([ix, iy, iz], dp))
434 :
435 : dist_vec = dist_vec_raw - cell_vector
436 :
437 : r2 = DOT_PRODUCT(dist_vec, dist_vec)
438 :
439 : IF (r2 > r2_threshold) CYCLE
440 :
441 : DO iset = 1, orb_basis_set%nset
442 : n_cart_total = ncoset(orb_basis_set%lmax(iset))
443 :
444 : DO ishell = 1, orb_basis_set%nshell(iset)
445 : l = orb_basis_set%l(ishell, iset)
446 : istart_co = ncoset(l - 1) + 1
447 : iend_co = ncoset(l)
448 :
449 : first_sgf = orb_basis_set%first_sgf(ishell, iset)
450 : last_sgf = orb_basis_set%last_sgf(ishell, iset)
451 :
452 : DO ipgf = 1, orb_basis_set%npgf(iset)
453 : alpha = orb_basis_set%zet(ipgf, iset)
454 : exp_val = EXP(-alpha*r2)
455 :
456 : DO isgf = first_sgf, last_sgf
457 : DO ico = istart_co, iend_co
458 : row_idx = (ipgf - 1)*n_cart_total + ico
459 : weight = orb_basis_set%sphi(row_idx, isgf)
460 : lx = indco(1, ico)
461 : ly = indco(2, ico)
462 : lz = indco(3, ico)
463 : poly = (dist_vec(1)**lx)*(dist_vec(2)**ly)*(dist_vec(3)**lz)
464 :
465 : phi_val(i_pt, isgf) = phi_val(i_pt, isgf) + (weight*poly*exp_val)
466 :
467 : END DO
468 : END DO
469 : END DO
470 : END DO
471 : END DO
472 : END DO
473 : END DO
474 : END DO
475 : END DO
476 : !$OMP END PARALLEL DO
477 :
478 0 : CALL timestop(handle)
479 :
480 : END SUBROUTINE fill_phi_for_atom
481 :
482 : ! **************************************************************************************************
483 : !> \brief Compute RI-RS Coefficients (Z_lP)
484 : !> \param qs_env ...
485 : !> \param bs_env ...
486 : !> \param ri_rs_grid_points ...
487 : !> \param mat_phi_mu_l ...
488 : !> \param mat_Z_lP ...
489 : ! **************************************************************************************************
490 :
491 0 : SUBROUTINE compute_coeff_Z_lP(qs_env, bs_env, ri_rs_grid_points, mat_phi_mu_l, mat_Z_lP)
492 :
493 : ! Arguments
494 : TYPE(qs_environment_type), POINTER :: qs_env
495 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
496 : REAL(KIND=dp), ALLOCATABLE, INTENT(INOUT) :: ri_rs_grid_points(:, :)
497 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l
498 : TYPE(dbcsr_type), INTENT(OUT) :: mat_Z_lP
499 :
500 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
501 : routineN = 'compute_coeff_Z_lP'
502 :
503 : INTEGER :: atom_j_mepos, atom_j_stride, atom_P, atom_P_start, atom_P_stride, col_end, &
504 : col_start, current_chunk_size, g, group_handle, handle, handle_dpotrf, handle_dpotrs, &
505 : handle_dsyrk, i, i_blk, ikind, info, j, j_ri, l, loc_idx, loc_ptr, max_ao_size, &
506 : max_loc_ri, my_group, n_ao_total, n_grid_total, n_groups, n_loc_ri, n_local_grid, &
507 : n_procs_per_atom, natom, nkind, npcol_phi, num_grid_chunks, P_loop_atom, r_end, r_start, &
508 : source_atom
509 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: local_grid_idx, row_offset
510 0 : INTEGER, DIMENSION(:), POINTER :: col_dist_ri, r_blk_sizes, ri_blk_sizes, &
511 0 : row_dist_grid
512 : REAL(KIND=dp) :: cutoff_ri, cutoff_ri_2, d_sP, dist2_min, &
513 : r2_threshold, r_c, t1, t2, t3
514 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_atom, cutoff_ri_per_kind, &
515 0 : d_vec_local
516 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_local, phi_local, &
517 0 : sphere_grid, Z_blk
518 : REAL(KIND=dp), DIMENSION(3) :: dist_vec_raw, pos_P
519 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
520 : TYPE(cell_type), POINTER :: cell
521 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
522 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_b, fm_struct_D
523 : TYPE(cp_fm_type) :: fm_b, fm_D
524 : TYPE(cp_logger_type), POINTER :: logger
525 : TYPE(dbcsr_distribution_type) :: dist_phi, dist_Z
526 0 : TYPE(gw_3c_ctx_type) :: ctx_3c
527 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
528 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
529 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
530 : TYPE(section_vals_type), POINTER :: input
531 :
532 0 : CALL timeset(routineN, handle)
533 :
534 0 : t1 = m_walltime()
535 :
536 : CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, input=input, &
537 0 : cell=cell, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
538 :
539 : ! ---------------------------------------------------------------------
540 : ! Subgroup setup. Default G=1 keeps the single-rank BLAS path; G>1 splits
541 : ! ranks into atom-groups so the Cholesky on D_local distributes across G
542 : ! ranks (memory ~1/G) and the compute_d_lp build also splits across the
543 : ! subgroup. G=1 leaves para_env_sub / blacs_env_sub NULL — no subgroup
544 : ! comms created, atom_P loop uses per-rank round-robin, compute_d_lp runs
545 : ! its full atom_j range on each rank, no allreduce.
546 : ! ---------------------------------------------------------------------
547 0 : n_procs_per_atom = MIN(bs_env%ri_rs%n_procs_per_atom_z_lp, para_env%num_pe)
548 0 : IF (n_procs_per_atom < 1) n_procs_per_atom = 1
549 :
550 0 : NULLIFY (para_env_sub, blacs_env_sub)
551 0 : IF (n_procs_per_atom > 1) THEN
552 0 : n_groups = para_env%num_pe/n_procs_per_atom
553 0 : my_group = MIN(para_env%mepos/n_procs_per_atom, n_groups - 1)
554 0 : ALLOCATE (para_env_sub)
555 0 : CALL para_env_sub%from_split(para_env, my_group)
556 0 : CALL cp_blacs_env_create(blacs_env=blacs_env_sub, para_env=para_env_sub)
557 0 : atom_P_start = my_group + 1
558 0 : atom_P_stride = n_groups
559 0 : atom_j_mepos = para_env_sub%mepos
560 0 : atom_j_stride = para_env_sub%num_pe
561 : ELSE
562 0 : atom_P_start = para_env%mepos + 1
563 0 : atom_P_stride = para_env%num_pe
564 0 : atom_j_mepos = 0
565 0 : atom_j_stride = 1
566 : END IF
567 :
568 0 : natom = SIZE(bs_env%i_RI_start_from_atom)
569 0 : n_ao_total = bs_env%i_ao_end_from_atom(natom)
570 0 : n_grid_total = SIZE(ri_rs_grid_points, 2)
571 :
572 : ! =========================================================================
573 : ! 1. SETUP DBCSR TOPOLOGY & EXACT OFFSETS
574 : ! =========================================================================
575 0 : CALL dbcsr_get_info(mat_phi_mu_l, row_blk_size=r_blk_sizes, distribution=dist_phi)
576 : CALL dbcsr_distribution_get(dist_phi, row_dist=row_dist_grid, &
577 0 : group=group_handle, npcols=npcol_phi)
578 :
579 0 : num_grid_chunks = SIZE(r_blk_sizes)
580 :
581 0 : ALLOCATE (row_offset(num_grid_chunks))
582 0 : row_offset(1) = 0
583 0 : DO i_blk = 2, num_grid_chunks
584 0 : row_offset(i_blk) = row_offset(i_blk - 1) + r_blk_sizes(i_blk - 1)
585 : END DO
586 :
587 0 : ALLOCATE (ri_blk_sizes(natom), col_dist_ri(natom))
588 0 : DO atom_P = 1, natom
589 0 : ri_blk_sizes(atom_P) = bs_env%i_RI_end_from_atom(atom_P) - bs_env%i_RI_start_from_atom(atom_P) + 1
590 0 : col_dist_ri(atom_P) = MOD(atom_P - 1, npcol_phi)
591 : END DO
592 :
593 0 : CALL dbcsr_distribution_new(dist_Z, template=dist_phi, row_dist=row_dist_grid, col_dist=col_dist_ri)
594 :
595 0 : IF (bs_env%ri_rs%Z_lP_exists) THEN
596 : CALL dbcsr_binary_read(filepath=TRIM(bs_env%prefix)//"Z_lP.matrix", &
597 : distribution=dist_Z, &
598 0 : matrix_new=mat_Z_lP)
599 0 : IF (bs_env%unit_nr > 0) THEN
600 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
601 0 : 'Read Z_lP from file ', ' Execution time', m_walltime() - t1, ' s'
602 0 : WRITE (bs_env%unit_nr, '(A)') ' '
603 : END IF
604 : ELSE
605 :
606 : CALL dbcsr_create(mat_Z_lP, name="mat_Z_lP", dist=dist_Z, &
607 : matrix_type=dbcsr_type_no_symmetry, &
608 0 : row_blk_size=r_blk_sizes, col_blk_size=ri_blk_sizes)
609 :
610 0 : max_ao_size = 0
611 0 : DO j = 1, SIZE(bs_env%i_ao_start_from_atom)
612 0 : max_ao_size = MAX(max_ao_size, bs_env%i_ao_end_from_atom(j) - bs_env%i_ao_start_from_atom(j) + 1)
613 : END DO
614 0 : max_loc_ri = MAXVAL(ri_blk_sizes)
615 :
616 : ! Per-atom RI-RS integration sphere:
617 : ! cutoff_ri(P) = r_c + r_AO(P)
618 : ! where r_c is the truncated-Coulomb cutoff of the RI metric. The
619 : ! CUTOFF_RADIUS_RI_RS keyword (when > 0) overrides the entire cutoff calculation.
620 0 : nkind = SIZE(atomic_kind_set)
621 0 : ALLOCATE (cutoff_ri_per_atom(natom))
622 :
623 0 : IF (bs_env%ri_rs%cutoff_radius_ri_rs > 0.0_dp) THEN
624 0 : cutoff_ri_per_atom(:) = bs_env%ri_rs%cutoff_radius_ri_rs
625 : ELSE
626 0 : r_c = bs_env%ri_metric%cutoff_radius
627 0 : DO P_loop_atom = 1, natom
628 : cutoff_ri_per_atom(P_loop_atom) = &
629 0 : r_c + bs_env%ri_rs%radius_ao_per_atom(P_loop_atom)
630 : END DO
631 : END IF
632 :
633 0 : ALLOCATE (cutoff_ri_per_kind(nkind))
634 0 : cutoff_ri_per_kind(:) = 0.0_dp
635 0 : IF (bs_env%unit_nr > 0) THEN
636 0 : DO P_loop_atom = 1, natom
637 0 : ikind = particle_set(P_loop_atom)%atomic_kind%kind_number
638 : cutoff_ri_per_kind(ikind) = MAX(cutoff_ri_per_kind(ikind), &
639 0 : cutoff_ri_per_atom(P_loop_atom))
640 : END DO
641 0 : WRITE (bs_env%unit_nr, '(T2,A)') 'Per-kind maximum RI-RS sphere cutoff (Bohr):'
642 0 : WRITE (bs_env%unit_nr, '(T4,A4,A14)') 'Kind', 'max cutoff_ri'
643 0 : DO ikind = 1, nkind
644 : WRITE (bs_env%unit_nr, '(T4,A4,F14.4)') &
645 0 : atomic_kind_set(ikind)%element_symbol, &
646 0 : cutoff_ri_per_kind(ikind)
647 : END DO
648 0 : WRITE (bs_env%unit_nr, '(A)') ' '
649 0 : DEALLOCATE (cutoff_ri_per_kind)
650 : END IF
651 :
652 : ! Shared 3c-integral context: hoists libint / t_c_g0 / md_ftable / contracted
653 : ! sphi tables out of the per-triple call so compute_d_lp threads only allocate
654 : ! a lightweight per-thread workspace. MPI-collective; must be outside any
655 : ! OMP region.
656 : CALL gw_3c_ctx_create(ctx_3c, qs_env, bs_env%ri_metric, &
657 : basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, &
658 0 : basis_i=bs_env%basis_set_RI)
659 :
660 : ! =========================================================================
661 : ! 2. MPI LOOP OVER ATOMS (Fully independent, no MPI barriers inside)
662 : ! phi_local for each atom_P's cutoff sphere is built on the fly via
663 : ! fill_phi_for_atom — no dense replicated phi_global, no allreduce.
664 : ! =========================================================================
665 0 : DO atom_P = atom_P_start, natom, atom_P_stride
666 :
667 0 : n_loc_ri = ri_blk_sizes(atom_P)
668 0 : pos_P(:) = particle_set(atom_P)%r(:)
669 :
670 0 : cutoff_ri = cutoff_ri_per_atom(atom_P)
671 0 : cutoff_ri_2 = cutoff_ri**2
672 :
673 : ! ---------------------------------------------------------------------
674 : ! A. Determine Local Grid Domain based on cutoff_ri (PBC distance)
675 : ! ---------------------------------------------------------------------
676 0 : n_local_grid = 0
677 0 : DO l = 1, n_grid_total
678 0 : dist_vec_raw = pbc(ri_rs_grid_points(1:3, l), pos_P(1:3), cell)
679 0 : dist2_min = DOT_PRODUCT(dist_vec_raw, dist_vec_raw)
680 0 : IF (dist2_min <= cutoff_ri_2) n_local_grid = n_local_grid + 1
681 : END DO
682 :
683 0 : ALLOCATE (local_grid_idx(n_local_grid))
684 :
685 0 : n_local_grid = 0
686 0 : DO l = 1, n_grid_total
687 0 : dist_vec_raw = pbc(ri_rs_grid_points(1:3, l), pos_P(1:3), cell)
688 0 : dist2_min = DOT_PRODUCT(dist_vec_raw, dist_vec_raw)
689 0 : IF (dist2_min <= cutoff_ri_2) THEN
690 0 : n_local_grid = n_local_grid + 1
691 0 : local_grid_idx(n_local_grid) = l
692 : END IF
693 : END DO
694 :
695 : ! ---------------------------------------------------------------------
696 : ! B. Build phi_local on the fly via fill_phi_for_atom.
697 : ! Only source atoms whose AO basis can reach the cutoff sphere of
698 : ! atom_P (MIC distance) contribute; the rest are pruned. The
699 : ! periodic fill_phi_for_atom sums over (ix, iy, iz) images of
700 : ! source_atom internally.
701 : ! ---------------------------------------------------------------------
702 0 : ALLOCATE (sphere_grid(3, n_local_grid))
703 0 : DO loc_idx = 1, n_local_grid
704 0 : sphere_grid(:, loc_idx) = ri_rs_grid_points(:, local_grid_idx(loc_idx))
705 : END DO
706 :
707 0 : ALLOCATE (phi_local(n_local_grid, n_ao_total))
708 0 : phi_local = 0.0_dp
709 :
710 0 : DO source_atom = 1, natom
711 0 : dist_vec_raw = pbc(particle_set(source_atom)%r(:), pos_P(:), cell)
712 0 : d_sP = NORM2(dist_vec_raw)
713 0 : IF (d_sP > bs_env%ri_rs%radius_ao_per_atom(source_atom) + cutoff_ri) CYCLE
714 :
715 0 : col_start = bs_env%i_ao_start_from_atom(source_atom)
716 0 : col_end = bs_env%i_ao_end_from_atom(source_atom)
717 0 : r2_threshold = bs_env%ri_rs%radius_ao_per_atom(source_atom)**2
718 :
719 : CALL fill_phi_for_atom(phi_local(:, col_start:col_end), sphere_grid, &
720 : n_local_grid, source_atom, particle_set, qs_kind_set, &
721 0 : cell, r2_threshold)
722 : END DO
723 :
724 0 : DEALLOCATE (sphere_grid)
725 :
726 : ! ---------------------------------------------------------------------
727 : ! C. Build Local RHS Matrix (d_lp_local) first so the subgroup-
728 : ! distributed compute_d_lp + allreduce is not entangled with the LHS
729 : ! build. compute_d_lp does not depend on D_local or d_vec_local.
730 : ! ---------------------------------------------------------------------
731 0 : ALLOCATE (d_lp_local(n_local_grid, n_loc_ri))
732 0 : d_lp_local = 0.0_dp
733 :
734 0 : t2 = m_walltime()
735 :
736 : CALL compute_d_lp(bs_env, ctx_3c, cell, phi_local, d_lp_local, n_local_grid, &
737 0 : n_loc_ri, atom_P, max_ao_size, atom_j_mepos, atom_j_stride)
738 :
739 : ! Reduce per-subgroup-rank partials into the replicated d_lp_local.
740 : ! Skipped for G=1 (BLAS path): each rank has the full sum locally.
741 0 : IF (n_procs_per_atom > 1) THEN
742 0 : CALL para_env_sub%sum(d_lp_local)
743 : END IF
744 :
745 0 : t3 = m_walltime()
746 :
747 : ! ---------------------------------------------------------------------
748 : ! D. Build d_vec_local (Jacobi diagonal) + LHS — BLAS or ScaLAPACK
749 : ! ---------------------------------------------------------------------
750 0 : ALLOCATE (d_vec_local(n_local_grid))
751 :
752 0 : IF (n_procs_per_atom == 1) THEN
753 : ! BLAS path: build D_local densely, compute d_vec as side-effect
754 : ! of the Jacobi step (preserves bit-identical arithmetic with the
755 : ! previous branch).
756 0 : ALLOCATE (D_local(n_local_grid, n_local_grid))
757 0 : D_local = 0.0_dp
758 :
759 0 : CALL timeset(routineN//"_dsyrk", handle_dsyrk)
760 : CALL dsyrk("L", "N", n_local_grid, n_ao_total, 1.0_dp, phi_local, &
761 0 : n_local_grid, 0.0_dp, D_local, n_local_grid)
762 0 : CALL timestop(handle_dsyrk)
763 :
764 : !$OMP PARALLEL DO DEFAULT(NONE) &
765 : !$OMP SHARED(n_local_grid, D_local, d_vec_local, bs_env) &
766 : !$OMP PRIVATE(i) &
767 0 : !$OMP SCHEDULE(STATIC)
768 : DO i = 1, n_local_grid
769 : D_local(i, i) = D_local(i, i)**2
770 : d_vec_local(i) = 1.0_dp/SQRT(MAX(D_local(i, i), 1.0E-16_dp))
771 : D_local(i, i) = (D_local(i, i)*d_vec_local(i)**2) + bs_env%ri_rs%tikhonov
772 : END DO
773 : !$OMP END PARALLEL DO
774 :
775 : !$OMP PARALLEL DO DEFAULT(NONE) &
776 : !$OMP SHARED(n_local_grid, D_local, d_vec_local) &
777 : !$OMP PRIVATE(j, i) &
778 0 : !$OMP SCHEDULE(DYNAMIC)
779 : DO j = 1, n_local_grid
780 : DO i = j + 1, n_local_grid
781 : D_local(i, j) = D_local(i, j)**2
782 : D_local(i, j) = D_local(i, j)*d_vec_local(i)*d_vec_local(j)
783 : D_local(j, i) = D_local(i, j)
784 : END DO
785 : END DO
786 : !$OMP END PARALLEL DO
787 : ELSE
788 : ! ScaLAPACK path: d_vec computed directly from phi (= 1/||phi_i||^2);
789 : ! solve_D_lp_distributed builds D block-cyclic internally with
790 : ! the squared+scaled values, so no dense D_local on this rank.
791 : !$OMP PARALLEL DO DEFAULT(NONE) &
792 : !$OMP SHARED(n_local_grid, n_ao_total, phi_local, d_vec_local) &
793 : !$OMP PRIVATE(i, j) &
794 0 : !$OMP SCHEDULE(STATIC)
795 : DO i = 1, n_local_grid
796 : d_vec_local(i) = 0.0_dp
797 : DO j = 1, n_ao_total
798 : d_vec_local(i) = d_vec_local(i) + phi_local(i, j)*phi_local(i, j)
799 : END DO
800 : d_vec_local(i) = 1.0_dp/MAX(d_vec_local(i), 1.0E-16_dp)
801 : END DO
802 : !$OMP END PARALLEL DO
803 : END IF
804 :
805 : ! ---------------------------------------------------------------------
806 : ! E. Pre-scale d_lp by d_vec
807 : ! ---------------------------------------------------------------------
808 : !$OMP PARALLEL DO DEFAULT(NONE) &
809 : !$OMP SHARED(n_loc_ri, n_local_grid, d_lp_local, d_vec_local) &
810 : !$OMP PRIVATE(j_ri, i) &
811 0 : !$OMP SCHEDULE(STATIC)
812 : DO j_ri = 1, n_loc_ri
813 : DO i = 1, n_local_grid
814 : d_lp_local(i, j_ri) = d_lp_local(i, j_ri)*d_vec_local(i)
815 : END DO
816 : END DO
817 : !$OMP END PARALLEL DO
818 :
819 : ! ---------------------------------------------------------------------
820 : ! F. Solve — BLAS dpotrf/dpotrs or ScaLAPACK pdpotrf/pdpotrs
821 : ! ---------------------------------------------------------------------
822 0 : IF (n_procs_per_atom == 1) THEN
823 0 : CALL timeset(routineN//"_dpotrf", handle_dpotrf)
824 0 : CALL dpotrf('L', n_local_grid, D_local, n_local_grid, info)
825 0 : CALL timestop(handle_dpotrf)
826 0 : CALL timeset(routineN//"_dpotrs", handle_dpotrs)
827 : CALL dpotrs('L', n_local_grid, n_loc_ri, D_local, n_local_grid, &
828 0 : d_lp_local, n_local_grid, info)
829 0 : CALL timestop(handle_dpotrs)
830 0 : DEALLOCATE (D_local)
831 : ELSE
832 : CALL solve_D_lp_distributed(phi_local, d_vec_local, d_lp_local, &
833 : n_local_grid, n_ao_total, n_loc_ri, &
834 : bs_env%ri_rs%tikhonov, &
835 : para_env_sub, blacs_env_sub, &
836 0 : fm_struct_D, fm_struct_b, fm_D, fm_b, info)
837 : END IF
838 :
839 : ! ---------------------------------------------------------------------
840 : ! G. Post-scale solution by d_vec (common to both paths)
841 : ! ---------------------------------------------------------------------
842 : !$OMP PARALLEL DO DEFAULT(NONE) &
843 : !$OMP SHARED(n_loc_ri, n_local_grid, d_lp_local, d_vec_local) &
844 : !$OMP PRIVATE(j_ri, i) &
845 0 : !$OMP SCHEDULE(STATIC)
846 : DO j_ri = 1, n_loc_ri
847 : DO i = 1, n_local_grid
848 : d_lp_local(i, j_ri) = d_lp_local(i, j_ri)*d_vec_local(i)
849 : END DO
850 : END DO
851 : !$OMP END PARALLEL DO
852 :
853 : ! ---------------------------------------------------------------------
854 : ! H. Scatter Local Solution Back to Global DBCSR Matrix.
855 : ! Under ScaLAPACK (G>1) the d_lp_local solution is identical on all
856 : ! G subgroup ranks (gathered via cp_fm_get_submatrix); only the
857 : ! subgroup root writes to mat_Z_lP so each atom column is emitted
858 : ! exactly once. DBCSR routes blocks to their global owner on finalize.
859 : ! local_grid_idx is ascending (built by the ordered scan above), so
860 : ! a single walking pointer over chunks works.
861 : ! ---------------------------------------------------------------------
862 0 : IF (n_procs_per_atom == 1 .OR. para_env_sub%mepos == 0) THEN
863 0 : ALLOCATE (Z_blk(MAXVAL(r_blk_sizes), n_loc_ri))
864 0 : loc_ptr = 1
865 :
866 0 : DO i_blk = 1, num_grid_chunks
867 0 : r_start = row_offset(i_blk) + 1
868 0 : r_end = row_offset(i_blk) + r_blk_sizes(i_blk)
869 0 : current_chunk_size = r_blk_sizes(i_blk)
870 :
871 0 : Z_blk = 0.0_dp
872 :
873 0 : DO WHILE (loc_ptr <= n_local_grid)
874 0 : g = local_grid_idx(loc_ptr)
875 0 : IF (g > r_end) EXIT
876 0 : Z_blk(g - r_start + 1, 1:n_loc_ri) = d_lp_local(loc_ptr, 1:n_loc_ri)
877 0 : loc_ptr = loc_ptr + 1
878 : END DO
879 :
880 0 : IF (MAXVAL(ABS(Z_blk(1:current_chunk_size, 1:n_loc_ri))) > bs_env%eps_filter) THEN
881 : CALL dbcsr_put_block(mat_Z_lP, row=i_blk, col=atom_P, &
882 0 : block=Z_blk(1:current_chunk_size, 1:n_loc_ri))
883 : END IF
884 : END DO
885 :
886 0 : DEALLOCATE (Z_blk)
887 : END IF
888 :
889 0 : DEALLOCATE (d_vec_local, d_lp_local)
890 0 : DEALLOCATE (local_grid_idx, phi_local)
891 :
892 : END DO
893 :
894 0 : DEALLOCATE (cutoff_ri_per_atom)
895 0 : CALL gw_3c_ctx_release(ctx_3c)
896 :
897 0 : CALL dbcsr_finalize(mat_Z_lP)
898 :
899 0 : IF (bs_env%unit_nr > 0) THEN
900 : WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
901 0 : 'Computed Z_lP ', ' Execution time', m_walltime() - t1, ' s'
902 0 : WRITE (bs_env%unit_nr, '(A)') ' '
903 : END IF
904 :
905 0 : logger => cp_get_default_logger()
906 :
907 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
908 0 : CALL dbcsr_binary_write(matrix=mat_Z_lP, filepath=TRIM(bs_env%prefix)//"Z_lP.matrix")
909 : END IF
910 :
911 : END IF
912 :
913 0 : DEALLOCATE (row_offset, ri_blk_sizes, col_dist_ri)
914 0 : CALL dbcsr_distribution_release(dist_Z)
915 :
916 0 : IF (n_procs_per_atom > 1) THEN
917 0 : CALL cp_blacs_env_release(blacs_env_sub)
918 0 : CALL para_env_sub%free()
919 0 : DEALLOCATE (para_env_sub)
920 : END IF
921 :
922 0 : DEALLOCATE (ri_rs_grid_points)
923 :
924 0 : CALL timestop(handle)
925 :
926 0 : END SUBROUTINE compute_coeff_Z_lP
927 :
928 : ! **************************************************************************************************
929 : !> \brief Computes the dense localized RHS d_lp(l,P) = Σ_{μν,R,S} Φ_μ(r_l)·Φ_ν(r_l)·(μν|P) for one
930 : !> RI atom P. OMP-threaded over (atom_j, atom_k) AO-pair blocks: per thread, sweep all
931 : !> (cell_R, cell_S) periodic images of (atom_j, atom_k) about atom_P at cell (0,0,0); each
932 : !> 3c block is built by build_3c_integral_block_ctx (cached libint / sphi tables in ctx,
933 : !> kind-radius triangle screen → `screened` short-circuits negligible image triples), and
934 : !> grid-chunked pair densities are contracted into a private d_lp partial that is reduced
935 : !> into d_lp at the end of the parallel region.
936 : !> \param bs_env ...
937 : !> \param ctx shared 3c-integral context (gw_3c_ctx_create)
938 : !> \param cell ...
939 : !> \param phi_val Φ_μ(r_l) on the local-sphere grid (n_grid_total × n_ao)
940 : !> \param d_lp output (n_grid_total × n_loc_ri), zeroed by the caller, accumulated here
941 : !> \param n_grid_total number of local-sphere grid rows
942 : !> \param n_loc_ri number of RI functions of atom_P
943 : !> \param atom_P RI atom (pinned to cell (0,0,0))
944 : !> \param max_ao_size ...
945 : !> \param atom_j_mepos ...
946 : !> \param atom_j_stride ...
947 : ! **************************************************************************************************
948 :
949 0 : SUBROUTINE compute_d_lp(bs_env, ctx, cell, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, &
950 : max_ao_size, atom_j_mepos, atom_j_stride)
951 :
952 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
953 : TYPE(gw_3c_ctx_type), INTENT(IN) :: ctx
954 : TYPE(cell_type), POINTER :: cell
955 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val
956 : INTEGER, INTENT(IN) :: n_grid_total, n_loc_ri
957 : REAL(KIND=dp), INTENT(INOUT) :: d_lp(n_grid_total, n_loc_ri)
958 : INTEGER, INTENT(IN) :: atom_P, max_ao_size, atom_j_mepos, &
959 : atom_j_stride
960 :
961 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp'
962 : INTEGER, PARAMETER :: grid_chunk = 1024
963 :
964 : INTEGER :: atom_j, atom_k, c, handle, handle_dgemm, ix_max, ix_min, ix_R, ix_S, iy_max, &
965 : iy_min, iy_R, iy_S, iz_max, iz_min, iz_R, iz_S, j, jk_idx, jsize, jstart, k, ksize, &
966 : kstart, l, l0, natom, ri
967 : INTEGER, DIMENSION(3) :: cell_R_vec, cell_S_vec
968 : LOGICAL :: any_kept, screened
969 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: d_lp_prv, int_2d_prv, rho_chunk
970 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_prv, int_3c_sum
971 0 : TYPE(gw_3c_ws_type) :: ws
972 :
973 0 : CALL timeset(routineN, handle)
974 :
975 0 : natom = SIZE(bs_env%i_ao_start_from_atom)
976 :
977 0 : IF (cell%perd(1) == 1) THEN; ix_min = -1; ix_max = 1; ELSE; ix_min = 0; ix_max = 0
978 : END IF
979 0 : IF (cell%perd(2) == 1) THEN; iy_min = -1; iy_max = 1; ELSE; iy_min = 0; iy_max = 0
980 : END IF
981 0 : IF (cell%perd(3) == 1) THEN; iz_min = -1; iz_max = 1; ELSE; iz_min = 0; iz_max = 0
982 : END IF
983 :
984 : !$OMP PARALLEL DEFAULT(NONE) &
985 : !$OMP SHARED(bs_env, ctx, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, max_ao_size, &
986 : !$OMP natom, ix_min, ix_max, iy_min, iy_max, iz_min, iz_max, &
987 : !$OMP atom_j_mepos, atom_j_stride) &
988 : !$OMP PRIVATE(any_kept, atom_j, atom_k, c, handle_dgemm, j, jk_idx, jsize, jstart, k, &
989 : !$OMP ksize, kstart, l, l0, ri, ix_R, iy_R, iz_R, ix_S, iy_S, iz_S, cell_R_vec, &
990 0 : !$OMP cell_S_vec, screened, d_lp_prv, int_2d_prv, rho_chunk, int_3c_prv, int_3c_sum, ws)
991 :
992 : CALL gw_3c_ws_create(ws, ctx)
993 : ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, n_loc_ri))
994 : ALLOCATE (int_3c_sum(max_ao_size, max_ao_size, n_loc_ri))
995 : ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, n_loc_ri))
996 : ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size))
997 : ALLOCATE (d_lp_prv(n_grid_total, n_loc_ri))
998 : d_lp_prv(:, :) = 0.0_dp
999 :
1000 : ! atom_P pinned at cell (0,0,0); enumerate (atom_j, cell_R) × (atom_k, cell_S). The ctx
1001 : ! integral builder's kind_radius triangle screen sets screened=.TRUE. for the bulk of
1002 : ! image triples (one or both AO atoms beyond the truncated-Coulomb reach of atom_P),
1003 : ! so the 27 × 27 = 729 candidate cells collapse to "adjacent cells" in practice.
1004 : ! MPI-stride atom_j over the subgroup (atom_j_stride = 1 for the BLAS path, > 1 for the
1005 : ! ScaLAPACK path). COLLAPSE(2) dropped because the outer stride is non-unit under
1006 : ! ScaLAPACK; the inner atom_k loop carries enough work for DYNAMIC.
1007 : !$OMP DO SCHEDULE(DYNAMIC)
1008 : DO atom_j = atom_j_mepos + 1, natom, atom_j_stride
1009 : DO atom_k = 1, natom
1010 : jstart = bs_env%i_ao_start_from_atom(atom_j)
1011 : jsize = bs_env%i_ao_end_from_atom(atom_j) - jstart + 1
1012 : kstart = bs_env%i_ao_start_from_atom(atom_k)
1013 : ksize = bs_env%i_ao_end_from_atom(atom_k) - kstart + 1
1014 :
1015 : int_3c_sum(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp
1016 : any_kept = .FALSE.
1017 :
1018 : DO ix_R = ix_min, ix_max
1019 : DO iy_R = iy_min, iy_max
1020 : DO iz_R = iz_min, iz_max
1021 : cell_R_vec = [ix_R, iy_R, iz_R]
1022 : DO ix_S = ix_min, ix_max
1023 : DO iy_S = iy_min, iy_max
1024 : DO iz_S = iz_min, iz_max
1025 : cell_S_vec = [ix_S, iy_S, iz_S]
1026 :
1027 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp
1028 :
1029 : CALL build_3c_integral_block_ctx( &
1030 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri), ctx, ws, &
1031 : atom_j=atom_j, atom_k=atom_k, atom_i=atom_P, &
1032 : cell_j=cell_R_vec, cell_k=cell_S_vec, cell_i=[0, 0, 0], &
1033 : screened=screened)
1034 : IF (screened) CYCLE
1035 :
1036 : any_kept = .TRUE.
1037 : int_3c_sum(1:jsize, 1:ksize, 1:n_loc_ri) = &
1038 : int_3c_sum(1:jsize, 1:ksize, 1:n_loc_ri) + &
1039 : int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri)
1040 : END DO
1041 : END DO
1042 : END DO
1043 : END DO
1044 : END DO
1045 : END DO
1046 :
1047 : IF (.NOT. any_kept) CYCLE
1048 :
1049 : ! Flatten 3D B_{μν,P} → 2D B_{(μν),P}
1050 : DO ri = 1, n_loc_ri
1051 : DO k = 1, ksize
1052 : DO j = 1, jsize
1053 : jk_idx = (k - 1)*jsize + j
1054 : int_2d_prv(jk_idx, ri) = int_3c_sum(j, k, ri)
1055 : END DO
1056 : END DO
1057 : END DO
1058 :
1059 : ! Pair density ρ(l,μν) = Φ_μ(r_l)Φ_ν(r_l) in grid chunks, contracted on the fly:
1060 : ! d_{l,P} += ρ(l,μν) B_{(μν),P} (dgemm runs serially inside the parallel region)
1061 : DO l0 = 1, n_grid_total, grid_chunk
1062 : c = MIN(grid_chunk, n_grid_total - l0 + 1)
1063 : DO k = 1, ksize
1064 : DO j = 1, jsize
1065 : jk_idx = (k - 1)*jsize + j
1066 : DO l = 1, c
1067 : rho_chunk(l, jk_idx) = phi_val(l0 + l - 1, jstart + j - 1)* &
1068 : phi_val(l0 + l - 1, kstart + k - 1)
1069 : END DO
1070 : END DO
1071 : END DO
1072 : CALL timeset(routineN//"_dgemm", handle_dgemm)
1073 : CALL dgemm("N", "N", c, n_loc_ri, jsize*ksize, &
1074 : 1.0_dp, rho_chunk, grid_chunk, &
1075 : int_2d_prv, max_ao_size*max_ao_size, &
1076 : 1.0_dp, d_lp_prv(l0, 1), n_grid_total)
1077 : CALL timestop(handle_dgemm)
1078 : END DO
1079 : END DO
1080 : END DO
1081 : !$OMP END DO
1082 :
1083 : !$OMP CRITICAL (compute_d_lp_reduce)
1084 : d_lp(1:n_grid_total, 1:n_loc_ri) = d_lp(1:n_grid_total, 1:n_loc_ri) + &
1085 : d_lp_prv(1:n_grid_total, 1:n_loc_ri)
1086 : !$OMP END CRITICAL (compute_d_lp_reduce)
1087 :
1088 : DEALLOCATE (int_3c_prv, int_3c_sum, int_2d_prv, rho_chunk, d_lp_prv)
1089 : CALL gw_3c_ws_release(ws)
1090 :
1091 : !$OMP END PARALLEL
1092 :
1093 0 : CALL timestop(handle)
1094 :
1095 0 : END SUBROUTINE compute_d_lp
1096 :
1097 : ! **************************************************************************************************
1098 : !> \brief Computes the χ(iτ, k=0) matrix
1099 : !> \param bs_env ...
1100 : !> \param mat_chi_Gamma_tau ...
1101 : !> \param mat_phi_mu_l ...
1102 : !> \param mat_Z_lP ...
1103 : ! **************************************************************************************************
1104 :
1105 0 : SUBROUTINE get_mat_chi_Gamma_tau(bs_env, mat_chi_Gamma_tau, mat_phi_mu_l, mat_Z_lP)
1106 :
1107 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1108 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1109 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
1110 :
1111 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
1112 :
1113 : INTEGER :: handle, i, i_t, ispin, npcol
1114 0 : INTEGER, DIMENSION(:), POINTER :: blk_ao, blk_grid, dist_col_grid, &
1115 0 : dist_row_grid
1116 : REAL(KIND=dp) :: t1, tau
1117 : TYPE(dbcsr_distribution_type) :: dist_grid_grid, dist_phi
1118 : TYPE(dbcsr_type) :: matrix_chi_grid, matrix_chi_grid_spin, &
1119 : matrix_G_occ_grid, matrix_G_vir_grid
1120 :
1121 0 : CALL timeset(routineN, handle)
1122 :
1123 : ! =========================================================================
1124 : ! 1. SETUP CORE TOPOLOGIES
1125 : ! =========================================================================
1126 0 : CALL dbcsr_get_info(mat_phi_mu_l, distribution=dist_phi, row_blk_size=blk_grid, col_blk_size=blk_ao)
1127 0 : CALL dbcsr_distribution_get(dist_phi, row_dist=dist_row_grid, npcols=npcol)
1128 :
1129 : ! Build a perfectly safe column distribution for the Grid dimension
1130 0 : ALLOCATE (dist_col_grid(SIZE(blk_grid)))
1131 0 : DO i = 1, SIZE(blk_grid)
1132 0 : dist_col_grid(i) = MOD(i - 1, npcol)
1133 : END DO
1134 :
1135 : CALL dbcsr_distribution_new(dist_grid_grid, template=dist_phi, &
1136 0 : row_dist=dist_row_grid, col_dist=dist_col_grid)
1137 :
1138 0 : CALL dbcsr_create(matrix_G_occ_grid, "G_occ_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1139 0 : CALL dbcsr_create(matrix_G_vir_grid, "G_vir_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1140 0 : CALL dbcsr_create(matrix_chi_grid, "chi_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1141 0 : CALL dbcsr_create(matrix_chi_grid_spin, "chi_grid_spin", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1142 :
1143 : ! =========================================================================
1144 : ! 2. MAIN IMAGINARY TIME LOOP
1145 : ! =========================================================================
1146 0 : DO i_t = 1, bs_env%num_time_freq_points
1147 0 : t1 = m_walltime()
1148 :
1149 0 : tau = bs_env%imag_time_points(i_t)
1150 0 : CALL dbcsr_set(matrix_chi_grid, 0.0_dp)
1151 :
1152 : ! ----------------------------------------------------------------------
1153 : ! A. SPIN LOOP (Allocations safely encapsulated in wrappers)
1154 : ! ----------------------------------------------------------------------
1155 0 : DO ispin = 1, bs_env%n_spin
1156 :
1157 : ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1158 : ! G^occ_ll'(i|τ|,k=0) = sum_µν Φ_µ(r_l) G^occ_µν Φ_ν(r_l')
1159 : CALL build_G_grid(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, &
1160 0 : matrix_G_occ_grid, bs_env%eps_filter)
1161 :
1162 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1163 : ! G^vir_ll'(i|τ|,k=0) = sum_µν Φ_µ(r_l) G^vir_µν Φ_ν(r_l')
1164 : CALL build_G_grid(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, &
1165 0 : matrix_G_vir_grid, bs_env%eps_filter)
1166 :
1167 : ! -------------------------------------------------------------------
1168 : ! B. ELEMENT-WISE HADAMARD PRODUCT
1169 : ! -------------------------------------------------------------------
1170 : ! χ_ll'(iτ,k=0) = G^occ_ll'(i|τ|,k=0) * G^vir_ll'(i|τ|,k=0)
1171 0 : CALL hadamard_product(matrix_G_occ_grid, matrix_G_vir_grid, matrix_chi_grid_spin, bs_env%spin_degeneracy)
1172 :
1173 : ! Accumulate spin contributions
1174 0 : CALL dbcsr_add(matrix_chi_grid, matrix_chi_grid_spin, 1.0_dp, 1.0_dp)
1175 :
1176 : END DO ! ispin
1177 :
1178 : ! ----------------------------------------------------------------------
1179 : ! C. TRANSFORM TO AUXILIARY BASIS & EXPORT DIRECTLY
1180 : ! χ_aux = Z^T * χ_grid * Z
1181 : ! χ_PQ(iτ,k=0) = sum_ll' Z_lP χ_ll'(iτ,k=0) Z_l'Q
1182 : ! Result is dumped directly into the final array mat_chi_Gamma_tau!
1183 : ! ----------------------------------------------------------------------
1184 : CALL contract_A_B_A("T", "N", mat_Z_lP, matrix_chi_grid, &
1185 0 : mat_chi_Gamma_tau(i_t)%matrix, bs_env%eps_filter)
1186 :
1187 0 : IF (bs_env%unit_nr > 0) THEN
1188 : WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
1189 0 : 'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
1190 0 : ', Execution time', m_walltime() - t1, ' s'
1191 : END IF
1192 :
1193 : END DO ! i_t
1194 :
1195 : ! =========================================================================
1196 : ! 3. FINAL CLEANUP
1197 : ! =========================================================================
1198 0 : CALL dbcsr_release(matrix_G_occ_grid)
1199 0 : CALL dbcsr_release(matrix_G_vir_grid)
1200 0 : CALL dbcsr_release(matrix_chi_grid)
1201 0 : CALL dbcsr_release(matrix_chi_grid_spin)
1202 0 : CALL dbcsr_distribution_release(dist_grid_grid)
1203 0 : DEALLOCATE (dist_col_grid)
1204 :
1205 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1206 :
1207 0 : CALL timestop(handle)
1208 :
1209 0 : END SUBROUTINE get_mat_chi_Gamma_tau
1210 :
1211 : ! **************************************************************************************************
1212 : !> \brief Computes Green's Function in grid basis
1213 : !> \param bs_env ...
1214 : !> \param tau ...
1215 : !> \param ispin ...
1216 : !> \param occ ...
1217 : !> \param vir ...
1218 : !> \param mat_phi_mu_l ...
1219 : !> \param matrix_G_grid ...
1220 : !> \param eps_filter ...
1221 : ! **************************************************************************************************
1222 :
1223 0 : SUBROUTINE build_G_grid(bs_env, tau, ispin, occ, vir, mat_phi_mu_l, matrix_G_grid, eps_filter)
1224 :
1225 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1226 : REAL(KIND=dp), INTENT(IN) :: tau
1227 : INTEGER, INTENT(IN) :: ispin
1228 : LOGICAL, INTENT(IN) :: occ, vir
1229 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, matrix_G_grid
1230 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1231 :
1232 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_G_grid'
1233 :
1234 : INTEGER :: handle
1235 0 : INTEGER, DIMENSION(:), POINTER :: blk_ao, dist_row_ao
1236 : TYPE(cp_fm_type), POINTER :: fm_G
1237 : TYPE(dbcsr_distribution_type) :: dist_ao_ao
1238 : TYPE(dbcsr_type) :: matrix_G_ao
1239 :
1240 0 : CALL timeset(routineN, handle)
1241 :
1242 : ! 1. Select the correct FM matrix based on occ/vir flags
1243 0 : IF (occ) THEN
1244 0 : fm_G => bs_env%fm_Gocc
1245 : ELSE
1246 0 : fm_G => bs_env%fm_Gvir
1247 : END IF
1248 :
1249 : ! 2. Compute Dense FM Green's Function
1250 : ! G^occ/vir_µλ(i|τ|,k=0) = sum_G^occ/vir_µλn^occ/vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1251 0 : CALL G_occ_vir(bs_env, tau, fm_G, ispin, occ=occ, vir=vir)
1252 :
1253 : ! 3. Setup AO DBCSR Topology and Create Matrix dynamically
1254 0 : CALL setup_square_topology(mat_phi_mu_l, 'COL', dist_ao_ao, blk_ao, dist_row_ao)
1255 :
1256 : CALL dbcsr_create(matrix_G_ao, name="G_ao", dist=dist_ao_ao, &
1257 : matrix_type=dbcsr_type_no_symmetry, &
1258 0 : row_blk_size=blk_ao, col_blk_size=blk_ao)
1259 :
1260 : ! 4. Convert FM to Sparse DBCSR
1261 0 : CALL copy_fm_to_dbcsr(fm_G, matrix_G_ao, keep_sparsity=.FALSE.)
1262 :
1263 : ! 5. Transform to Grid Basis: G_grid = phi * G_ao * phi^T
1264 : ! G^occ/vir_ll'(i|τ|,k=0) = sum_µν Φ_µ(r_l) G^occ/vir_µν Φ_ν(r_l')
1265 0 : CALL contract_A_B_A("N", "T", mat_phi_mu_l, matrix_G_ao, matrix_G_grid, eps_filter)
1266 :
1267 : ! 6. Release AO matrix and topology
1268 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_ao_ao, mapped_dist=dist_row_ao, m1=matrix_G_ao)
1269 :
1270 0 : CALL timestop(handle)
1271 :
1272 0 : END SUBROUTINE build_G_grid
1273 :
1274 : ! **************************************************************************************************
1275 : !> \brief Generalized routine to compute OUT = A * B * A^T OR OUT = A^T * B * A using DBCSR
1276 : !> \param transA_left ...
1277 : !> \param transA_right ...
1278 : !> \param matrix_A ...
1279 : !> \param matrix_B ...
1280 : !> \param matrix_out ...
1281 : !> \param eps_filter ...
1282 : ! **************************************************************************************************
1283 :
1284 7468 : SUBROUTINE contract_A_B_A(transA_left, transA_right, matrix_A, matrix_B, matrix_out, eps_filter)
1285 :
1286 : CHARACTER(LEN=1), INTENT(IN) :: transA_left, transA_right
1287 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B, matrix_out
1288 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1289 :
1290 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_A_B_A'
1291 :
1292 : INTEGER :: handle
1293 : TYPE(dbcsr_type) :: matrix_tmp
1294 :
1295 7468 : CALL timeset(routineN, handle)
1296 :
1297 7468 : CALL dbcsr_create(matrix_tmp, template=matrix_A)
1298 :
1299 7468 : IF (transA_left == "N" .AND. transA_right == "T") THEN
1300 : ! Path 1: Out = A * B * A^T
1301 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_A, matrix_B, &
1302 3964 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
1303 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, matrix_A, &
1304 3964 : 0.0_dp, matrix_out, filter_eps=eps_filter)
1305 :
1306 3504 : ELSE IF (transA_left == "T" .AND. transA_right == "N") THEN
1307 : ! Path 2: Out = A^T * B * A
1308 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_B, matrix_A, &
1309 3504 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
1310 : CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_A, matrix_tmp, &
1311 3504 : 0.0_dp, matrix_out, filter_eps=eps_filter)
1312 : ELSE
1313 0 : CPABORT("Unsupported transposition pair in contract_A_B_A")
1314 : END IF
1315 :
1316 7468 : CALL dbcsr_release(matrix_tmp)
1317 :
1318 7468 : CALL timestop(handle)
1319 :
1320 7468 : END SUBROUTINE contract_A_B_A
1321 :
1322 : ! **************************************************************************************************
1323 : !> \brief Computes C = A ◦ B (Element-wise Hadamard product) for sparse DBCSR matrices.
1324 : !> \param matrix_A ...
1325 : !> \param matrix_B ...
1326 : !> \param matrix_C ...
1327 : !> \param fac (Scaling factor applied to the product)
1328 : ! **************************************************************************************************
1329 :
1330 0 : SUBROUTINE hadamard_product(matrix_A, matrix_B, matrix_C, fac)
1331 :
1332 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B, matrix_C
1333 : REAL(KIND=dp), INTENT(IN) :: fac
1334 :
1335 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hadamard_product'
1336 :
1337 : INTEGER :: col, handle, row
1338 : LOGICAL :: found
1339 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk_B, blk_C
1340 : TYPE(dbcsr_iterator_type) :: iter
1341 :
1342 0 : CALL timeset(routineN, handle)
1343 :
1344 0 : CALL dbcsr_copy(matrix_C, matrix_A)
1345 :
1346 0 : CALL dbcsr_iterator_start(iter, matrix_C)
1347 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1348 0 : CALL dbcsr_iterator_next_block(iter, row, col, blk_C)
1349 :
1350 0 : CALL dbcsr_get_block_p(matrix_B, row, col, blk_B, found)
1351 :
1352 0 : IF (found) THEN
1353 0 : blk_C(:, :) = fac*blk_C(:, :)*blk_B(:, :)
1354 : ELSE
1355 : ! If B is sparse here, the product is zero
1356 0 : blk_C(:, :) = 0.0_dp
1357 : END IF
1358 : END DO
1359 0 : CALL dbcsr_iterator_stop(iter)
1360 :
1361 0 : CALL timestop(handle)
1362 :
1363 0 : END SUBROUTINE hadamard_product
1364 :
1365 : ! **************************************************************************************************
1366 : !> \brief In-place Hadamard A <- fac * (A ◦ B). Value mutation only (no block insert/delete),
1367 : !> so iterating A while writing through the block pointer is safe.
1368 : !> \param matrix_A in/out factor (overwritten by the product)
1369 : !> \param matrix_B second factor (looked up; blocks absent in B zero the A block)
1370 : !> \param fac (Scaling factor applied to the product)
1371 : ! **************************************************************************************************
1372 7008 : SUBROUTINE hadamard_product_inplace(matrix_A, matrix_B, fac)
1373 :
1374 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_A, matrix_B
1375 : REAL(KIND=dp), INTENT(IN) :: fac
1376 :
1377 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hadamard_product_inplace'
1378 :
1379 : INTEGER :: col, handle, row
1380 : LOGICAL :: found
1381 3504 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: blk_A, blk_B
1382 : TYPE(dbcsr_iterator_type) :: iter
1383 :
1384 3504 : CALL timeset(routineN, handle)
1385 :
1386 3504 : CALL dbcsr_iterator_start(iter, matrix_A)
1387 214608 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1388 211104 : CALL dbcsr_iterator_next_block(iter, row, col, blk_A)
1389 :
1390 211104 : CALL dbcsr_get_block_p(matrix_B, row, col, blk_B, found)
1391 :
1392 214608 : IF (found) THEN
1393 332612544 : blk_A(:, :) = fac*blk_A(:, :)*blk_B(:, :)
1394 : ELSE
1395 : ! If B is sparse here, the product is zero
1396 0 : blk_A(:, :) = 0.0_dp
1397 : END IF
1398 : END DO
1399 3504 : CALL dbcsr_iterator_stop(iter)
1400 :
1401 3504 : CALL timestop(handle)
1402 :
1403 3504 : END SUBROUTINE hadamard_product_inplace
1404 :
1405 : ! **************************************************************************************************
1406 : !> \brief Compute screened Coulomb interaction matrix
1407 : !> \param bs_env ...
1408 : !> \param qs_env ...
1409 : !> \param mat_chi_Gamma_tau ...
1410 : !> \param fm_W_time ...
1411 : ! **************************************************************************************************
1412 :
1413 0 : SUBROUTINE compute_W(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_time)
1414 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1415 : TYPE(qs_environment_type), POINTER :: qs_env
1416 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
1417 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
1418 :
1419 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_W'
1420 :
1421 : INTEGER :: handle, i_t, j_w
1422 : REAL(KIND=dp) :: t1
1423 : TYPE(cp_fm_type) :: fm_M_inv_V_sqrt, fm_V, fm_V_sqrt
1424 :
1425 0 : CALL timeset(routineN, handle)
1426 :
1427 0 : t1 = m_walltime()
1428 :
1429 0 : CALL create_fm_W_MIC_time(bs_env, fm_W_time)
1430 :
1431 : ! 1. Allocate V and M matrices
1432 0 : CALL cp_fm_create(fm_V, bs_env%fm_RI_RI%matrix_struct)
1433 0 : CALL cp_fm_create(fm_V_sqrt, bs_env%fm_RI_RI%matrix_struct)
1434 0 : CALL cp_fm_create(fm_M_inv_V_sqrt, bs_env%fm_RI_RI%matrix_struct)
1435 :
1436 : ! Compute V and M^-1 * V^0.5
1437 0 : CALL compute_V_MinvVsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_M_inv_V_sqrt)
1438 :
1439 : ! 2. Loop over frequencies
1440 0 : DO j_w = 1, bs_env%num_time_freq_points
1441 : ! Fourier transformation of χ_PQ(iτ) to χ_PQ(iω_j)
1442 0 : CALL compute_fm_chi_Gamma_freq(bs_env, bs_env%fm_chi_Gamma_freq, j_w, mat_chi_Gamma_tau)
1443 :
1444 : ! ε(iω_j) = Id - V^0.5*M^-1*χ(iω_j)*M^-1*V^0.5
1445 : ! W(iω_j) = V^0.5*(ε^-1(iω_j)-Id)*V^0.5
1446 : CALL compute_fm_W_freq(bs_env, bs_env%fm_chi_Gamma_freq, fm_V_sqrt, &
1447 0 : fm_M_inv_V_sqrt, bs_env%fm_W_MIC_freq)
1448 :
1449 : ! Fourier transform from W_PQ^MIC(iω_j) to W_PQ^MIC(iτ)
1450 0 : CALL Fourier_transform_w_to_t(bs_env, fm_W_time, bs_env%fm_W_MIC_freq, j_w)
1451 : END DO
1452 :
1453 : ! M^-1*W^MIC(iτ)*M^-1
1454 0 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_W_time)
1455 :
1456 0 : IF (bs_env%unit_nr > 0) THEN
1457 : WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
1458 0 : 'Computed W(iτ),', ' Execution time', m_walltime() - t1, ' s'
1459 : END IF
1460 :
1461 0 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
1462 :
1463 : ! Cleanup
1464 0 : CALL cp_fm_release(fm_V)
1465 0 : CALL cp_fm_release(fm_V_sqrt)
1466 0 : CALL cp_fm_release(fm_M_inv_V_sqrt)
1467 :
1468 : ! Marek : Fourier transform W^MIC(itau) back to get it at a specific im.frequency point - iomega = 0
1469 0 : IF (bs_env%rtp_method == rtp_method_bse) THEN
1470 0 : t1 = m_walltime()
1471 0 : CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1472 : ! Set to zero
1473 0 : CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
1474 : ! Sum over all times
1475 0 : DO i_t = 1, bs_env%num_time_freq_points
1476 : ! Add the relevant structure with correct weight
1477 : CALL cp_fm_scale_and_add(1.0_dp, bs_env%fm_W_MIC_freq_zero, &
1478 0 : bs_env%imag_time_weights_freq_zero(i_t), fm_W_time(i_t))
1479 : END DO
1480 : ! Done, save to file
1481 0 : CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
1482 : ! Report calculation
1483 0 : IF (bs_env%unit_nr > 0) THEN
1484 : WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
1485 0 : 'Computed W(0),', ' Execution time', m_walltime() - t1, ' s'
1486 : END IF
1487 : END IF
1488 :
1489 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1490 :
1491 0 : CALL timestop(handle)
1492 :
1493 0 : END SUBROUTINE compute_W
1494 :
1495 : ! **************************************************************************************************
1496 : !> \brief Computes V, V^0.5, and M^-1 * V^0.5
1497 : !> \param bs_env ...
1498 : !> \param qs_env ...
1499 : !> \param fm_V ...
1500 : !> \param fm_V_sqrt ...
1501 : !> \param fm_Minv_Vsqrt ...
1502 : ! **************************************************************************************************
1503 :
1504 0 : SUBROUTINE compute_V_MinvVsqrt(bs_env, qs_env, fm_V, fm_V_sqrt, fm_Minv_Vsqrt)
1505 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1506 : TYPE(qs_environment_type), POINTER :: qs_env
1507 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_V, fm_V_sqrt, fm_Minv_Vsqrt
1508 :
1509 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_MinvVsqrt'
1510 :
1511 : INTEGER :: handle, info, n_RI, ndep
1512 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1513 : TYPE(cell_type), POINTER :: cell
1514 : TYPE(cp_fm_type) :: fm_work
1515 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_M
1516 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_V_kp
1517 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1518 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1519 :
1520 0 : CALL timeset(routineN, handle)
1521 :
1522 0 : n_RI = bs_env%n_RI
1523 0 : CALL cp_fm_create(fm_work, fm_V%matrix_struct)
1524 :
1525 : ! -----------------------------------------------------------------------
1526 : ! 1. Build Coulomb Matrix V(k=0) using the kp-routine but only for ikp=1
1527 : ! -----------------------------------------------------------------------
1528 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, cell=cell, &
1529 0 : qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
1530 :
1531 0 : ALLOCATE (mat_V_kp(1:1, 1:2))
1532 0 : NULLIFY (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
1533 0 : ALLOCATE (mat_V_kp(1, 1)%matrix, mat_V_kp(1, 2)%matrix)
1534 :
1535 0 : CALL dbcsr_create(mat_V_kp(1, 1)%matrix, template=bs_env%mat_RI_RI%matrix)
1536 0 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 1)%matrix)
1537 0 : CALL dbcsr_set(mat_V_kp(1, 1)%matrix, 0.0_dp)
1538 :
1539 : ! Dummy imaginary part just to satisfy the routine
1540 0 : CALL dbcsr_create(mat_V_kp(1, 2)%matrix, template=bs_env%mat_RI_RI%matrix)
1541 0 : CALL dbcsr_reserve_all_blocks(mat_V_kp(1, 2)%matrix)
1542 0 : CALL dbcsr_set(mat_V_kp(1, 2)%matrix, 0.0_dp)
1543 :
1544 0 : bs_env%kpoints_chi_eps_W%nkp_grid = bs_env%nkp_grid_chi_eps_W_orig
1545 :
1546 : CALL build_2c_coulomb_matrix_kp(mat_V_kp, bs_env%kpoints_chi_eps_W, "RI_AUX", cell, &
1547 : particle_set, qs_kind_set, atomic_kind_set, &
1548 0 : bs_env%size_lattice_sum_V, operator_coulomb, 1, 1)
1549 :
1550 : ! Copy real part to fm_V
1551 0 : CALL copy_dbcsr_to_fm(mat_V_kp(1, 1)%matrix, fm_V)
1552 :
1553 0 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 1)%matrix)
1554 0 : CALL dbcsr_deallocate_matrix(mat_V_kp(1, 2)%matrix)
1555 0 : DEALLOCATE (mat_V_kp)
1556 :
1557 : ! -----------------------------------------------------------------------
1558 : ! 2. Get RI-Metric Matrix M(k=0)
1559 : ! -----------------------------------------------------------------------
1560 : CALL RI_2c_integral_mat(qs_env, fm_M, fm_V, n_RI, bs_env%ri_metric, &
1561 0 : do_kpoints=.FALSE., regularization_RI=bs_env%regularization_RI)
1562 :
1563 : ! -----------------------------------------------------------------------
1564 : ! 3. M -> M^-1
1565 : ! -----------------------------------------------------------------------
1566 0 : CALL cp_fm_cholesky_decompose(fm_M(1, 1), info_out=info)
1567 0 : IF (info == 0) THEN
1568 0 : CALL cp_fm_cholesky_invert(fm_M(1, 1))
1569 0 : CALL cp_fm_uplo_to_full(fm_M(1, 1), fm_work)
1570 : ELSE
1571 : ! Fallback if Cholesky fails due to conditioning
1572 0 : CALL cp_fm_power(fm_M(1, 1), fm_work, -1.0_dp, bs_env%eps_eigval_mat_RI, ndep)
1573 0 : CALL cp_fm_to_fm(fm_work, fm_M(1, 1))
1574 : END IF
1575 :
1576 : ! -----------------------------------------------------------------------
1577 : ! 4. V -> V^0.5
1578 : ! -----------------------------------------------------------------------
1579 0 : CALL cp_fm_to_fm(fm_V, fm_V_sqrt)
1580 0 : CALL cp_fm_cholesky_decompose(fm_V_sqrt, info_out=info)
1581 0 : IF (info == 0) THEN
1582 0 : CALL clean_lower_part(fm_V_sqrt)
1583 : ELSE
1584 0 : CALL cp_fm_power(fm_V, fm_V_sqrt, 0.5_dp, bs_env%eps_eigval_mat_RI, ndep)
1585 : END IF
1586 :
1587 : ! -----------------------------------------------------------------------
1588 : ! 5. M^-1 * V^0.5
1589 : ! -----------------------------------------------------------------------
1590 : CALL parallel_gemm("N", "T", n_RI, n_RI, n_RI, 1.0_dp, fm_M(1, 1), fm_V_sqrt, &
1591 0 : 0.0_dp, fm_Minv_Vsqrt)
1592 :
1593 0 : CALL cp_fm_release(fm_M)
1594 0 : CALL cp_fm_release(fm_work)
1595 :
1596 0 : CALL timestop(handle)
1597 :
1598 0 : END SUBROUTINE compute_V_MinvVsqrt
1599 :
1600 : ! **************************************************************************************************
1601 : !> \brief Computes W(iω) from χ_PQ(iω_j)
1602 : !> \param bs_env ...
1603 : !> \param fm_chi_freq_j ...
1604 : !> \param fm_V_sqrt ...
1605 : !> \param fm_Minv_Vsqrt ...
1606 : !> \param fm_W_freq_j ...
1607 : ! **************************************************************************************************
1608 :
1609 0 : SUBROUTINE compute_fm_W_freq(bs_env, fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt, fm_W_freq_j)
1610 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1611 : TYPE(cp_fm_type), INTENT(IN) :: fm_chi_freq_j, fm_V_sqrt, fm_Minv_Vsqrt
1612 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_W_freq_j
1613 :
1614 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_fm_W_freq'
1615 :
1616 : INTEGER :: handle, info, n_RI, ndep
1617 : TYPE(cp_fm_type) :: fm_eps_freq_j, fm_work
1618 :
1619 0 : CALL timeset(routineN, handle)
1620 :
1621 0 : n_RI = bs_env%n_RI
1622 :
1623 0 : CALL cp_fm_create(fm_eps_freq_j, fm_chi_freq_j%matrix_struct)
1624 0 : CALL cp_fm_create(fm_work, fm_chi_freq_j%matrix_struct)
1625 :
1626 : ! -----------------------------------------------------------------------
1627 : ! 1. ε(iω_j) = Id - (M^-1 * V^0.5)^T * χ(iω_j) * (M^-1 * V^0.5)
1628 : ! -----------------------------------------------------------------------
1629 : ! work = χ(iω_j) * (M^-1 * V^0.5)
1630 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
1631 0 : fm_chi_freq_j, fm_Minv_Vsqrt, 0.0_dp, fm_work)
1632 :
1633 : ! eps_work = (M^-1 * V^0.5)^T * work
1634 : CALL parallel_gemm('T', 'N', n_RI, n_RI, n_RI, 1.0_dp, &
1635 0 : fm_Minv_Vsqrt, fm_work, 0.0_dp, fm_eps_freq_j)
1636 :
1637 : ! ε(iω_j) = Id - eps_work --> -eps_work + Id
1638 0 : CALL fm_add_on_diag(fm_eps_freq_j, 1.0_dp)
1639 :
1640 : ! Force perfect symmetry before Cholesky to avoid info != 0 due to GEMM noise
1641 0 : CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
1642 :
1643 : ! -----------------------------------------------------------------------
1644 : ! 2. W(iω_j) = V^0.5^T * (ε^-1(iω_j) - Id) * V^0.5
1645 : ! -----------------------------------------------------------------------
1646 :
1647 : ! a) Cholesky decomposition of ε(iω_j)
1648 0 : CALL cp_fm_cholesky_decompose(fm_eps_freq_j, info_out=info)
1649 :
1650 : ! b) Inversion
1651 0 : IF (info == 0) THEN
1652 0 : CALL cp_fm_cholesky_invert(fm_eps_freq_j)
1653 0 : CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
1654 : ELSE
1655 : ! Fallback to expensive diagonalization if Cholesky fails
1656 0 : CALL cp_fm_power(fm_eps_freq_j, fm_work, -1.0_dp, bs_env%eps_eigval_mat_RI, ndep)
1657 0 : CALL cp_fm_to_fm(fm_work, fm_eps_freq_j)
1658 : END IF
1659 :
1660 : ! c) ε^-1(iω_j) - Id
1661 0 : CALL fm_add_on_diag(fm_eps_freq_j, -1.0_dp)
1662 :
1663 : ! d) work = (ε^-1(iω_j) - Id) * V^0.5
1664 : CALL parallel_gemm('N', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_eps_freq_j, fm_V_sqrt, &
1665 0 : 0.0_dp, fm_work)
1666 :
1667 : ! e) W(iw) = V^0.5^T * work
1668 : CALL parallel_gemm('T', 'N', n_RI, n_RI, n_RI, 1.0_dp, fm_V_sqrt, fm_work, &
1669 0 : 0.0_dp, fm_W_freq_j)
1670 :
1671 : ! Cleanup
1672 0 : CALL cp_fm_release(fm_work)
1673 0 : CALL cp_fm_release(fm_eps_freq_j)
1674 :
1675 0 : CALL timestop(handle)
1676 :
1677 0 : END SUBROUTINE compute_fm_W_freq
1678 :
1679 : ! **************************************************************************************************
1680 : !> \brief Adds a real scalar value to the diagonal of a real full matrix (fm)
1681 : !> \param fm ...
1682 : !> \param alpha ...
1683 : ! **************************************************************************************************
1684 :
1685 0 : SUBROUTINE fm_add_on_diag(fm, alpha)
1686 : TYPE(cp_fm_type), INTENT(INOUT) :: fm
1687 : REAL(KIND=dp), INTENT(IN) :: alpha
1688 :
1689 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_add_on_diag'
1690 :
1691 : INTEGER :: handle, i_global, i_row, j_col, &
1692 : j_global, ncol_local, nrow_local
1693 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1694 :
1695 0 : CALL timeset(routineN, handle)
1696 :
1697 : CALL cp_fm_get_info(matrix=fm, &
1698 : nrow_local=nrow_local, &
1699 : ncol_local=ncol_local, &
1700 : row_indices=row_indices, &
1701 0 : col_indices=col_indices)
1702 :
1703 0 : DO j_col = 1, ncol_local
1704 0 : j_global = col_indices(j_col)
1705 0 : DO i_row = 1, nrow_local
1706 0 : i_global = row_indices(i_row)
1707 0 : IF (j_global == i_global) THEN
1708 0 : fm%local_data(i_row, j_col) = fm%local_data(i_row, j_col) + alpha
1709 : END IF
1710 : END DO
1711 : END DO
1712 :
1713 0 : CALL timestop(handle)
1714 :
1715 0 : END SUBROUTINE fm_add_on_diag
1716 :
1717 : ! **************************************************************************************************
1718 : !> \brief Zeroes out the strictly lower triangular part of a real matrix
1719 : !> \param fm_mat ...
1720 : ! **************************************************************************************************
1721 0 : SUBROUTINE clean_lower_part(fm_mat)
1722 : TYPE(cp_fm_type) :: fm_mat
1723 :
1724 : CHARACTER(LEN=*), PARAMETER :: routineN = 'clean_lower_part'
1725 :
1726 : INTEGER :: handle, i_row, j_col, j_global, &
1727 : ncol_local, nrow_local
1728 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1729 :
1730 0 : CALL timeset(routineN, handle)
1731 :
1732 : CALL cp_fm_get_info(matrix=fm_mat, &
1733 : nrow_local=nrow_local, ncol_local=ncol_local, &
1734 0 : row_indices=row_indices, col_indices=col_indices)
1735 :
1736 0 : DO j_col = 1, ncol_local
1737 0 : j_global = col_indices(j_col)
1738 0 : DO i_row = 1, nrow_local
1739 0 : IF (j_global < row_indices(i_row)) fm_mat%local_data(i_row, j_col) = 0.0_dp
1740 : END DO
1741 : END DO
1742 :
1743 0 : CALL timestop(handle)
1744 :
1745 0 : END SUBROUTINE clean_lower_part
1746 :
1747 : ! **************************************************************************************************
1748 : !> \brief Computes the exact exchange part of the GW self-energy
1749 : !> \param bs_env ...
1750 : !> \param qs_env ...
1751 : !> \param mat_phi_mu_l ...
1752 : !> \param mat_Z_lP ...
1753 : !> \param fm_Sigma_x_Gamma ...
1754 : ! **************************************************************************************************
1755 :
1756 0 : SUBROUTINE compute_Sigma_x(bs_env, qs_env, mat_phi_mu_l, mat_Z_lP, fm_Sigma_x_Gamma)
1757 :
1758 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1759 : TYPE(qs_environment_type), POINTER :: qs_env
1760 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
1761 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma
1762 :
1763 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_x'
1764 :
1765 : INTEGER :: handle, ispin
1766 0 : INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
1767 0 : dist_row_aux
1768 : REAL(KIND=dp) :: t1
1769 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_Vtr_Gamma
1770 : TYPE(dbcsr_distribution_type) :: dist_aux_aux, dist_grid_grid
1771 : TYPE(dbcsr_type) :: mat_Sigma_x_Gamma, matrix_D_grid, &
1772 : matrix_Sigma_x_grid, matrix_V_aux, &
1773 : matrix_V_grid
1774 :
1775 0 : CALL timeset(routineN, handle)
1776 :
1777 0 : t1 = m_walltime()
1778 :
1779 0 : ALLOCATE (fm_Sigma_x_Gamma(bs_env%n_spin))
1780 0 : DO ispin = 1, bs_env%n_spin
1781 0 : CALL cp_fm_create(fm_Sigma_x_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1782 : END DO
1783 :
1784 0 : CALL dbcsr_create(mat_Sigma_x_Gamma, template=bs_env%mat_ao_ao%matrix)
1785 :
1786 : ! =========================================================================
1787 : ! 1. SETUP CORE TOPOLOGIES
1788 : ! =========================================================================
1789 0 : CALL setup_square_topology(mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, dist_col_grid)
1790 0 : CALL setup_square_topology(mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
1791 :
1792 : ! =========================================================================
1793 : ! 2. COMPUTE V^tr_ll'
1794 : ! =========================================================================
1795 : CALL RI_2c_integral_mat(qs_env, fm_Vtr_Gamma, bs_env%fm_RI_RI, bs_env%n_RI, &
1796 0 : bs_env%trunc_coulomb, do_kpoints=.FALSE.)
1797 :
1798 : ! M^-1 * V^tr * M^-1 directly modifies fm_Vtr_Gamma(:, 1)
1799 0 : CALL multiply_fm_W_MIC_time_with_Minv_Gamma(bs_env, qs_env, fm_Vtr_Gamma(:, 1))
1800 :
1801 0 : CALL dbcsr_create(matrix_V_aux, "V_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1802 0 : CALL dbcsr_create(matrix_V_grid, "V_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1803 :
1804 0 : CALL copy_fm_to_dbcsr(fm_Vtr_Gamma(1, 1), matrix_V_aux, keep_sparsity=.FALSE.)
1805 :
1806 : ! V^tr_ll' = sum_PQ Z_lP V^trunc_PQ Z_l'Q
1807 0 : CALL contract_A_B_A("N", "T", mat_Z_lP, matrix_V_aux, matrix_V_grid, bs_env%eps_filter)
1808 0 : CALL dbcsr_release(matrix_V_aux)
1809 :
1810 : ! =========================================================================
1811 : ! 3. SPIN LOOP FOR EXACT EXCHANGE
1812 : ! =========================================================================
1813 0 : DO ispin = 1, bs_env%n_spin
1814 :
1815 : ! Density matrix on grid is essentially G_occ at tau = 0.0
1816 : ! D_µν = sum_n^occ C_µn(k=0) C_νn(k=0)
1817 : ! D_ll' = sum_µν Φ_µ(r_l) D_µν Φ_ν(r_l')
1818 0 : CALL dbcsr_create(matrix_D_grid, "D_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1819 0 : CALL build_G_grid(bs_env, 0.0_dp, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_D_grid, bs_env%eps_filter)
1820 :
1821 : ! Element-wise Hadamard product: Σ^x_grid = D_grid ◦ V_grid
1822 : ! Σ^x_ll' = D_ll' * V^tr_ll'
1823 0 : CALL dbcsr_create(matrix_Sigma_x_grid, template=matrix_V_grid)
1824 0 : CALL hadamard_product(matrix_D_grid, matrix_V_grid, matrix_Sigma_x_grid, 1.0_dp)
1825 :
1826 0 : CALL dbcsr_release(matrix_D_grid)
1827 :
1828 : ! Transform back to AO basis: Σ^x_ao = -1.0 * phi^T * Σ^x_grid * phi
1829 : ! Σ^x_λσ(k=0) = -sum_ll' Φ_λ(r_l) Σ^x_ll' Φ_σ(r_l')
1830 0 : CALL contract_A_B_A("T", "N", mat_phi_mu_l, matrix_Sigma_x_grid, mat_Sigma_x_Gamma, bs_env%eps_filter)
1831 0 : CALL dbcsr_scale(mat_Sigma_x_Gamma, -1.0_dp)
1832 :
1833 0 : CALL dbcsr_release(matrix_Sigma_x_grid)
1834 :
1835 : ! Data I/O and Export to CP2K Full Matrices
1836 0 : CALL copy_dbcsr_to_fm(mat_Sigma_x_Gamma, fm_Sigma_x_Gamma(ispin))
1837 :
1838 : END DO ! ispin
1839 :
1840 0 : IF (bs_env%unit_nr > 0) THEN
1841 : WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
1842 0 : 'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
1843 0 : WRITE (bs_env%unit_nr, '(A)') ' '
1844 : END IF
1845 :
1846 : ! =========================================================================
1847 : ! 4. CLEANUP
1848 : ! =========================================================================
1849 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid, &
1850 0 : m1=mat_Sigma_x_Gamma, m2=matrix_V_grid)
1851 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
1852 :
1853 0 : CALL cp_fm_release(fm_Vtr_Gamma)
1854 :
1855 0 : CALL timestop(handle)
1856 :
1857 0 : END SUBROUTINE compute_Sigma_x
1858 :
1859 : ! **************************************************************************************************
1860 : !> \brief Computes the correlation part of the GW self-energy
1861 : !> \param bs_env ...
1862 : !> \param fm_W_time ...
1863 : !> \param mat_phi_mu_l ...
1864 : !> \param mat_Z_lP ...
1865 : !> \param fm_Sigma_c_Gamma_time ...
1866 : ! **************************************************************************************************
1867 :
1868 0 : SUBROUTINE compute_Sigma_c(bs_env, fm_W_time, mat_phi_mu_l, mat_Z_lP, fm_Sigma_c_Gamma_time)
1869 :
1870 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
1871 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_time
1872 : TYPE(dbcsr_type), INTENT(INOUT) :: mat_phi_mu_l, mat_Z_lP
1873 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_Sigma_c_Gamma_time
1874 :
1875 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_Sigma_c'
1876 :
1877 : INTEGER :: handle, i_t, ispin
1878 0 : INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
1879 0 : dist_row_aux
1880 : REAL(KIND=dp) :: t1, tau
1881 : TYPE(dbcsr_distribution_type) :: dist_aux_aux, dist_grid_grid
1882 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg_tau, mat_Sigma_pos_tau
1883 : TYPE(dbcsr_type) :: matrix_G_occ_grid, matrix_G_vir_grid, matrix_Sigma_neg_grid, &
1884 : matrix_Sigma_pos_grid, matrix_W_aux, matrix_W_grid
1885 :
1886 0 : CALL timeset(routineN, handle)
1887 :
1888 : ! =========================================================================
1889 : ! 1. SETUP CORE TOPOLOGIES AND PRE-ALLOCATE OUTPUT ARRAYS
1890 : ! =========================================================================
1891 0 : CALL setup_square_topology(mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, dist_col_grid)
1892 0 : CALL setup_square_topology(mat_Z_lP, 'COL', dist_aux_aux, blk_aux, dist_row_aux)
1893 :
1894 : ! Pre-allocate local DBCSR matrices to act as targets for final output
1895 0 : NULLIFY (mat_Sigma_neg_tau, mat_Sigma_pos_tau)
1896 0 : ALLOCATE (mat_Sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1897 0 : ALLOCATE (mat_Sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1898 :
1899 0 : DO i_t = 1, bs_env%num_time_freq_points
1900 0 : DO ispin = 1, bs_env%n_spin
1901 0 : ALLOCATE (mat_Sigma_neg_tau(i_t, ispin)%matrix)
1902 0 : ALLOCATE (mat_Sigma_pos_tau(i_t, ispin)%matrix)
1903 0 : CALL dbcsr_create(mat_Sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1904 0 : CALL dbcsr_create(mat_Sigma_pos_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1905 : END DO
1906 : END DO
1907 :
1908 : ! =========================================================================
1909 : ! 2. MAIN IMAGINARY TIME LOOP
1910 : ! =========================================================================
1911 0 : DO i_t = 1, bs_env%num_time_freq_points
1912 0 : tau = bs_env%imag_time_points(i_t)
1913 :
1914 : ! -------------------------------------------------------------------
1915 : ! Compute W_grid = Z * W_aux * Z^T
1916 : ! -------------------------------------------------------------------
1917 0 : CALL dbcsr_create(matrix_W_aux, "W_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1918 0 : CALL dbcsr_create(matrix_W_grid, "W_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1919 :
1920 0 : CALL copy_fm_to_dbcsr(fm_W_time(i_t), matrix_W_aux, keep_sparsity=.FALSE.)
1921 :
1922 : ! W^MIC_ll'(iτ,k=0) = sum_PQ Z_lP W^MIC_PQ(iτ) Z_l'Q
1923 0 : CALL contract_A_B_A("N", "T", mat_Z_lP, matrix_W_aux, matrix_W_grid, bs_env%eps_filter)
1924 :
1925 0 : CALL dbcsr_release(matrix_W_aux) ! Clean up aux basis immediately
1926 :
1927 0 : DO ispin = 1, bs_env%n_spin
1928 0 : t1 = m_walltime()
1929 :
1930 : ! -------------------------------------------------------------------
1931 : ! A. Transform Green's Functions to the Grid
1932 : ! -------------------------------------------------------------------
1933 0 : CALL dbcsr_create(matrix_G_occ_grid, "G_occ_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1934 0 : CALL dbcsr_create(matrix_G_vir_grid, "G_vir_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1935 :
1936 : ! G^occ_µλ(i|τ|,k=0) = sum_G^occ_µλn^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1937 : ! G^occ_ll'(i|τ|,k=0) = sum_µν Φ_µ(r_l) G^occ_µν Φ_ν(r_l')
1938 0 : CALL build_G_grid(bs_env, tau, ispin, .TRUE., .FALSE., mat_phi_mu_l, matrix_G_occ_grid, bs_env%eps_filter)
1939 :
1940 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
1941 : ! G^vir_ll'(i|τ|,k=0) = sum_µν Φ_µ(r_l) G^vir_µν Φ_ν(r_l')
1942 0 : CALL build_G_grid(bs_env, tau, ispin, .FALSE., .TRUE., mat_phi_mu_l, matrix_G_vir_grid, bs_env%eps_filter)
1943 :
1944 : ! -------------------------------------------------------------------
1945 : ! B. Element-wise Hadamard Products for Sigma_c on Grid
1946 : ! Σ_neg_grid = G_occ_grid ◦ W_grid
1947 : ! Σ_pos_grid = G_vir_grid ◦ W_grid
1948 : ! -------------------------------------------------------------------
1949 0 : CALL dbcsr_create(matrix_Sigma_neg_grid, template=matrix_W_grid)
1950 0 : CALL dbcsr_create(matrix_Sigma_pos_grid, template=matrix_W_grid)
1951 :
1952 : ! Σ^c_ll'(iτ,k=0) = -G^occ_ll'(i|τ|,k=0) * W^MIC_ll'(iτ,k=0), for τ < 0
1953 0 : CALL hadamard_product(matrix_G_occ_grid, matrix_W_grid, matrix_Sigma_neg_grid, 1.0_dp)
1954 :
1955 : ! Σ^c_ll'(iτ,k=0) = G^vir_ll'(i|τ|,k=0) * W^MIC_ll'(iτ,k=0), for τ > 0
1956 0 : CALL hadamard_product(matrix_G_vir_grid, matrix_W_grid, matrix_Sigma_pos_grid, 1.0_dp)
1957 :
1958 : ! Instantly purge massive G_grid arrays to save memory
1959 0 : CALL dbcsr_release(matrix_G_occ_grid)
1960 0 : CALL dbcsr_release(matrix_G_vir_grid)
1961 :
1962 : ! -------------------------------------------------------------------
1963 : ! C. Transform Sigma back to AO Basis
1964 : ! Σ_AO = phi^T * Σ_grid * phi
1965 : ! -------------------------------------------------------------------
1966 :
1967 : ! Σ^c_λσ(iτ,k=0) = sum_ll' Φ_λ(r_l) Σ^c_ll'(iτ,k=0) Φ_σ(r_l'), for τ < 0
1968 : CALL contract_A_B_A("T", "N", mat_phi_mu_l, matrix_Sigma_neg_grid, &
1969 0 : mat_Sigma_neg_tau(i_t, ispin)%matrix, bs_env%eps_filter)
1970 0 : CALL dbcsr_scale(mat_Sigma_neg_tau(i_t, ispin)%matrix, -1.0_dp)
1971 :
1972 : ! Σ^c_λσ(iτ,k=0) = sum_ll' Φ_λ(r_l) Σ^c_ll'(iτ,k=0) Φ_σ(r_l'), for τ > 0
1973 : CALL contract_A_B_A("T", "N", mat_phi_mu_l, matrix_Sigma_pos_grid, &
1974 0 : mat_Sigma_pos_tau(i_t, ispin)%matrix, bs_env%eps_filter)
1975 :
1976 : ! Purge Grid Sigma arrays
1977 0 : CALL dbcsr_release(matrix_Sigma_neg_grid)
1978 0 : CALL dbcsr_release(matrix_Sigma_pos_grid)
1979 :
1980 0 : IF (bs_env%unit_nr > 0) THEN
1981 : WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
1982 0 : 'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1983 0 : ', Execution time', m_walltime() - t1, ' s'
1984 : END IF
1985 :
1986 : END DO ! ispin
1987 :
1988 : ! Release the W_grid for this time point
1989 0 : CALL dbcsr_release(matrix_W_grid)
1990 :
1991 : END DO ! i_t
1992 :
1993 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1994 :
1995 : ! -------------------------------------------------------------------------
1996 : ! 3. FINALIZE AND CLEANUP
1997 : ! -------------------------------------------------------------------------
1998 : CALL fill_fm_Sigma_c_Gamma_time(fm_Sigma_c_Gamma_time, bs_env, &
1999 0 : mat_Sigma_pos_tau, mat_Sigma_neg_tau)
2000 :
2001 0 : CALL cp_fm_release(fm_W_time)
2002 :
2003 0 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_neg_tau)
2004 0 : CALL dbcsr_deallocate_matrix_set(mat_Sigma_pos_tau)
2005 :
2006 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid)
2007 0 : CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
2008 :
2009 0 : CALL delete_unnecessary_files(bs_env)
2010 0 : CALL timestop(handle)
2011 :
2012 0 : END SUBROUTINE compute_Sigma_c
2013 :
2014 : ! **************************************************************************************************
2015 : !> \brief DBCSR Topology Generation
2016 : !> \param matrix_template ...
2017 : !> \param dim_type ...
2018 : !> \param square_dist ...
2019 : !> \param blk_sizes ...
2020 : !> \param mapped_dist ...
2021 : ! **************************************************************************************************
2022 :
2023 30444 : SUBROUTINE setup_square_topology(matrix_template, dim_type, square_dist, blk_sizes, mapped_dist)
2024 :
2025 : TYPE(dbcsr_type), INTENT(IN) :: matrix_template
2026 : CHARACTER(LEN=*), INTENT(IN) :: dim_type
2027 : TYPE(dbcsr_distribution_type), INTENT(OUT) :: square_dist
2028 : INTEGER, DIMENSION(:), INTENT(OUT), POINTER :: blk_sizes, mapped_dist
2029 :
2030 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_square_topology'
2031 :
2032 : INTEGER :: handle, i, np, npcols, nprows
2033 10148 : INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
2034 : TYPE(dbcsr_distribution_type) :: dist_template
2035 :
2036 10148 : CALL timeset(routineN, handle)
2037 :
2038 : CALL dbcsr_get_info(matrix_template, distribution=dist_template, &
2039 10148 : row_blk_size=row_blk, col_blk_size=col_blk)
2040 : CALL dbcsr_distribution_get(dist_template, row_dist=row_dist, col_dist=col_dist, &
2041 10148 : nprows=nprows, npcols=npcols)
2042 :
2043 10148 : IF (TRIM(dim_type) == 'ROW') THEN
2044 : ! Creates ROW x ROW (e.g., Grid x Grid from mat_phi_mu_l)
2045 3964 : blk_sizes => row_blk
2046 3964 : np = npcols
2047 11892 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
2048 47228 : DO i = 1, SIZE(blk_sizes)
2049 47228 : mapped_dist(i) = MOD(i - 1, np)
2050 : END DO
2051 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
2052 3964 : row_dist=row_dist, col_dist=mapped_dist)
2053 :
2054 6184 : ELSE IF (TRIM(dim_type) == 'COL') THEN
2055 : ! Creates COL x COL (e.g., Aux x Aux from mat_Z_lP)
2056 6184 : blk_sizes => col_blk
2057 6184 : np = nprows
2058 18552 : ALLOCATE (mapped_dist(SIZE(blk_sizes)))
2059 18552 : DO i = 1, SIZE(blk_sizes)
2060 18552 : mapped_dist(i) = MOD(i - 1, np)
2061 : END DO
2062 : CALL dbcsr_distribution_new(square_dist, template=dist_template, &
2063 6184 : row_dist=mapped_dist, col_dist=col_dist)
2064 : END IF
2065 :
2066 10148 : CALL timestop(handle)
2067 :
2068 10148 : END SUBROUTINE setup_square_topology
2069 :
2070 : ! **************************************************************************************************
2071 : !> \brief DBCSR matrices deallocation
2072 : !> \param dist ...
2073 : !> \param mapped_dist ...
2074 : !> \param m1 ...
2075 : !> \param m2 ...
2076 : !> \param m3 ...
2077 : !> \param m4 ...
2078 : ! **************************************************************************************************
2079 :
2080 10148 : SUBROUTINE release_dbcsr_topology_and_matrices(dist, mapped_dist, m1, m2, m3, m4)
2081 :
2082 : TYPE(dbcsr_distribution_type), INTENT(INOUT), &
2083 : OPTIONAL :: dist
2084 : INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL, &
2085 : POINTER :: mapped_dist
2086 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: m1, m2, m3, m4
2087 :
2088 : CHARACTER(LEN=*), PARAMETER :: routineN = 'release_dbcsr_topology_and_matrices'
2089 :
2090 : INTEGER :: handle
2091 :
2092 10148 : CALL timeset(routineN, handle)
2093 :
2094 10148 : IF (PRESENT(dist)) CALL dbcsr_distribution_release(dist)
2095 10148 : IF (PRESENT(mapped_dist)) THEN
2096 10148 : IF (ASSOCIATED(mapped_dist)) THEN
2097 10148 : DEALLOCATE (mapped_dist)
2098 : NULLIFY (mapped_dist)
2099 : END IF
2100 : END IF
2101 10148 : IF (PRESENT(m1)) CALL dbcsr_release(m1)
2102 10148 : IF (PRESENT(m2)) CALL dbcsr_release(m2)
2103 10148 : IF (PRESENT(m3)) CALL dbcsr_release(m3)
2104 10148 : IF (PRESENT(m4)) CALL dbcsr_release(m4)
2105 :
2106 10148 : CALL timestop(handle)
2107 :
2108 10148 : END SUBROUTINE release_dbcsr_topology_and_matrices
2109 :
2110 : END MODULE gw_large_cell_Gamma_ri_rs
|
