Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Framework for 2c-integrals for RI
10 : !> \par History
11 : !> 06.2012 created [Mauro Del Ben]
12 : !> 03.2019 separated from mp2_ri_gpw [Frederick Stein]
13 : ! **************************************************************************************************
14 : MODULE mp2_ri_2c
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind_set
17 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
18 : gto_basis_set_type
19 : USE cell_types, ONLY: cell_type,&
20 : get_cell,&
21 : plane_distance
22 : USE constants_operator, ONLY: operator_coulomb
23 : USE cp_blacs_env, ONLY: cp_blacs_env_create,&
24 : cp_blacs_env_release,&
25 : cp_blacs_env_type
26 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_scale_and_add_fm
27 : USE cp_cfm_types, ONLY: cp_cfm_create,&
28 : cp_cfm_release,&
29 : cp_cfm_to_fm,&
30 : cp_cfm_type
31 : USE cp_control_types, ONLY: dft_control_type
32 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34 : cp_dbcsr_dist2d_to_dist,&
35 : dbcsr_allocate_matrix_set,&
36 : dbcsr_deallocate_matrix_set
37 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_param
38 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
39 : cp_fm_triangular_invert
40 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
41 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
42 : cp_fm_power,&
43 : cp_fm_syevx
44 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
45 : cp_fm_struct_release,&
46 : cp_fm_struct_type
47 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
48 : cp_fm_create,&
49 : cp_fm_get_info,&
50 : cp_fm_indxg2p,&
51 : cp_fm_release,&
52 : cp_fm_set_all,&
53 : cp_fm_to_fm,&
54 : cp_fm_type
55 : USE dbcsr_api, ONLY: &
56 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
57 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, dbcsr_release, &
58 : dbcsr_reserve_all_blocks, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
59 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
60 : USE distribution_2d_types, ONLY: distribution_2d_type
61 : USE group_dist_types, ONLY: create_group_dist,&
62 : get_group_dist,&
63 : group_dist_d1_type,&
64 : release_group_dist
65 : USE input_constants, ONLY: do_eri_gpw,&
66 : do_eri_mme,&
67 : do_eri_os,&
68 : do_potential_coulomb,&
69 : do_potential_id,&
70 : do_potential_long,&
71 : do_potential_truncated
72 : USE kinds, ONLY: dp
73 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
74 : USE kpoint_methods, ONLY: rskp_transform
75 : USE kpoint_types, ONLY: get_kpoint_info,&
76 : kpoint_type
77 : USE libint_2c_3c, ONLY: compare_potential_types,&
78 : libint_potential_type
79 : USE machine, ONLY: m_flush
80 : USE message_passing, ONLY: mp_comm_type,&
81 : mp_para_env_release,&
82 : mp_para_env_type,&
83 : mp_request_type
84 : USE mp2_eri, ONLY: mp2_eri_2c_integrate
85 : USE mp2_eri_gpw, ONLY: mp2_eri_2c_integrate_gpw
86 : USE mp2_types, ONLY: integ_mat_buffer_type
87 : USE parallel_gemm_api, ONLY: parallel_gemm
88 : USE particle_methods, ONLY: get_particle_set
89 : USE particle_types, ONLY: particle_type
90 : USE qs_environment_types, ONLY: get_qs_env,&
91 : qs_environment_type
92 : USE qs_integral_utils, ONLY: basis_set_list_setup
93 : USE qs_interactions, ONLY: init_interaction_radii
94 : USE qs_kind_types, ONLY: get_qs_kind,&
95 : qs_kind_type
96 : USE qs_ks_types, ONLY: get_ks_env,&
97 : set_ks_env
98 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
99 : release_neighbor_list_sets
100 : USE qs_tensors, ONLY: build_2c_integrals,&
101 : build_2c_neighbor_lists
102 : USE rpa_communication, ONLY: communicate_buffer
103 : USE rpa_gw_kpoints_util, ONLY: cp_cfm_power
104 :
105 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
106 : #include "./base/base_uses.f90"
107 :
108 : IMPLICIT NONE
109 :
110 : PRIVATE
111 :
112 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_2c'
113 :
114 : PUBLIC :: get_2c_integrals, setup_trunc_coulomb_pot_for_exchange_self_energy, RI_2c_integral_mat, &
115 : inversion_of_M_and_mult_with_chol_dec_of_V
116 :
117 : CONTAINS
118 :
119 : ! **************************************************************************************************
120 : !> \brief ...
121 : !> \param qs_env ...
122 : !> \param eri_method ...
123 : !> \param eri_param ...
124 : !> \param para_env ...
125 : !> \param para_env_sub ...
126 : !> \param mp2_memory ...
127 : !> \param my_Lrows ...
128 : !> \param my_Vrows ...
129 : !> \param fm_matrix_PQ ...
130 : !> \param ngroup ...
131 : !> \param color_sub ...
132 : !> \param dimen_RI ...
133 : !> \param dimen_RI_red reduced RI dimension, needed if we perform SVD instead of Cholesky
134 : !> \param kpoints ...
135 : !> \param my_group_L_size ...
136 : !> \param my_group_L_start ...
137 : !> \param my_group_L_end ...
138 : !> \param gd_array ...
139 : !> \param calc_PQ_cond_num ...
140 : !> \param do_svd ...
141 : !> \param potential ...
142 : !> \param ri_metric ...
143 : !> \param fm_matrix_L_kpoints ...
144 : !> \param fm_matrix_Minv_L_kpoints ...
145 : !> \param fm_matrix_Minv ...
146 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
147 : !> \param do_im_time ...
148 : !> \param do_kpoints ...
149 : !> \param mp2_eps_pgf_orb_S ...
150 : !> \param qs_kind_set ...
151 : !> \param sab_orb_sub ...
152 : !> \param calc_forces ...
153 : !> \param unit_nr ...
154 : ! **************************************************************************************************
155 592 : SUBROUTINE get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, &
156 : my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, &
157 : kpoints, my_group_L_size, my_group_L_start, my_group_L_end, &
158 : gd_array, calc_PQ_cond_num, do_svd, potential, ri_metric, &
159 : fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
160 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
161 : do_im_time, do_kpoints, mp2_eps_pgf_orb_S, qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
162 :
163 : TYPE(qs_environment_type), POINTER :: qs_env
164 : INTEGER, INTENT(IN) :: eri_method
165 : TYPE(cp_eri_mme_param), POINTER :: eri_param
166 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
167 : REAL(KIND=dp), INTENT(IN) :: mp2_memory
168 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
169 : INTENT(OUT) :: my_Lrows, my_Vrows
170 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_PQ
171 : INTEGER, INTENT(IN) :: ngroup, color_sub
172 : INTEGER, INTENT(OUT) :: dimen_RI, dimen_RI_red
173 : TYPE(kpoint_type), POINTER :: kpoints
174 : INTEGER, INTENT(OUT) :: my_group_L_size, my_group_L_start, &
175 : my_group_L_end
176 : TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
177 : LOGICAL, INTENT(IN) :: calc_PQ_cond_num, do_svd
178 : TYPE(libint_potential_type) :: potential, ri_metric
179 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_L_kpoints, &
180 : fm_matrix_Minv_L_kpoints, &
181 : fm_matrix_Minv, &
182 : fm_matrix_Minv_Vtrunc_Minv
183 : LOGICAL, INTENT(IN) :: do_im_time, do_kpoints
184 : REAL(KIND=dp), INTENT(IN) :: mp2_eps_pgf_orb_S
185 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
186 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
187 : POINTER :: sab_orb_sub
188 : LOGICAL, INTENT(IN) :: calc_forces
189 : INTEGER, INTENT(IN) :: unit_nr
190 :
191 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_2c_integrals'
192 :
193 : INTEGER :: handle, num_small_eigen
194 : REAL(KIND=dp) :: cond_num, eps_pgf_orb_old
195 : TYPE(cp_fm_type) :: fm_matrix_L_work, fm_matrix_M_inv_work, &
196 : fm_matrix_V
197 : TYPE(dft_control_type), POINTER :: dft_control
198 : TYPE(libint_potential_type) :: trunc_coulomb
199 : TYPE(mp_para_env_type), POINTER :: para_env_L
200 :
201 592 : CALL timeset(routineN, handle)
202 :
203 : ! calculate V and store it in fm_matrix_L_work
204 : CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
205 : fm_matrix_L_work, ngroup, color_sub, dimen_RI, &
206 : my_group_L_size, my_group_L_start, my_group_L_end, &
207 : gd_array, calc_PQ_cond_num, cond_num, &
208 592 : num_small_eigen, potential, sab_orb_sub, do_im_time=do_im_time)
209 :
210 592 : IF (do_im_time .AND. calc_forces) THEN
211 : !need a copy of the (P|Q) integrals
212 50 : CALL cp_fm_create(fm_matrix_PQ, fm_matrix_L_work%matrix_struct)
213 50 : CALL cp_fm_to_fm(fm_matrix_L_work, fm_matrix_PQ)
214 : END IF
215 :
216 592 : dimen_RI_red = dimen_RI
217 :
218 592 : IF (compare_potential_types(ri_metric, potential) .AND. .NOT. do_im_time) THEN
219 : CALL decomp_mat_L(fm_matrix_L_work, do_svd, num_small_eigen, cond_num, .TRUE., gd_array, ngroup, &
220 422 : dimen_RI, dimen_RI_red, para_env)
221 : ELSE
222 :
223 : ! RI-metric matrix (depending on RI metric: Coulomb, overlap or truncated Coulomb)
224 170 : IF (do_im_time) THEN
225 134 : CALL get_qs_env(qs_env, dft_control=dft_control)
226 :
227 : ! re-init the radii to be able to generate pair lists with appropriate screening for overlap matrix
228 134 : eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
229 134 : dft_control%qs_control%eps_pgf_orb = mp2_eps_pgf_orb_S
230 134 : CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
231 :
232 : CALL RI_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L_work, dimen_RI, ri_metric, &
233 : do_kpoints, kpoints, put_mat_KS_env=.TRUE., &
234 134 : regularization_RI=qs_env%mp2_env%ri_rpa_im_time%regularization_RI)
235 :
236 : ! re-init the radii to previous values
237 134 : dft_control%qs_control%eps_pgf_orb = eps_pgf_orb_old
238 134 : CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
239 :
240 : ! GPW overlap matrix
241 : ELSE
242 :
243 36 : CALL mp_para_env_release(para_env_L)
244 36 : CALL release_group_dist(gd_array)
245 :
246 108 : ALLOCATE (fm_matrix_Minv_L_kpoints(1, 1))
247 :
248 : ! Calculate matrix of RI operator (for overlap metric: S), store it in fm_matrix_Minv_L_kpoints
249 : CALL compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
250 : fm_matrix_Minv_L_kpoints(1, 1), ngroup, color_sub, dimen_RI, &
251 : my_group_L_size, my_group_L_start, my_group_L_end, &
252 : gd_array, calc_PQ_cond_num, cond_num, &
253 : num_small_eigen, ri_metric, sab_orb_sub, &
254 36 : fm_matrix_L_extern=fm_matrix_L_work)
255 :
256 : END IF
257 :
258 170 : IF (do_kpoints) THEN
259 :
260 : ! k-dependent 1/r Coulomb matrix V_PQ(k)
261 22 : CALL compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, kpoints)
262 :
263 : CALL inversion_of_M_and_mult_with_chol_dec_of_V(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, dimen_RI, &
264 22 : kpoints, qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S)
265 :
266 22 : CALL setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, trunc_coulomb)
267 :
268 : ! Gamma-only truncated Coulomb matrix V^tr with cutoff radius = half the unit cell size; for exchange self-energy
269 : CALL RI_2c_integral_mat(qs_env, fm_matrix_Minv_Vtrunc_Minv, fm_matrix_L_work, dimen_RI, trunc_coulomb, &
270 22 : do_kpoints=.FALSE., kpoints=kpoints, put_mat_KS_env=.FALSE., regularization_RI=0.0_dp)
271 :
272 : ! Gamma-only RI-metric matrix; for computing Gamma-only/MIC self-energy
273 : CALL RI_2c_integral_mat(qs_env, fm_matrix_Minv, fm_matrix_L_work, dimen_RI, ri_metric, &
274 22 : do_kpoints=.FALSE., kpoints=kpoints, put_mat_KS_env=.FALSE., regularization_RI=0.0_dp)
275 :
276 : CALL Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
277 22 : fm_matrix_Minv, qs_env)
278 : ELSE
279 148 : IF (calc_forces .AND. (.NOT. do_im_time)) THEN
280 : ! For the gradients, we make a copy of the inverse root of L
281 12 : CALL cp_fm_create(fm_matrix_V, fm_matrix_L_work%matrix_struct)
282 12 : CALL cp_fm_to_fm(fm_matrix_L_work, fm_matrix_V)
283 :
284 : CALL decomp_mat_L(fm_matrix_V, do_svd, num_small_eigen, cond_num, .TRUE., gd_array, ngroup, &
285 12 : dimen_RI, dimen_RI_red, para_env)
286 : END IF
287 :
288 : CALL decomp_mat_L(fm_matrix_L_work, do_svd, num_small_eigen, cond_num, .FALSE., gd_array, ngroup, &
289 148 : dimen_RI, dimen_RI_red, para_env)
290 :
291 : CALL decomp_mat_L(fm_matrix_Minv_L_kpoints(1, 1), .FALSE., num_small_eigen, cond_num, .TRUE., &
292 148 : gd_array, ngroup, dimen_RI, dimen_RI_red, para_env)
293 :
294 148 : CALL cp_fm_create(fm_matrix_M_inv_work, fm_matrix_Minv_L_kpoints(1, 1)%matrix_struct)
295 148 : CALL cp_fm_set_all(fm_matrix_M_inv_work, 0.0_dp)
296 :
297 : CALL parallel_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_Minv_L_kpoints(1, 1), &
298 148 : fm_matrix_Minv_L_kpoints(1, 1), 0.0_dp, fm_matrix_M_inv_work)
299 :
300 148 : IF (do_svd) THEN
301 : ! We have to reset the size of fm_matrix_Minv_L_kpoints
302 32 : CALL reset_size_matrix(fm_matrix_Minv_L_kpoints(1, 1), dimen_RI_red, fm_matrix_L_work%matrix_struct)
303 :
304 : ! L (=fm_matrix_Minv_L_kpoints) = S^(-1)*chol_dec(V)
305 : CALL parallel_gemm('T', 'N', dimen_RI, dimen_RI_red, dimen_RI, 1.0_dp, fm_matrix_M_inv_work, &
306 32 : fm_matrix_L_work, 0.0_dp, fm_matrix_Minv_L_kpoints(1, 1))
307 : ELSE
308 : ! L (=fm_matrix_Minv_L_kpoints) = S^(-1)*chol_dec(V)
309 : CALL parallel_gemm('T', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_M_inv_work, &
310 116 : fm_matrix_L_work, 0.0_dp, fm_matrix_Minv_L_kpoints(1, 1))
311 : END IF
312 :
313 : ! clean the S_inv matrix
314 148 : CALL cp_fm_release(fm_matrix_M_inv_work)
315 : END IF
316 :
317 170 : IF (.NOT. do_im_time) THEN
318 :
319 36 : CALL cp_fm_to_fm(fm_matrix_Minv_L_kpoints(1, 1), fm_matrix_L_work)
320 36 : CALL cp_fm_release(fm_matrix_Minv_L_kpoints)
321 :
322 : END IF
323 :
324 : END IF
325 :
326 : ! If the group distribution changed because of an SVD, we get the new values here
327 592 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
328 :
329 1050 : IF (.NOT. do_im_time) THEN
330 458 : IF (unit_nr > 0) THEN
331 229 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Cholesky decomposition group size:", para_env_L%num_pe
332 229 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Number of groups for auxiliary basis functions", ngroup
333 229 : IF (calc_PQ_cond_num .OR. do_svd) THEN
334 : WRITE (UNIT=unit_nr, FMT="(T3,A,T67,ES14.5)") &
335 16 : "RI_INFO| Condition number of the (P|Q):", cond_num
336 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
337 16 : "RI_INFO| Number of non-positive Eigenvalues of (P|Q):", num_small_eigen
338 : END IF
339 229 : CALL m_flush(unit_nr)
340 : END IF
341 :
342 : ! replicate the necessary row of the L^{-1} matrix on each proc
343 458 : CALL grep_Lcols(fm_matrix_L_work, my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
344 470 : IF (calc_forces .AND. .NOT. compare_potential_types(qs_env%mp2_env%ri_metric, &
345 : qs_env%mp2_env%potential_parameter)) THEN
346 12 : CALL grep_Lcols(fm_matrix_V, my_group_L_start, my_group_L_end, my_group_L_size, my_Vrows)
347 : END IF
348 : END IF
349 592 : CALL cp_fm_release(fm_matrix_L_work)
350 592 : IF (calc_forces .AND. .NOT. (do_im_time .OR. compare_potential_types(qs_env%mp2_env%ri_metric, &
351 12 : qs_env%mp2_env%potential_parameter))) CALL cp_fm_release(fm_matrix_V)
352 592 : CALL mp_para_env_release(para_env_L)
353 :
354 592 : CALL timestop(handle)
355 :
356 592 : END SUBROUTINE get_2c_integrals
357 :
358 : ! **************************************************************************************************
359 : !> \brief ...
360 : !> \param fm_matrix_L ...
361 : !> \param do_svd ...
362 : !> \param num_small_eigen ...
363 : !> \param cond_num ...
364 : !> \param do_inversion ...
365 : !> \param gd_array ...
366 : !> \param ngroup ...
367 : !> \param dimen_RI ...
368 : !> \param dimen_RI_red ...
369 : !> \param para_env ...
370 : ! **************************************************************************************************
371 730 : SUBROUTINE decomp_mat_L(fm_matrix_L, do_svd, num_small_eigen, cond_num, do_inversion, gd_array, ngroup, &
372 : dimen_RI, dimen_RI_red, para_env)
373 :
374 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_matrix_L
375 : LOGICAL, INTENT(IN) :: do_svd
376 : INTEGER, INTENT(INOUT) :: num_small_eigen
377 : REAL(KIND=dp), INTENT(INOUT) :: cond_num
378 : LOGICAL, INTENT(IN) :: do_inversion
379 : TYPE(group_dist_d1_type), INTENT(INOUT) :: gd_array
380 : INTEGER, INTENT(IN) :: ngroup, dimen_RI
381 : INTEGER, INTENT(INOUT) :: dimen_RI_red
382 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
383 :
384 730 : IF (do_svd) THEN
385 44 : CALL matrix_root_with_svd(fm_matrix_L, num_small_eigen, cond_num, do_inversion, para_env)
386 :
387 44 : dimen_RI_red = dimen_RI - num_small_eigen
388 :
389 : ! We changed the size of fm_matrix_L in matrix_root_with_svd.
390 : ! So, we have to get new group sizes
391 44 : CALL release_group_dist(gd_array)
392 :
393 44 : CALL create_group_dist(gd_array, ngroup, dimen_RI_red)
394 : ELSE
395 :
396 : ! calculate Cholesky decomposition L (V=LL^T)
397 686 : CALL cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion=do_inversion)
398 :
399 686 : IF (do_inversion) CALL invert_mat(fm_matrix_L)
400 : END IF
401 :
402 730 : END SUBROUTINE decomp_mat_L
403 :
404 : ! **************************************************************************************************
405 : !> \brief ...
406 : !> \param qs_env ...
407 : !> \param fm_matrix_L_kpoints ...
408 : !> \param fm_matrix_Minv_L_kpoints ...
409 : !> \param kpoints ...
410 : ! **************************************************************************************************
411 22 : SUBROUTINE compute_V_by_lattice_sum(qs_env, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, kpoints)
412 : TYPE(qs_environment_type), POINTER :: qs_env
413 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_L_kpoints, &
414 : fm_matrix_Minv_L_kpoints
415 : TYPE(kpoint_type), POINTER :: kpoints
416 :
417 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_V_by_lattice_sum'
418 :
419 : INTEGER :: handle, i_dim, i_real_imag, ikp, nkp, &
420 : nkp_extra, nkp_orig
421 : INTEGER, DIMENSION(3) :: nkp_grid_orig, periodic
422 22 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
423 : TYPE(cell_type), POINTER :: cell
424 22 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_RI_aux_transl, matrix_v_RI_kp
425 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
426 22 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
427 :
428 22 : CALL timeset(routineN, handle)
429 :
430 22 : NULLIFY (matrix_s_RI_aux_transl, particle_set, cell, qs_kind_set)
431 :
432 : CALL get_qs_env(qs_env=qs_env, &
433 : matrix_s_RI_aux_kp=matrix_s_RI_aux_transl, &
434 : particle_set=particle_set, &
435 : cell=cell, &
436 : qs_kind_set=qs_kind_set, &
437 22 : atomic_kind_set=atomic_kind_set)
438 :
439 : ! check that we have a 2n x 2m x 2k mesh to guarantee good convergence
440 22 : CALL get_cell(cell=cell, periodic=periodic)
441 88 : DO i_dim = 1, 3
442 88 : IF (periodic(i_dim) == 1) THEN
443 40 : CPASSERT(MODULO(kpoints%nkp_grid(i_dim), 2) == 0)
444 : END IF
445 : END DO
446 :
447 22 : nkp = kpoints%nkp
448 :
449 1062 : ALLOCATE (fm_matrix_L_kpoints(nkp, 2))
450 498 : DO ikp = 1, nkp
451 1450 : DO i_real_imag = 1, 2
452 : CALL cp_fm_create(fm_matrix_L_kpoints(ikp, i_real_imag), &
453 952 : fm_matrix_Minv_L_kpoints(1, i_real_imag)%matrix_struct)
454 1428 : CALL cp_fm_set_all(fm_matrix_L_kpoints(ikp, i_real_imag), 0.0_dp)
455 : END DO
456 : END DO
457 :
458 22 : CALL allocate_matrix_v_RI_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
459 :
460 22 : IF (qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints) THEN
461 :
462 18 : nkp_orig = qs_env%mp2_env%ri_rpa_im_time%nkp_orig
463 18 : nkp_extra = qs_env%mp2_env%ri_rpa_im_time%nkp_extra
464 :
465 : CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", &
466 : cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
467 : atomic_kind_set=atomic_kind_set, &
468 : size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
469 18 : operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp_orig)
470 :
471 72 : nkp_grid_orig = kpoints%nkp_grid
472 144 : kpoints%nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid_extra(1:3)
473 :
474 : CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", &
475 : cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
476 : atomic_kind_set=atomic_kind_set, &
477 : size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
478 18 : operator_type=operator_coulomb, ikp_start=nkp_orig + 1, ikp_end=nkp)
479 :
480 72 : kpoints%nkp_grid = nkp_grid_orig
481 :
482 : ELSE
483 :
484 : CALL build_2c_coulomb_matrix_kp(matrix_v_RI_kp, kpoints, basis_type="RI_AUX", &
485 : cell=cell, particle_set=particle_set, qs_kind_set=qs_kind_set, &
486 : atomic_kind_set=atomic_kind_set, &
487 : size_lattice_sum=qs_env%mp2_env%mp2_gpw%size_lattice_sum, &
488 4 : operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp)
489 :
490 : END IF
491 :
492 498 : DO ikp = 1, nkp
493 :
494 476 : CALL copy_dbcsr_to_fm(matrix_v_RI_kp(ikp, 1)%matrix, fm_matrix_L_kpoints(ikp, 1))
495 498 : CALL copy_dbcsr_to_fm(matrix_v_RI_kp(ikp, 2)%matrix, fm_matrix_L_kpoints(ikp, 2))
496 :
497 : END DO
498 :
499 22 : CALL dbcsr_deallocate_matrix_set(matrix_v_RI_kp)
500 :
501 22 : CALL timestop(handle)
502 :
503 22 : END SUBROUTINE compute_V_by_lattice_sum
504 :
505 : ! **************************************************************************************************
506 : !> \brief ...
507 : !> \param matrix_v_RI_kp ...
508 : !> \param matrix_s_RI_aux_transl ...
509 : !> \param nkp ...
510 : ! **************************************************************************************************
511 22 : SUBROUTINE allocate_matrix_v_RI_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
512 :
513 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_v_RI_kp, matrix_s_RI_aux_transl
514 : INTEGER :: nkp
515 :
516 : INTEGER :: ikp
517 :
518 22 : NULLIFY (matrix_v_RI_kp)
519 22 : CALL dbcsr_allocate_matrix_set(matrix_v_RI_kp, nkp, 2)
520 :
521 498 : DO ikp = 1, nkp
522 :
523 476 : ALLOCATE (matrix_v_RI_kp(ikp, 1)%matrix)
524 : CALL dbcsr_create(matrix_v_RI_kp(ikp, 1)%matrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
525 476 : matrix_type=dbcsr_type_no_symmetry)
526 476 : CALL dbcsr_reserve_all_blocks(matrix_v_RI_kp(ikp, 1)%matrix)
527 476 : CALL dbcsr_set(matrix_v_RI_kp(ikp, 1)%matrix, 0.0_dp)
528 :
529 476 : ALLOCATE (matrix_v_RI_kp(ikp, 2)%matrix)
530 : CALL dbcsr_create(matrix_v_RI_kp(ikp, 2)%matrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
531 476 : matrix_type=dbcsr_type_no_symmetry)
532 476 : CALL dbcsr_reserve_all_blocks(matrix_v_RI_kp(ikp, 2)%matrix)
533 498 : CALL dbcsr_set(matrix_v_RI_kp(ikp, 2)%matrix, 0.0_dp)
534 :
535 : END DO
536 :
537 22 : END SUBROUTINE allocate_matrix_v_RI_kp
538 :
539 : ! **************************************************************************************************
540 : !> \brief ...
541 : !> \param qs_env ...
542 : !> \param fm_matrix_Minv_L_kpoints ...
543 : !> \param fm_matrix_L ...
544 : !> \param dimen_RI ...
545 : !> \param ri_metric ...
546 : !> \param do_kpoints ...
547 : !> \param kpoints ...
548 : !> \param put_mat_KS_env ...
549 : !> \param regularization_RI ...
550 : !> \param ikp_ext ...
551 : ! **************************************************************************************************
552 584 : SUBROUTINE RI_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L, dimen_RI, ri_metric, &
553 : do_kpoints, kpoints, put_mat_KS_env, regularization_RI, ikp_ext)
554 :
555 : TYPE(qs_environment_type), POINTER :: qs_env
556 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_kpoints
557 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
558 : INTEGER, INTENT(IN) :: dimen_RI
559 : TYPE(libint_potential_type) :: ri_metric
560 : LOGICAL, INTENT(IN) :: do_kpoints
561 : TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
562 : LOGICAL, OPTIONAL :: put_mat_KS_env
563 : REAL(KIND=dp), OPTIONAL :: regularization_RI
564 : INTEGER, OPTIONAL :: ikp_ext
565 :
566 : CHARACTER(LEN=*), PARAMETER :: routineN = 'RI_2c_integral_mat'
567 :
568 : INTEGER :: handle, i_real_imag, i_size, ikp, &
569 : ikp_for_xkp, img, n_real_imag, natom, &
570 : nimg, nkind, nkp
571 584 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI
572 584 : INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
573 584 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
574 : LOGICAL :: my_put_mat_KS_env
575 : REAL(KIND=dp) :: my_regularization_RI
576 584 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
577 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
578 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
579 : TYPE(cp_fm_type) :: fm_matrix_S_global
580 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
581 584 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_RI_aux_transl
582 584 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_2c
583 : TYPE(dbcsr_type), POINTER :: cmatrix, matrix_s_RI_aux_desymm, &
584 : rmatrix, tmpmat
585 : TYPE(dft_control_type), POINTER :: dft_control
586 : TYPE(distribution_2d_type), POINTER :: dist_2d
587 584 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI
588 : TYPE(mp_para_env_type), POINTER :: para_env
589 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
590 584 : POINTER :: sab_RI
591 584 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
592 584 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
593 :
594 584 : CALL timeset(routineN, handle)
595 :
596 584 : NULLIFY (sab_RI, matrix_s_RI_aux_transl, dist_2d)
597 :
598 584 : IF (PRESENT(regularization_RI)) THEN
599 584 : my_regularization_RI = regularization_RI
600 : ELSE
601 0 : my_regularization_RI = 0.0_dp
602 : END IF
603 :
604 584 : IF (PRESENT(put_mat_KS_env)) THEN
605 178 : my_put_mat_KS_env = put_mat_KS_env
606 : ELSE
607 : my_put_mat_KS_env = .FALSE.
608 : END IF
609 :
610 : CALL get_qs_env(qs_env=qs_env, &
611 : para_env=para_env, &
612 : blacs_env=blacs_env, &
613 : nkind=nkind, &
614 : distribution_2d=dist_2d, &
615 : qs_kind_set=qs_kind_set, &
616 : particle_set=particle_set, &
617 : dft_control=dft_control, &
618 584 : natom=natom)
619 :
620 1752 : ALLOCATE (sizes_RI(natom))
621 2908 : ALLOCATE (basis_set_RI(nkind))
622 :
623 584 : CALL basis_set_list_setup(basis_set_RI, 'RI_AUX', qs_kind_set)
624 584 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_RI)
625 :
626 : CALL build_2c_neighbor_lists(sab_RI, basis_set_RI, basis_set_RI, ri_metric, "RPA_2c_nl_RI", qs_env, &
627 584 : sym_ij=.TRUE., dist_2d=dist_2d)
628 :
629 584 : CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
630 1752 : ALLOCATE (row_bsize(SIZE(sizes_RI)))
631 1752 : ALLOCATE (col_bsize(SIZE(sizes_RI)))
632 2818 : row_bsize(:) = sizes_RI
633 2818 : col_bsize(:) = sizes_RI
634 :
635 584 : IF (do_kpoints) THEN
636 386 : nimg = dft_control%nimages
637 : ELSE
638 198 : nimg = 1
639 : END IF
640 :
641 12762 : ALLOCATE (mat_2c(nimg))
642 : CALL dbcsr_create(mat_2c(1), "(RI|RI)", dbcsr_dist, dbcsr_type_symmetric, &
643 584 : row_bsize, col_bsize)
644 584 : DEALLOCATE (row_bsize, col_bsize)
645 :
646 11010 : DO img = 2, nimg
647 11010 : CALL dbcsr_create(mat_2c(img), template=mat_2c(1))
648 : END DO
649 :
650 : CALL build_2c_integrals(mat_2c, 0.0_dp, qs_env, sab_RI, basis_set_RI, basis_set_RI, &
651 : ri_metric, do_kpoints=do_kpoints, ext_kpoints=kpoints, &
652 584 : regularization_RI=my_regularization_RI)
653 :
654 584 : CALL dbcsr_distribution_release(dbcsr_dist)
655 584 : DEALLOCATE (basis_set_RI)
656 :
657 584 : IF (my_put_mat_KS_env) THEN
658 134 : CALL get_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
659 134 : IF (ASSOCIATED(matrix_s_RI_aux_transl)) CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
660 : END IF
661 584 : CALL dbcsr_allocate_matrix_set(matrix_s_RI_aux_transl, 1, nimg)
662 11594 : DO img = 1, nimg
663 11010 : ALLOCATE (matrix_s_RI_aux_transl(1, img)%matrix)
664 11010 : CALL dbcsr_copy(matrix_s_RI_aux_transl(1, img)%matrix, mat_2c(img))
665 11594 : CALL dbcsr_release(mat_2c(img))
666 : END DO
667 :
668 584 : IF (my_put_mat_KS_env) THEN
669 134 : CALL set_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
670 : END IF
671 :
672 584 : IF (do_kpoints) THEN
673 386 : CPASSERT(PRESENT(kpoints))
674 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, &
675 386 : cell_to_index=cell_to_index)
676 386 : n_real_imag = 2
677 : ELSE
678 198 : nkp = 1
679 198 : n_real_imag = 1
680 : END IF
681 :
682 584 : IF (PRESENT(ikp_ext)) nkp = 1
683 :
684 5184 : ALLOCATE (fm_matrix_Minv_L_kpoints(nkp, n_real_imag))
685 1622 : DO i_size = 1, nkp
686 3500 : DO i_real_imag = 1, n_real_imag
687 1878 : CALL cp_fm_create(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), fm_matrix_L%matrix_struct)
688 2916 : CALL cp_fm_set_all(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), 0.0_dp)
689 : END DO
690 : END DO
691 :
692 584 : NULLIFY (fm_struct)
693 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
694 584 : ncol_global=dimen_RI, para_env=para_env)
695 :
696 584 : CALL cp_fm_create(fm_matrix_S_global, fm_struct)
697 584 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
698 :
699 584 : IF (do_kpoints) THEN
700 :
701 386 : ALLOCATE (rmatrix, cmatrix, tmpmat)
702 : CALL dbcsr_create(rmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
703 386 : matrix_type=dbcsr_type_symmetric)
704 : CALL dbcsr_create(cmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
705 386 : matrix_type=dbcsr_type_antisymmetric)
706 : CALL dbcsr_create(tmpmat, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
707 386 : matrix_type=dbcsr_type_no_symmetry)
708 386 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_RI)
709 386 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_RI)
710 :
711 1226 : DO ikp = 1, nkp
712 :
713 : ! s matrix is not spin dependent, double the work
714 840 : CALL dbcsr_set(rmatrix, 0.0_dp)
715 840 : CALL dbcsr_set(cmatrix, 0.0_dp)
716 :
717 840 : IF (PRESENT(ikp_ext)) THEN
718 364 : ikp_for_xkp = ikp_ext
719 : ELSE
720 : ikp_for_xkp = ikp
721 : END IF
722 :
723 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s_RI_aux_transl, ispin=1, &
724 840 : xkp=xkp(1:3, ikp_for_xkp), cell_to_index=cell_to_index, sab_nl=sab_RI)
725 :
726 840 : CALL dbcsr_set(tmpmat, 0.0_dp)
727 840 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
728 :
729 840 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
730 840 : CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
731 840 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 1), para_env)
732 :
733 840 : CALL dbcsr_set(tmpmat, 0.0_dp)
734 840 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
735 :
736 840 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
737 840 : CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
738 1226 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 2), para_env)
739 :
740 : END DO
741 :
742 386 : CALL dbcsr_deallocate_matrix(rmatrix)
743 386 : CALL dbcsr_deallocate_matrix(cmatrix)
744 386 : CALL dbcsr_deallocate_matrix(tmpmat)
745 :
746 : ELSE
747 :
748 : NULLIFY (matrix_s_RI_aux_desymm)
749 198 : ALLOCATE (matrix_s_RI_aux_desymm)
750 : CALL dbcsr_create(matrix_s_RI_aux_desymm, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
751 198 : name='S_RI non_symm', matrix_type=dbcsr_type_no_symmetry)
752 :
753 198 : CALL dbcsr_desymmetrize(matrix_s_RI_aux_transl(1, 1)%matrix, matrix_s_RI_aux_desymm)
754 :
755 198 : CALL copy_dbcsr_to_fm(matrix_s_RI_aux_desymm, fm_matrix_S_global)
756 :
757 198 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(1, 1), para_env)
758 :
759 198 : CALL dbcsr_deallocate_matrix(matrix_s_RI_aux_desymm)
760 :
761 : END IF
762 :
763 584 : CALL release_neighbor_list_sets(sab_RI)
764 :
765 584 : CALL cp_fm_struct_release(fm_struct)
766 :
767 584 : CALL cp_fm_release(fm_matrix_S_global)
768 :
769 584 : IF (.NOT. my_put_mat_KS_env) THEN
770 450 : CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
771 : END IF
772 :
773 584 : CALL timestop(handle)
774 :
775 1168 : END SUBROUTINE RI_2c_integral_mat
776 :
777 : ! **************************************************************************************************
778 : !> \brief ...
779 : !> \param qs_env ...
780 : !> \param eri_method ...
781 : !> \param eri_param ...
782 : !> \param para_env ...
783 : !> \param para_env_sub ...
784 : !> \param para_env_L ...
785 : !> \param mp2_memory ...
786 : !> \param fm_matrix_L ...
787 : !> \param ngroup ...
788 : !> \param color_sub ...
789 : !> \param dimen_RI ...
790 : !> \param my_group_L_size ...
791 : !> \param my_group_L_start ...
792 : !> \param my_group_L_end ...
793 : !> \param gd_array ...
794 : !> \param calc_PQ_cond_num ...
795 : !> \param cond_num ...
796 : !> \param num_small_eigen ...
797 : !> \param potential ...
798 : !> \param sab_orb_sub ...
799 : !> \param do_im_time ...
800 : !> \param fm_matrix_L_extern ...
801 : ! **************************************************************************************************
802 628 : SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
803 : fm_matrix_L, ngroup, color_sub, dimen_RI, &
804 : my_group_L_size, my_group_L_start, my_group_L_end, &
805 : gd_array, calc_PQ_cond_num, cond_num, num_small_eigen, potential, &
806 : sab_orb_sub, do_im_time, fm_matrix_L_extern)
807 :
808 : TYPE(qs_environment_type), POINTER :: qs_env
809 : INTEGER, INTENT(IN) :: eri_method
810 : TYPE(cp_eri_mme_param), POINTER :: eri_param
811 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
812 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
813 : TYPE(mp_para_env_type), INTENT(OUT), POINTER :: para_env_L
814 : REAL(KIND=dp), INTENT(IN) :: mp2_memory
815 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_L
816 : INTEGER, INTENT(IN) :: ngroup, color_sub
817 : INTEGER, INTENT(OUT) :: dimen_RI, my_group_L_size, &
818 : my_group_L_start, my_group_L_end
819 : TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
820 : LOGICAL, INTENT(IN) :: calc_PQ_cond_num
821 : REAL(KIND=dp), INTENT(OUT) :: cond_num
822 : INTEGER, INTENT(OUT) :: num_small_eigen
823 : TYPE(libint_potential_type), INTENT(IN) :: potential
824 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
825 : POINTER :: sab_orb_sub
826 : LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
827 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_L_extern
828 :
829 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_2c_integrals'
830 :
831 : INTEGER :: best_group_size, color_L, group_size, handle, handle2, i_global, iatom, iiB, &
832 : ikind, iproc, j_global, jjB, natom, ncol_local, nkind, nrow_local, nsgf, potential_type, &
833 : proc_receive, proc_receive_static, proc_send_static, proc_shift, rec_L_end, rec_L_size, &
834 : rec_L_start, strat_group_size
835 628 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
836 628 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
837 : LOGICAL :: my_do_im_time
838 : REAL(KIND=dp) :: min_mem_for_QK
839 628 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egen_L
840 628 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: L_external_col, L_local_col
841 628 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
842 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_L
843 : TYPE(cp_fm_type) :: fm_matrix_L_diag
844 628 : TYPE(group_dist_d1_type) :: gd_sub_array
845 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
846 628 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
847 628 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
848 :
849 628 : CALL timeset(routineN, handle)
850 :
851 628 : my_do_im_time = .FALSE.
852 628 : IF (PRESENT(do_im_time)) THEN
853 592 : my_do_im_time = do_im_time
854 : END IF
855 :
856 : ! get stuff
857 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
858 628 : particle_set=particle_set)
859 :
860 628 : nkind = SIZE(qs_kind_set)
861 628 : natom = SIZE(particle_set)
862 :
863 628 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
864 :
865 1780 : DO ikind = 1, nkind
866 1152 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
867 1780 : CPASSERT(ASSOCIATED(basis_set_a))
868 : END DO
869 :
870 628 : dimen_RI = 0
871 2474 : DO iatom = 1, natom
872 1846 : ikind = kind_of(iatom)
873 1846 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
874 2474 : dimen_RI = dimen_RI + nsgf
875 : END DO
876 :
877 : ! check that very small systems are not running on too many processes
878 628 : IF (dimen_RI < ngroup) THEN
879 : CALL cp_abort(__LOCATION__, "Product of block size and number "// &
880 0 : "of RI functions should not exceed total number of processes")
881 : END IF
882 :
883 628 : CALL create_group_dist(gd_array, ngroup, dimen_RI)
884 :
885 628 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
886 :
887 628 : CALL timeset(routineN//"_loop_lm", handle2)
888 :
889 2512 : ALLOCATE (L_local_col(dimen_RI, my_group_L_size))
890 2163416 : L_local_col = 0.0_dp
891 :
892 628 : potential_type = potential%potential_type
893 :
894 : IF ((eri_method .EQ. do_eri_mme .OR. eri_method .EQ. do_eri_os) &
895 628 : .AND. potential_type .EQ. do_potential_coulomb .OR. &
896 : (eri_method .EQ. do_eri_mme .AND. potential_type .EQ. do_potential_long)) THEN
897 :
898 : CALL mp2_eri_2c_integrate(eri_param, potential, para_env_sub, qs_env, &
899 : basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
900 : hab=L_local_col, first_b=my_group_L_start, last_b=my_group_L_end, &
901 260 : eri_method=eri_method)
902 :
903 : ELSEIF (eri_method .EQ. do_eri_gpw .OR. &
904 : (potential_type == do_potential_long .AND. qs_env%mp2_env%eri_method == do_eri_os) &
905 368 : .OR. (potential_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)) THEN
906 :
907 : CALL mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, &
908 368 : natom, potential, sab_orb_sub, L_local_col, kind_of)
909 :
910 : ELSE
911 0 : CPABORT("unknown ERI method")
912 : END IF
913 :
914 628 : CALL timestop(handle2)
915 :
916 : ! split the total number of proc in a subgroup of the size of ~1/10 of the
917 : ! total num of proc
918 628 : best_group_size = para_env%num_pe
919 :
920 628 : strat_group_size = MAX(1, para_env%num_pe/10)
921 :
922 628 : min_mem_for_QK = REAL(dimen_RI, KIND=dp)*dimen_RI*3.0_dp*8.0_dp/1024_dp/1024_dp
923 :
924 628 : group_size = strat_group_size - 1
925 662 : DO iproc = strat_group_size, para_env%num_pe
926 662 : group_size = group_size + 1
927 : ! check that group_size is a multiple of sub_group_size and a divisor of
928 : ! the total num of proc
929 662 : IF (MOD(para_env%num_pe, group_size) /= 0 .OR. MOD(group_size, para_env_sub%num_pe) /= 0) CYCLE
930 :
931 : ! check for memory
932 628 : IF (REAL(group_size, KIND=dp)*mp2_memory < min_mem_for_QK) CYCLE
933 :
934 : best_group_size = group_size
935 34 : EXIT
936 : END DO
937 :
938 628 : IF (my_do_im_time) THEN
939 : ! para_env_L is the para_env for im. time to avoid bug
940 134 : best_group_size = para_env%num_pe
941 : END IF
942 :
943 : ! create the L group
944 628 : color_L = para_env%mepos/best_group_size
945 628 : ALLOCATE (para_env_L)
946 628 : CALL para_env_L%from_split(para_env, color_L)
947 :
948 : ! create the blacs_L
949 628 : NULLIFY (blacs_env_L)
950 628 : CALL cp_blacs_env_create(blacs_env=blacs_env_L, para_env=para_env_L)
951 :
952 : ! create the full matrix L defined in the L group
953 628 : CALL create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, "fm_matrix_L", fm_matrix_L_extern)
954 :
955 628 : IF (my_do_im_time .AND. para_env%num_pe > 1) THEN
956 :
957 : CALL fill_fm_L_from_L_loc_non_blocking(fm_matrix_L, L_local_col, para_env, &
958 : my_group_L_start, my_group_L_end, &
959 134 : dimen_RI)
960 :
961 : ELSE
962 494 : BLOCK
963 : TYPE(mp_comm_type) :: comm_exchange
964 :
965 494 : CALL comm_exchange%from_split(para_env_L, para_env_sub%mepos)
966 :
967 : CALL create_group_dist(gd_sub_array, my_group_L_start, &
968 494 : my_group_L_end, my_group_L_size, comm_exchange)
969 :
970 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
971 : nrow_local=nrow_local, &
972 : ncol_local=ncol_local, &
973 : row_indices=row_indices, &
974 494 : col_indices=col_indices)
975 :
976 37180 : DO jjB = 1, ncol_local
977 36686 : j_global = col_indices(jjB)
978 37180 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
979 1508996 : DO iiB = 1, nrow_local
980 1489533 : i_global = row_indices(iiB)
981 1508996 : fm_matrix_L%local_data(iiB, jjB) = L_local_col(i_global, j_global - my_group_L_start + 1)
982 : END DO
983 : END IF
984 : END DO
985 :
986 494 : proc_send_static = MODULO(comm_exchange%mepos + 1, comm_exchange%num_pe)
987 494 : proc_receive_static = MODULO(comm_exchange%mepos - 1, comm_exchange%num_pe)
988 :
989 494 : DO proc_shift = 1, comm_exchange%num_pe - 1
990 0 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
991 :
992 0 : CALL get_group_dist(gd_sub_array, proc_receive, rec_L_start, rec_L_end, rec_L_size)
993 :
994 0 : ALLOCATE (L_external_col(dimen_RI, rec_L_size))
995 0 : L_external_col = 0.0_dp
996 :
997 0 : CALL comm_exchange%sendrecv(L_local_col, proc_send_static, L_external_col, proc_receive_static)
998 :
999 0 : DO jjB = 1, ncol_local
1000 0 : j_global = col_indices(jjB)
1001 0 : IF (j_global >= rec_L_start .AND. j_global <= rec_L_end) THEN
1002 0 : DO iiB = 1, nrow_local
1003 0 : i_global = row_indices(iiB)
1004 0 : fm_matrix_L%local_data(iiB, jjB) = L_external_col(i_global, j_global - rec_L_start + 1)
1005 : END DO
1006 : END IF
1007 : END DO
1008 :
1009 494 : CALL MOVE_ALLOC(L_external_col, L_local_col)
1010 : END DO
1011 :
1012 494 : CALL release_group_dist(gd_sub_array)
1013 988 : CALL comm_exchange%free()
1014 : END BLOCK
1015 : END IF
1016 :
1017 628 : DEALLOCATE (L_local_col)
1018 :
1019 : ! create the new group for the sum of the local data over all processes
1020 : BLOCK
1021 : TYPE(mp_comm_type) :: comm_exchange
1022 628 : comm_exchange = fm_matrix_L%matrix_struct%context%interconnect(para_env)
1023 628 : CALL comm_exchange%sum(fm_matrix_L%local_data)
1024 1256 : CALL comm_exchange%free()
1025 : END BLOCK
1026 :
1027 628 : cond_num = 1.0_dp
1028 628 : num_small_eigen = 0
1029 628 : IF (calc_PQ_cond_num) THEN
1030 : ! calculate the condition number of the (P|Q) matrix
1031 : ! create a copy of the matrix
1032 :
1033 0 : CALL create_matrix_L(fm_matrix_L_diag, blacs_env_L, dimen_RI, para_env_L, "fm_matrix_L_diag", fm_matrix_L_extern)
1034 :
1035 0 : CALL cp_fm_to_fm(source=fm_matrix_L, destination=fm_matrix_L_diag)
1036 :
1037 0 : ALLOCATE (egen_L(dimen_RI))
1038 :
1039 0 : egen_L = 0.0_dp
1040 0 : CALL cp_fm_syevx(matrix=fm_matrix_L_diag, eigenvalues=egen_L)
1041 :
1042 0 : num_small_eigen = 0
1043 0 : DO iiB = 1, dimen_RI
1044 0 : IF (ABS(egen_L(iiB)) < 0.001_dp) num_small_eigen = num_small_eigen + 1
1045 : END DO
1046 :
1047 0 : cond_num = MAXVAL(ABS(egen_L))/MINVAL(ABS(egen_L))
1048 :
1049 0 : CALL cp_fm_release(fm_matrix_L_diag)
1050 :
1051 0 : DEALLOCATE (egen_L)
1052 : END IF
1053 :
1054 : ! release blacs_env
1055 628 : CALL cp_blacs_env_release(blacs_env_L)
1056 :
1057 628 : CALL timestop(handle)
1058 :
1059 2512 : END SUBROUTINE compute_2c_integrals
1060 :
1061 : ! **************************************************************************************************
1062 : !> \brief ...
1063 : !> \param matrix ...
1064 : !> \param num_small_evals ...
1065 : !> \param cond_num ...
1066 : !> \param do_inversion ...
1067 : !> \param para_env ...
1068 : ! **************************************************************************************************
1069 44 : SUBROUTINE matrix_root_with_svd(matrix, num_small_evals, cond_num, do_inversion, para_env)
1070 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1071 : INTEGER, INTENT(OUT) :: num_small_evals
1072 : REAL(KIND=dp), INTENT(OUT) :: cond_num
1073 : LOGICAL, INTENT(IN) :: do_inversion
1074 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1075 :
1076 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_root_with_svd'
1077 :
1078 : INTEGER :: group_size_L, handle, ii, needed_evals, &
1079 : nrow, pos_max(1)
1080 44 : INTEGER, ALLOCATABLE, DIMENSION(:) :: num_eval
1081 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
1082 : TYPE(cp_fm_type) :: evecs
1083 : TYPE(mp_comm_type) :: comm_exchange
1084 :
1085 44 : CALL timeset(routineN, handle)
1086 :
1087 : ! Create arrays carrying eigenvalues and eigenvectors later
1088 44 : CALL cp_fm_get_info(matrix=matrix, nrow_global=nrow)
1089 :
1090 132 : ALLOCATE (evals(nrow))
1091 3724 : evals = 0
1092 :
1093 44 : CALL cp_fm_create(evecs, matrix%matrix_struct)
1094 :
1095 : ! Perform the EVD (=SVD of a positive semidefinite matrix)
1096 44 : CALL choose_eigv_solver(matrix, evecs, evals)
1097 :
1098 : ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
1099 44 : num_small_evals = 0
1100 210 : DO ii = 1, nrow
1101 210 : IF (evals(ii) > 0.0_dp) THEN
1102 44 : num_small_evals = ii - 1
1103 44 : EXIT
1104 : END IF
1105 : END DO
1106 44 : needed_evals = nrow - num_small_evals
1107 :
1108 : ! Get the condition number w.r.t. considered eigenvalues
1109 44 : cond_num = evals(nrow)/evals(num_small_evals + 1)
1110 :
1111 : ! Determine the eigenvalues of the request matrix root or its inverse
1112 210 : evals(1:num_small_evals) = 0.0_dp
1113 44 : IF (do_inversion) THEN
1114 1008 : evals(num_small_evals + 1:nrow) = 1.0_dp/SQRT(evals(num_small_evals + 1:nrow))
1115 : ELSE
1116 2550 : evals(num_small_evals + 1:nrow) = SQRT(evals(num_small_evals + 1:nrow))
1117 : END IF
1118 :
1119 44 : CALL cp_fm_column_scale(evecs, evals)
1120 :
1121 : ! As it turns out, the results in the subgroups differ.
1122 : ! Thus, we have to equilibrate the results.
1123 : ! Step 1: Get a communicator connecting processes with same id within subgroups
1124 : ! use that group_size_L is divisible by the total number of processes (see above)
1125 44 : group_size_L = para_env%num_pe/matrix%matrix_struct%para_env%num_pe
1126 44 : comm_exchange = matrix%matrix_struct%context%interconnect(para_env)
1127 :
1128 132 : ALLOCATE (num_eval(0:group_size_L - 1))
1129 44 : CALL comm_exchange%allgather(num_small_evals, num_eval)
1130 :
1131 118 : num_small_evals = MINVAL(num_eval)
1132 :
1133 118 : IF (num_small_evals /= MAXVAL(num_eval)) THEN
1134 : ! Step 2: Get position of maximum value
1135 0 : pos_max = MAXLOC(num_eval)
1136 0 : num_small_evals = num_eval(pos_max(1))
1137 0 : needed_evals = nrow - num_small_evals
1138 :
1139 : ! Step 3: Broadcast your local data to all other processes
1140 0 : CALL comm_exchange%bcast(evecs%local_data, pos_max(1))
1141 0 : CALL comm_exchange%bcast(cond_num, pos_max(1))
1142 : END IF
1143 :
1144 44 : DEALLOCATE (num_eval)
1145 :
1146 44 : CALL comm_exchange%free()
1147 :
1148 44 : CALL reset_size_matrix(matrix, needed_evals, matrix%matrix_struct)
1149 :
1150 : ! Copy the needed eigenvectors
1151 44 : CALL cp_fm_to_fm(evecs, matrix, needed_evals, num_small_evals + 1)
1152 :
1153 44 : CALL cp_fm_release(evecs)
1154 :
1155 44 : CALL timestop(handle)
1156 :
1157 88 : END SUBROUTINE matrix_root_with_svd
1158 :
1159 : ! **************************************************************************************************
1160 : !> \brief ...
1161 : !> \param matrix ...
1162 : !> \param new_size ...
1163 : !> \param fm_struct_template ...
1164 : ! **************************************************************************************************
1165 76 : SUBROUTINE reset_size_matrix(matrix, new_size, fm_struct_template)
1166 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1167 : INTEGER, INTENT(IN) :: new_size
1168 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_template
1169 :
1170 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_size_matrix'
1171 :
1172 : INTEGER :: handle
1173 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1174 :
1175 76 : CALL timeset(routineN, handle)
1176 :
1177 : ! Choose the block sizes as large as in the template for the later steps
1178 76 : NULLIFY (fm_struct)
1179 76 : CALL cp_fm_struct_create(fm_struct, ncol_global=new_size, template_fmstruct=fm_struct_template, force_block=.TRUE.)
1180 :
1181 76 : CALL cp_fm_release(matrix)
1182 :
1183 76 : CALL cp_fm_create(matrix, fm_struct)
1184 76 : CALL cp_fm_set_all(matrix, 0.0_dp)
1185 :
1186 76 : CALL cp_fm_struct_release(fm_struct)
1187 :
1188 76 : CALL timestop(handle)
1189 :
1190 76 : END SUBROUTINE reset_size_matrix
1191 :
1192 : ! **************************************************************************************************
1193 : !> \brief ...
1194 : !> \param fm_matrix_L ...
1195 : !> \param L_local_col ...
1196 : !> \param para_env ...
1197 : !> \param my_group_L_start ...
1198 : !> \param my_group_L_end ...
1199 : !> \param dimen_RI ...
1200 : ! **************************************************************************************************
1201 134 : SUBROUTINE fill_fm_L_from_L_loc_non_blocking(fm_matrix_L, L_local_col, para_env, my_group_L_start, &
1202 : my_group_L_end, dimen_RI)
1203 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1204 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1205 : INTENT(IN) :: L_local_col
1206 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1207 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, &
1208 : dimen_RI
1209 :
1210 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_fm_L_from_L_loc_non_blocking'
1211 :
1212 : INTEGER :: dummy_proc, handle, handle2, i_entry_rec, i_row, i_row_global, iproc, j_col, &
1213 : j_col_global, LLL, MMM, ncol_block, ncol_local, npcol, nprow, nrow_block, nrow_local, &
1214 : proc_send, send_pcol, send_prow
1215 134 : INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
1216 : num_entries_send
1217 134 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1218 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1219 134 : DIMENSION(:) :: buffer_rec, buffer_send
1220 134 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
1221 :
1222 134 : CALL timeset(routineN, handle)
1223 :
1224 134 : CALL timeset(routineN//"_1", handle2)
1225 :
1226 : ! get info for the dest
1227 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1228 : nrow_local=nrow_local, &
1229 : ncol_local=ncol_local, &
1230 : row_indices=row_indices, &
1231 : col_indices=col_indices, &
1232 : nrow_block=nrow_block, &
1233 134 : ncol_block=ncol_block)
1234 :
1235 134 : nprow = fm_matrix_L%matrix_struct%context%num_pe(1)
1236 134 : npcol = fm_matrix_L%matrix_struct%context%num_pe(2)
1237 :
1238 402 : ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
1239 402 : ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
1240 :
1241 402 : num_entries_rec(:) = 0
1242 402 : num_entries_send(:) = 0
1243 :
1244 134 : dummy_proc = 0
1245 :
1246 : ! get the process, where the elements from L_local_col have to be sent
1247 11380 : DO LLL = 1, dimen_RI
1248 :
1249 : send_prow = cp_fm_indxg2p(LLL, nrow_block, dummy_proc, &
1250 11246 : fm_matrix_L%matrix_struct%first_p_pos(1), nprow)
1251 :
1252 561929 : DO MMM = my_group_L_start, my_group_L_end
1253 :
1254 : send_pcol = cp_fm_indxg2p(MMM, ncol_block, dummy_proc, &
1255 550549 : fm_matrix_L%matrix_struct%first_p_pos(2), npcol)
1256 :
1257 550549 : proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1258 :
1259 561795 : num_entries_send(proc_send) = num_entries_send(proc_send) + 1
1260 :
1261 : END DO
1262 :
1263 : END DO
1264 :
1265 134 : CALL timestop(handle2)
1266 :
1267 134 : CALL timeset(routineN//"_2", handle2)
1268 :
1269 134 : CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
1270 :
1271 134 : CALL timestop(handle2)
1272 :
1273 134 : CALL timeset(routineN//"_3", handle2)
1274 :
1275 : ! allocate buffers to send the entries and the information of the entries
1276 670 : ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
1277 670 : ALLOCATE (buffer_send(0:para_env%num_pe - 1))
1278 :
1279 : ! allocate data message and corresponding indices
1280 402 : DO iproc = 0, para_env%num_pe - 1
1281 :
1282 804 : ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
1283 550951 : buffer_rec(iproc)%msg = 0.0_dp
1284 :
1285 : END DO
1286 :
1287 134 : CALL timestop(handle2)
1288 :
1289 134 : CALL timeset(routineN//"_4", handle2)
1290 :
1291 402 : DO iproc = 0, para_env%num_pe - 1
1292 :
1293 804 : ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
1294 550951 : buffer_send(iproc)%msg = 0.0_dp
1295 :
1296 : END DO
1297 :
1298 134 : CALL timestop(handle2)
1299 :
1300 134 : CALL timeset(routineN//"_5", handle2)
1301 :
1302 402 : DO iproc = 0, para_env%num_pe - 1
1303 :
1304 804 : ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
1305 1102036 : buffer_rec(iproc)%indx = 0
1306 :
1307 : END DO
1308 :
1309 134 : CALL timestop(handle2)
1310 :
1311 134 : CALL timeset(routineN//"_6", handle2)
1312 :
1313 402 : DO iproc = 0, para_env%num_pe - 1
1314 :
1315 804 : ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
1316 1102036 : buffer_send(iproc)%indx = 0
1317 :
1318 : END DO
1319 :
1320 134 : CALL timestop(handle2)
1321 :
1322 134 : CALL timeset(routineN//"_7", handle2)
1323 :
1324 402 : ALLOCATE (entry_counter(0:para_env%num_pe - 1))
1325 402 : entry_counter(:) = 0
1326 :
1327 : ! get the process, where the elements from L_local_col have to be sent and
1328 : ! write the data and the info to buffer_send
1329 11380 : DO LLL = 1, dimen_RI
1330 :
1331 : send_prow = cp_fm_indxg2p(LLL, nrow_block, dummy_proc, &
1332 11246 : fm_matrix_L%matrix_struct%first_p_pos(1), nprow)
1333 :
1334 561929 : DO MMM = my_group_L_start, my_group_L_end
1335 :
1336 : send_pcol = cp_fm_indxg2p(MMM, ncol_block, dummy_proc, &
1337 550549 : fm_matrix_L%matrix_struct%first_p_pos(2), npcol)
1338 :
1339 550549 : proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1340 :
1341 550549 : entry_counter(proc_send) = entry_counter(proc_send) + 1
1342 :
1343 : buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
1344 550549 : L_local_col(LLL, MMM - my_group_L_start + 1)
1345 :
1346 550549 : buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = LLL
1347 561795 : buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = MMM
1348 :
1349 : END DO
1350 :
1351 : END DO
1352 :
1353 2010 : ALLOCATE (req_array(1:para_env%num_pe, 4))
1354 :
1355 134 : CALL timestop(handle2)
1356 :
1357 134 : CALL timeset(routineN//"_8", handle2)
1358 :
1359 : ! communicate the buffer
1360 : CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, &
1361 134 : buffer_send, req_array)
1362 :
1363 534373 : fm_matrix_L%local_data = 0.0_dp
1364 :
1365 134 : CALL timestop(handle2)
1366 :
1367 134 : CALL timeset(routineN//"_9", handle2)
1368 :
1369 : ! fill fm_matrix_L with the entries from buffer_rec
1370 402 : DO iproc = 0, para_env%num_pe - 1
1371 :
1372 550951 : DO i_entry_rec = 1, num_entries_rec(iproc)
1373 :
1374 31219430 : DO i_row = 1, nrow_local
1375 :
1376 30668613 : i_row_global = row_indices(i_row)
1377 :
1378 4667031343 : DO j_col = 1, ncol_local
1379 :
1380 4635812181 : j_col_global = col_indices(j_col)
1381 :
1382 4635812181 : IF (i_row_global == buffer_rec(iproc)%indx(i_entry_rec, 1) .AND. &
1383 30668613 : j_col_global == buffer_rec(iproc)%indx(i_entry_rec, 2)) THEN
1384 :
1385 550549 : fm_matrix_L%local_data(i_row, j_col) = buffer_rec(iproc)%msg(i_entry_rec)
1386 :
1387 : END IF
1388 :
1389 : END DO
1390 :
1391 : END DO
1392 :
1393 : END DO
1394 :
1395 : END DO
1396 :
1397 134 : CALL timestop(handle2)
1398 :
1399 134 : CALL timeset(routineN//"_10", handle2)
1400 :
1401 402 : DO iproc = 0, para_env%num_pe - 1
1402 268 : DEALLOCATE (buffer_rec(iproc)%msg)
1403 268 : DEALLOCATE (buffer_rec(iproc)%indx)
1404 268 : DEALLOCATE (buffer_send(iproc)%msg)
1405 402 : DEALLOCATE (buffer_send(iproc)%indx)
1406 : END DO
1407 :
1408 670 : DEALLOCATE (buffer_rec, buffer_send)
1409 134 : DEALLOCATE (req_array)
1410 134 : DEALLOCATE (entry_counter)
1411 134 : DEALLOCATE (num_entries_rec, num_entries_send)
1412 :
1413 134 : CALL timestop(handle2)
1414 :
1415 134 : CALL timestop(handle)
1416 :
1417 1340 : END SUBROUTINE fill_fm_L_from_L_loc_non_blocking
1418 :
1419 : ! **************************************************************************************************
1420 : !> \brief ...
1421 : !> \param fm_matrix_L ...
1422 : !> \param dimen_RI ...
1423 : !> \param do_inversion ...
1424 : ! **************************************************************************************************
1425 1416 : SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion)
1426 :
1427 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1428 : INTEGER, INTENT(IN) :: dimen_RI
1429 : LOGICAL, INTENT(IN) :: do_inversion
1430 :
1431 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholesky_decomp'
1432 :
1433 : INTEGER :: handle, i_global, iiB, info_chol, &
1434 : j_global, jjB, ncol_local, nrow_local
1435 708 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1436 :
1437 708 : CALL timeset(routineN, handle)
1438 :
1439 : ! do cholesky decomposition
1440 708 : CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=dimen_RI, info_out=info_chol)
1441 708 : CPASSERT(info_chol == 0)
1442 :
1443 708 : IF (.NOT. do_inversion) THEN
1444 : ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1445 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1446 : nrow_local=nrow_local, &
1447 : ncol_local=ncol_local, &
1448 : row_indices=row_indices, &
1449 116 : col_indices=col_indices)
1450 5803 : DO iiB = 1, nrow_local
1451 5687 : i_global = row_indices(iiB)
1452 538330 : DO jjB = 1, ncol_local
1453 532527 : j_global = col_indices(jjB)
1454 538214 : IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
1455 : END DO
1456 : END DO
1457 :
1458 : END IF
1459 :
1460 708 : CALL timestop(handle)
1461 :
1462 708 : END SUBROUTINE cholesky_decomp
1463 :
1464 : ! **************************************************************************************************
1465 : !> \brief ...
1466 : !> \param fm_matrix_Minv_L_kpoints ...
1467 : !> \param fm_matrix_L_kpoints ...
1468 : !> \param dimen_RI ...
1469 : !> \param kpoints ...
1470 : !> \param eps_eigval_S ...
1471 : ! **************************************************************************************************
1472 22 : SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, &
1473 : dimen_RI, kpoints, eps_eigval_S)
1474 :
1475 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: fm_matrix_Minv_L_kpoints, &
1476 : fm_matrix_L_kpoints
1477 : INTEGER, INTENT(IN) :: dimen_RI
1478 : TYPE(kpoint_type), POINTER :: kpoints
1479 : REAL(KIND=dp), INTENT(IN) :: eps_eigval_S
1480 :
1481 : CHARACTER(LEN=*), PARAMETER :: routineN = 'inversion_of_M_and_mult_with_chol_dec_of_V'
1482 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
1483 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp), ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
1484 :
1485 : INTEGER :: handle, ikp, nkp
1486 : TYPE(cp_cfm_type) :: cfm_matrix_K_tmp, cfm_matrix_M_tmp, &
1487 : cfm_matrix_V_tmp, cfm_matrix_Vtrunc_tmp
1488 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1489 :
1490 22 : CALL timeset(routineN, handle)
1491 :
1492 22 : CALL cp_fm_get_info(fm_matrix_Minv_L_kpoints(1, 1), matrix_struct=matrix_struct)
1493 :
1494 22 : CALL cp_cfm_create(cfm_matrix_M_tmp, matrix_struct)
1495 22 : CALL cp_cfm_create(cfm_matrix_V_tmp, matrix_struct)
1496 22 : CALL cp_cfm_create(cfm_matrix_K_tmp, matrix_struct)
1497 22 : CALL cp_cfm_create(cfm_matrix_Vtrunc_tmp, matrix_struct)
1498 :
1499 22 : CALL get_kpoint_info(kpoints, nkp=nkp)
1500 :
1501 498 : DO ikp = 1, nkp
1502 :
1503 476 : CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_M_tmp, cone, fm_matrix_Minv_L_kpoints(ikp, 1))
1504 476 : CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_M_tmp, ione, fm_matrix_Minv_L_kpoints(ikp, 2))
1505 :
1506 476 : CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_V_tmp, cone, fm_matrix_L_kpoints(ikp, 1))
1507 476 : CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_V_tmp, ione, fm_matrix_L_kpoints(ikp, 2))
1508 :
1509 476 : CALL cp_cfm_power(cfm_matrix_M_tmp, threshold=eps_eigval_S, exponent=-1.0_dp)
1510 :
1511 476 : CALL cp_cfm_power(cfm_matrix_V_tmp, threshold=0.0_dp, exponent=0.5_dp)
1512 :
1513 : ! save L(k) = SQRT(V(k))
1514 476 : CALL cp_cfm_to_fm(cfm_matrix_V_tmp, fm_matrix_L_kpoints(ikp, 1), fm_matrix_L_kpoints(ikp, 2))
1515 :
1516 : ! get K = M^(-1)*L, where L is the Cholesky decomposition of V = L^T*L
1517 : CALL parallel_gemm("N", "C", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_M_tmp, cfm_matrix_V_tmp, &
1518 476 : czero, cfm_matrix_K_tmp)
1519 :
1520 : ! move
1521 498 : CALL cp_cfm_to_fm(cfm_matrix_K_tmp, fm_matrix_Minv_L_kpoints(ikp, 1), fm_matrix_Minv_L_kpoints(ikp, 2))
1522 :
1523 : END DO
1524 :
1525 22 : CALL cp_cfm_release(cfm_matrix_M_tmp)
1526 22 : CALL cp_cfm_release(cfm_matrix_V_tmp)
1527 22 : CALL cp_cfm_release(cfm_matrix_K_tmp)
1528 22 : CALL cp_cfm_release(cfm_matrix_Vtrunc_tmp)
1529 :
1530 22 : CALL timestop(handle)
1531 :
1532 22 : END SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V
1533 :
1534 : ! **************************************************************************************************
1535 : !> \brief ...
1536 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
1537 : !> \param fm_matrix_Minv ...
1538 : !> \param qs_env ...
1539 : ! **************************************************************************************************
1540 44 : SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
1541 22 : fm_matrix_Minv, qs_env)
1542 :
1543 : TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_Minv_Vtrunc_Minv, &
1544 : fm_matrix_Minv
1545 : TYPE(qs_environment_type), POINTER :: qs_env
1546 :
1547 : CHARACTER(LEN=*), PARAMETER :: &
1548 : routineN = 'Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc'
1549 :
1550 : INTEGER :: dimen_RI, handle, ndep
1551 : REAL(KIND=dp) :: eps_eigval_S_Gamma
1552 : TYPE(cp_fm_type) :: fm_matrix_RI_metric_inv_work, fm_work
1553 :
1554 22 : CALL timeset(routineN, handle)
1555 :
1556 22 : CALL cp_fm_create(fm_work, fm_matrix_Minv(1, 1)%matrix_struct)
1557 22 : CALL cp_fm_set_all(fm_work, 0.0_dp)
1558 :
1559 22 : CALL cp_fm_create(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1)%matrix_struct)
1560 22 : CALL cp_fm_set_all(fm_matrix_RI_metric_inv_work, 0.0_dp)
1561 :
1562 22 : CALL cp_fm_get_info(matrix=fm_matrix_Minv(1, 1), nrow_global=dimen_RI)
1563 :
1564 22 : eps_eigval_S_Gamma = qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S_Gamma
1565 :
1566 22 : IF (eps_eigval_S_Gamma > 1.0E-18) THEN
1567 :
1568 : ! remove small eigenvalues
1569 : CALL cp_fm_power(fm_matrix_Minv(1, 1), fm_matrix_RI_metric_inv_work, -0.5_dp, &
1570 0 : eps_eigval_S_Gamma, ndep)
1571 :
1572 : ELSE
1573 :
1574 22 : CALL cholesky_decomp(fm_matrix_Minv(1, 1), dimen_RI, do_inversion=.TRUE.)
1575 :
1576 22 : CALL invert_mat(fm_matrix_Minv(1, 1))
1577 :
1578 : END IF
1579 :
1580 : CALL parallel_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_Minv(1, 1), &
1581 22 : fm_matrix_Minv(1, 1), 0.0_dp, fm_matrix_RI_metric_inv_work)
1582 :
1583 22 : CALL cp_fm_to_fm(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1))
1584 :
1585 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_RI_metric_inv_work, &
1586 22 : fm_matrix_Minv_Vtrunc_Minv(1, 1), 0.0_dp, fm_work)
1587 :
1588 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work, &
1589 22 : fm_matrix_RI_metric_inv_work, 0.0_dp, fm_matrix_Minv_Vtrunc_Minv(1, 1))
1590 :
1591 22 : CALL cp_fm_release(fm_work)
1592 22 : CALL cp_fm_release(fm_matrix_RI_metric_inv_work)
1593 :
1594 22 : CALL timestop(handle)
1595 :
1596 22 : END SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc
1597 :
1598 : ! **************************************************************************************************
1599 : !> \brief ...
1600 : !> \param qs_env ...
1601 : !> \param trunc_coulomb ...
1602 : !> \param rel_cutoff_trunc_coulomb_ri_x ...
1603 : ! **************************************************************************************************
1604 40 : SUBROUTINE setup_trunc_coulomb_pot_for_exchange_self_energy(qs_env, trunc_coulomb, &
1605 : rel_cutoff_trunc_coulomb_ri_x)
1606 : TYPE(qs_environment_type), POINTER :: qs_env
1607 : TYPE(libint_potential_type), OPTIONAL :: trunc_coulomb
1608 : REAL(KIND=dp), OPTIONAL :: rel_cutoff_trunc_coulomb_ri_x
1609 :
1610 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_trunc_coulomb_pot_for_exchange_self_energy'
1611 :
1612 : INTEGER :: handle
1613 : INTEGER, DIMENSION(3) :: periodic
1614 : REAL(KIND=dp) :: my_rel_cutoff_trunc_coulomb_ri_x, shortest_dist_cell_planes
1615 : TYPE(cell_type), POINTER :: cell
1616 :
1617 40 : CALL timeset(routineN, handle)
1618 :
1619 40 : NULLIFY (cell)
1620 40 : CALL get_qs_env(qs_env, cell=cell)
1621 :
1622 40 : CALL get_cell(cell=cell, periodic=periodic)
1623 :
1624 40 : shortest_dist_cell_planes = 1.0E4_dp
1625 40 : IF (periodic(1) == 1) THEN
1626 16 : IF (shortest_dist_cell_planes > plane_distance(1, 0, 0, cell)) THEN
1627 16 : shortest_dist_cell_planes = plane_distance(1, 0, 0, cell)
1628 : END IF
1629 : END IF
1630 40 : IF (periodic(2) == 1) THEN
1631 32 : IF (shortest_dist_cell_planes > plane_distance(0, 1, 0, cell)) THEN
1632 20 : shortest_dist_cell_planes = plane_distance(0, 1, 0, cell)
1633 : END IF
1634 : END IF
1635 40 : IF (periodic(3) == 1) THEN
1636 28 : IF (shortest_dist_cell_planes > plane_distance(0, 0, 1, cell)) THEN
1637 8 : shortest_dist_cell_planes = plane_distance(0, 0, 1, cell)
1638 : END IF
1639 : END IF
1640 :
1641 40 : IF (PRESENT(rel_cutoff_trunc_coulomb_ri_x)) THEN
1642 18 : my_rel_cutoff_trunc_coulomb_ri_x = rel_cutoff_trunc_coulomb_ri_x
1643 : ELSE
1644 22 : my_rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
1645 : END IF
1646 :
1647 40 : IF (PRESENT(trunc_coulomb)) THEN
1648 40 : trunc_coulomb%potential_type = do_potential_truncated
1649 40 : trunc_coulomb%cutoff_radius = shortest_dist_cell_planes*my_rel_cutoff_trunc_coulomb_ri_x
1650 40 : trunc_coulomb%filename = "t_c_g.dat"
1651 : ! dummy
1652 40 : trunc_coulomb%omega = 0.0_dp
1653 : END IF
1654 :
1655 40 : CALL timestop(handle)
1656 :
1657 40 : END SUBROUTINE setup_trunc_coulomb_pot_for_exchange_self_energy
1658 :
1659 : ! **************************************************************************************************
1660 : !> \brief ...
1661 : !> \param fm_matrix_L ...
1662 : ! **************************************************************************************************
1663 1184 : SUBROUTINE invert_mat(fm_matrix_L)
1664 :
1665 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1666 :
1667 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_mat'
1668 :
1669 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
1670 : ncol_local, nrow_local
1671 592 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1672 :
1673 592 : CALL timeset(routineN, handle)
1674 :
1675 592 : CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U')
1676 :
1677 : ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1678 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1679 : nrow_local=nrow_local, &
1680 : ncol_local=ncol_local, &
1681 : row_indices=row_indices, &
1682 592 : col_indices=col_indices)
1683 38848 : DO iiB = 1, nrow_local
1684 38256 : i_global = row_indices(iiB)
1685 3196680 : DO jjB = 1, ncol_local
1686 3157832 : j_global = col_indices(jjB)
1687 3196088 : IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
1688 : END DO
1689 : END DO
1690 :
1691 592 : CALL timestop(handle)
1692 :
1693 592 : END SUBROUTINE invert_mat
1694 :
1695 : ! **************************************************************************************************
1696 : !> \brief ...
1697 : !> \param fm_matrix_L ...
1698 : !> \param blacs_env_L ...
1699 : !> \param dimen_RI ...
1700 : !> \param para_env_L ...
1701 : !> \param name ...
1702 : !> \param fm_matrix_L_extern ...
1703 : ! **************************************************************************************************
1704 628 : SUBROUTINE create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, fm_matrix_L_extern)
1705 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_L
1706 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_L
1707 : INTEGER, INTENT(IN) :: dimen_RI
1708 : TYPE(mp_para_env_type), POINTER :: para_env_L
1709 : CHARACTER(LEN=*), INTENT(IN) :: name
1710 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_L_extern
1711 :
1712 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_matrix_L'
1713 :
1714 : INTEGER :: handle
1715 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1716 :
1717 628 : CALL timeset(routineN, handle)
1718 :
1719 : ! create the full matrix L defined in the L group
1720 628 : IF (PRESENT(fm_matrix_L_extern)) THEN
1721 36 : CALL cp_fm_create(fm_matrix_L, fm_matrix_L_extern%matrix_struct, name=name)
1722 : ELSE
1723 592 : NULLIFY (fm_struct)
1724 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_L, nrow_global=dimen_RI, &
1725 592 : ncol_global=dimen_RI, para_env=para_env_L)
1726 :
1727 592 : CALL cp_fm_create(fm_matrix_L, fm_struct, name=name)
1728 :
1729 592 : CALL cp_fm_struct_release(fm_struct)
1730 : END IF
1731 :
1732 628 : CALL cp_fm_set_all(matrix=fm_matrix_L, alpha=0.0_dp)
1733 :
1734 628 : CALL timestop(handle)
1735 :
1736 628 : END SUBROUTINE create_matrix_L
1737 :
1738 : ! **************************************************************************************************
1739 : !> \brief ...
1740 : !> \param fm_matrix_L ...
1741 : !> \param my_group_L_start ...
1742 : !> \param my_group_L_end ...
1743 : !> \param my_group_L_size ...
1744 : !> \param my_Lrows ...
1745 : ! **************************************************************************************************
1746 470 : SUBROUTINE grep_Lcols(fm_matrix_L, &
1747 : my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
1748 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1749 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, &
1750 : my_group_L_size
1751 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1752 : INTENT(OUT) :: my_Lrows
1753 :
1754 : CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_Lcols'
1755 :
1756 : INTEGER :: dimen_RI, handle, handle2, i_global, iiB, j_global, jjB, max_row_col_local, &
1757 : ncol_local, ncol_rec, nrow_local, nrow_rec, proc_receive_static, proc_send_static, &
1758 : proc_shift
1759 470 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info
1760 470 : INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, &
1761 470 : row_indices, row_indices_rec
1762 470 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: local_L, rec_L
1763 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1764 470 : POINTER :: local_L_internal
1765 : TYPE(mp_para_env_type), POINTER :: para_env
1766 :
1767 470 : CALL timeset(routineN, handle)
1768 :
1769 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1770 : nrow_local=nrow_local, &
1771 : ncol_local=ncol_local, &
1772 : row_indices=row_indices, &
1773 : col_indices=col_indices, &
1774 : nrow_global=dimen_RI, &
1775 : local_data=local_L_internal, &
1776 470 : para_env=para_env)
1777 :
1778 1880 : ALLOCATE (my_Lrows(dimen_RI, my_group_L_size))
1779 1505436 : my_Lrows = 0.0_dp
1780 :
1781 1880 : ALLOCATE (local_L(nrow_local, ncol_local))
1782 2741398 : local_L(:, :) = local_L_internal(1:nrow_local, 1:ncol_local)
1783 :
1784 470 : max_row_col_local = MAX(nrow_local, ncol_local)
1785 470 : CALL para_env%max(max_row_col_local)
1786 :
1787 1880 : ALLOCATE (local_col_row_info(0:max_row_col_local, 2))
1788 71614 : local_col_row_info = 0
1789 : ! 0,1 nrows
1790 470 : local_col_row_info(0, 1) = nrow_local
1791 34099 : local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
1792 : ! 0,2 ncols
1793 470 : local_col_row_info(0, 2) = ncol_local
1794 34984 : local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
1795 :
1796 1880 : ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
1797 :
1798 : ! accumulate data on my_Lrows starting from myself
1799 34984 : DO jjB = 1, ncol_local
1800 34514 : j_global = col_indices(jjB)
1801 34984 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
1802 1415978 : DO iiB = 1, nrow_local
1803 1397701 : i_global = row_indices(iiB)
1804 1415978 : my_Lrows(i_global, j_global - my_group_L_start + 1) = local_L(iiB, jjB)
1805 : END DO
1806 : END IF
1807 : END DO
1808 :
1809 470 : proc_send_static = MODULO(para_env%mepos + 1, para_env%num_pe)
1810 470 : proc_receive_static = MODULO(para_env%mepos - 1, para_env%num_pe)
1811 :
1812 470 : CALL timeset(routineN//"_comm", handle2)
1813 :
1814 494 : DO proc_shift = 1, para_env%num_pe - 1
1815 : ! first exchange information on the local data
1816 4200 : rec_col_row_info = 0
1817 24 : CALL para_env%sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static)
1818 24 : nrow_rec = rec_col_row_info(0, 1)
1819 24 : ncol_rec = rec_col_row_info(0, 2)
1820 :
1821 72 : ALLOCATE (row_indices_rec(nrow_rec))
1822 1061 : row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1823 :
1824 72 : ALLOCATE (col_indices_rec(ncol_rec))
1825 2064 : col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1826 :
1827 96 : ALLOCATE (rec_L(nrow_rec, ncol_rec))
1828 91052 : rec_L = 0.0_dp
1829 :
1830 : ! then send and receive the real data
1831 24 : CALL para_env%sendrecv(local_L, proc_send_static, rec_L, proc_receive_static)
1832 :
1833 : ! accumulate the received data on my_Lrows
1834 2064 : DO jjB = 1, ncol_rec
1835 2040 : j_global = col_indices_rec(jjB)
1836 2064 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
1837 91028 : DO iiB = 1, nrow_rec
1838 88988 : i_global = row_indices_rec(iiB)
1839 91028 : my_Lrows(i_global, j_global - my_group_L_start + 1) = rec_L(iiB, jjB)
1840 : END DO
1841 : END IF
1842 : END DO
1843 :
1844 4200 : local_col_row_info(:, :) = rec_col_row_info
1845 24 : CALL MOVE_ALLOC(rec_L, local_L)
1846 :
1847 24 : DEALLOCATE (col_indices_rec)
1848 494 : DEALLOCATE (row_indices_rec)
1849 : END DO
1850 470 : CALL timestop(handle2)
1851 :
1852 470 : DEALLOCATE (local_col_row_info)
1853 470 : DEALLOCATE (rec_col_row_info)
1854 470 : DEALLOCATE (local_L)
1855 :
1856 470 : CALL timestop(handle)
1857 :
1858 1410 : END SUBROUTINE grep_Lcols
1859 :
1860 : END MODULE mp2_ri_2c
|
