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_api, ONLY: &
33 : dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
34 : dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_p_type, dbcsr_release, &
35 : dbcsr_reserve_all_blocks, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
36 : dbcsr_type_no_symmetry, dbcsr_type_symmetric
37 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
38 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
39 : cp_dbcsr_dist2d_to_dist,&
40 : dbcsr_allocate_matrix_set,&
41 : dbcsr_deallocate_matrix_set
42 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_param
43 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
44 : cp_fm_triangular_invert
45 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
46 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
47 : cp_fm_power,&
48 : cp_fm_syevx
49 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
50 : cp_fm_struct_release,&
51 : cp_fm_struct_type
52 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
53 : cp_fm_create,&
54 : cp_fm_get_info,&
55 : cp_fm_release,&
56 : cp_fm_set_all,&
57 : cp_fm_to_fm,&
58 : cp_fm_type
59 : USE distribution_2d_types, ONLY: distribution_2d_type
60 : USE group_dist_types, ONLY: create_group_dist,&
61 : get_group_dist,&
62 : group_dist_d1_type,&
63 : release_group_dist
64 : USE input_constants, ONLY: do_eri_gpw,&
65 : do_eri_mme,&
66 : do_eri_os,&
67 : do_potential_coulomb,&
68 : do_potential_id,&
69 : do_potential_long,&
70 : do_potential_truncated
71 : USE kinds, ONLY: dp
72 : USE kpoint_coulomb_2c, ONLY: build_2c_coulomb_matrix_kp
73 : USE kpoint_methods, ONLY: kpoint_init_cell_index,&
74 : 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, trunc_coulomb_for_exchange, 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 618 : 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 618 : 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 618 : num_small_eigen, potential, sab_orb_sub, do_im_time=do_im_time)
209 :
210 618 : 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 618 : dimen_RI_red = dimen_RI
217 :
218 618 : 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 448 : 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 trunc_coulomb_for_exchange(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 618 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
328 :
329 1102 : IF (.NOT. do_im_time) THEN
330 484 : IF (unit_nr > 0) THEN
331 242 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Cholesky decomposition group size:", para_env_L%num_pe
332 242 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Number of groups for auxiliary basis functions", ngroup
333 242 : 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 242 : CALL m_flush(unit_nr)
340 : END IF
341 :
342 : ! replicate the necessary row of the L^{-1} matrix on each proc
343 484 : CALL grep_Lcols(fm_matrix_L_work, my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
344 496 : 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 618 : CALL cp_fm_release(fm_matrix_L_work)
350 618 : 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 618 : CALL mp_para_env_release(para_env_L)
353 :
354 618 : CALL timestop(handle)
355 :
356 618 : 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 756 : 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 756 : 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 712 : CALL cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion=do_inversion)
398 :
399 712 : IF (do_inversion) CALL invert_mat(fm_matrix_L)
400 : END IF
401 :
402 756 : 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 : !> \param do_build_cell_index ...
552 : ! **************************************************************************************************
553 468 : SUBROUTINE RI_2c_integral_mat(qs_env, fm_matrix_Minv_L_kpoints, fm_matrix_L, dimen_RI, ri_metric, &
554 : do_kpoints, kpoints, put_mat_KS_env, regularization_RI, ikp_ext, &
555 : do_build_cell_index)
556 :
557 : TYPE(qs_environment_type), POINTER :: qs_env
558 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_Minv_L_kpoints
559 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
560 : INTEGER, INTENT(IN) :: dimen_RI
561 : TYPE(libint_potential_type) :: ri_metric
562 : LOGICAL, INTENT(IN) :: do_kpoints
563 : TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
564 : LOGICAL, OPTIONAL :: put_mat_KS_env
565 : REAL(KIND=dp), OPTIONAL :: regularization_RI
566 : INTEGER, OPTIONAL :: ikp_ext
567 : LOGICAL, OPTIONAL :: do_build_cell_index
568 :
569 : CHARACTER(LEN=*), PARAMETER :: routineN = 'RI_2c_integral_mat'
570 :
571 : INTEGER :: handle, i_real_imag, i_size, ikp, &
572 : ikp_for_xkp, img, n_real_imag, natom, &
573 : nimg, nkind, nkp
574 468 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI
575 468 : INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
576 468 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
577 : LOGICAL :: my_do_build_cell_index, my_put_mat_KS_env
578 : REAL(KIND=dp) :: my_regularization_RI
579 468 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xkp
580 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
581 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
582 : TYPE(cp_fm_type) :: fm_matrix_S_global
583 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
584 468 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_RI_aux_transl
585 468 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: mat_2c
586 : TYPE(dbcsr_type), POINTER :: cmatrix, matrix_s_RI_aux_desymm, &
587 : rmatrix, tmpmat
588 : TYPE(dft_control_type), POINTER :: dft_control
589 : TYPE(distribution_2d_type), POINTER :: dist_2d
590 468 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_RI
591 : TYPE(mp_para_env_type), POINTER :: para_env
592 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
593 468 : POINTER :: sab_RI
594 468 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
595 468 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
596 :
597 468 : CALL timeset(routineN, handle)
598 :
599 468 : NULLIFY (sab_RI, matrix_s_RI_aux_transl, dist_2d)
600 :
601 468 : IF (PRESENT(regularization_RI)) THEN
602 438 : my_regularization_RI = regularization_RI
603 : ELSE
604 30 : my_regularization_RI = 0.0_dp
605 : END IF
606 :
607 468 : IF (PRESENT(put_mat_KS_env)) THEN
608 178 : my_put_mat_KS_env = put_mat_KS_env
609 : ELSE
610 : my_put_mat_KS_env = .FALSE.
611 : END IF
612 :
613 468 : IF (PRESENT(do_build_cell_index)) THEN
614 260 : my_do_build_cell_index = do_build_cell_index
615 : ELSE
616 : my_do_build_cell_index = .FALSE.
617 : END IF
618 :
619 : CALL get_qs_env(qs_env=qs_env, &
620 : para_env=para_env, &
621 : blacs_env=blacs_env, &
622 : nkind=nkind, &
623 : distribution_2d=dist_2d, &
624 : qs_kind_set=qs_kind_set, &
625 : particle_set=particle_set, &
626 : dft_control=dft_control, &
627 468 : natom=natom)
628 :
629 1404 : ALLOCATE (sizes_RI(natom))
630 2328 : ALLOCATE (basis_set_RI(nkind))
631 :
632 468 : CALL basis_set_list_setup(basis_set_RI, 'RI_AUX', qs_kind_set)
633 468 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_RI)
634 :
635 : CALL build_2c_neighbor_lists(sab_RI, basis_set_RI, basis_set_RI, ri_metric, "RPA_2c_nl_RI", qs_env, &
636 468 : sym_ij=.TRUE., dist_2d=dist_2d)
637 :
638 468 : CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
639 1404 : ALLOCATE (row_bsize(SIZE(sizes_RI)))
640 1404 : ALLOCATE (col_bsize(SIZE(sizes_RI)))
641 1716 : row_bsize(:) = sizes_RI
642 1716 : col_bsize(:) = sizes_RI
643 :
644 468 : IF (do_kpoints) THEN
645 282 : CPASSERT(PRESENT(kpoints))
646 282 : IF (my_do_build_cell_index) THEN
647 10 : CALL kpoint_init_cell_index(kpoints, sab_RI, para_env, dft_control)
648 : END IF
649 : CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, &
650 282 : cell_to_index=cell_to_index)
651 282 : n_real_imag = 2
652 282 : nimg = dft_control%nimages
653 : ELSE
654 186 : nkp = 1
655 186 : n_real_imag = 1
656 186 : nimg = 1
657 : END IF
658 :
659 2990 : ALLOCATE (mat_2c(nimg))
660 : CALL dbcsr_create(mat_2c(1), "(RI|RI)", dbcsr_dist, dbcsr_type_symmetric, &
661 468 : row_bsize, col_bsize)
662 468 : DEALLOCATE (row_bsize, col_bsize)
663 :
664 1586 : DO img = 2, nimg
665 1586 : CALL dbcsr_create(mat_2c(img), template=mat_2c(1))
666 : END DO
667 :
668 : CALL build_2c_integrals(mat_2c, 0.0_dp, qs_env, sab_RI, basis_set_RI, basis_set_RI, &
669 : ri_metric, do_kpoints=do_kpoints, ext_kpoints=kpoints, &
670 468 : regularization_RI=my_regularization_RI)
671 :
672 468 : CALL dbcsr_distribution_release(dbcsr_dist)
673 468 : DEALLOCATE (basis_set_RI)
674 :
675 468 : IF (my_put_mat_KS_env) THEN
676 134 : CALL get_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
677 134 : IF (ASSOCIATED(matrix_s_RI_aux_transl)) CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
678 : END IF
679 468 : CALL dbcsr_allocate_matrix_set(matrix_s_RI_aux_transl, 1, nimg)
680 2054 : DO img = 1, nimg
681 1586 : ALLOCATE (matrix_s_RI_aux_transl(1, img)%matrix)
682 1586 : CALL dbcsr_copy(matrix_s_RI_aux_transl(1, img)%matrix, mat_2c(img))
683 2054 : CALL dbcsr_release(mat_2c(img))
684 : END DO
685 :
686 468 : IF (my_put_mat_KS_env) THEN
687 134 : CALL set_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
688 : END IF
689 :
690 468 : IF (PRESENT(ikp_ext)) nkp = 1
691 :
692 4280 : ALLOCATE (fm_matrix_Minv_L_kpoints(nkp, n_real_imag))
693 1390 : DO i_size = 1, nkp
694 3048 : DO i_real_imag = 1, n_real_imag
695 1658 : CALL cp_fm_create(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), fm_matrix_L%matrix_struct)
696 2580 : CALL cp_fm_set_all(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), 0.0_dp)
697 : END DO
698 : END DO
699 :
700 468 : NULLIFY (fm_struct)
701 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
702 468 : ncol_global=dimen_RI, para_env=para_env)
703 :
704 468 : CALL cp_fm_create(fm_matrix_S_global, fm_struct)
705 468 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
706 :
707 468 : IF (do_kpoints) THEN
708 :
709 282 : ALLOCATE (rmatrix, cmatrix, tmpmat)
710 : CALL dbcsr_create(rmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
711 282 : matrix_type=dbcsr_type_symmetric)
712 : CALL dbcsr_create(cmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
713 282 : matrix_type=dbcsr_type_antisymmetric)
714 : CALL dbcsr_create(tmpmat, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
715 282 : matrix_type=dbcsr_type_no_symmetry)
716 282 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_RI)
717 282 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_RI)
718 :
719 1018 : DO ikp = 1, nkp
720 :
721 : ! s matrix is not spin dependent, double the work
722 736 : CALL dbcsr_set(rmatrix, 0.0_dp)
723 736 : CALL dbcsr_set(cmatrix, 0.0_dp)
724 :
725 736 : IF (PRESENT(ikp_ext)) THEN
726 260 : ikp_for_xkp = ikp_ext
727 : ELSE
728 : ikp_for_xkp = ikp
729 : END IF
730 :
731 : CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s_RI_aux_transl, ispin=1, &
732 736 : xkp=xkp(1:3, ikp_for_xkp), cell_to_index=cell_to_index, sab_nl=sab_RI)
733 :
734 736 : CALL dbcsr_set(tmpmat, 0.0_dp)
735 736 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
736 :
737 736 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
738 736 : CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
739 736 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 1), para_env)
740 :
741 736 : CALL dbcsr_set(tmpmat, 0.0_dp)
742 736 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
743 :
744 736 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
745 736 : CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
746 1018 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 2), para_env)
747 :
748 : END DO
749 :
750 282 : CALL dbcsr_deallocate_matrix(rmatrix)
751 282 : CALL dbcsr_deallocate_matrix(cmatrix)
752 282 : CALL dbcsr_deallocate_matrix(tmpmat)
753 :
754 : ELSE
755 :
756 : NULLIFY (matrix_s_RI_aux_desymm)
757 186 : ALLOCATE (matrix_s_RI_aux_desymm)
758 : CALL dbcsr_create(matrix_s_RI_aux_desymm, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
759 186 : name='S_RI non_symm', matrix_type=dbcsr_type_no_symmetry)
760 :
761 186 : CALL dbcsr_desymmetrize(matrix_s_RI_aux_transl(1, 1)%matrix, matrix_s_RI_aux_desymm)
762 :
763 186 : CALL copy_dbcsr_to_fm(matrix_s_RI_aux_desymm, fm_matrix_S_global)
764 :
765 186 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(1, 1), para_env)
766 :
767 186 : CALL dbcsr_deallocate_matrix(matrix_s_RI_aux_desymm)
768 :
769 : END IF
770 :
771 468 : CALL release_neighbor_list_sets(sab_RI)
772 :
773 468 : CALL cp_fm_struct_release(fm_struct)
774 :
775 468 : CALL cp_fm_release(fm_matrix_S_global)
776 :
777 468 : IF (.NOT. my_put_mat_KS_env) THEN
778 334 : CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
779 : END IF
780 :
781 468 : CALL timestop(handle)
782 :
783 1872 : END SUBROUTINE RI_2c_integral_mat
784 :
785 : ! **************************************************************************************************
786 : !> \brief ...
787 : !> \param qs_env ...
788 : !> \param eri_method ...
789 : !> \param eri_param ...
790 : !> \param para_env ...
791 : !> \param para_env_sub ...
792 : !> \param para_env_L ...
793 : !> \param mp2_memory ...
794 : !> \param fm_matrix_L ...
795 : !> \param ngroup ...
796 : !> \param color_sub ...
797 : !> \param dimen_RI ...
798 : !> \param my_group_L_size ...
799 : !> \param my_group_L_start ...
800 : !> \param my_group_L_end ...
801 : !> \param gd_array ...
802 : !> \param calc_PQ_cond_num ...
803 : !> \param cond_num ...
804 : !> \param num_small_eigen ...
805 : !> \param potential ...
806 : !> \param sab_orb_sub ...
807 : !> \param do_im_time ...
808 : !> \param fm_matrix_L_extern ...
809 : ! **************************************************************************************************
810 654 : SUBROUTINE compute_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, para_env_L, mp2_memory, &
811 : fm_matrix_L, ngroup, color_sub, dimen_RI, &
812 : my_group_L_size, my_group_L_start, my_group_L_end, &
813 : gd_array, calc_PQ_cond_num, cond_num, num_small_eigen, potential, &
814 : sab_orb_sub, do_im_time, fm_matrix_L_extern)
815 :
816 : TYPE(qs_environment_type), POINTER :: qs_env
817 : INTEGER, INTENT(IN) :: eri_method
818 : TYPE(cp_eri_mme_param), POINTER :: eri_param
819 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
820 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env_sub
821 : TYPE(mp_para_env_type), INTENT(OUT), POINTER :: para_env_L
822 : REAL(KIND=dp), INTENT(IN) :: mp2_memory
823 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_L
824 : INTEGER, INTENT(IN) :: ngroup, color_sub
825 : INTEGER, INTENT(OUT) :: dimen_RI, my_group_L_size, &
826 : my_group_L_start, my_group_L_end
827 : TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
828 : LOGICAL, INTENT(IN) :: calc_PQ_cond_num
829 : REAL(KIND=dp), INTENT(OUT) :: cond_num
830 : INTEGER, INTENT(OUT) :: num_small_eigen
831 : TYPE(libint_potential_type), INTENT(IN) :: potential
832 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
833 : POINTER :: sab_orb_sub
834 : LOGICAL, INTENT(IN), OPTIONAL :: do_im_time
835 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_L_extern
836 :
837 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_2c_integrals'
838 :
839 : INTEGER :: best_group_size, color_L, group_size, handle, handle2, i_global, iatom, iiB, &
840 : ikind, iproc, j_global, jjB, natom, ncol_local, nkind, nrow_local, nsgf, potential_type, &
841 : proc_receive, proc_receive_static, proc_send_static, proc_shift, rec_L_end, rec_L_size, &
842 : rec_L_start, strat_group_size
843 654 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
844 654 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
845 : LOGICAL :: my_do_im_time
846 : REAL(KIND=dp) :: min_mem_for_QK
847 654 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egen_L
848 654 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: L_external_col, L_local_col
849 654 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
850 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_L
851 : TYPE(cp_fm_type) :: fm_matrix_L_diag
852 654 : TYPE(group_dist_d1_type) :: gd_sub_array
853 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
854 654 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
855 654 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
856 :
857 654 : CALL timeset(routineN, handle)
858 :
859 654 : my_do_im_time = .FALSE.
860 654 : IF (PRESENT(do_im_time)) THEN
861 618 : my_do_im_time = do_im_time
862 : END IF
863 :
864 : ! get stuff
865 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
866 654 : particle_set=particle_set)
867 :
868 654 : nkind = SIZE(qs_kind_set)
869 654 : natom = SIZE(particle_set)
870 :
871 654 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
872 :
873 1858 : DO ikind = 1, nkind
874 1204 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
875 1858 : CPASSERT(ASSOCIATED(basis_set_a))
876 : END DO
877 :
878 654 : dimen_RI = 0
879 2580 : DO iatom = 1, natom
880 1926 : ikind = kind_of(iatom)
881 1926 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
882 2580 : dimen_RI = dimen_RI + nsgf
883 : END DO
884 :
885 : ! check that very small systems are not running on too many processes
886 654 : IF (dimen_RI < ngroup) THEN
887 : CALL cp_abort(__LOCATION__, "Product of block size and number "// &
888 0 : "of RI functions should not exceed total number of processes")
889 : END IF
890 :
891 654 : CALL create_group_dist(gd_array, ngroup, dimen_RI)
892 :
893 654 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
894 :
895 654 : CALL timeset(routineN//"_loop_lm", handle2)
896 :
897 2616 : ALLOCATE (L_local_col(dimen_RI, my_group_L_size))
898 2256612 : L_local_col = 0.0_dp
899 :
900 654 : potential_type = potential%potential_type
901 :
902 : IF ((eri_method .EQ. do_eri_mme .OR. eri_method .EQ. do_eri_os) &
903 654 : .AND. potential_type .EQ. do_potential_coulomb .OR. &
904 : (eri_method .EQ. do_eri_mme .AND. potential_type .EQ. do_potential_long)) THEN
905 :
906 : CALL mp2_eri_2c_integrate(eri_param, potential, para_env_sub, qs_env, &
907 : basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
908 : hab=L_local_col, first_b=my_group_L_start, last_b=my_group_L_end, &
909 280 : eri_method=eri_method)
910 :
911 : ELSEIF (eri_method .EQ. do_eri_gpw .OR. &
912 : (potential_type == do_potential_long .AND. qs_env%mp2_env%eri_method == do_eri_os) &
913 374 : .OR. (potential_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)) THEN
914 :
915 : CALL mp2_eri_2c_integrate_gpw(qs_env, para_env_sub, my_group_L_start, my_group_L_end, &
916 374 : natom, potential, sab_orb_sub, L_local_col, kind_of)
917 :
918 : ELSE
919 0 : CPABORT("unknown ERI method")
920 : END IF
921 :
922 654 : CALL timestop(handle2)
923 :
924 : ! split the total number of proc in a subgroup of the size of ~1/10 of the
925 : ! total num of proc
926 654 : best_group_size = para_env%num_pe
927 :
928 654 : strat_group_size = MAX(1, para_env%num_pe/10)
929 :
930 654 : min_mem_for_QK = REAL(dimen_RI, KIND=dp)*dimen_RI*3.0_dp*8.0_dp/1024_dp/1024_dp
931 :
932 654 : group_size = strat_group_size - 1
933 688 : DO iproc = strat_group_size, para_env%num_pe
934 688 : group_size = group_size + 1
935 : ! check that group_size is a multiple of sub_group_size and a divisor of
936 : ! the total num of proc
937 688 : IF (MOD(para_env%num_pe, group_size) /= 0 .OR. MOD(group_size, para_env_sub%num_pe) /= 0) CYCLE
938 :
939 : ! check for memory
940 654 : IF (REAL(group_size, KIND=dp)*mp2_memory < min_mem_for_QK) CYCLE
941 :
942 : best_group_size = group_size
943 34 : EXIT
944 : END DO
945 :
946 654 : IF (my_do_im_time) THEN
947 : ! para_env_L is the para_env for im. time to avoid bug
948 134 : best_group_size = para_env%num_pe
949 : END IF
950 :
951 : ! create the L group
952 654 : color_L = para_env%mepos/best_group_size
953 654 : ALLOCATE (para_env_L)
954 654 : CALL para_env_L%from_split(para_env, color_L)
955 :
956 : ! create the blacs_L
957 654 : NULLIFY (blacs_env_L)
958 654 : CALL cp_blacs_env_create(blacs_env=blacs_env_L, para_env=para_env_L)
959 :
960 : ! create the full matrix L defined in the L group
961 654 : CALL create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, "fm_matrix_L", fm_matrix_L_extern)
962 :
963 654 : IF (my_do_im_time .AND. para_env%num_pe > 1) THEN
964 :
965 : CALL fill_fm_L_from_L_loc_non_blocking(fm_matrix_L, L_local_col, para_env, &
966 : my_group_L_start, my_group_L_end, &
967 134 : dimen_RI)
968 :
969 : ELSE
970 520 : BLOCK
971 : TYPE(mp_comm_type) :: comm_exchange
972 :
973 520 : CALL comm_exchange%from_split(para_env_L, para_env_sub%mepos)
974 :
975 : CALL create_group_dist(gd_sub_array, my_group_L_start, &
976 520 : my_group_L_end, my_group_L_size, comm_exchange)
977 :
978 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
979 : nrow_local=nrow_local, &
980 : ncol_local=ncol_local, &
981 : row_indices=row_indices, &
982 520 : col_indices=col_indices)
983 :
984 39392 : DO jjB = 1, ncol_local
985 38872 : j_global = col_indices(jjB)
986 39392 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
987 1602166 : DO iiB = 1, nrow_local
988 1581610 : i_global = row_indices(iiB)
989 1602166 : fm_matrix_L%local_data(iiB, jjB) = L_local_col(i_global, j_global - my_group_L_start + 1)
990 : END DO
991 : END IF
992 : END DO
993 :
994 520 : proc_send_static = MODULO(comm_exchange%mepos + 1, comm_exchange%num_pe)
995 520 : proc_receive_static = MODULO(comm_exchange%mepos - 1, comm_exchange%num_pe)
996 :
997 520 : DO proc_shift = 1, comm_exchange%num_pe - 1
998 0 : proc_receive = MODULO(comm_exchange%mepos - proc_shift, comm_exchange%num_pe)
999 :
1000 0 : CALL get_group_dist(gd_sub_array, proc_receive, rec_L_start, rec_L_end, rec_L_size)
1001 :
1002 0 : ALLOCATE (L_external_col(dimen_RI, rec_L_size))
1003 0 : L_external_col = 0.0_dp
1004 :
1005 0 : CALL comm_exchange%sendrecv(L_local_col, proc_send_static, L_external_col, proc_receive_static)
1006 :
1007 0 : DO jjB = 1, ncol_local
1008 0 : j_global = col_indices(jjB)
1009 0 : IF (j_global >= rec_L_start .AND. j_global <= rec_L_end) THEN
1010 0 : DO iiB = 1, nrow_local
1011 0 : i_global = row_indices(iiB)
1012 0 : fm_matrix_L%local_data(iiB, jjB) = L_external_col(i_global, j_global - rec_L_start + 1)
1013 : END DO
1014 : END IF
1015 : END DO
1016 :
1017 520 : CALL MOVE_ALLOC(L_external_col, L_local_col)
1018 : END DO
1019 :
1020 520 : CALL release_group_dist(gd_sub_array)
1021 1040 : CALL comm_exchange%free()
1022 : END BLOCK
1023 : END IF
1024 :
1025 654 : DEALLOCATE (L_local_col)
1026 :
1027 : ! create the new group for the sum of the local data over all processes
1028 : BLOCK
1029 : TYPE(mp_comm_type) :: comm_exchange
1030 654 : comm_exchange = fm_matrix_L%matrix_struct%context%interconnect(para_env)
1031 654 : CALL comm_exchange%sum(fm_matrix_L%local_data)
1032 1308 : CALL comm_exchange%free()
1033 : END BLOCK
1034 :
1035 654 : cond_num = 1.0_dp
1036 654 : num_small_eigen = 0
1037 654 : IF (calc_PQ_cond_num) THEN
1038 : ! calculate the condition number of the (P|Q) matrix
1039 : ! create a copy of the matrix
1040 :
1041 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)
1042 :
1043 0 : CALL cp_fm_to_fm(source=fm_matrix_L, destination=fm_matrix_L_diag)
1044 :
1045 0 : ALLOCATE (egen_L(dimen_RI))
1046 :
1047 0 : egen_L = 0.0_dp
1048 0 : CALL cp_fm_syevx(matrix=fm_matrix_L_diag, eigenvalues=egen_L)
1049 :
1050 0 : num_small_eigen = 0
1051 0 : DO iiB = 1, dimen_RI
1052 0 : IF (ABS(egen_L(iiB)) < 0.001_dp) num_small_eigen = num_small_eigen + 1
1053 : END DO
1054 :
1055 0 : cond_num = MAXVAL(ABS(egen_L))/MINVAL(ABS(egen_L))
1056 :
1057 0 : CALL cp_fm_release(fm_matrix_L_diag)
1058 :
1059 0 : DEALLOCATE (egen_L)
1060 : END IF
1061 :
1062 : ! release blacs_env
1063 654 : CALL cp_blacs_env_release(blacs_env_L)
1064 :
1065 654 : CALL timestop(handle)
1066 :
1067 3270 : END SUBROUTINE compute_2c_integrals
1068 :
1069 : ! **************************************************************************************************
1070 : !> \brief ...
1071 : !> \param matrix ...
1072 : !> \param num_small_evals ...
1073 : !> \param cond_num ...
1074 : !> \param do_inversion ...
1075 : !> \param para_env ...
1076 : ! **************************************************************************************************
1077 44 : SUBROUTINE matrix_root_with_svd(matrix, num_small_evals, cond_num, do_inversion, para_env)
1078 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1079 : INTEGER, INTENT(OUT) :: num_small_evals
1080 : REAL(KIND=dp), INTENT(OUT) :: cond_num
1081 : LOGICAL, INTENT(IN) :: do_inversion
1082 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1083 :
1084 : CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_root_with_svd'
1085 :
1086 : INTEGER :: group_size_L, handle, ii, needed_evals, &
1087 : nrow, pos_max(1)
1088 44 : INTEGER, ALLOCATABLE, DIMENSION(:) :: num_eval
1089 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
1090 : TYPE(cp_fm_type) :: evecs
1091 : TYPE(mp_comm_type) :: comm_exchange
1092 :
1093 44 : CALL timeset(routineN, handle)
1094 :
1095 : ! Create arrays carrying eigenvalues and eigenvectors later
1096 44 : CALL cp_fm_get_info(matrix=matrix, nrow_global=nrow)
1097 :
1098 132 : ALLOCATE (evals(nrow))
1099 3724 : evals = 0
1100 :
1101 44 : CALL cp_fm_create(evecs, matrix%matrix_struct)
1102 :
1103 : ! Perform the EVD (=SVD of a positive semidefinite matrix)
1104 44 : CALL choose_eigv_solver(matrix, evecs, evals)
1105 :
1106 : ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
1107 44 : num_small_evals = 0
1108 202 : DO ii = 1, nrow
1109 202 : IF (evals(ii) > 0.0_dp) THEN
1110 44 : num_small_evals = ii - 1
1111 44 : EXIT
1112 : END IF
1113 : END DO
1114 44 : needed_evals = nrow - num_small_evals
1115 :
1116 : ! Get the condition number w.r.t. considered eigenvalues
1117 44 : cond_num = evals(nrow)/evals(num_small_evals + 1)
1118 :
1119 : ! Determine the eigenvalues of the request matrix root or its inverse
1120 202 : evals(1:num_small_evals) = 0.0_dp
1121 44 : IF (do_inversion) THEN
1122 1008 : evals(num_small_evals + 1:nrow) = 1.0_dp/SQRT(evals(num_small_evals + 1:nrow))
1123 : ELSE
1124 2558 : evals(num_small_evals + 1:nrow) = SQRT(evals(num_small_evals + 1:nrow))
1125 : END IF
1126 :
1127 44 : CALL cp_fm_column_scale(evecs, evals)
1128 :
1129 : ! As it turns out, the results in the subgroups differ.
1130 : ! Thus, we have to equilibrate the results.
1131 : ! Step 1: Get a communicator connecting processes with same id within subgroups
1132 : ! use that group_size_L is divisible by the total number of processes (see above)
1133 44 : group_size_L = para_env%num_pe/matrix%matrix_struct%para_env%num_pe
1134 44 : comm_exchange = matrix%matrix_struct%context%interconnect(para_env)
1135 :
1136 132 : ALLOCATE (num_eval(0:group_size_L - 1))
1137 44 : CALL comm_exchange%allgather(num_small_evals, num_eval)
1138 :
1139 118 : num_small_evals = MINVAL(num_eval)
1140 :
1141 118 : IF (num_small_evals /= MAXVAL(num_eval)) THEN
1142 : ! Step 2: Get position of maximum value
1143 0 : pos_max = MAXLOC(num_eval)
1144 0 : num_small_evals = num_eval(pos_max(1))
1145 0 : needed_evals = nrow - num_small_evals
1146 :
1147 : ! Step 3: Broadcast your local data to all other processes
1148 0 : CALL comm_exchange%bcast(evecs%local_data, pos_max(1))
1149 0 : CALL comm_exchange%bcast(cond_num, pos_max(1))
1150 : END IF
1151 :
1152 44 : DEALLOCATE (num_eval)
1153 :
1154 44 : CALL comm_exchange%free()
1155 :
1156 44 : CALL reset_size_matrix(matrix, needed_evals, matrix%matrix_struct)
1157 :
1158 : ! Copy the needed eigenvectors
1159 44 : CALL cp_fm_to_fm(evecs, matrix, needed_evals, num_small_evals + 1)
1160 :
1161 44 : CALL cp_fm_release(evecs)
1162 :
1163 44 : CALL timestop(handle)
1164 :
1165 132 : END SUBROUTINE matrix_root_with_svd
1166 :
1167 : ! **************************************************************************************************
1168 : !> \brief ...
1169 : !> \param matrix ...
1170 : !> \param new_size ...
1171 : !> \param fm_struct_template ...
1172 : ! **************************************************************************************************
1173 152 : SUBROUTINE reset_size_matrix(matrix, new_size, fm_struct_template)
1174 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix
1175 : INTEGER, INTENT(IN) :: new_size
1176 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_template
1177 :
1178 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reset_size_matrix'
1179 :
1180 : INTEGER :: handle
1181 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1182 :
1183 76 : CALL timeset(routineN, handle)
1184 :
1185 : ! Choose the block sizes as large as in the template for the later steps
1186 76 : NULLIFY (fm_struct)
1187 76 : CALL cp_fm_struct_create(fm_struct, ncol_global=new_size, template_fmstruct=fm_struct_template, force_block=.TRUE.)
1188 :
1189 76 : CALL cp_fm_release(matrix)
1190 :
1191 76 : CALL cp_fm_create(matrix, fm_struct)
1192 76 : CALL cp_fm_set_all(matrix, 0.0_dp)
1193 :
1194 76 : CALL cp_fm_struct_release(fm_struct)
1195 :
1196 76 : CALL timestop(handle)
1197 :
1198 76 : END SUBROUTINE reset_size_matrix
1199 :
1200 : ! **************************************************************************************************
1201 : !> \brief ...
1202 : !> \param fm_matrix_L ...
1203 : !> \param L_local_col ...
1204 : !> \param para_env ...
1205 : !> \param my_group_L_start ...
1206 : !> \param my_group_L_end ...
1207 : !> \param dimen_RI ...
1208 : ! **************************************************************************************************
1209 134 : SUBROUTINE fill_fm_L_from_L_loc_non_blocking(fm_matrix_L, L_local_col, para_env, my_group_L_start, &
1210 : my_group_L_end, dimen_RI)
1211 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1212 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1213 : INTENT(IN) :: L_local_col
1214 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1215 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, &
1216 : dimen_RI
1217 :
1218 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fill_fm_L_from_L_loc_non_blocking'
1219 :
1220 : INTEGER :: dummy_proc, handle, handle2, i_entry_rec, i_row, i_row_global, iproc, j_col, &
1221 : j_col_global, LLL, MMM, ncol_local, nrow_local, proc_send, send_pcol, send_prow
1222 134 : INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
1223 : num_entries_send
1224 134 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1225 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1226 134 : DIMENSION(:) :: buffer_rec, buffer_send
1227 134 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
1228 :
1229 134 : CALL timeset(routineN, handle)
1230 :
1231 134 : CALL timeset(routineN//"_1", handle2)
1232 :
1233 : ! get info for the dest
1234 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1235 : nrow_local=nrow_local, &
1236 : ncol_local=ncol_local, &
1237 : row_indices=row_indices, &
1238 134 : col_indices=col_indices)
1239 :
1240 402 : ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
1241 268 : ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
1242 :
1243 402 : num_entries_rec(:) = 0
1244 402 : num_entries_send(:) = 0
1245 :
1246 134 : dummy_proc = 0
1247 :
1248 : ! get the process, where the elements from L_local_col have to be sent
1249 11380 : DO LLL = 1, dimen_RI
1250 :
1251 11246 : send_prow = fm_matrix_L%matrix_struct%g2p_row(LLL)
1252 :
1253 561929 : DO MMM = my_group_L_start, my_group_L_end
1254 :
1255 550549 : send_pcol = fm_matrix_L%matrix_struct%g2p_col(MMM)
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 536 : 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 268 : 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 11246 : send_prow = fm_matrix_L%matrix_struct%g2p_row(LLL)
1332 :
1333 561929 : DO MMM = my_group_L_start, my_group_L_end
1334 :
1335 550549 : send_pcol = fm_matrix_L%matrix_struct%g2p_col(MMM)
1336 :
1337 550549 : proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1338 :
1339 550549 : entry_counter(proc_send) = entry_counter(proc_send) + 1
1340 :
1341 : buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
1342 550549 : L_local_col(LLL, MMM - my_group_L_start + 1)
1343 :
1344 550549 : buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = LLL
1345 561795 : buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = MMM
1346 :
1347 : END DO
1348 :
1349 : END DO
1350 :
1351 2010 : ALLOCATE (req_array(1:para_env%num_pe, 4))
1352 :
1353 134 : CALL timestop(handle2)
1354 :
1355 134 : CALL timeset(routineN//"_8", handle2)
1356 :
1357 : ! communicate the buffer
1358 : CALL communicate_buffer(para_env, num_entries_rec, num_entries_send, buffer_rec, &
1359 134 : buffer_send, req_array)
1360 :
1361 534373 : fm_matrix_L%local_data = 0.0_dp
1362 :
1363 134 : CALL timestop(handle2)
1364 :
1365 134 : CALL timeset(routineN//"_9", handle2)
1366 :
1367 : ! fill fm_matrix_L with the entries from buffer_rec
1368 402 : DO iproc = 0, para_env%num_pe - 1
1369 :
1370 550951 : DO i_entry_rec = 1, num_entries_rec(iproc)
1371 :
1372 31219430 : DO i_row = 1, nrow_local
1373 :
1374 30668613 : i_row_global = row_indices(i_row)
1375 :
1376 4667031343 : DO j_col = 1, ncol_local
1377 :
1378 4635812181 : j_col_global = col_indices(j_col)
1379 :
1380 4635812181 : IF (i_row_global == buffer_rec(iproc)%indx(i_entry_rec, 1) .AND. &
1381 30668613 : j_col_global == buffer_rec(iproc)%indx(i_entry_rec, 2)) THEN
1382 :
1383 550549 : fm_matrix_L%local_data(i_row, j_col) = buffer_rec(iproc)%msg(i_entry_rec)
1384 :
1385 : END IF
1386 :
1387 : END DO
1388 :
1389 : END DO
1390 :
1391 : END DO
1392 :
1393 : END DO
1394 :
1395 134 : CALL timestop(handle2)
1396 :
1397 134 : CALL timeset(routineN//"_10", handle2)
1398 :
1399 402 : DO iproc = 0, para_env%num_pe - 1
1400 268 : DEALLOCATE (buffer_rec(iproc)%msg)
1401 268 : DEALLOCATE (buffer_rec(iproc)%indx)
1402 268 : DEALLOCATE (buffer_send(iproc)%msg)
1403 402 : DEALLOCATE (buffer_send(iproc)%indx)
1404 : END DO
1405 :
1406 670 : DEALLOCATE (buffer_rec, buffer_send)
1407 134 : DEALLOCATE (req_array)
1408 134 : DEALLOCATE (entry_counter)
1409 134 : DEALLOCATE (num_entries_rec, num_entries_send)
1410 :
1411 134 : CALL timestop(handle2)
1412 :
1413 134 : CALL timestop(handle)
1414 :
1415 1340 : END SUBROUTINE fill_fm_L_from_L_loc_non_blocking
1416 :
1417 : ! **************************************************************************************************
1418 : !> \brief ...
1419 : !> \param fm_matrix_L ...
1420 : !> \param dimen_RI ...
1421 : !> \param do_inversion ...
1422 : ! **************************************************************************************************
1423 1468 : SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion)
1424 :
1425 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1426 : INTEGER, INTENT(IN) :: dimen_RI
1427 : LOGICAL, INTENT(IN) :: do_inversion
1428 :
1429 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholesky_decomp'
1430 :
1431 : INTEGER :: handle, i_global, iiB, info_chol, &
1432 : j_global, jjB, ncol_local, nrow_local
1433 734 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1434 :
1435 734 : CALL timeset(routineN, handle)
1436 :
1437 : ! do cholesky decomposition
1438 734 : CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=dimen_RI, info_out=info_chol)
1439 734 : CPASSERT(info_chol == 0)
1440 :
1441 734 : IF (.NOT. do_inversion) THEN
1442 : ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1443 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1444 : nrow_local=nrow_local, &
1445 : ncol_local=ncol_local, &
1446 : row_indices=row_indices, &
1447 116 : col_indices=col_indices)
1448 5803 : DO iiB = 1, nrow_local
1449 5687 : i_global = row_indices(iiB)
1450 538330 : DO jjB = 1, ncol_local
1451 532527 : j_global = col_indices(jjB)
1452 538214 : IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
1453 : END DO
1454 : END DO
1455 :
1456 : END IF
1457 :
1458 734 : CALL timestop(handle)
1459 :
1460 734 : END SUBROUTINE cholesky_decomp
1461 :
1462 : ! **************************************************************************************************
1463 : !> \brief ...
1464 : !> \param fm_matrix_Minv_L_kpoints ...
1465 : !> \param fm_matrix_L_kpoints ...
1466 : !> \param dimen_RI ...
1467 : !> \param kpoints ...
1468 : !> \param eps_eigval_S ...
1469 : ! **************************************************************************************************
1470 22 : SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, &
1471 : dimen_RI, kpoints, eps_eigval_S)
1472 :
1473 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: fm_matrix_Minv_L_kpoints, &
1474 : fm_matrix_L_kpoints
1475 : INTEGER, INTENT(IN) :: dimen_RI
1476 : TYPE(kpoint_type), POINTER :: kpoints
1477 : REAL(KIND=dp), INTENT(IN) :: eps_eigval_S
1478 :
1479 : CHARACTER(LEN=*), PARAMETER :: routineN = 'inversion_of_M_and_mult_with_chol_dec_of_V'
1480 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
1481 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp), ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
1482 :
1483 : INTEGER :: handle, ikp, nkp
1484 : TYPE(cp_cfm_type) :: cfm_matrix_K_tmp, cfm_matrix_M_tmp, &
1485 : cfm_matrix_V_tmp, cfm_matrix_Vtrunc_tmp
1486 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1487 :
1488 22 : CALL timeset(routineN, handle)
1489 :
1490 22 : CALL cp_fm_get_info(fm_matrix_Minv_L_kpoints(1, 1), matrix_struct=matrix_struct)
1491 :
1492 22 : CALL cp_cfm_create(cfm_matrix_M_tmp, matrix_struct)
1493 22 : CALL cp_cfm_create(cfm_matrix_V_tmp, matrix_struct)
1494 22 : CALL cp_cfm_create(cfm_matrix_K_tmp, matrix_struct)
1495 22 : CALL cp_cfm_create(cfm_matrix_Vtrunc_tmp, matrix_struct)
1496 :
1497 22 : CALL get_kpoint_info(kpoints, nkp=nkp)
1498 :
1499 498 : DO ikp = 1, nkp
1500 :
1501 476 : CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_M_tmp, cone, fm_matrix_Minv_L_kpoints(ikp, 1))
1502 476 : CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_M_tmp, ione, fm_matrix_Minv_L_kpoints(ikp, 2))
1503 :
1504 476 : CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_V_tmp, cone, fm_matrix_L_kpoints(ikp, 1))
1505 476 : CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_V_tmp, ione, fm_matrix_L_kpoints(ikp, 2))
1506 :
1507 476 : CALL cp_cfm_power(cfm_matrix_M_tmp, threshold=eps_eigval_S, exponent=-1.0_dp)
1508 :
1509 476 : CALL cp_cfm_power(cfm_matrix_V_tmp, threshold=0.0_dp, exponent=0.5_dp)
1510 :
1511 : ! save L(k) = SQRT(V(k))
1512 476 : CALL cp_cfm_to_fm(cfm_matrix_V_tmp, fm_matrix_L_kpoints(ikp, 1), fm_matrix_L_kpoints(ikp, 2))
1513 :
1514 : ! get K = M^(-1)*L, where L is the Cholesky decomposition of V = L^T*L
1515 : CALL parallel_gemm("N", "C", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_M_tmp, cfm_matrix_V_tmp, &
1516 476 : czero, cfm_matrix_K_tmp)
1517 :
1518 : ! move
1519 498 : CALL cp_cfm_to_fm(cfm_matrix_K_tmp, fm_matrix_Minv_L_kpoints(ikp, 1), fm_matrix_Minv_L_kpoints(ikp, 2))
1520 :
1521 : END DO
1522 :
1523 22 : CALL cp_cfm_release(cfm_matrix_M_tmp)
1524 22 : CALL cp_cfm_release(cfm_matrix_V_tmp)
1525 22 : CALL cp_cfm_release(cfm_matrix_K_tmp)
1526 22 : CALL cp_cfm_release(cfm_matrix_Vtrunc_tmp)
1527 :
1528 22 : CALL timestop(handle)
1529 :
1530 22 : END SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V
1531 :
1532 : ! **************************************************************************************************
1533 : !> \brief ...
1534 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
1535 : !> \param fm_matrix_Minv ...
1536 : !> \param qs_env ...
1537 : ! **************************************************************************************************
1538 88 : SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
1539 22 : fm_matrix_Minv, qs_env)
1540 :
1541 : TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_Minv_Vtrunc_Minv, &
1542 : fm_matrix_Minv
1543 : TYPE(qs_environment_type), POINTER :: qs_env
1544 :
1545 : CHARACTER(LEN=*), PARAMETER :: &
1546 : routineN = 'Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc'
1547 :
1548 : INTEGER :: dimen_RI, handle, ndep
1549 : REAL(KIND=dp) :: eps_eigval_S_Gamma
1550 : TYPE(cp_fm_type) :: fm_matrix_RI_metric_inv_work, fm_work
1551 :
1552 22 : CALL timeset(routineN, handle)
1553 :
1554 22 : CALL cp_fm_create(fm_work, fm_matrix_Minv(1, 1)%matrix_struct)
1555 22 : CALL cp_fm_set_all(fm_work, 0.0_dp)
1556 :
1557 22 : CALL cp_fm_create(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1)%matrix_struct)
1558 22 : CALL cp_fm_set_all(fm_matrix_RI_metric_inv_work, 0.0_dp)
1559 :
1560 22 : CALL cp_fm_get_info(matrix=fm_matrix_Minv(1, 1), nrow_global=dimen_RI)
1561 :
1562 22 : eps_eigval_S_Gamma = qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S_Gamma
1563 :
1564 22 : IF (eps_eigval_S_Gamma > 1.0E-18) THEN
1565 :
1566 : ! remove small eigenvalues
1567 : CALL cp_fm_power(fm_matrix_Minv(1, 1), fm_matrix_RI_metric_inv_work, -0.5_dp, &
1568 0 : eps_eigval_S_Gamma, ndep)
1569 :
1570 : ELSE
1571 :
1572 22 : CALL cholesky_decomp(fm_matrix_Minv(1, 1), dimen_RI, do_inversion=.TRUE.)
1573 :
1574 22 : CALL invert_mat(fm_matrix_Minv(1, 1))
1575 :
1576 : END IF
1577 :
1578 : CALL parallel_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_Minv(1, 1), &
1579 22 : fm_matrix_Minv(1, 1), 0.0_dp, fm_matrix_RI_metric_inv_work)
1580 :
1581 22 : CALL cp_fm_to_fm(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1))
1582 :
1583 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_RI_metric_inv_work, &
1584 22 : fm_matrix_Minv_Vtrunc_Minv(1, 1), 0.0_dp, fm_work)
1585 :
1586 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work, &
1587 22 : fm_matrix_RI_metric_inv_work, 0.0_dp, fm_matrix_Minv_Vtrunc_Minv(1, 1))
1588 :
1589 22 : CALL cp_fm_release(fm_work)
1590 22 : CALL cp_fm_release(fm_matrix_RI_metric_inv_work)
1591 :
1592 22 : CALL timestop(handle)
1593 :
1594 22 : END SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc
1595 :
1596 : ! **************************************************************************************************
1597 : !> \brief ...
1598 : !> \param qs_env ...
1599 : !> \param trunc_coulomb ...
1600 : !> \param rel_cutoff_trunc_coulomb_ri_x ...
1601 : !> \param cell_grid ...
1602 : !> \param do_BvK_cell ...
1603 : ! **************************************************************************************************
1604 44 : SUBROUTINE trunc_coulomb_for_exchange(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x, &
1605 : cell_grid, do_BvK_cell)
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 : INTEGER, DIMENSION(3), OPTIONAL :: cell_grid
1610 : LOGICAL, OPTIONAL :: do_BvK_cell
1611 :
1612 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trunc_coulomb_for_exchange'
1613 :
1614 : INTEGER :: handle, i_dim
1615 : INTEGER, DIMENSION(3) :: periodic
1616 : LOGICAL :: my_do_BvK_cell
1617 : REAL(KIND=dp) :: kp_fac, kp_fac_idim, my_rel_cutoff_trunc_coulomb_ri_x, &
1618 : shortest_dist_cell_planes
1619 : TYPE(cell_type), POINTER :: cell
1620 : TYPE(kpoint_type), POINTER :: kpoints_scf
1621 :
1622 44 : CALL timeset(routineN, handle)
1623 :
1624 44 : NULLIFY (cell)
1625 44 : CALL get_qs_env(qs_env, cell=cell, kpoints=kpoints_scf)
1626 44 : CALL get_cell(cell=cell, periodic=periodic)
1627 :
1628 44 : my_do_BvK_cell = .FALSE.
1629 44 : IF (PRESENT(do_BvK_cell)) my_do_BvK_cell = do_BvK_cell
1630 22 : IF (my_do_BvK_cell) THEN
1631 : kp_fac = 1.0E10_dp
1632 24 : DO i_dim = 1, 3
1633 : ! look for smallest k-point mesh in periodic direction
1634 24 : IF (periodic(i_dim) == 1) THEN
1635 12 : IF (PRESENT(cell_grid)) THEN
1636 12 : kp_fac_idim = REAL(cell_grid(i_dim), KIND=dp)
1637 : ELSE
1638 0 : kp_fac_idim = REAL(kpoints_scf%nkp_grid(i_dim), KIND=dp)
1639 : END IF
1640 12 : IF (kp_fac > kp_fac_idim) kp_fac = kp_fac_idim
1641 : END IF
1642 : END DO
1643 : ELSE
1644 : kp_fac = 1.0_dp
1645 : END IF
1646 :
1647 44 : shortest_dist_cell_planes = 1.0E4_dp
1648 44 : IF (periodic(1) == 1) THEN
1649 16 : IF (shortest_dist_cell_planes > plane_distance(1, 0, 0, cell)) THEN
1650 16 : shortest_dist_cell_planes = plane_distance(1, 0, 0, cell)
1651 : END IF
1652 : END IF
1653 44 : IF (periodic(2) == 1) THEN
1654 36 : IF (shortest_dist_cell_planes > plane_distance(0, 1, 0, cell)) THEN
1655 24 : shortest_dist_cell_planes = plane_distance(0, 1, 0, cell)
1656 : END IF
1657 : END IF
1658 44 : IF (periodic(3) == 1) THEN
1659 32 : IF (shortest_dist_cell_planes > plane_distance(0, 0, 1, cell)) THEN
1660 8 : shortest_dist_cell_planes = plane_distance(0, 0, 1, cell)
1661 : END IF
1662 : END IF
1663 :
1664 44 : IF (PRESENT(rel_cutoff_trunc_coulomb_ri_x)) THEN
1665 22 : my_rel_cutoff_trunc_coulomb_ri_x = rel_cutoff_trunc_coulomb_ri_x
1666 : ELSE
1667 22 : my_rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
1668 : END IF
1669 :
1670 44 : IF (PRESENT(trunc_coulomb)) THEN
1671 44 : trunc_coulomb%potential_type = do_potential_truncated
1672 : trunc_coulomb%cutoff_radius = shortest_dist_cell_planes* &
1673 : my_rel_cutoff_trunc_coulomb_ri_x* &
1674 44 : kp_fac
1675 44 : trunc_coulomb%filename = "t_c_g.dat"
1676 : ! dummy
1677 44 : trunc_coulomb%omega = 0.0_dp
1678 : END IF
1679 :
1680 44 : CALL timestop(handle)
1681 :
1682 44 : END SUBROUTINE trunc_coulomb_for_exchange
1683 :
1684 : ! **************************************************************************************************
1685 : !> \brief ...
1686 : !> \param fm_matrix_L ...
1687 : ! **************************************************************************************************
1688 1236 : SUBROUTINE invert_mat(fm_matrix_L)
1689 :
1690 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1691 :
1692 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_mat'
1693 :
1694 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
1695 : ncol_local, nrow_local
1696 618 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1697 :
1698 618 : CALL timeset(routineN, handle)
1699 :
1700 618 : CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U')
1701 :
1702 : ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1703 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1704 : nrow_local=nrow_local, &
1705 : ncol_local=ncol_local, &
1706 : row_indices=row_indices, &
1707 618 : col_indices=col_indices)
1708 41060 : DO iiB = 1, nrow_local
1709 40442 : i_global = row_indices(iiB)
1710 3383046 : DO jjB = 1, ncol_local
1711 3341986 : j_global = col_indices(jjB)
1712 3382428 : IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
1713 : END DO
1714 : END DO
1715 :
1716 618 : CALL timestop(handle)
1717 :
1718 618 : END SUBROUTINE invert_mat
1719 :
1720 : ! **************************************************************************************************
1721 : !> \brief ...
1722 : !> \param fm_matrix_L ...
1723 : !> \param blacs_env_L ...
1724 : !> \param dimen_RI ...
1725 : !> \param para_env_L ...
1726 : !> \param name ...
1727 : !> \param fm_matrix_L_extern ...
1728 : ! **************************************************************************************************
1729 654 : SUBROUTINE create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, fm_matrix_L_extern)
1730 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_L
1731 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_L
1732 : INTEGER, INTENT(IN) :: dimen_RI
1733 : TYPE(mp_para_env_type), POINTER :: para_env_L
1734 : CHARACTER(LEN=*), INTENT(IN) :: name
1735 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_L_extern
1736 :
1737 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_matrix_L'
1738 :
1739 : INTEGER :: handle
1740 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1741 :
1742 654 : CALL timeset(routineN, handle)
1743 :
1744 : ! create the full matrix L defined in the L group
1745 654 : IF (PRESENT(fm_matrix_L_extern)) THEN
1746 36 : CALL cp_fm_create(fm_matrix_L, fm_matrix_L_extern%matrix_struct, name=name)
1747 : ELSE
1748 618 : NULLIFY (fm_struct)
1749 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_L, nrow_global=dimen_RI, &
1750 618 : ncol_global=dimen_RI, para_env=para_env_L)
1751 :
1752 618 : CALL cp_fm_create(fm_matrix_L, fm_struct, name=name)
1753 :
1754 618 : CALL cp_fm_struct_release(fm_struct)
1755 : END IF
1756 :
1757 654 : CALL cp_fm_set_all(matrix=fm_matrix_L, alpha=0.0_dp)
1758 :
1759 654 : CALL timestop(handle)
1760 :
1761 654 : END SUBROUTINE create_matrix_L
1762 :
1763 : ! **************************************************************************************************
1764 : !> \brief ...
1765 : !> \param fm_matrix_L ...
1766 : !> \param my_group_L_start ...
1767 : !> \param my_group_L_end ...
1768 : !> \param my_group_L_size ...
1769 : !> \param my_Lrows ...
1770 : ! **************************************************************************************************
1771 496 : SUBROUTINE grep_Lcols(fm_matrix_L, &
1772 : my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
1773 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1774 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, &
1775 : my_group_L_size
1776 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1777 : INTENT(OUT) :: my_Lrows
1778 :
1779 : CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_Lcols'
1780 :
1781 : INTEGER :: dimen_RI, handle, handle2, i_global, iiB, j_global, jjB, max_row_col_local, &
1782 : ncol_local, ncol_rec, nrow_local, nrow_rec, proc_receive_static, proc_send_static, &
1783 : proc_shift
1784 496 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info
1785 496 : INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, &
1786 496 : row_indices, row_indices_rec
1787 496 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: local_L, rec_L
1788 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1789 496 : POINTER :: local_L_internal
1790 : TYPE(mp_para_env_type), POINTER :: para_env
1791 :
1792 496 : CALL timeset(routineN, handle)
1793 :
1794 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1795 : nrow_local=nrow_local, &
1796 : ncol_local=ncol_local, &
1797 : row_indices=row_indices, &
1798 : col_indices=col_indices, &
1799 : nrow_global=dimen_RI, &
1800 : local_data=local_L_internal, &
1801 496 : para_env=para_env)
1802 :
1803 1984 : ALLOCATE (my_Lrows(dimen_RI, my_group_L_size))
1804 1598968 : my_Lrows = 0.0_dp
1805 :
1806 1984 : ALLOCATE (local_L(nrow_local, ncol_local))
1807 2928436 : local_L(:, :) = local_L_internal(1:nrow_local, 1:ncol_local)
1808 :
1809 496 : max_row_col_local = MAX(nrow_local, ncol_local)
1810 496 : CALL para_env%max(max_row_col_local)
1811 :
1812 1984 : ALLOCATE (local_col_row_info(0:max_row_col_local, 2))
1813 76116 : local_col_row_info = 0
1814 : ! 0,1 nrows
1815 496 : local_col_row_info(0, 1) = nrow_local
1816 36311 : local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
1817 : ! 0,2 ncols
1818 496 : local_col_row_info(0, 2) = ncol_local
1819 37204 : local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
1820 :
1821 992 : ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
1822 :
1823 : ! accumulate data on my_Lrows starting from myself
1824 37204 : DO jjB = 1, ncol_local
1825 36708 : j_global = col_indices(jjB)
1826 37204 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
1827 1509484 : DO iiB = 1, nrow_local
1828 1490110 : i_global = row_indices(iiB)
1829 1509484 : my_Lrows(i_global, j_global - my_group_L_start + 1) = local_L(iiB, jjB)
1830 : END DO
1831 : END IF
1832 : END DO
1833 :
1834 496 : proc_send_static = MODULO(para_env%mepos + 1, para_env%num_pe)
1835 496 : proc_receive_static = MODULO(para_env%mepos - 1, para_env%num_pe)
1836 :
1837 496 : CALL timeset(routineN//"_comm", handle2)
1838 :
1839 520 : DO proc_shift = 1, para_env%num_pe - 1
1840 : ! first exchange information on the local data
1841 4200 : rec_col_row_info = 0
1842 24 : CALL para_env%sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static)
1843 24 : nrow_rec = rec_col_row_info(0, 1)
1844 24 : ncol_rec = rec_col_row_info(0, 2)
1845 :
1846 72 : ALLOCATE (row_indices_rec(nrow_rec))
1847 1061 : row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1848 :
1849 72 : ALLOCATE (col_indices_rec(ncol_rec))
1850 2064 : col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1851 :
1852 96 : ALLOCATE (rec_L(nrow_rec, ncol_rec))
1853 91052 : rec_L = 0.0_dp
1854 :
1855 : ! then send and receive the real data
1856 24 : CALL para_env%sendrecv(local_L, proc_send_static, rec_L, proc_receive_static)
1857 :
1858 : ! accumulate the received data on my_Lrows
1859 2064 : DO jjB = 1, ncol_rec
1860 2040 : j_global = col_indices_rec(jjB)
1861 2064 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
1862 91028 : DO iiB = 1, nrow_rec
1863 88988 : i_global = row_indices_rec(iiB)
1864 91028 : my_Lrows(i_global, j_global - my_group_L_start + 1) = rec_L(iiB, jjB)
1865 : END DO
1866 : END IF
1867 : END DO
1868 :
1869 4200 : local_col_row_info(:, :) = rec_col_row_info
1870 24 : CALL MOVE_ALLOC(rec_L, local_L)
1871 :
1872 24 : DEALLOCATE (col_indices_rec)
1873 520 : DEALLOCATE (row_indices_rec)
1874 : END DO
1875 496 : CALL timestop(handle2)
1876 :
1877 496 : DEALLOCATE (local_col_row_info)
1878 496 : DEALLOCATE (rec_col_row_info)
1879 496 : DEALLOCATE (local_L)
1880 :
1881 496 : CALL timestop(handle)
1882 :
1883 1488 : END SUBROUTINE grep_Lcols
1884 :
1885 : END MODULE mp2_ri_2c
|
