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