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