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 654 : 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 654 : 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 654 : num_small_eigen, potential, sab_orb_sub, do_im_time=do_im_time)
212 :
213 654 : 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 654 : dimen_RI_red = dimen_RI
220 :
221 654 : 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 484 : dimen_RI, dimen_RI_red, para_env)
224 : ELSE
225 :
226 : ! RI-metric matrix (depending on RI metric: Coulomb, overlap or truncated Coulomb)
227 170 : IF (do_im_time) THEN
228 134 : 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 134 : eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
232 134 : dft_control%qs_control%eps_pgf_orb = mp2_eps_pgf_orb_S
233 134 : 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 134 : regularization_RI=qs_env%mp2_env%ri_rpa_im_time%regularization_RI)
238 :
239 : ! re-init the radii to previous values
240 134 : dft_control%qs_control%eps_pgf_orb = eps_pgf_orb_old
241 134 : 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 170 : IF (do_kpoints) THEN
262 :
263 : ! k-dependent 1/r Coulomb matrix V_PQ(k)
264 20 : 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 20 : kpoints, qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S)
268 :
269 20 : 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 20 : 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 20 : 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 20 : 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 170 : 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 654 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
331 :
332 1174 : IF (.NOT. do_im_time) THEN
333 520 : IF (unit_nr > 0) THEN
334 260 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Cholesky decomposition group size:", para_env_L%num_pe
335 260 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") "RI_INFO| Number of groups for auxiliary basis functions", ngroup
336 260 : 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 260 : CALL m_flush(unit_nr)
343 : END IF
344 :
345 : ! replicate the necessary row of the L^{-1} matrix on each proc
346 520 : CALL grep_Lcols(fm_matrix_L_work, my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
347 532 : 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 654 : CALL cp_fm_release(fm_matrix_L_work)
353 654 : 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 654 : CALL mp_para_env_release(para_env_L)
356 :
357 654 : CALL timestop(handle)
358 :
359 654 : 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 796 : 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 796 : 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 752 : CALL cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion=do_inversion)
403 :
404 752 : IF (do_inversion) CALL invert_mat(fm_matrix_L)
405 : END IF
406 :
407 796 : 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 20 : 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 20 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
428 : TYPE(cell_type), POINTER :: cell
429 20 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_RI_aux_transl, matrix_v_RI_kp
430 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
431 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
432 :
433 20 : CALL timeset(routineN, handle)
434 :
435 20 : 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 20 : atomic_kind_set=atomic_kind_set)
443 :
444 : ! check that we have a 2n x 2m x 2k mesh to guarantee good convergence
445 20 : CALL get_cell(cell=cell, periodic=periodic)
446 80 : DO i_dim = 1, 3
447 80 : IF (periodic(i_dim) == 1) THEN
448 36 : CPASSERT(MODULO(kpoints%nkp_grid(i_dim), 2) == 0)
449 : END IF
450 : END DO
451 :
452 20 : nkp = kpoints%nkp
453 :
454 948 : ALLOCATE (fm_matrix_L_kpoints(nkp, 2))
455 60 : DO i_real_imag = 1, 2
456 908 : DO ikp = 1, nkp
457 : CALL cp_fm_create(fm_matrix_L_kpoints(ikp, i_real_imag), &
458 848 : fm_matrix_Minv_L_kpoints(1, i_real_imag)%matrix_struct)
459 888 : CALL cp_fm_set_all(fm_matrix_L_kpoints(ikp, i_real_imag), 0.0_dp)
460 : END DO
461 : END DO
462 :
463 20 : CALL allocate_matrix_v_RI_kp(matrix_v_RI_kp, matrix_s_RI_aux_transl, nkp)
464 :
465 20 : IF (qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints) THEN
466 :
467 16 : nkp_orig = qs_env%mp2_env%ri_rpa_im_time%nkp_orig
468 16 : 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 16 : operator_type=operator_coulomb, ikp_start=1, ikp_end=nkp_orig)
475 :
476 64 : nkp_grid_orig = kpoints%nkp_grid
477 128 : 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 16 : operator_type=operator_coulomb, ikp_start=nkp_orig + 1, ikp_end=nkp)
484 :
485 64 : 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 444 : DO ikp = 1, nkp
498 :
499 424 : CALL copy_dbcsr_to_fm(matrix_v_RI_kp(ikp, 1)%matrix, fm_matrix_L_kpoints(ikp, 1))
500 444 : CALL copy_dbcsr_to_fm(matrix_v_RI_kp(ikp, 2)%matrix, fm_matrix_L_kpoints(ikp, 2))
501 :
502 : END DO
503 :
504 20 : CALL dbcsr_deallocate_matrix_set(matrix_v_RI_kp)
505 :
506 20 : CALL timestop(handle)
507 :
508 20 : 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 20 : 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 20 : NULLIFY (matrix_v_RI_kp)
524 20 : CALL dbcsr_allocate_matrix_set(matrix_v_RI_kp, nkp, 2)
525 :
526 444 : DO ikp = 1, nkp
527 :
528 424 : 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 424 : matrix_type=dbcsr_type_no_symmetry)
531 424 : CALL dbcsr_reserve_all_blocks(matrix_v_RI_kp(ikp, 1)%matrix)
532 424 : CALL dbcsr_set(matrix_v_RI_kp(ikp, 1)%matrix, 0.0_dp)
533 :
534 424 : 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 424 : matrix_type=dbcsr_type_no_symmetry)
537 424 : CALL dbcsr_reserve_all_blocks(matrix_v_RI_kp(ikp, 2)%matrix)
538 444 : CALL dbcsr_set(matrix_v_RI_kp(ikp, 2)%matrix, 0.0_dp)
539 :
540 : END DO
541 :
542 20 : 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 438 : 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 438 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sizes_RI
580 438 : INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
581 438 : 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 438 : 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 438 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_RI_aux_transl
590 438 : 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 438 : 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 438 : POINTER :: sab_RI
599 438 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
600 438 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
601 :
602 438 : CALL timeset(routineN, handle)
603 :
604 438 : NULLIFY (sab_RI, matrix_s_RI_aux_transl, dist_2d)
605 :
606 438 : IF (PRESENT(regularization_RI)) THEN
607 390 : my_regularization_RI = regularization_RI
608 : ELSE
609 48 : my_regularization_RI = 0.0_dp
610 : END IF
611 :
612 438 : IF (PRESENT(put_mat_KS_env)) THEN
613 174 : my_put_mat_KS_env = put_mat_KS_env
614 : ELSE
615 : my_put_mat_KS_env = .FALSE.
616 : END IF
617 :
618 438 : IF (PRESENT(do_build_cell_index)) THEN
619 216 : 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 438 : natom=natom)
633 :
634 1314 : ALLOCATE (sizes_RI(natom))
635 2152 : ALLOCATE (basis_set_RI(nkind))
636 :
637 438 : CALL basis_set_list_setup(basis_set_RI, 'RI_AUX', qs_kind_set)
638 438 : 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 438 : sym_ij=.TRUE., dist_2d=dist_2d)
642 :
643 438 : CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
644 1314 : ALLOCATE (row_bsize(SIZE(sizes_RI)))
645 1314 : ALLOCATE (col_bsize(SIZE(sizes_RI)))
646 1622 : row_bsize(:) = sizes_RI
647 1622 : col_bsize(:) = sizes_RI
648 :
649 438 : IF (do_kpoints) THEN
650 236 : CPASSERT(PRESENT(kpoints))
651 236 : 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 236 : cell_to_index=cell_to_index)
656 236 : n_real_imag = 2
657 236 : nimg = dft_control%nimages
658 : ELSE
659 202 : nkp = 1
660 202 : n_real_imag = 1
661 202 : nimg = 1
662 : END IF
663 :
664 2646 : ALLOCATE (mat_2c(nimg))
665 : CALL dbcsr_create(mat_2c(1), "(RI|RI)", dbcsr_dist, dbcsr_type_symmetric, &
666 438 : row_bsize, col_bsize)
667 438 : DEALLOCATE (row_bsize, col_bsize)
668 :
669 1332 : DO img = 2, nimg
670 1332 : 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 438 : regularization_RI=my_regularization_RI)
676 :
677 438 : CALL dbcsr_distribution_release(dbcsr_dist)
678 438 : DEALLOCATE (basis_set_RI)
679 :
680 438 : IF (my_put_mat_KS_env) THEN
681 134 : CALL get_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
682 134 : IF (ASSOCIATED(matrix_s_RI_aux_transl)) CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
683 : END IF
684 438 : CALL dbcsr_allocate_matrix_set(matrix_s_RI_aux_transl, 1, nimg)
685 1770 : DO img = 1, nimg
686 1332 : ALLOCATE (matrix_s_RI_aux_transl(1, img)%matrix)
687 1332 : CALL dbcsr_copy(matrix_s_RI_aux_transl(1, img)%matrix, mat_2c(img))
688 1770 : CALL dbcsr_release(mat_2c(img))
689 : END DO
690 :
691 438 : IF (my_put_mat_KS_env) THEN
692 134 : CALL set_ks_env(qs_env%ks_env, matrix_s_RI_aux_kp=matrix_s_RI_aux_transl)
693 : END IF
694 :
695 438 : IF (PRESENT(ikp_ext)) nkp = 1
696 :
697 3908 : ALLOCATE (fm_matrix_Minv_L_kpoints(nkp, n_real_imag))
698 1112 : DO i_real_imag = 1, n_real_imag
699 2594 : DO i_size = 1, nkp
700 1482 : CALL cp_fm_create(fm_matrix_Minv_L_kpoints(i_size, i_real_imag), fm_matrix_L%matrix_struct)
701 2156 : 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 438 : NULLIFY (fm_struct)
706 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=dimen_RI, &
707 438 : ncol_global=dimen_RI, para_env=para_env)
708 :
709 438 : CALL cp_fm_create(fm_matrix_S_global, fm_struct)
710 438 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
711 :
712 438 : IF (do_kpoints) THEN
713 :
714 236 : ALLOCATE (rmatrix, cmatrix, tmpmat)
715 : CALL dbcsr_create(rmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
716 236 : matrix_type=dbcsr_type_symmetric)
717 : CALL dbcsr_create(cmatrix, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
718 236 : matrix_type=dbcsr_type_antisymmetric)
719 : CALL dbcsr_create(tmpmat, template=matrix_s_RI_aux_transl(1, 1)%matrix, &
720 236 : matrix_type=dbcsr_type_no_symmetry)
721 236 : CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_RI)
722 236 : CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_RI)
723 :
724 876 : DO ikp = 1, nkp
725 :
726 : ! s matrix is not spin dependent, double the work
727 640 : CALL dbcsr_set(rmatrix, 0.0_dp)
728 640 : CALL dbcsr_set(cmatrix, 0.0_dp)
729 :
730 640 : IF (PRESENT(ikp_ext)) THEN
731 216 : 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 640 : xkp=xkp(1:3, ikp_for_xkp), cell_to_index=cell_to_index, sab_nl=sab_RI)
738 :
739 640 : CALL dbcsr_set(tmpmat, 0.0_dp)
740 640 : CALL dbcsr_desymmetrize(rmatrix, tmpmat)
741 :
742 640 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
743 640 : CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
744 640 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 1), para_env)
745 :
746 640 : CALL dbcsr_set(tmpmat, 0.0_dp)
747 640 : CALL dbcsr_desymmetrize(cmatrix, tmpmat)
748 :
749 640 : CALL cp_fm_set_all(fm_matrix_S_global, 0.0_dp)
750 640 : CALL copy_dbcsr_to_fm(tmpmat, fm_matrix_S_global)
751 876 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(ikp, 2), para_env)
752 :
753 : END DO
754 :
755 236 : CALL dbcsr_deallocate_matrix(rmatrix)
756 236 : CALL dbcsr_deallocate_matrix(cmatrix)
757 236 : CALL dbcsr_deallocate_matrix(tmpmat)
758 :
759 : ELSE
760 :
761 : NULLIFY (matrix_s_RI_aux_desymm)
762 202 : 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 202 : name='S_RI non_symm', matrix_type=dbcsr_type_no_symmetry)
765 :
766 202 : CALL dbcsr_desymmetrize(matrix_s_RI_aux_transl(1, 1)%matrix, matrix_s_RI_aux_desymm)
767 :
768 202 : CALL copy_dbcsr_to_fm(matrix_s_RI_aux_desymm, fm_matrix_S_global)
769 :
770 202 : CALL cp_fm_copy_general(fm_matrix_S_global, fm_matrix_Minv_L_kpoints(1, 1), para_env)
771 :
772 202 : CALL dbcsr_deallocate_matrix(matrix_s_RI_aux_desymm)
773 :
774 : END IF
775 :
776 438 : CALL release_neighbor_list_sets(sab_RI)
777 :
778 438 : CALL cp_fm_struct_release(fm_struct)
779 :
780 438 : CALL cp_fm_release(fm_matrix_S_global)
781 :
782 438 : IF (.NOT. my_put_mat_KS_env) THEN
783 304 : CALL dbcsr_deallocate_matrix_set(matrix_s_RI_aux_transl)
784 : END IF
785 :
786 438 : CALL timestop(handle)
787 :
788 1752 : 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 690 : 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 690 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
849 690 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
850 : LOGICAL :: my_do_im_time
851 : REAL(KIND=dp) :: min_mem_for_QK
852 690 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: egen_L
853 690 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: L_external_col, L_local_col
854 690 : 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 690 : TYPE(group_dist_d1_type) :: gd_sub_array
858 : TYPE(gto_basis_set_type), POINTER :: basis_set_a
859 690 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
860 690 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
861 :
862 690 : CALL timeset(routineN, handle)
863 :
864 690 : my_do_im_time = .FALSE.
865 690 : IF (PRESENT(do_im_time)) THEN
866 654 : 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 690 : particle_set=particle_set)
872 :
873 690 : nkind = SIZE(qs_kind_set)
874 690 : natom = SIZE(particle_set)
875 :
876 690 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
877 :
878 1960 : DO ikind = 1, nkind
879 1270 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=basis_set_a, basis_type="RI_AUX")
880 1960 : CPASSERT(ASSOCIATED(basis_set_a))
881 : END DO
882 :
883 690 : dimen_RI = 0
884 2714 : DO iatom = 1, natom
885 2024 : ikind = kind_of(iatom)
886 2024 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), nsgf=nsgf, basis_type="RI_AUX")
887 2714 : dimen_RI = dimen_RI + nsgf
888 : END DO
889 :
890 : ! check that very small systems are not running on too many processes
891 690 : 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 690 : CALL create_group_dist(gd_array, ngroup, dimen_RI)
897 :
898 690 : CALL get_group_dist(gd_array, color_sub, my_group_L_start, my_group_L_end, my_group_L_size)
899 :
900 690 : CALL timeset(routineN//"_loop_lm", handle2)
901 :
902 2760 : ALLOCATE (L_local_col(dimen_RI, my_group_L_size))
903 2365350 : L_local_col = 0.0_dp
904 :
905 690 : potential_type = potential%potential_type
906 :
907 : IF ((eri_method .EQ. do_eri_mme .OR. eri_method .EQ. do_eri_os) &
908 690 : .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 292 : 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 398 : .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 398 : natom, potential, sab_orb_sub, L_local_col, kind_of)
922 :
923 : ELSE
924 0 : CPABORT("unknown ERI method")
925 : END IF
926 :
927 690 : 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 690 : best_group_size = para_env%num_pe
932 :
933 690 : strat_group_size = MAX(1, para_env%num_pe/10)
934 :
935 690 : min_mem_for_QK = REAL(dimen_RI, KIND=dp)*dimen_RI*3.0_dp*8.0_dp/1024_dp/1024_dp
936 :
937 690 : group_size = strat_group_size - 1
938 724 : DO iproc = strat_group_size, para_env%num_pe
939 724 : 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 724 : 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 690 : 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 690 : IF (my_do_im_time) THEN
952 : ! para_env_L is the para_env for im. time to avoid bug
953 134 : best_group_size = para_env%num_pe
954 : END IF
955 :
956 : ! create the L group
957 690 : color_L = para_env%mepos/best_group_size
958 690 : ALLOCATE (para_env_L)
959 690 : CALL para_env_L%from_split(para_env, color_L)
960 :
961 : ! create the blacs_L
962 690 : NULLIFY (blacs_env_L)
963 690 : 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 690 : CALL create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, "fm_matrix_L", fm_matrix_L_extern)
967 :
968 690 : 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 134 : dimen_RI)
973 :
974 : ELSE
975 556 : BLOCK
976 : TYPE(mp_comm_type) :: comm_exchange
977 :
978 556 : 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 556 : 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 556 : col_indices=col_indices)
988 :
989 41984 : DO jjB = 1, ncol_local
990 41428 : j_global = col_indices(jjB)
991 41984 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
992 1705702 : DO iiB = 1, nrow_local
993 1683868 : i_global = row_indices(iiB)
994 1705702 : 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 556 : proc_send_static = MODULO(comm_exchange%mepos + 1, comm_exchange%num_pe)
1000 556 : proc_receive_static = MODULO(comm_exchange%mepos - 1, comm_exchange%num_pe)
1001 :
1002 556 : 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 556 : CALL MOVE_ALLOC(L_external_col, L_local_col)
1023 : END DO
1024 :
1025 556 : CALL release_group_dist(gd_sub_array)
1026 1112 : CALL comm_exchange%free()
1027 : END BLOCK
1028 : END IF
1029 :
1030 690 : 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 690 : comm_exchange = fm_matrix_L%matrix_struct%context%interconnect(para_env)
1036 690 : CALL comm_exchange%sum(fm_matrix_L%local_data)
1037 1380 : CALL comm_exchange%free()
1038 : END BLOCK
1039 :
1040 690 : cond_num = 1.0_dp
1041 690 : num_small_eigen = 0
1042 690 : 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 690 : CALL cp_blacs_env_release(blacs_env_L)
1069 :
1070 690 : CALL timestop(handle)
1071 :
1072 3450 : 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 134 : 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, iproc, j_col, LLL, MMM, &
1228 : ncol_local, nrow_local, proc_send, send_pcol, send_prow
1229 134 : INTEGER, ALLOCATABLE, DIMENSION(:) :: entry_counter, num_entries_rec, &
1230 : num_entries_send
1231 134 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1232 : TYPE(integ_mat_buffer_type), ALLOCATABLE, &
1233 134 : DIMENSION(:) :: buffer_rec, buffer_send
1234 134 : TYPE(mp_request_type), DIMENSION(:, :), POINTER :: req_array
1235 :
1236 134 : CALL timeset(routineN, handle)
1237 :
1238 134 : 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 134 : col_indices=col_indices)
1246 :
1247 402 : ALLOCATE (num_entries_rec(0:para_env%num_pe - 1))
1248 268 : ALLOCATE (num_entries_send(0:para_env%num_pe - 1))
1249 :
1250 402 : num_entries_rec(:) = 0
1251 402 : num_entries_send(:) = 0
1252 :
1253 134 : dummy_proc = 0
1254 :
1255 : ! get the process, where the elements from L_local_col have to be sent
1256 11462 : DO LLL = 1, dimen_RI
1257 :
1258 11328 : send_prow = fm_matrix_L%matrix_struct%g2p_row(LLL)
1259 :
1260 567136 : DO MMM = my_group_L_start, my_group_L_end
1261 :
1262 555674 : send_pcol = fm_matrix_L%matrix_struct%g2p_col(MMM)
1263 :
1264 555674 : proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1265 :
1266 567002 : num_entries_send(proc_send) = num_entries_send(proc_send) + 1
1267 :
1268 : END DO
1269 :
1270 : END DO
1271 :
1272 134 : CALL timestop(handle2)
1273 :
1274 134 : CALL timeset(routineN//"_2", handle2)
1275 :
1276 134 : CALL para_env%alltoall(num_entries_send, num_entries_rec, 1)
1277 :
1278 134 : CALL timestop(handle2)
1279 :
1280 134 : CALL timeset(routineN//"_3", handle2)
1281 :
1282 : ! allocate buffers to send the entries and the information of the entries
1283 670 : ALLOCATE (buffer_rec(0:para_env%num_pe - 1))
1284 536 : ALLOCATE (buffer_send(0:para_env%num_pe - 1))
1285 :
1286 : ! allocate data message and corresponding indices
1287 402 : DO iproc = 0, para_env%num_pe - 1
1288 :
1289 804 : ALLOCATE (buffer_rec(iproc)%msg(num_entries_rec(iproc)))
1290 556076 : buffer_rec(iproc)%msg = 0.0_dp
1291 :
1292 : END DO
1293 :
1294 134 : CALL timestop(handle2)
1295 :
1296 134 : CALL timeset(routineN//"_4", handle2)
1297 :
1298 402 : DO iproc = 0, para_env%num_pe - 1
1299 :
1300 804 : ALLOCATE (buffer_send(iproc)%msg(num_entries_send(iproc)))
1301 556076 : buffer_send(iproc)%msg = 0.0_dp
1302 :
1303 : END DO
1304 :
1305 134 : CALL timestop(handle2)
1306 :
1307 134 : CALL timeset(routineN//"_5", handle2)
1308 :
1309 402 : DO iproc = 0, para_env%num_pe - 1
1310 :
1311 804 : ALLOCATE (buffer_rec(iproc)%indx(num_entries_rec(iproc), 2))
1312 1112286 : buffer_rec(iproc)%indx = 0
1313 :
1314 : END DO
1315 :
1316 134 : CALL timestop(handle2)
1317 :
1318 134 : CALL timeset(routineN//"_6", handle2)
1319 :
1320 402 : DO iproc = 0, para_env%num_pe - 1
1321 :
1322 804 : ALLOCATE (buffer_send(iproc)%indx(num_entries_send(iproc), 2))
1323 1112286 : buffer_send(iproc)%indx = 0
1324 :
1325 : END DO
1326 :
1327 134 : CALL timestop(handle2)
1328 :
1329 134 : CALL timeset(routineN//"_7", handle2)
1330 :
1331 268 : ALLOCATE (entry_counter(0:para_env%num_pe - 1))
1332 402 : 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 6130 : DO MMM = my_group_L_start, my_group_L_end
1337 :
1338 5996 : send_pcol = fm_matrix_L%matrix_struct%g2p_col(MMM)
1339 :
1340 561804 : DO LLL = 1, dimen_RI
1341 :
1342 555674 : send_prow = fm_matrix_L%matrix_struct%g2p_row(LLL)
1343 :
1344 555674 : proc_send = fm_matrix_L%matrix_struct%context%blacs2mpi(send_prow, send_pcol)
1345 :
1346 555674 : entry_counter(proc_send) = entry_counter(proc_send) + 1
1347 :
1348 : buffer_send(proc_send)%msg(entry_counter(proc_send)) = &
1349 555674 : L_local_col(LLL, MMM - my_group_L_start + 1)
1350 :
1351 555674 : buffer_send(proc_send)%indx(entry_counter(proc_send), 1) = LLL
1352 561670 : buffer_send(proc_send)%indx(entry_counter(proc_send), 2) = MMM
1353 :
1354 : END DO
1355 :
1356 : END DO
1357 :
1358 2010 : ALLOCATE (req_array(1:para_env%num_pe, 4))
1359 :
1360 134 : CALL timestop(handle2)
1361 :
1362 134 : 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 134 : buffer_send, req_array)
1367 :
1368 539580 : fm_matrix_L%local_data = 0.0_dp
1369 :
1370 134 : CALL timestop(handle2)
1371 :
1372 134 : CALL timeset(routineN//"_9", handle2)
1373 :
1374 : ! fill fm_matrix_L with the entries from buffer_rec
1375 402 : DO iproc = 0, para_env%num_pe - 1
1376 556076 : DO i_entry_rec = 1, num_entries_rec(iproc)
1377 62274974 : DO j_col = 1, ncol_local
1378 62274706 : IF (col_indices(j_col) == buffer_rec(iproc)%indx(i_entry_rec, 2)) THEN
1379 31493602 : DO i_row = 1, nrow_local
1380 31493602 : IF (row_indices(i_row) == buffer_rec(iproc)%indx(i_entry_rec, 1)) THEN
1381 555674 : fm_matrix_L%local_data(i_row, j_col) = buffer_rec(iproc)%msg(i_entry_rec)
1382 : END IF
1383 : END DO
1384 : END IF
1385 : END DO
1386 : END DO
1387 : END DO
1388 :
1389 134 : CALL timestop(handle2)
1390 :
1391 134 : CALL timeset(routineN//"_10", handle2)
1392 :
1393 402 : DO iproc = 0, para_env%num_pe - 1
1394 268 : DEALLOCATE (buffer_rec(iproc)%msg)
1395 268 : DEALLOCATE (buffer_rec(iproc)%indx)
1396 268 : DEALLOCATE (buffer_send(iproc)%msg)
1397 402 : DEALLOCATE (buffer_send(iproc)%indx)
1398 : END DO
1399 :
1400 670 : DEALLOCATE (buffer_rec, buffer_send)
1401 134 : DEALLOCATE (req_array)
1402 134 : DEALLOCATE (entry_counter)
1403 134 : DEALLOCATE (num_entries_rec, num_entries_send)
1404 :
1405 134 : CALL timestop(handle2)
1406 :
1407 134 : CALL timestop(handle)
1408 :
1409 1340 : END SUBROUTINE fill_fm_L_from_L_loc_non_blocking
1410 :
1411 : ! **************************************************************************************************
1412 : !> \brief ...
1413 : !> \param fm_matrix_L ...
1414 : !> \param dimen_RI ...
1415 : !> \param do_inversion ...
1416 : ! **************************************************************************************************
1417 1544 : SUBROUTINE cholesky_decomp(fm_matrix_L, dimen_RI, do_inversion)
1418 :
1419 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1420 : INTEGER, INTENT(IN) :: dimen_RI
1421 : LOGICAL, INTENT(IN) :: do_inversion
1422 :
1423 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cholesky_decomp'
1424 :
1425 : INTEGER :: handle, i_global, iiB, info_chol, &
1426 : j_global, jjB, ncol_local, nrow_local
1427 772 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1428 :
1429 772 : CALL timeset(routineN, handle)
1430 :
1431 : ! do cholesky decomposition
1432 772 : CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=dimen_RI, info_out=info_chol)
1433 772 : CPASSERT(info_chol == 0)
1434 :
1435 772 : IF (.NOT. do_inversion) THEN
1436 : ! clean the lower part of the L^{-1} matrix (just to not have surprises afterwards)
1437 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1438 : nrow_local=nrow_local, &
1439 : ncol_local=ncol_local, &
1440 : row_indices=row_indices, &
1441 118 : col_indices=col_indices)
1442 10330 : DO jjB = 1, ncol_local
1443 10212 : j_global = col_indices(jjB)
1444 549746 : DO iiB = 1, nrow_local
1445 539416 : i_global = row_indices(iiB)
1446 549628 : IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
1447 : END DO
1448 : END DO
1449 :
1450 : END IF
1451 :
1452 772 : CALL timestop(handle)
1453 :
1454 772 : END SUBROUTINE cholesky_decomp
1455 :
1456 : ! **************************************************************************************************
1457 : !> \brief ...
1458 : !> \param fm_matrix_Minv_L_kpoints ...
1459 : !> \param fm_matrix_L_kpoints ...
1460 : !> \param dimen_RI ...
1461 : !> \param kpoints ...
1462 : !> \param eps_eigval_S ...
1463 : ! **************************************************************************************************
1464 20 : SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V(fm_matrix_Minv_L_kpoints, fm_matrix_L_kpoints, &
1465 : dimen_RI, kpoints, eps_eigval_S)
1466 :
1467 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: fm_matrix_Minv_L_kpoints, &
1468 : fm_matrix_L_kpoints
1469 : INTEGER, INTENT(IN) :: dimen_RI
1470 : TYPE(kpoint_type), POINTER :: kpoints
1471 : REAL(KIND=dp), INTENT(IN) :: eps_eigval_S
1472 :
1473 : CHARACTER(LEN=*), PARAMETER :: routineN = 'inversion_of_M_and_mult_with_chol_dec_of_V'
1474 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
1475 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp), ione = CMPLX(0.0_dp, 1.0_dp, KIND=dp)
1476 :
1477 : INTEGER :: handle, ikp, nkp
1478 : TYPE(cp_cfm_type) :: cfm_matrix_K_tmp, cfm_matrix_M_tmp, &
1479 : cfm_matrix_V_tmp, cfm_matrix_Vtrunc_tmp
1480 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
1481 :
1482 20 : CALL timeset(routineN, handle)
1483 :
1484 20 : CALL cp_fm_get_info(fm_matrix_Minv_L_kpoints(1, 1), matrix_struct=matrix_struct)
1485 :
1486 20 : CALL cp_cfm_create(cfm_matrix_M_tmp, matrix_struct)
1487 20 : CALL cp_cfm_create(cfm_matrix_V_tmp, matrix_struct)
1488 20 : CALL cp_cfm_create(cfm_matrix_K_tmp, matrix_struct)
1489 20 : CALL cp_cfm_create(cfm_matrix_Vtrunc_tmp, matrix_struct)
1490 :
1491 20 : CALL get_kpoint_info(kpoints, nkp=nkp)
1492 :
1493 444 : DO ikp = 1, nkp
1494 :
1495 424 : CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_M_tmp, cone, fm_matrix_Minv_L_kpoints(ikp, 1))
1496 424 : CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_M_tmp, ione, fm_matrix_Minv_L_kpoints(ikp, 2))
1497 :
1498 424 : CALL cp_cfm_scale_and_add_fm(czero, cfm_matrix_V_tmp, cone, fm_matrix_L_kpoints(ikp, 1))
1499 424 : CALL cp_cfm_scale_and_add_fm(cone, cfm_matrix_V_tmp, ione, fm_matrix_L_kpoints(ikp, 2))
1500 :
1501 424 : CALL cp_cfm_power(cfm_matrix_M_tmp, threshold=eps_eigval_S, exponent=-1.0_dp)
1502 :
1503 424 : CALL cp_cfm_power(cfm_matrix_V_tmp, threshold=0.0_dp, exponent=0.5_dp)
1504 :
1505 : ! save L(k) = SQRT(V(k))
1506 424 : CALL cp_cfm_to_fm(cfm_matrix_V_tmp, fm_matrix_L_kpoints(ikp, 1), fm_matrix_L_kpoints(ikp, 2))
1507 :
1508 : ! get K = M^(-1)*L, where L is the Cholesky decomposition of V = L^T*L
1509 : CALL parallel_gemm("N", "C", dimen_RI, dimen_RI, dimen_RI, cone, cfm_matrix_M_tmp, cfm_matrix_V_tmp, &
1510 424 : czero, cfm_matrix_K_tmp)
1511 :
1512 : ! move
1513 444 : CALL cp_cfm_to_fm(cfm_matrix_K_tmp, fm_matrix_Minv_L_kpoints(ikp, 1), fm_matrix_Minv_L_kpoints(ikp, 2))
1514 :
1515 : END DO
1516 :
1517 20 : CALL cp_cfm_release(cfm_matrix_M_tmp)
1518 20 : CALL cp_cfm_release(cfm_matrix_V_tmp)
1519 20 : CALL cp_cfm_release(cfm_matrix_K_tmp)
1520 20 : CALL cp_cfm_release(cfm_matrix_Vtrunc_tmp)
1521 :
1522 20 : CALL timestop(handle)
1523 :
1524 20 : END SUBROUTINE inversion_of_M_and_mult_with_chol_dec_of_V
1525 :
1526 : ! **************************************************************************************************
1527 : !> \brief ...
1528 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
1529 : !> \param fm_matrix_Minv ...
1530 : !> \param qs_env ...
1531 : ! **************************************************************************************************
1532 80 : SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc(fm_matrix_Minv_Vtrunc_Minv, &
1533 20 : fm_matrix_Minv, qs_env)
1534 :
1535 : TYPE(cp_fm_type), DIMENSION(:, :) :: fm_matrix_Minv_Vtrunc_Minv, &
1536 : fm_matrix_Minv
1537 : TYPE(qs_environment_type), POINTER :: qs_env
1538 :
1539 : CHARACTER(LEN=*), PARAMETER :: &
1540 : routineN = 'Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc'
1541 :
1542 : INTEGER :: dimen_RI, handle, ndep
1543 : REAL(KIND=dp) :: eps_eigval_S_Gamma
1544 : TYPE(cp_fm_type) :: fm_matrix_RI_metric_inv_work, fm_work
1545 :
1546 20 : CALL timeset(routineN, handle)
1547 :
1548 20 : CALL cp_fm_create(fm_work, fm_matrix_Minv(1, 1)%matrix_struct)
1549 20 : CALL cp_fm_set_all(fm_work, 0.0_dp)
1550 :
1551 20 : CALL cp_fm_create(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1)%matrix_struct)
1552 20 : CALL cp_fm_set_all(fm_matrix_RI_metric_inv_work, 0.0_dp)
1553 :
1554 20 : CALL cp_fm_get_info(matrix=fm_matrix_Minv(1, 1), nrow_global=dimen_RI)
1555 :
1556 20 : eps_eigval_S_Gamma = qs_env%mp2_env%ri_rpa_im_time%eps_eigval_S_Gamma
1557 :
1558 20 : IF (eps_eigval_S_Gamma > 1.0E-18) THEN
1559 :
1560 : ! remove small eigenvalues
1561 : CALL cp_fm_power(fm_matrix_Minv(1, 1), fm_matrix_RI_metric_inv_work, -0.5_dp, &
1562 0 : eps_eigval_S_Gamma, ndep)
1563 :
1564 : ELSE
1565 :
1566 20 : CALL cholesky_decomp(fm_matrix_Minv(1, 1), dimen_RI, do_inversion=.TRUE.)
1567 :
1568 20 : CALL invert_mat(fm_matrix_Minv(1, 1))
1569 :
1570 : END IF
1571 :
1572 : CALL parallel_gemm('N', 'T', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_Minv(1, 1), &
1573 20 : fm_matrix_Minv(1, 1), 0.0_dp, fm_matrix_RI_metric_inv_work)
1574 :
1575 20 : CALL cp_fm_to_fm(fm_matrix_RI_metric_inv_work, fm_matrix_Minv(1, 1))
1576 :
1577 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_matrix_RI_metric_inv_work, &
1578 20 : fm_matrix_Minv_Vtrunc_Minv(1, 1), 0.0_dp, fm_work)
1579 :
1580 : CALL parallel_gemm('N', 'N', dimen_RI, dimen_RI, dimen_RI, 1.0_dp, fm_work, &
1581 20 : fm_matrix_RI_metric_inv_work, 0.0_dp, fm_matrix_Minv_Vtrunc_Minv(1, 1))
1582 :
1583 20 : CALL cp_fm_release(fm_work)
1584 20 : CALL cp_fm_release(fm_matrix_RI_metric_inv_work)
1585 :
1586 20 : CALL timestop(handle)
1587 :
1588 20 : END SUBROUTINE Gamma_only_inversion_of_M_and_mult_with_chol_dec_of_Vtrunc
1589 :
1590 : ! **************************************************************************************************
1591 : !> \brief ...
1592 : !> \param qs_env ...
1593 : !> \param trunc_coulomb ...
1594 : !> \param rel_cutoff_trunc_coulomb_ri_x ...
1595 : !> \param cell_grid ...
1596 : !> \param do_BvK_cell ...
1597 : ! **************************************************************************************************
1598 48 : SUBROUTINE trunc_coulomb_for_exchange(qs_env, trunc_coulomb, rel_cutoff_trunc_coulomb_ri_x, &
1599 : cell_grid, do_BvK_cell)
1600 : TYPE(qs_environment_type), POINTER :: qs_env
1601 : TYPE(libint_potential_type), OPTIONAL :: trunc_coulomb
1602 : REAL(KIND=dp), OPTIONAL :: rel_cutoff_trunc_coulomb_ri_x
1603 : INTEGER, DIMENSION(3), OPTIONAL :: cell_grid
1604 : LOGICAL, OPTIONAL :: do_BvK_cell
1605 :
1606 : CHARACTER(LEN=*), PARAMETER :: routineN = 'trunc_coulomb_for_exchange'
1607 :
1608 : INTEGER :: handle, i_dim
1609 : INTEGER, DIMENSION(3) :: periodic
1610 : LOGICAL :: my_do_BvK_cell
1611 : REAL(KIND=dp) :: kp_fac, kp_fac_idim, my_rel_cutoff_trunc_coulomb_ri_x, &
1612 : shortest_dist_cell_planes
1613 : TYPE(cell_type), POINTER :: cell
1614 : TYPE(kpoint_type), POINTER :: kpoints_scf
1615 :
1616 48 : CALL timeset(routineN, handle)
1617 :
1618 48 : NULLIFY (cell)
1619 48 : CALL get_qs_env(qs_env, cell=cell, kpoints=kpoints_scf)
1620 48 : CALL get_cell(cell=cell, periodic=periodic)
1621 :
1622 48 : my_do_BvK_cell = .FALSE.
1623 48 : IF (PRESENT(do_BvK_cell)) my_do_BvK_cell = do_BvK_cell
1624 28 : IF (my_do_BvK_cell) THEN
1625 : kp_fac = 1.0E10_dp
1626 24 : DO i_dim = 1, 3
1627 : ! look for smallest k-point mesh in periodic direction
1628 24 : IF (periodic(i_dim) == 1) THEN
1629 12 : IF (PRESENT(cell_grid)) THEN
1630 12 : kp_fac_idim = REAL(cell_grid(i_dim), KIND=dp)
1631 : ELSE
1632 0 : kp_fac_idim = REAL(kpoints_scf%nkp_grid(i_dim), KIND=dp)
1633 : END IF
1634 12 : IF (kp_fac > kp_fac_idim) kp_fac = kp_fac_idim
1635 : END IF
1636 : END DO
1637 : ELSE
1638 : kp_fac = 1.0_dp
1639 : END IF
1640 :
1641 48 : shortest_dist_cell_planes = 1.0E4_dp
1642 48 : IF (periodic(1) == 1) THEN
1643 10 : IF (shortest_dist_cell_planes > plane_distance(1, 0, 0, cell)) THEN
1644 10 : shortest_dist_cell_planes = plane_distance(1, 0, 0, cell)
1645 : END IF
1646 : END IF
1647 48 : IF (periodic(2) == 1) THEN
1648 28 : IF (shortest_dist_cell_planes > plane_distance(0, 1, 0, cell)) THEN
1649 22 : shortest_dist_cell_planes = plane_distance(0, 1, 0, cell)
1650 : END IF
1651 : END IF
1652 48 : IF (periodic(3) == 1) THEN
1653 30 : IF (shortest_dist_cell_planes > plane_distance(0, 0, 1, cell)) THEN
1654 6 : shortest_dist_cell_planes = plane_distance(0, 0, 1, cell)
1655 : END IF
1656 : END IF
1657 :
1658 48 : IF (PRESENT(rel_cutoff_trunc_coulomb_ri_x)) THEN
1659 28 : my_rel_cutoff_trunc_coulomb_ri_x = rel_cutoff_trunc_coulomb_ri_x
1660 : ELSE
1661 20 : my_rel_cutoff_trunc_coulomb_ri_x = qs_env%mp2_env%ri_rpa_im_time%rel_cutoff_trunc_coulomb_ri_x
1662 : END IF
1663 :
1664 48 : IF (PRESENT(trunc_coulomb)) THEN
1665 48 : trunc_coulomb%potential_type = do_potential_truncated
1666 : trunc_coulomb%cutoff_radius = shortest_dist_cell_planes* &
1667 : my_rel_cutoff_trunc_coulomb_ri_x* &
1668 48 : kp_fac
1669 48 : trunc_coulomb%filename = "t_c_g.dat"
1670 : ! dummy
1671 48 : trunc_coulomb%omega = 0.0_dp
1672 : END IF
1673 :
1674 48 : CALL timestop(handle)
1675 :
1676 48 : END SUBROUTINE trunc_coulomb_for_exchange
1677 :
1678 : ! **************************************************************************************************
1679 : !> \brief ...
1680 : !> \param fm_matrix_L ...
1681 : ! **************************************************************************************************
1682 1308 : SUBROUTINE invert_mat(fm_matrix_L)
1683 :
1684 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1685 :
1686 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invert_mat'
1687 :
1688 : INTEGER :: handle, i_global, iiB, j_global, jjB, &
1689 : ncol_local, nrow_local
1690 654 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1691 :
1692 654 : CALL timeset(routineN, handle)
1693 :
1694 654 : CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U')
1695 :
1696 : ! clear the lower part of the L^{-1} matrix (just to have no surprises afterwards)
1697 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1698 : nrow_local=nrow_local, &
1699 : ncol_local=ncol_local, &
1700 : row_indices=row_indices, &
1701 654 : col_indices=col_indices)
1702 : ! assume upper triangular format
1703 50394 : DO jjB = 1, ncol_local
1704 49740 : j_global = col_indices(jjB)
1705 3602021 : DO iiB = 1, nrow_local
1706 3551627 : i_global = row_indices(iiB)
1707 3601367 : IF (j_global < i_global) fm_matrix_L%local_data(iiB, jjB) = 0.0_dp
1708 : END DO
1709 : END DO
1710 :
1711 654 : CALL timestop(handle)
1712 :
1713 654 : END SUBROUTINE invert_mat
1714 :
1715 : ! **************************************************************************************************
1716 : !> \brief ...
1717 : !> \param fm_matrix_L ...
1718 : !> \param blacs_env_L ...
1719 : !> \param dimen_RI ...
1720 : !> \param para_env_L ...
1721 : !> \param name ...
1722 : !> \param fm_matrix_L_extern ...
1723 : ! **************************************************************************************************
1724 690 : SUBROUTINE create_matrix_L(fm_matrix_L, blacs_env_L, dimen_RI, para_env_L, name, fm_matrix_L_extern)
1725 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_L
1726 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_L
1727 : INTEGER, INTENT(IN) :: dimen_RI
1728 : TYPE(mp_para_env_type), POINTER :: para_env_L
1729 : CHARACTER(LEN=*), INTENT(IN) :: name
1730 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: fm_matrix_L_extern
1731 :
1732 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_matrix_L'
1733 :
1734 : INTEGER :: handle
1735 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1736 :
1737 690 : CALL timeset(routineN, handle)
1738 :
1739 : ! create the full matrix L defined in the L group
1740 690 : IF (PRESENT(fm_matrix_L_extern)) THEN
1741 36 : CALL cp_fm_create(fm_matrix_L, fm_matrix_L_extern%matrix_struct, name=name)
1742 : ELSE
1743 654 : NULLIFY (fm_struct)
1744 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_L, nrow_global=dimen_RI, &
1745 654 : ncol_global=dimen_RI, para_env=para_env_L)
1746 :
1747 654 : CALL cp_fm_create(fm_matrix_L, fm_struct, name=name)
1748 :
1749 654 : CALL cp_fm_struct_release(fm_struct)
1750 : END IF
1751 :
1752 690 : CALL cp_fm_set_all(matrix=fm_matrix_L, alpha=0.0_dp)
1753 :
1754 690 : CALL timestop(handle)
1755 :
1756 690 : END SUBROUTINE create_matrix_L
1757 :
1758 : ! **************************************************************************************************
1759 : !> \brief ...
1760 : !> \param fm_matrix_L ...
1761 : !> \param my_group_L_start ...
1762 : !> \param my_group_L_end ...
1763 : !> \param my_group_L_size ...
1764 : !> \param my_Lrows ...
1765 : ! **************************************************************************************************
1766 532 : SUBROUTINE grep_Lcols(fm_matrix_L, &
1767 : my_group_L_start, my_group_L_end, my_group_L_size, my_Lrows)
1768 : TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_L
1769 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end, &
1770 : my_group_L_size
1771 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
1772 : INTENT(OUT) :: my_Lrows
1773 :
1774 : CHARACTER(LEN=*), PARAMETER :: routineN = 'grep_Lcols'
1775 :
1776 : INTEGER :: dimen_RI, handle, handle2, i_global, iiB, j_global, jjB, max_row_col_local, &
1777 : ncol_local, ncol_rec, nrow_local, nrow_rec, proc_receive_static, proc_send_static, &
1778 : proc_shift
1779 532 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, rec_col_row_info
1780 532 : INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_rec, &
1781 532 : row_indices, row_indices_rec
1782 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: local_L, rec_L
1783 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
1784 532 : POINTER :: local_L_internal
1785 : TYPE(mp_para_env_type), POINTER :: para_env
1786 :
1787 532 : CALL timeset(routineN, handle)
1788 :
1789 : CALL cp_fm_get_info(matrix=fm_matrix_L, &
1790 : nrow_local=nrow_local, &
1791 : ncol_local=ncol_local, &
1792 : row_indices=row_indices, &
1793 : col_indices=col_indices, &
1794 : nrow_global=dimen_RI, &
1795 : local_data=local_L_internal, &
1796 532 : para_env=para_env)
1797 :
1798 2128 : ALLOCATE (my_Lrows(dimen_RI, my_group_L_size))
1799 1702540 : my_Lrows = 0.0_dp
1800 :
1801 2128 : ALLOCATE (local_L(nrow_local, ncol_local))
1802 3135544 : local_L(:, :) = local_L_internal(1:nrow_local, 1:ncol_local)
1803 :
1804 532 : max_row_col_local = MAX(nrow_local, ncol_local)
1805 532 : CALL para_env%max(max_row_col_local)
1806 :
1807 2128 : ALLOCATE (local_col_row_info(0:max_row_col_local, 2))
1808 81408 : local_col_row_info = 0
1809 : ! 0,1 nrows
1810 532 : local_col_row_info(0, 1) = nrow_local
1811 38903 : local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
1812 : ! 0,2 ncols
1813 532 : local_col_row_info(0, 2) = ncol_local
1814 39796 : local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
1815 :
1816 1064 : ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
1817 :
1818 : ! accumulate data on my_Lrows starting from myself
1819 39796 : DO jjB = 1, ncol_local
1820 39264 : j_global = col_indices(jjB)
1821 39796 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
1822 1613020 : DO iiB = 1, nrow_local
1823 1592368 : i_global = row_indices(iiB)
1824 1613020 : my_Lrows(i_global, j_global - my_group_L_start + 1) = local_L(iiB, jjB)
1825 : END DO
1826 : END IF
1827 : END DO
1828 :
1829 532 : proc_send_static = MODULO(para_env%mepos + 1, para_env%num_pe)
1830 532 : proc_receive_static = MODULO(para_env%mepos - 1, para_env%num_pe)
1831 :
1832 532 : CALL timeset(routineN//"_comm", handle2)
1833 :
1834 556 : DO proc_shift = 1, para_env%num_pe - 1
1835 : ! first exchange information on the local data
1836 4200 : rec_col_row_info = 0
1837 24 : CALL para_env%sendrecv(local_col_row_info, proc_send_static, rec_col_row_info, proc_receive_static)
1838 24 : nrow_rec = rec_col_row_info(0, 1)
1839 24 : ncol_rec = rec_col_row_info(0, 2)
1840 :
1841 72 : ALLOCATE (row_indices_rec(nrow_rec))
1842 1061 : row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1843 :
1844 72 : ALLOCATE (col_indices_rec(ncol_rec))
1845 2064 : col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1846 :
1847 96 : ALLOCATE (rec_L(nrow_rec, ncol_rec))
1848 91052 : rec_L = 0.0_dp
1849 :
1850 : ! then send and receive the real data
1851 24 : CALL para_env%sendrecv(local_L, proc_send_static, rec_L, proc_receive_static)
1852 :
1853 : ! accumulate the received data on my_Lrows
1854 2064 : DO jjB = 1, ncol_rec
1855 2040 : j_global = col_indices_rec(jjB)
1856 2064 : IF (j_global >= my_group_L_start .AND. j_global <= my_group_L_end) THEN
1857 91028 : DO iiB = 1, nrow_rec
1858 88988 : i_global = row_indices_rec(iiB)
1859 91028 : my_Lrows(i_global, j_global - my_group_L_start + 1) = rec_L(iiB, jjB)
1860 : END DO
1861 : END IF
1862 : END DO
1863 :
1864 4200 : local_col_row_info(:, :) = rec_col_row_info
1865 24 : CALL MOVE_ALLOC(rec_L, local_L)
1866 :
1867 24 : DEALLOCATE (col_indices_rec)
1868 556 : DEALLOCATE (row_indices_rec)
1869 : END DO
1870 532 : CALL timestop(handle2)
1871 :
1872 532 : DEALLOCATE (local_col_row_info)
1873 532 : DEALLOCATE (rec_col_row_info)
1874 532 : DEALLOCATE (local_L)
1875 :
1876 532 : CALL timestop(handle)
1877 :
1878 1596 : END SUBROUTINE grep_Lcols
1879 :
1880 : END MODULE mp2_ri_2c
|
