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