(git:d3d49ac)
Loading...
Searching...
No Matches
gw_large_cell_gamma_ri_rs.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE 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
98CONTAINS
99
100! **************************************************************************************************
101!> \brief GW calculation using RI-RS formalism for molecules
102!> \param qs_env ...
103!> \param bs_env ...
104! **************************************************************************************************
105
106 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_sigma_x_gamma, fm_w_time
115 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: fm_sigma_c_gamma_time
116
117 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 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 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 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 bs_env%ri_rs%mat_phi_mu_l, bs_env%ri_rs%mat_Z_lP)
159 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 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 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 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 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 CALL compute_qp_energies(bs_env, qs_env, fm_sigma_x_gamma, fm_sigma_c_gamma_time)
206
207 CALL de_init_bs_env(qs_env, bs_env)
208
209 CALL timestop(handle)
210
211 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 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 INTEGER, DIMENSION(:), POINTER :: c_blk_sizes, col_dist, r_blk_sizes, &
234 row_dist
235 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atom_col_buffer
236 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
242 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
243
244 CALL timeset(routinen, handle)
245
246 ! Setup Grid Blocking
247 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 para_env=para_env)
253
254 natom = SIZE(particle_set)
255 total_grid_npts = SIZE(ri_rs_grid_points, 2)
256
257 ! Map the starting indices of spherical gaussian functions (SGF) for each atom
258 ALLOCATE (first_sgf(natom + 1))
259 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 ALLOCATE (c_blk_sizes(natom))
267 DO iatom = 1, natom
268 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 num_grid_chunks = ceiling(real(total_grid_npts, kind=dp)/real(chunk_size, kind=dp))
273 ALLOCATE (r_blk_sizes(num_grid_chunks))
274 r_blk_sizes = chunk_size
275 IF (mod(total_grid_npts, chunk_size) /= 0) THEN
276 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 CALL get_qs_env(qs_env, dbcsr_dist=dbcsr_dist_ks)
281 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 ALLOCATE (row_dist(num_grid_chunks))
286 DO i = 1, num_grid_chunks
287 row_dist(i) = mod(i - 1, nprow)
288 END DO
289
290 ALLOCATE (col_dist(natom))
291 DO i = 1, natom
292 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 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 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 DO iatom = para_env%mepos + 1, natom, para_env%num_pe
308
309 c_size = c_blk_sizes(iatom)
310
311 ! Allocate a temporary dense buffer just for this specific atom
312 ALLOCATE (atom_col_buffer(total_grid_npts, c_size))
313 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 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 DO i_blk = 1, num_grid_chunks
324 r_start = (i_blk - 1)*chunk_size + 1
325 r_end = min(i_blk*chunk_size, total_grid_npts)
326
327 ! Apply dynamic sparsity filtering: Only store blocks with physical significance
328 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 block=atom_col_buffer(r_start:r_end, 1:c_size))
331 END IF
332 END DO
333
334 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 CALL dbcsr_finalize(mat_phi_mu_l)
340
341 IF (bs_env%unit_nr > 0) THEN
342 WRITE (bs_env%unit_nr, *) "Done with evaluation of phi"
343 END IF
344
345 ! -------------------------------------------------------------------------
346 ! CLEANUP
347 ! -------------------------------------------------------------------------
348 DEALLOCATE (first_sgf, r_blk_sizes, c_blk_sizes, row_dist, col_dist)
349 CALL dbcsr_distribution_release(dist)
350
351 CALL timestop(handle)
352
353 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 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 CALL timeset(routinen, handle)
396
397 ! Get Atom Info
398 ikind = particle_set(iatom)%atomic_kind%kind_number
399 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
400 CALL get_cell(cell=cell, h=hmat)
401
402 IF (.NOT. ASSOCIATED(orb_basis_set)) THEN
403 CALL timestop(handle)
404 RETURN
405 END IF
406
407 IF (cell%perd(1) == 1) then; ix_min = -1; ix_max = 1; else; ix_min = 0; ix_max = 0
408 END IF
409 IF (cell%perd(2) == 1) then; iy_min = -1; iy_max = 1; else; iy_min = 0; iy_max = 0
410 END IF
411 IF (cell%perd(3) == 1) then; iz_min = -1; iz_max = 1; else; iz_min = 0; iz_max = 0
412 END IF
413
414 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 !$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 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: local_grid_idx, row_offset
510 INTEGER, DIMENSION(:), POINTER :: col_dist_ri, r_blk_sizes, ri_blk_sizes, &
511 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cutoff_ri_per_atom, cutoff_ri_per_kind, &
515 d_vec_local
516 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_local, d_lp_local, phi_local, &
517 sphere_grid, z_blk
518 REAL(kind=dp), DIMENSION(3) :: dist_vec_raw, pos_p
519 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 TYPE(gw_3c_ctx_type) :: ctx_3c
527 TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
528 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
529 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
530 TYPE(section_vals_type), POINTER :: input
531
532 CALL timeset(routinen, handle)
533
534 t1 = m_walltime()
535
536 CALL get_qs_env(qs_env, para_env=para_env, particle_set=particle_set, input=input, &
537 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 n_procs_per_atom = min(bs_env%ri_rs%n_procs_per_atom_z_lp, para_env%num_pe)
548 IF (n_procs_per_atom < 1) n_procs_per_atom = 1
549
550 NULLIFY (para_env_sub, blacs_env_sub)
551 IF (n_procs_per_atom > 1) THEN
552 n_groups = para_env%num_pe/n_procs_per_atom
553 my_group = min(para_env%mepos/n_procs_per_atom, n_groups - 1)
554 ALLOCATE (para_env_sub)
555 CALL para_env_sub%from_split(para_env, my_group)
556 CALL cp_blacs_env_create(blacs_env=blacs_env_sub, para_env=para_env_sub)
557 atom_p_start = my_group + 1
558 atom_p_stride = n_groups
559 atom_j_mepos = para_env_sub%mepos
560 atom_j_stride = para_env_sub%num_pe
561 ELSE
562 atom_p_start = para_env%mepos + 1
563 atom_p_stride = para_env%num_pe
564 atom_j_mepos = 0
565 atom_j_stride = 1
566 END IF
567
568 natom = SIZE(bs_env%i_RI_start_from_atom)
569 n_ao_total = bs_env%i_ao_end_from_atom(natom)
570 n_grid_total = SIZE(ri_rs_grid_points, 2)
571
572 ! =========================================================================
573 ! 1. SETUP DBCSR TOPOLOGY & EXACT OFFSETS
574 ! =========================================================================
575 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 group=group_handle, npcols=npcol_phi)
578
579 num_grid_chunks = SIZE(r_blk_sizes)
580
581 ALLOCATE (row_offset(num_grid_chunks))
582 row_offset(1) = 0
583 DO i_blk = 2, num_grid_chunks
584 row_offset(i_blk) = row_offset(i_blk - 1) + r_blk_sizes(i_blk - 1)
585 END DO
586
587 ALLOCATE (ri_blk_sizes(natom), col_dist_ri(natom))
588 DO atom_p = 1, natom
589 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 col_dist_ri(atom_p) = mod(atom_p - 1, npcol_phi)
591 END DO
592
593 CALL dbcsr_distribution_new(dist_z, template=dist_phi, row_dist=row_dist_grid, col_dist=col_dist_ri)
594
595 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 matrix_new=mat_z_lp)
599 IF (bs_env%unit_nr > 0) THEN
600 WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
601 'Read Z_lP from file ', ' Execution time', m_walltime() - t1, ' s'
602 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 row_blk_size=r_blk_sizes, col_blk_size=ri_blk_sizes)
609
610 max_ao_size = 0
611 DO j = 1, SIZE(bs_env%i_ao_start_from_atom)
612 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 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 nkind = SIZE(atomic_kind_set)
621 ALLOCATE (cutoff_ri_per_atom(natom))
622
623 IF (bs_env%ri_rs%cutoff_radius_ri_rs > 0.0_dp) THEN
624 cutoff_ri_per_atom(:) = bs_env%ri_rs%cutoff_radius_ri_rs
625 ELSE
626 r_c = bs_env%ri_metric%cutoff_radius
627 DO p_loop_atom = 1, natom
628 cutoff_ri_per_atom(p_loop_atom) = &
629 r_c + bs_env%ri_rs%radius_ao_per_atom(p_loop_atom)
630 END DO
631 END IF
632
633 ALLOCATE (cutoff_ri_per_kind(nkind))
634 cutoff_ri_per_kind(:) = 0.0_dp
635 IF (bs_env%unit_nr > 0) THEN
636 DO p_loop_atom = 1, natom
637 ikind = particle_set(p_loop_atom)%atomic_kind%kind_number
638 cutoff_ri_per_kind(ikind) = max(cutoff_ri_per_kind(ikind), &
639 cutoff_ri_per_atom(p_loop_atom))
640 END DO
641 WRITE (bs_env%unit_nr, '(T2,A)') 'Per-kind maximum RI-RS sphere cutoff (Bohr):'
642 WRITE (bs_env%unit_nr, '(T4,A4,A14)') 'Kind', 'max cutoff_ri'
643 DO ikind = 1, nkind
644 WRITE (bs_env%unit_nr, '(T4,A4,F14.4)') &
645 atomic_kind_set(ikind)%element_symbol, &
646 cutoff_ri_per_kind(ikind)
647 END DO
648 WRITE (bs_env%unit_nr, '(A)') ' '
649 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 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 DO atom_p = atom_p_start, natom, atom_p_stride
666
667 n_loc_ri = ri_blk_sizes(atom_p)
668 pos_p(:) = particle_set(atom_p)%r(:)
669
670 cutoff_ri = cutoff_ri_per_atom(atom_p)
671 cutoff_ri_2 = cutoff_ri**2
672
673 ! ---------------------------------------------------------------------
674 ! A. Determine Local Grid Domain based on cutoff_ri (PBC distance)
675 ! ---------------------------------------------------------------------
676 n_local_grid = 0
677 DO l = 1, n_grid_total
678 dist_vec_raw = pbc(ri_rs_grid_points(1:3, l), pos_p(1:3), cell)
679 dist2_min = dot_product(dist_vec_raw, dist_vec_raw)
680 IF (dist2_min <= cutoff_ri_2) n_local_grid = n_local_grid + 1
681 END DO
682
683 ALLOCATE (local_grid_idx(n_local_grid))
684
685 n_local_grid = 0
686 DO l = 1, n_grid_total
687 dist_vec_raw = pbc(ri_rs_grid_points(1:3, l), pos_p(1:3), cell)
688 dist2_min = dot_product(dist_vec_raw, dist_vec_raw)
689 IF (dist2_min <= cutoff_ri_2) THEN
690 n_local_grid = n_local_grid + 1
691 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 ALLOCATE (sphere_grid(3, n_local_grid))
703 DO loc_idx = 1, n_local_grid
704 sphere_grid(:, loc_idx) = ri_rs_grid_points(:, local_grid_idx(loc_idx))
705 END DO
706
707 ALLOCATE (phi_local(n_local_grid, n_ao_total))
708 phi_local = 0.0_dp
709
710 DO source_atom = 1, natom
711 dist_vec_raw = pbc(particle_set(source_atom)%r(:), pos_p(:), cell)
712 d_sp = norm2(dist_vec_raw)
713 IF (d_sp > bs_env%ri_rs%radius_ao_per_atom(source_atom) + cutoff_ri) cycle
714
715 col_start = bs_env%i_ao_start_from_atom(source_atom)
716 col_end = bs_env%i_ao_end_from_atom(source_atom)
717 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 cell, r2_threshold)
722 END DO
723
724 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 ALLOCATE (d_lp_local(n_local_grid, n_loc_ri))
732 d_lp_local = 0.0_dp
733
734 t2 = m_walltime()
735
736 CALL compute_d_lp(bs_env, ctx_3c, cell, phi_local, d_lp_local, n_local_grid, &
737 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 IF (n_procs_per_atom > 1) THEN
742 CALL para_env_sub%sum(d_lp_local)
743 END IF
744
745 t3 = m_walltime()
746
747 ! ---------------------------------------------------------------------
748 ! D. Build d_vec_local (Jacobi diagonal) + LHS — BLAS or ScaLAPACK
749 ! ---------------------------------------------------------------------
750 ALLOCATE (d_vec_local(n_local_grid))
751
752 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 ALLOCATE (d_local(n_local_grid, n_local_grid))
757 d_local = 0.0_dp
758
759 CALL timeset(routinen//"_dsyrk", handle_dsyrk)
760 CALL dsyrk("L", "N", n_local_grid, n_ao_total, 1.0_dp, phi_local, &
761 n_local_grid, 0.0_dp, d_local, n_local_grid)
762 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 !$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 !$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 !$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 !$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 IF (n_procs_per_atom == 1) THEN
823 CALL timeset(routinen//"_dpotrf", handle_dpotrf)
824 CALL dpotrf('L', n_local_grid, d_local, n_local_grid, info)
825 CALL timestop(handle_dpotrf)
826 CALL timeset(routinen//"_dpotrs", handle_dpotrs)
827 CALL dpotrs('L', n_local_grid, n_loc_ri, d_local, n_local_grid, &
828 d_lp_local, n_local_grid, info)
829 CALL timestop(handle_dpotrs)
830 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 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 !$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 IF (n_procs_per_atom == 1 .OR. para_env_sub%mepos == 0) THEN
863 ALLOCATE (z_blk(maxval(r_blk_sizes), n_loc_ri))
864 loc_ptr = 1
865
866 DO i_blk = 1, num_grid_chunks
867 r_start = row_offset(i_blk) + 1
868 r_end = row_offset(i_blk) + r_blk_sizes(i_blk)
869 current_chunk_size = r_blk_sizes(i_blk)
870
871 z_blk = 0.0_dp
872
873 DO WHILE (loc_ptr <= n_local_grid)
874 g = local_grid_idx(loc_ptr)
875 IF (g > r_end) EXIT
876 z_blk(g - r_start + 1, 1:n_loc_ri) = d_lp_local(loc_ptr, 1:n_loc_ri)
877 loc_ptr = loc_ptr + 1
878 END DO
879
880 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 block=z_blk(1:current_chunk_size, 1:n_loc_ri))
883 END IF
884 END DO
885
886 DEALLOCATE (z_blk)
887 END IF
888
889 DEALLOCATE (d_vec_local, d_lp_local)
890 DEALLOCATE (local_grid_idx, phi_local)
891
892 END DO
893
894 DEALLOCATE (cutoff_ri_per_atom)
895 CALL gw_3c_ctx_release(ctx_3c)
896
897 CALL dbcsr_finalize(mat_z_lp)
898
899 IF (bs_env%unit_nr > 0) THEN
900 WRITE (bs_env%unit_nr, '(T2,A,T57,A,F7.1,A)') &
901 'Computed Z_lP ', ' Execution time', m_walltime() - t1, ' s'
902 WRITE (bs_env%unit_nr, '(A)') ' '
903 END IF
904
905 logger => cp_get_default_logger()
906
907 IF (btest(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
908 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 DEALLOCATE (row_offset, ri_blk_sizes, col_dist_ri)
914 CALL dbcsr_distribution_release(dist_z)
915
916 IF (n_procs_per_atom > 1) THEN
917 CALL cp_blacs_env_release(blacs_env_sub)
918 CALL para_env_sub%free()
919 DEALLOCATE (para_env_sub)
920 END IF
921
922 DEALLOCATE (ri_rs_grid_points)
923
924 CALL timestop(handle)
925
926 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: d_lp_prv, int_2d_prv, rho_chunk
970 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_prv, int_3c_sum
971 TYPE(gw_3c_ws_type) :: ws
972
973 CALL timeset(routinen, handle)
974
975 natom = SIZE(bs_env%i_ao_start_from_atom)
976
977 IF (cell%perd(1) == 1) then; ix_min = -1; ix_max = 1; else; ix_min = 0; ix_max = 0
978 END IF
979 IF (cell%perd(2) == 1) then; iy_min = -1; iy_max = 1; else; iy_min = 0; iy_max = 0
980 END IF
981 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 !$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 CALL timestop(handle)
1094
1095 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 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 INTEGER, DIMENSION(:), POINTER :: blk_ao, blk_grid, dist_col_grid, &
1115 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 CALL timeset(routinen, handle)
1122
1123 ! =========================================================================
1124 ! 1. SETUP CORE TOPOLOGIES
1125 ! =========================================================================
1126 CALL dbcsr_get_info(mat_phi_mu_l, distribution=dist_phi, row_blk_size=blk_grid, col_blk_size=blk_ao)
1127 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 ALLOCATE (dist_col_grid(SIZE(blk_grid)))
1131 DO i = 1, SIZE(blk_grid)
1132 dist_col_grid(i) = mod(i - 1, npcol)
1133 END DO
1134
1135 CALL dbcsr_distribution_new(dist_grid_grid, template=dist_phi, &
1136 row_dist=dist_row_grid, col_dist=dist_col_grid)
1137
1138 CALL dbcsr_create(matrix_g_occ_grid, "G_occ_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1139 CALL dbcsr_create(matrix_g_vir_grid, "G_vir_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1140 CALL dbcsr_create(matrix_chi_grid, "chi_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1141 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 DO i_t = 1, bs_env%num_time_freq_points
1147 t1 = m_walltime()
1148
1149 tau = bs_env%imag_time_points(i_t)
1150 CALL dbcsr_set(matrix_chi_grid, 0.0_dp)
1151
1152 ! ----------------------------------------------------------------------
1153 ! A. SPIN LOOP (Allocations safely encapsulated in wrappers)
1154 ! ----------------------------------------------------------------------
1155 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 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 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 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 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 mat_chi_gamma_tau(i_t)%matrix, bs_env%eps_filter)
1186
1187 IF (bs_env%unit_nr > 0) THEN
1188 WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
1189 χτ'Computed (i,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
1190 ', Execution time', m_walltime() - t1, ' s'
1191 END IF
1192
1193 END DO ! i_t
1194
1195 ! =========================================================================
1196 ! 3. FINAL CLEANUP
1197 ! =========================================================================
1198 CALL dbcsr_release(matrix_g_occ_grid)
1199 CALL dbcsr_release(matrix_g_vir_grid)
1200 CALL dbcsr_release(matrix_chi_grid)
1201 CALL dbcsr_release(matrix_chi_grid_spin)
1202 CALL dbcsr_distribution_release(dist_grid_grid)
1203 DEALLOCATE (dist_col_grid)
1204
1205 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1206
1207 CALL timestop(handle)
1208
1209 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 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 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 CALL timeset(routinen, handle)
1241
1242 ! 1. Select the correct FM matrix based on occ/vir flags
1243 IF (occ) THEN
1244 fm_g => bs_env%fm_Gocc
1245 ELSE
1246 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 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 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 row_blk_size=blk_ao, col_blk_size=blk_ao)
1259
1260 ! 4. Convert FM to Sparse DBCSR
1261 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 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 CALL release_dbcsr_topology_and_matrices(dist=dist_ao_ao, mapped_dist=dist_row_ao, m1=matrix_g_ao)
1269
1270 CALL timestop(handle)
1271
1272 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 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 CALL timeset(routinen, handle)
1296
1297 CALL dbcsr_create(matrix_tmp, template=matrix_a)
1298
1299 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 0.0_dp, matrix_tmp, filter_eps=eps_filter)
1303 CALL dbcsr_multiply("N", "T", 1.0_dp, matrix_tmp, matrix_a, &
1304 0.0_dp, matrix_out, filter_eps=eps_filter)
1305
1306 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 0.0_dp, matrix_tmp, filter_eps=eps_filter)
1310 CALL dbcsr_multiply("T", "N", 1.0_dp, matrix_a, matrix_tmp, &
1311 0.0_dp, matrix_out, filter_eps=eps_filter)
1312 ELSE
1313 cpabort("Unsupported transposition pair in contract_A_B_A")
1314 END IF
1315
1316 CALL dbcsr_release(matrix_tmp)
1317
1318 CALL timestop(handle)
1319
1320 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: blk_b, blk_c
1340 TYPE(dbcsr_iterator_type) :: iter
1341
1342 CALL timeset(routinen, handle)
1343
1344 CALL dbcsr_copy(matrix_c, matrix_a)
1345
1346 CALL dbcsr_iterator_start(iter, matrix_c)
1347 DO WHILE (dbcsr_iterator_blocks_left(iter))
1348 CALL dbcsr_iterator_next_block(iter, row, col, blk_c)
1349
1350 CALL dbcsr_get_block_p(matrix_b, row, col, blk_b, found)
1351
1352 IF (found) THEN
1353 blk_c(:, :) = fac*blk_c(:, :)*blk_b(:, :)
1354 ELSE
1355 ! If B is sparse here, the product is zero
1356 blk_c(:, :) = 0.0_dp
1357 END IF
1358 END DO
1359 CALL dbcsr_iterator_stop(iter)
1360
1361 CALL timestop(handle)
1362
1363 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: blk_a, blk_b
1382 TYPE(dbcsr_iterator_type) :: iter
1383
1384 CALL timeset(routinen, handle)
1385
1386 CALL dbcsr_iterator_start(iter, matrix_a)
1387 DO WHILE (dbcsr_iterator_blocks_left(iter))
1388 CALL dbcsr_iterator_next_block(iter, row, col, blk_a)
1389
1390 CALL dbcsr_get_block_p(matrix_b, row, col, blk_b, found)
1391
1392 IF (found) THEN
1393 blk_a(:, :) = fac*blk_a(:, :)*blk_b(:, :)
1394 ELSE
1395 ! If B is sparse here, the product is zero
1396 blk_a(:, :) = 0.0_dp
1397 END IF
1398 END DO
1399 CALL dbcsr_iterator_stop(iter)
1400
1401 CALL timestop(handle)
1402
1403 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 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 CALL timeset(routinen, handle)
1426
1427 t1 = m_walltime()
1428
1429 CALL create_fm_w_mic_time(bs_env, fm_w_time)
1430
1431 ! 1. Allocate V and M matrices
1432 CALL cp_fm_create(fm_v, bs_env%fm_RI_RI%matrix_struct)
1433 CALL cp_fm_create(fm_v_sqrt, bs_env%fm_RI_RI%matrix_struct)
1434 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 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 DO j_w = 1, bs_env%num_time_freq_points
1441 ! Fourier transformation of χ_PQ(iτ) to χ_PQ(iω_j)
1442 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 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 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 CALL multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_w_time)
1455
1456 IF (bs_env%unit_nr > 0) THEN
1457 WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
1458 τ'Computed W(i),', ' Execution time', m_walltime() - t1, ' s'
1459 END IF
1460
1461 CALL dbcsr_deallocate_matrix_set(mat_chi_gamma_tau)
1462
1463 ! Cleanup
1464 CALL cp_fm_release(fm_v)
1465 CALL cp_fm_release(fm_v_sqrt)
1466 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 IF (bs_env%rtp_method == rtp_method_bse) THEN
1470 t1 = m_walltime()
1471 CALL cp_fm_create(bs_env%fm_W_MIC_freq_zero, bs_env%fm_W_MIC_freq%matrix_struct)
1472 ! Set to zero
1473 CALL cp_fm_set_all(bs_env%fm_W_MIC_freq_zero, 0.0_dp)
1474 ! Sum over all times
1475 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 bs_env%imag_time_weights_freq_zero(i_t), fm_w_time(i_t))
1479 END DO
1480 ! Done, save to file
1481 CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
1482 ! Report calculation
1483 IF (bs_env%unit_nr > 0) THEN
1484 WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
1485 'Computed W(0),', ' Execution time', m_walltime() - t1, ' s'
1486 END IF
1487 END IF
1488
1489 IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
1490
1491 CALL timestop(handle)
1492
1493 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 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 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1513 TYPE(cell_type), POINTER :: cell
1514 TYPE(cp_fm_type) :: fm_work
1515 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_m
1516 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_v_kp
1517 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1518 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1519
1520 CALL timeset(routinen, handle)
1521
1522 n_ri = bs_env%n_RI
1523 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 qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set)
1530
1531 ALLOCATE (mat_v_kp(1:1, 1:2))
1532 NULLIFY (mat_v_kp(1, 1)%matrix, mat_v_kp(1, 2)%matrix)
1533 ALLOCATE (mat_v_kp(1, 1)%matrix, mat_v_kp(1, 2)%matrix)
1534
1535 CALL dbcsr_create(mat_v_kp(1, 1)%matrix, template=bs_env%mat_RI_RI%matrix)
1536 CALL dbcsr_reserve_all_blocks(mat_v_kp(1, 1)%matrix)
1537 CALL dbcsr_set(mat_v_kp(1, 1)%matrix, 0.0_dp)
1538
1539 ! Dummy imaginary part just to satisfy the routine
1540 CALL dbcsr_create(mat_v_kp(1, 2)%matrix, template=bs_env%mat_RI_RI%matrix)
1541 CALL dbcsr_reserve_all_blocks(mat_v_kp(1, 2)%matrix)
1542 CALL dbcsr_set(mat_v_kp(1, 2)%matrix, 0.0_dp)
1543
1544 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 bs_env%size_lattice_sum_V, operator_coulomb, 1, 1)
1549
1550 ! Copy real part to fm_V
1551 CALL copy_dbcsr_to_fm(mat_v_kp(1, 1)%matrix, fm_v)
1552
1553 CALL dbcsr_deallocate_matrix(mat_v_kp(1, 1)%matrix)
1554 CALL dbcsr_deallocate_matrix(mat_v_kp(1, 2)%matrix)
1555 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 do_kpoints=.false., regularization_ri=bs_env%regularization_RI)
1562
1563 ! -----------------------------------------------------------------------
1564 ! 3. M -> M^-1
1565 ! -----------------------------------------------------------------------
1566 CALL cp_fm_cholesky_decompose(fm_m(1, 1), info_out=info)
1567 IF (info == 0) THEN
1568 CALL cp_fm_cholesky_invert(fm_m(1, 1))
1569 CALL cp_fm_uplo_to_full(fm_m(1, 1), fm_work)
1570 ELSE
1571 ! Fallback if Cholesky fails due to conditioning
1572 CALL cp_fm_power(fm_m(1, 1), fm_work, -1.0_dp, bs_env%eps_eigval_mat_RI, ndep)
1573 CALL cp_fm_to_fm(fm_work, fm_m(1, 1))
1574 END IF
1575
1576 ! -----------------------------------------------------------------------
1577 ! 4. V -> V^0.5
1578 ! -----------------------------------------------------------------------
1579 CALL cp_fm_to_fm(fm_v, fm_v_sqrt)
1580 CALL cp_fm_cholesky_decompose(fm_v_sqrt, info_out=info)
1581 IF (info == 0) THEN
1582 CALL clean_lower_part(fm_v_sqrt)
1583 ELSE
1584 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_dp, fm_minv_vsqrt)
1592
1593 CALL cp_fm_release(fm_m)
1594 CALL cp_fm_release(fm_work)
1595
1596 CALL timestop(handle)
1597
1598 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 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 CALL timeset(routinen, handle)
1620
1621 n_ri = bs_env%n_RI
1622
1623 CALL cp_fm_create(fm_eps_freq_j, fm_chi_freq_j%matrix_struct)
1624 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 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 fm_minv_vsqrt, fm_work, 0.0_dp, fm_eps_freq_j)
1636
1637 ! ε(iω_j) = Id - eps_work --> -eps_work + Id
1638 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 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 CALL cp_fm_cholesky_decompose(fm_eps_freq_j, info_out=info)
1649
1650 ! b) Inversion
1651 IF (info == 0) THEN
1652 CALL cp_fm_cholesky_invert(fm_eps_freq_j)
1653 CALL cp_fm_uplo_to_full(fm_eps_freq_j, fm_work)
1654 ELSE
1655 ! Fallback to expensive diagonalization if Cholesky fails
1656 CALL cp_fm_power(fm_eps_freq_j, fm_work, -1.0_dp, bs_env%eps_eigval_mat_RI, ndep)
1657 CALL cp_fm_to_fm(fm_work, fm_eps_freq_j)
1658 END IF
1659
1660 ! c) ε^-1(iω_j) - Id
1661 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_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_dp, fm_w_freq_j)
1670
1671 ! Cleanup
1672 CALL cp_fm_release(fm_work)
1673 CALL cp_fm_release(fm_eps_freq_j)
1674
1675 CALL timestop(handle)
1676
1677 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1694
1695 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 col_indices=col_indices)
1702
1703 DO j_col = 1, ncol_local
1704 j_global = col_indices(j_col)
1705 DO i_row = 1, nrow_local
1706 i_global = row_indices(i_row)
1707 IF (j_global == i_global) THEN
1708 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 CALL timestop(handle)
1714
1715 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1729
1730 CALL timeset(routinen, handle)
1731
1732 CALL cp_fm_get_info(matrix=fm_mat, &
1733 nrow_local=nrow_local, ncol_local=ncol_local, &
1734 row_indices=row_indices, col_indices=col_indices)
1735
1736 DO j_col = 1, ncol_local
1737 j_global = col_indices(j_col)
1738 DO i_row = 1, nrow_local
1739 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 CALL timestop(handle)
1744
1745 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 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 INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
1767 dist_row_aux
1768 REAL(kind=dp) :: t1
1769 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 CALL timeset(routinen, handle)
1776
1777 t1 = m_walltime()
1778
1779 ALLOCATE (fm_sigma_x_gamma(bs_env%n_spin))
1780 DO ispin = 1, bs_env%n_spin
1781 CALL cp_fm_create(fm_sigma_x_gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
1782 END DO
1783
1784 CALL dbcsr_create(mat_sigma_x_gamma, template=bs_env%mat_ao_ao%matrix)
1785
1786 ! =========================================================================
1787 ! 1. SETUP CORE TOPOLOGIES
1788 ! =========================================================================
1789 CALL setup_square_topology(mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, dist_col_grid)
1790 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 bs_env%trunc_coulomb, do_kpoints=.false.)
1797
1798 ! M^-1 * V^tr * M^-1 directly modifies fm_Vtr_Gamma(:, 1)
1799 CALL multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_vtr_gamma(:, 1))
1800
1801 CALL dbcsr_create(matrix_v_aux, "V_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1802 CALL dbcsr_create(matrix_v_grid, "V_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1803
1804 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 CALL contract_a_b_a("N", "T", mat_z_lp, matrix_v_aux, matrix_v_grid, bs_env%eps_filter)
1808 CALL dbcsr_release(matrix_v_aux)
1809
1810 ! =========================================================================
1811 ! 3. SPIN LOOP FOR EXACT EXCHANGE
1812 ! =========================================================================
1813 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 CALL dbcsr_create(matrix_d_grid, "D_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1819 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 CALL dbcsr_create(matrix_sigma_x_grid, template=matrix_v_grid)
1824 CALL hadamard_product(matrix_d_grid, matrix_v_grid, matrix_sigma_x_grid, 1.0_dp)
1825
1826 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 CALL contract_a_b_a("T", "N", mat_phi_mu_l, matrix_sigma_x_grid, mat_sigma_x_gamma, bs_env%eps_filter)
1831 CALL dbcsr_scale(mat_sigma_x_gamma, -1.0_dp)
1832
1833 CALL dbcsr_release(matrix_sigma_x_grid)
1834
1835 ! Data I/O and Export to CP2K Full Matrices
1836 CALL copy_dbcsr_to_fm(mat_sigma_x_gamma, fm_sigma_x_gamma(ispin))
1837
1838 END DO ! ispin
1839
1840 IF (bs_env%unit_nr > 0) THEN
1841 WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
1842 Σ'Computed ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
1843 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 m1=mat_sigma_x_gamma, m2=matrix_v_grid)
1851 CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
1852
1853 CALL cp_fm_release(fm_vtr_gamma)
1854
1855 CALL timestop(handle)
1856
1857 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 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 INTEGER, DIMENSION(:), POINTER :: blk_aux, blk_grid, dist_col_grid, &
1879 dist_row_aux
1880 REAL(kind=dp) :: t1, tau
1881 TYPE(dbcsr_distribution_type) :: dist_aux_aux, dist_grid_grid
1882 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 CALL timeset(routinen, handle)
1887
1888 ! =========================================================================
1889 ! 1. SETUP CORE TOPOLOGIES AND PRE-ALLOCATE OUTPUT ARRAYS
1890 ! =========================================================================
1891 CALL setup_square_topology(mat_phi_mu_l, 'ROW', dist_grid_grid, blk_grid, dist_col_grid)
1892 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 NULLIFY (mat_sigma_neg_tau, mat_sigma_pos_tau)
1896 ALLOCATE (mat_sigma_neg_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1897 ALLOCATE (mat_sigma_pos_tau(bs_env%num_time_freq_points, bs_env%n_spin))
1898
1899 DO i_t = 1, bs_env%num_time_freq_points
1900 DO ispin = 1, bs_env%n_spin
1901 ALLOCATE (mat_sigma_neg_tau(i_t, ispin)%matrix)
1902 ALLOCATE (mat_sigma_pos_tau(i_t, ispin)%matrix)
1903 CALL dbcsr_create(mat_sigma_neg_tau(i_t, ispin)%matrix, template=bs_env%mat_ao_ao%matrix)
1904 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 DO i_t = 1, bs_env%num_time_freq_points
1912 tau = bs_env%imag_time_points(i_t)
1913
1914 ! -------------------------------------------------------------------
1915 ! Compute W_grid = Z * W_aux * Z^T
1916 ! -------------------------------------------------------------------
1917 CALL dbcsr_create(matrix_w_aux, "W_aux", dist_aux_aux, dbcsr_type_no_symmetry, blk_aux, blk_aux)
1918 CALL dbcsr_create(matrix_w_grid, "W_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1919
1920 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 CALL contract_a_b_a("N", "T", mat_z_lp, matrix_w_aux, matrix_w_grid, bs_env%eps_filter)
1924
1925 CALL dbcsr_release(matrix_w_aux) ! Clean up aux basis immediately
1926
1927 DO ispin = 1, bs_env%n_spin
1928 t1 = m_walltime()
1929
1930 ! -------------------------------------------------------------------
1931 ! A. Transform Green's Functions to the Grid
1932 ! -------------------------------------------------------------------
1933 CALL dbcsr_create(matrix_g_occ_grid, "G_occ_grid", dist_grid_grid, dbcsr_type_no_symmetry, blk_grid, blk_grid)
1934 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 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 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 CALL dbcsr_create(matrix_sigma_neg_grid, template=matrix_w_grid)
1950 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 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 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 CALL dbcsr_release(matrix_g_occ_grid)
1960 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 mat_sigma_neg_tau(i_t, ispin)%matrix, bs_env%eps_filter)
1970 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 mat_sigma_pos_tau(i_t, ispin)%matrix, bs_env%eps_filter)
1975
1976 ! Purge Grid Sigma arrays
1977 CALL dbcsr_release(matrix_sigma_neg_grid)
1978 CALL dbcsr_release(matrix_sigma_pos_grid)
1979
1980 IF (bs_env%unit_nr > 0) THEN
1981 WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
1982 Στ'Computed ^c(i,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
1983 ', 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 CALL dbcsr_release(matrix_w_grid)
1990
1991 END DO ! i_t
1992
1993 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 mat_sigma_pos_tau, mat_sigma_neg_tau)
2000
2001 CALL cp_fm_release(fm_w_time)
2002
2003 CALL dbcsr_deallocate_matrix_set(mat_sigma_neg_tau)
2004 CALL dbcsr_deallocate_matrix_set(mat_sigma_pos_tau)
2005
2006 CALL release_dbcsr_topology_and_matrices(dist=dist_grid_grid, mapped_dist=dist_col_grid)
2007 CALL release_dbcsr_topology_and_matrices(dist=dist_aux_aux, mapped_dist=dist_row_aux)
2008
2009 CALL delete_unnecessary_files(bs_env)
2010 CALL timestop(handle)
2011
2012 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 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 INTEGER, DIMENSION(:), POINTER :: col_blk, col_dist, row_blk, row_dist
2034 TYPE(dbcsr_distribution_type) :: dist_template
2035
2036 CALL timeset(routinen, handle)
2037
2038 CALL dbcsr_get_info(matrix_template, distribution=dist_template, &
2039 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 nprows=nprows, npcols=npcols)
2042
2043 IF (trim(dim_type) == 'ROW') THEN
2044 ! Creates ROW x ROW (e.g., Grid x Grid from mat_phi_mu_l)
2045 blk_sizes => row_blk
2046 np = npcols
2047 ALLOCATE (mapped_dist(SIZE(blk_sizes)))
2048 DO i = 1, SIZE(blk_sizes)
2049 mapped_dist(i) = mod(i - 1, np)
2050 END DO
2051 CALL dbcsr_distribution_new(square_dist, template=dist_template, &
2052 row_dist=row_dist, col_dist=mapped_dist)
2053
2054 ELSE IF (trim(dim_type) == 'COL') THEN
2055 ! Creates COL x COL (e.g., Aux x Aux from mat_Z_lP)
2056 blk_sizes => col_blk
2057 np = nprows
2058 ALLOCATE (mapped_dist(SIZE(blk_sizes)))
2059 DO i = 1, SIZE(blk_sizes)
2060 mapped_dist(i) = mod(i - 1, np)
2061 END DO
2062 CALL dbcsr_distribution_new(square_dist, template=dist_template, &
2063 row_dist=mapped_dist, col_dist=col_dist)
2064 END IF
2065
2066 CALL timestop(handle)
2067
2068 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 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 CALL timeset(routinen, handle)
2093
2094 IF (PRESENT(dist)) CALL dbcsr_distribution_release(dist)
2095 IF (PRESENT(mapped_dist)) THEN
2096 IF (ASSOCIATED(mapped_dist)) THEN
2097 DEALLOCATE (mapped_dist)
2098 NULLIFY (mapped_dist)
2099 END IF
2100 END IF
2101 IF (PRESENT(m1)) CALL dbcsr_release(m1)
2102 IF (PRESENT(m2)) CALL dbcsr_release(m2)
2103 IF (PRESENT(m3)) CALL dbcsr_release(m3)
2104 IF (PRESENT(m4)) CALL dbcsr_release(m4)
2105
2106 CALL timestop(handle)
2107
2108 END SUBROUTINE release_dbcsr_topology_and_matrices
2109
2110END MODULE gw_large_cell_gamma_ri_rs
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:56
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cell_types::pbc
Definition cell_types.F:104
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
basis_set_types
Definition basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
constants_operator
constants for the different operators of the 2c-integrals
Definition constants_operator.F:14
constants_operator::operator_coulomb
integer, parameter, public operator_coulomb
Definition constants_operator.F:22
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_blacs_env::cp_blacs_env_release
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
Definition cp_blacs_env.F:274
cp_blacs_env::cp_blacs_env_create
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Definition cp_blacs_env.F:123
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:643
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:616
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_binary_write
subroutine, public dbcsr_binary_write(matrix, filepath)
...
Definition cp_dbcsr_api.F:290
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_binary_read
subroutine, public dbcsr_binary_read(filepath, distribution, matrix_new)
...
Definition cp_dbcsr_api.F:273
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_reserve_all_blocks
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
Definition cp_dbcsr_contrib.F:295
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::max_elements_per_block
integer, save, public max_elements_per_block
Definition cp_dbcsr_operations.F:57
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:212
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1584
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_invert
subroutine, public cp_fm_cholesky_invert(matrix, n, info_out)
used to replace the cholesky decomposition by the inverse
Definition cp_fm_cholesky.F:182
cp_fm_cholesky::cp_fm_cholesky_decompose
subroutine, public cp_fm_cholesky_decompose(matrix, n, info_out)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_fm_cholesky.F:61
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_power
subroutine, public cp_fm_power(matrix, work, exponent, threshold, n_dependent, verbose, eigvals)
...
Definition cp_fm_diag.F:1127
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
gw_integrals
Utility method to build 3-center integrals for small cell GW.
Definition gw_integrals.F:11
gw_integrals::build_3c_integral_block_ctx
subroutine, public build_3c_integral_block_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, cell_j, cell_k, cell_i, j_offset, k_offset, i_offset, screened)
Computes the 3c integral block (mu(atom_j) nu(atom_k) | P(atom_i)) for ONE atom triple,...
Definition gw_integrals.F:305
gw_integrals::gw_3c_ws_create
subroutine, public gw_3c_ws_create(ws, ctx)
Creates a per-thread 3c workspace: libint object + LIBXSMM contraction buffers.
Definition gw_integrals.F:259
gw_integrals::gw_3c_ws_release
subroutine, public gw_3c_ws_release(ws)
Releases a per-thread 3c workspace.
Definition gw_integrals.F:275
gw_integrals::gw_3c_ctx_release
subroutine, public gw_3c_ctx_release(ctx)
Releases the shared 3c-integral context.
Definition gw_integrals.F:232
gw_integrals::gw_3c_ctx_create
subroutine, public gw_3c_ctx_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i)
Builds the shared 3c-integral context: screening radii, basis maxima, contracted sphi tables,...
Definition gw_integrals.F:107
gw_large_cell_gamma_ri_rs
GW using RI-RS Approximation for molecules.
Definition gw_large_cell_gamma_ri_rs.F:13
gw_large_cell_gamma_ri_rs::release_dbcsr_topology_and_matrices
subroutine, public release_dbcsr_topology_and_matrices(dist, mapped_dist, m1, m2, m3, m4)
DBCSR matrices deallocation.
Definition gw_large_cell_gamma_ri_rs.F:2081
gw_large_cell_gamma_ri_rs::hadamard_product_inplace
subroutine, public hadamard_product_inplace(matrix_a, matrix_b, fac)
In-place Hadamard A <- fac * (A ◦ B). Value mutation only (no block insert/delete),...
Definition gw_large_cell_gamma_ri_rs.F:1373
gw_large_cell_gamma_ri_rs::contract_a_b_a
subroutine, public contract_a_b_a(transa_left, transa_right, matrix_a, matrix_b, matrix_out, eps_filter)
Generalized routine to compute OUT = A * B * A^T OR OUT = A^T * B * A using DBCSR.
Definition gw_large_cell_gamma_ri_rs.F:1285
gw_large_cell_gamma_ri_rs::gw_calc_large_cell_gamma_ri_rs
subroutine, public gw_calc_large_cell_gamma_ri_rs(qs_env, bs_env)
GW calculation using RI-RS formalism for molecules.
Definition gw_large_cell_gamma_ri_rs.F:107
gw_large_cell_gamma_ri_rs::setup_square_topology
subroutine, public setup_square_topology(matrix_template, dim_type, square_dist, blk_sizes, mapped_dist)
DBCSR Topology Generation.
Definition gw_large_cell_gamma_ri_rs.F:2024
gw_large_cell_gamma
Routines from paper [Graml2024].
Definition gw_large_cell_gamma.F:15
gw_large_cell_gamma::compute_fm_chi_gamma_freq
subroutine, public compute_fm_chi_gamma_freq(bs_env, fm_chi_gamma_freq, j_w, mat_chi_gamma_tau)
...
Definition gw_large_cell_gamma.F:1342
gw_large_cell_gamma::delete_unnecessary_files
subroutine, public delete_unnecessary_files(bs_env)
...
Definition gw_large_cell_gamma.F:2381
gw_large_cell_gamma::fill_fm_sigma_c_gamma_time
subroutine, public fill_fm_sigma_c_gamma_time(fm_sigma_c_gamma_time, bs_env, mat_sigma_pos_tau, mat_sigma_neg_tau)
...
Definition gw_large_cell_gamma.F:2272
gw_large_cell_gamma::fm_write
subroutine, public fm_write(fm, matrix_index, matrix_name, qs_env)
...
Definition gw_large_cell_gamma.F:500
gw_large_cell_gamma::compute_qp_energies
subroutine, public compute_qp_energies(bs_env, qs_env, fm_sigma_x_gamma, fm_sigma_c_gamma_time)
...
Definition gw_large_cell_gamma.F:2448
gw_large_cell_gamma::multiply_fm_w_mic_time_with_minv_gamma
subroutine, public multiply_fm_w_mic_time_with_minv_gamma(bs_env, qs_env, fm_w_mic_time)
...
Definition gw_large_cell_gamma.F:1645
gw_large_cell_gamma::get_w_mic
subroutine, public get_w_mic(bs_env, qs_env, mat_chi_gamma_tau, fm_w_mic_time)
...
Definition gw_large_cell_gamma.F:796
gw_large_cell_gamma::g_occ_vir
subroutine, public g_occ_vir(bs_env, tau, fm_g_gamma, ispin, occ, vir)
...
Definition gw_large_cell_gamma.F:553
gw_large_cell_gamma::create_fm_w_mic_time
subroutine, public create_fm_w_mic_time(bs_env, fm_w_mic_time)
...
Definition gw_large_cell_gamma.F:1583
gw_large_cell_gamma::fourier_transform_w_to_t
subroutine, public fourier_transform_w_to_t(bs_env, fm_w_mic_time, fm_w_mic_freq_j, j_w)
...
Definition gw_large_cell_gamma.F:1609
gw_non_periodic_ri_rs
GW using RI-RS Approximation for molecules.
Definition gw_non_periodic_ri_rs.F:14
gw_non_periodic_ri_rs::ri_rs_grid_assembler
subroutine, public ri_rs_grid_assembler(qs_env, bs_env, ri_rs_grid_points)
Compute grid points for RI-RS Right now based on Ivan and Xavier implementation JCP 150,...
Definition gw_non_periodic_ri_rs.F:406
gw_non_periodic_ri_rs::precompute_ri_rs_radii
subroutine, public precompute_ri_rs_radii(qs_env, bs_env)
Compute per-atom AO and RI basis radii from the most diffuse Gaussian primitive in the AO ("ORB") and...
Definition gw_non_periodic_ri_rs.F:242
gw_non_periodic_ri_rs::get_basis_offsets
subroutine, public get_basis_offsets(particle_set, qs_kind_set, first_sgf, total_sgf)
Computes the AO basis offsets: first_sgf(iatom) is the global index of the first spherical Gaussian f...
Definition gw_non_periodic_ri_rs.F:910
gw_non_periodic_ri_rs::solve_d_lp_distributed
subroutine, public solve_d_lp_distributed(phi_local, d_vec, d_lp, n_loc, n_ao, n_rhs, tikhonov, para_env_sub, blacs_env_sub, fm_struct_d, fm_struct_b, fm_d, fm_b, info)
Distributed pdpotrf/pdpotrs solve of D x = b for one atom of the RI-RS Z_lP build (Phase B,...
Definition gw_non_periodic_ri_rs.F:2088
gw_utils
Definition gw_utils.F:15
gw_utils::de_init_bs_env
subroutine, public de_init_bs_env(qs_env, bs_env)
Releases the memory-heavy GW intermediates that cannot be freed in bs_env_release,...
Definition gw_utils.F:240
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::rtp_method_bse
integer, parameter, public rtp_method_bse
Definition input_constants.F:986
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kpoint_coulomb_2c
Routines to compute the Coulomb integral V_(alpha beta)(k) for a k-point k using lattice summation in...
Definition kpoint_coulomb_2c.F:15
kpoint_coulomb_2c::build_2c_coulomb_matrix_kp
subroutine, public build_2c_coulomb_matrix_kp(matrix_v_kp, kpoints, basis_type, cell, particle_set, qs_kind_set, atomic_kind_set, size_lattice_sum, operator_type, ikp_start, ikp_end)
...
Definition kpoint_coulomb_2c.F:76
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mp2_ri_2c
Framework for 2c-integrals for RI.
Definition mp2_ri_2c.F:14
mp2_ri_2c::ri_2c_integral_mat
subroutine, public ri_2c_integral_mat(qs_env, fm_matrix_minv_l_kpoints, fm_matrix_l, dimen_ri, ri_metric, do_kpoints, kpoints, put_mat_ks_env, regularization_ri, ikp_ext, do_build_cell_index)
...
Definition mp2_ri_2c.F:564
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
post_scf_bandstructure_types
Definition post_scf_bandstructure_types.F:13
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
gw_integrals::gw_3c_ctx_type
Shared read-only context for repeated 3-center integral block builds: screening parameters,...
Definition gw_integrals.F:65
gw_integrals::gw_3c_ws_type
Per-thread workspace for 3-center integral block builds: libint object + contraction buffers....
Definition gw_integrals.F:89
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_types::particle_type
Definition particle_types.F:35
post_scf_bandstructure_types::post_scf_bandstructure_type
Definition post_scf_bandstructure_types.F:125
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190