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