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 Routines to calculate and distribute 2c- and 3c- 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_integrals
15 : USE OMP_LIB, ONLY: omp_get_num_threads,&
16 : omp_get_thread_num
17 : USE atomic_kind_types, ONLY: atomic_kind_type
18 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
19 : gto_basis_set_type
20 : USE bibliography, ONLY: DelBen2013,&
21 : cite_reference
22 : USE cell_types, ONLY: cell_type,&
23 : get_cell
24 : USE cp_blacs_env, ONLY: cp_blacs_env_type
25 : USE cp_control_types, ONLY: dft_control_type
26 : USE cp_dbcsr_api, ONLY: &
27 : dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
28 : dbcsr_release_p, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
29 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
30 : cp_dbcsr_m_by_n_from_template
31 : USE cp_eri_mme_interface, ONLY: cp_eri_mme_param,&
32 : cp_eri_mme_set_params
33 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 : cp_fm_struct_release,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_get_info,&
38 : cp_fm_release,&
39 : cp_fm_type
40 : USE cp_log_handling, ONLY: cp_to_string
41 : USE cp_units, ONLY: cp_unit_from_cp2k
42 : USE dbt_api, ONLY: &
43 : dbt_clear, dbt_contract, dbt_copy, dbt_create, dbt_destroy, dbt_distribution_destroy, &
44 : dbt_distribution_new, dbt_distribution_type, dbt_filter, dbt_get_block, dbt_get_info, &
45 : dbt_get_stored_coordinates, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
46 : dbt_pgrid_type, dbt_put_block, dbt_reserve_blocks, dbt_scale, dbt_split_blocks, dbt_type
47 : USE group_dist_types, ONLY: create_group_dist,&
48 : get_group_dist,&
49 : group_dist_d1_type
50 : USE hfx_types, ONLY: alloc_containers,&
51 : block_ind_type,&
52 : hfx_compression_type
53 : USE input_constants, ONLY: &
54 : do_eri_gpw, do_eri_mme, do_eri_os, do_potential_coulomb, do_potential_id, &
55 : do_potential_long, do_potential_short, do_potential_truncated, kp_weights_W_auto, &
56 : kp_weights_W_tailored, kp_weights_W_uniform
57 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
58 : section_vals_type,&
59 : section_vals_val_get
60 : USE kinds, ONLY: default_string_length,&
61 : dp,&
62 : int_8
63 : USE kpoint_methods, ONLY: kpoint_init_cell_index
64 : USE kpoint_types, ONLY: kpoint_type
65 : USE libint_2c_3c, ONLY: compare_potential_types,&
66 : libint_potential_type
67 : USE machine, ONLY: m_flush
68 : USE message_passing, ONLY: mp_cart_type,&
69 : mp_comm_type,&
70 : mp_para_env_type
71 : USE mp2_eri, ONLY: mp2_eri_3c_integrate
72 : USE mp2_eri_gpw, ONLY: cleanup_gpw,&
73 : mp2_eri_3c_integrate_gpw,&
74 : prepare_gpw
75 : USE mp2_ri_2c, ONLY: get_2c_integrals
76 : USE mp2_types, ONLY: three_dim_real_array
77 : USE particle_methods, ONLY: get_particle_set
78 : USE particle_types, ONLY: particle_type
79 : USE pw_env_types, ONLY: pw_env_type
80 : USE pw_poisson_types, ONLY: pw_poisson_type
81 : USE pw_pool_types, ONLY: pw_pool_type
82 : USE pw_types, ONLY: pw_c1d_gs_type,&
83 : pw_r3d_rs_type
84 : USE qs_environment_types, ONLY: get_qs_env,&
85 : qs_environment_type,&
86 : set_qs_env
87 : USE qs_integral_utils, ONLY: basis_set_list_setup
88 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
89 : USE qs_kind_types, ONLY: qs_kind_type
90 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
91 : USE qs_tensors, ONLY: build_3c_integrals,&
92 : build_3c_neighbor_lists,&
93 : compress_tensor,&
94 : get_tensor_occupancy,&
95 : neighbor_list_3c_destroy
96 : USE qs_tensors_types, ONLY: create_3c_tensor,&
97 : create_tensor_batches,&
98 : distribution_3d_create,&
99 : distribution_3d_type,&
100 : neighbor_list_3c_type,&
101 : pgf_block_sizes
102 : USE task_list_types, ONLY: task_list_type
103 : USE util, ONLY: get_limit
104 : #include "./base/base_uses.f90"
105 :
106 : IMPLICIT NONE
107 :
108 : PRIVATE
109 :
110 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_integrals'
111 :
112 : PUBLIC :: mp2_ri_gpw_compute_in, compute_kpoints
113 :
114 : TYPE intermediate_matrix_type
115 : TYPE(dbcsr_type) :: matrix_ia_jnu, matrix_ia_jb
116 : INTEGER :: max_row_col_local = 0
117 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info
118 : TYPE(cp_fm_type) :: fm_BIb_jb = cp_fm_type()
119 : CHARACTER(LEN=default_string_length) :: descr = ""
120 : END TYPE intermediate_matrix_type
121 :
122 : CONTAINS
123 :
124 : ! **************************************************************************************************
125 : !> \brief with ri mp2 gpw
126 : !> \param BIb_C ...
127 : !> \param BIb_C_gw ...
128 : !> \param BIb_C_bse_ij ...
129 : !> \param BIb_C_bse_ab ...
130 : !> \param gd_array ...
131 : !> \param gd_B_virtual ...
132 : !> \param dimen_RI ...
133 : !> \param dimen_RI_red ...
134 : !> \param qs_env ...
135 : !> \param para_env ...
136 : !> \param para_env_sub ...
137 : !> \param color_sub ...
138 : !> \param cell ...
139 : !> \param particle_set ...
140 : !> \param atomic_kind_set ...
141 : !> \param qs_kind_set ...
142 : !> \param fm_matrix_PQ ...
143 : !> \param fm_matrix_L_kpoints ...
144 : !> \param fm_matrix_Minv_L_kpoints ...
145 : !> \param fm_matrix_Minv ...
146 : !> \param fm_matrix_Minv_Vtrunc_Minv ...
147 : !> \param nmo ...
148 : !> \param homo ...
149 : !> \param mat_munu ...
150 : !> \param sab_orb_sub ...
151 : !> \param mo_coeff_o ...
152 : !> \param mo_coeff_v ...
153 : !> \param mo_coeff_all ...
154 : !> \param mo_coeff_gw ...
155 : !> \param mo_coeff_o_bse ...
156 : !> \param mo_coeff_v_bse ...
157 : !> \param eps_filter ...
158 : !> \param unit_nr ...
159 : !> \param mp2_memory ...
160 : !> \param calc_PQ_cond_num ...
161 : !> \param calc_forces ...
162 : !> \param blacs_env_sub ...
163 : !> \param my_do_gw ...
164 : !> \param do_bse ...
165 : !> \param gd_B_all ...
166 : !> \param starts_array_mc ...
167 : !> \param ends_array_mc ...
168 : !> \param starts_array_mc_block ...
169 : !> \param ends_array_mc_block ...
170 : !> \param gw_corr_lev_occ ...
171 : !> \param gw_corr_lev_virt ...
172 : !> \param bse_lev_virt ...
173 : !> \param do_im_time ...
174 : !> \param do_kpoints_cubic_RPA ...
175 : !> \param kpoints ...
176 : !> \param t_3c_M ...
177 : !> \param t_3c_O ...
178 : !> \param t_3c_O_compressed ...
179 : !> \param t_3c_O_ind ...
180 : !> \param ri_metric ...
181 : !> \param gd_B_occ_bse ...
182 : !> \param gd_B_virt_bse ...
183 : !> \author Mauro Del Ben
184 : ! **************************************************************************************************
185 5344 : SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd_array, gd_B_virtual, &
186 : dimen_RI, dimen_RI_red, qs_env, para_env, para_env_sub, color_sub, &
187 : cell, particle_set, atomic_kind_set, qs_kind_set, &
188 : fm_matrix_PQ, fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, &
189 : fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
190 668 : nmo, homo, mat_munu, &
191 668 : sab_orb_sub, mo_coeff_o, mo_coeff_v, mo_coeff_all, &
192 668 : mo_coeff_gw, mo_coeff_o_bse, mo_coeff_v_bse, eps_filter, unit_nr, &
193 : mp2_memory, calc_PQ_cond_num, calc_forces, blacs_env_sub, my_do_gw, do_bse, &
194 : gd_B_all, starts_array_mc, ends_array_mc, &
195 : starts_array_mc_block, ends_array_mc_block, &
196 : gw_corr_lev_occ, gw_corr_lev_virt, &
197 668 : bse_lev_virt, &
198 : do_im_time, do_kpoints_cubic_RPA, kpoints, &
199 : t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
200 : ri_metric, gd_B_occ_bse, gd_B_virt_bse)
201 :
202 : TYPE(three_dim_real_array), ALLOCATABLE, &
203 : DIMENSION(:), INTENT(OUT) :: BIb_C, BIb_C_gw, BIb_C_bse_ij, &
204 : BIb_C_bse_ab
205 : TYPE(group_dist_d1_type), INTENT(OUT) :: gd_array
206 : TYPE(group_dist_d1_type), ALLOCATABLE, &
207 : DIMENSION(:), INTENT(OUT) :: gd_B_virtual
208 : INTEGER, INTENT(OUT) :: dimen_RI, dimen_RI_red
209 : TYPE(qs_environment_type), POINTER :: qs_env
210 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
211 : INTEGER, INTENT(IN) :: color_sub
212 : TYPE(cell_type), POINTER :: cell
213 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
214 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
215 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
216 : TYPE(cp_fm_type), INTENT(OUT) :: fm_matrix_PQ
217 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_matrix_L_kpoints, &
218 : fm_matrix_Minv_L_kpoints, &
219 : fm_matrix_Minv, &
220 : fm_matrix_Minv_Vtrunc_Minv
221 : INTEGER, INTENT(IN) :: nmo
222 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
223 : TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu
224 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
225 : INTENT(IN), POINTER :: sab_orb_sub
226 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mo_coeff_o, mo_coeff_v, mo_coeff_all, &
227 : mo_coeff_gw, mo_coeff_o_bse, &
228 : mo_coeff_v_bse
229 : REAL(KIND=dp), INTENT(IN) :: eps_filter
230 : INTEGER, INTENT(IN) :: unit_nr
231 : REAL(KIND=dp), INTENT(IN) :: mp2_memory
232 : LOGICAL, INTENT(IN) :: calc_PQ_cond_num, calc_forces
233 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
234 : LOGICAL, INTENT(IN) :: my_do_gw, do_bse
235 : TYPE(group_dist_d1_type), INTENT(OUT) :: gd_B_all
236 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: starts_array_mc, ends_array_mc, &
237 : starts_array_mc_block, &
238 : ends_array_mc_block
239 : INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
240 : INTEGER, DIMENSION(:), INTENT(IN) :: bse_lev_virt
241 : LOGICAL, INTENT(IN) :: do_im_time, do_kpoints_cubic_RPA
242 : TYPE(kpoint_type), POINTER :: kpoints
243 : TYPE(dbt_type), INTENT(OUT) :: t_3c_M
244 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :), &
245 : INTENT(OUT) :: t_3c_O
246 : TYPE(hfx_compression_type), ALLOCATABLE, &
247 : DIMENSION(:, :, :), INTENT(INOUT) :: t_3c_O_compressed
248 : TYPE(block_ind_type), ALLOCATABLE, &
249 : DIMENSION(:, :, :) :: t_3c_O_ind
250 : TYPE(libint_potential_type), INTENT(IN) :: ri_metric
251 : TYPE(group_dist_d1_type), ALLOCATABLE, &
252 : DIMENSION(:), INTENT(OUT) :: gd_B_occ_bse, gd_B_virt_bse
253 :
254 : CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_ri_gpw_compute_in'
255 :
256 : INTEGER :: cm, cut_memory, cut_memory_int, eri_method, gw_corr_lev_total, handle, handle2, &
257 : handle4, i, i_counter, i_mem, ibasis, ispin, itmp(2), j, jcell, kcell, LLL, min_bsize, &
258 : my_B_all_end, my_B_all_size, my_B_all_start, my_group_L_end, my_group_L_size, &
259 : my_group_L_start, n_rep, natom, ngroup, nimg, nkind, nspins, potential_type, &
260 : ri_metric_type
261 : INTEGER(int_8) :: nze
262 668 : INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI, &
263 668 : ends_array_mc_block_int, ends_array_mc_int, my_B_occ_bse_end, my_B_occ_bse_size, &
264 668 : my_B_occ_bse_start, my_B_size, my_B_virt_bse_end, my_B_virt_bse_size, &
265 1336 : my_B_virt_bse_start, my_B_virtual_end, my_B_virtual_start, sizes_AO, sizes_AO_split, &
266 1336 : sizes_RI, sizes_RI_split, starts_array_mc_block_int, starts_array_mc_int, virtual
267 : INTEGER, DIMENSION(2, 3) :: bounds
268 : INTEGER, DIMENSION(3) :: bounds_3c, pcoord, pdims, pdims_t3c, &
269 : periodic
270 : LOGICAL :: do_gpw, do_kpoints_from_Gamma, do_svd, &
271 : memory_info
272 : REAL(KIND=dp) :: compression_factor, cutoff_old, eps_pgf_orb, eps_pgf_orb_old, eps_svd, &
273 : mem_for_abK, mem_for_iaK, mem_for_ijK, memory_3c, occ, omega_pot, rc_ang, &
274 : relative_cutoff_old
275 668 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old
276 668 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: my_Lrows, my_Vrows
277 : TYPE(cp_eri_mme_param), POINTER :: eri_param
278 668 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_munu_local_L
279 3340 : TYPE(dbt_pgrid_type) :: pgrid_t3c_M, pgrid_t3c_overl
280 8684 : TYPE(dbt_type) :: t_3c_overl_int_template, t_3c_tmp
281 668 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_overl_int
282 : TYPE(dft_control_type), POINTER :: dft_control
283 : TYPE(distribution_3d_type) :: dist_3d
284 668 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_ao, basis_set_ri_aux
285 : TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
286 : TYPE(intermediate_matrix_type), ALLOCATABLE, &
287 668 : DIMENSION(:) :: intermed_mat, intermed_mat_bse_ab, &
288 668 : intermed_mat_bse_ij, intermed_mat_gw
289 668 : TYPE(mp_cart_type) :: mp_comm_t3c_2
290 : TYPE(neighbor_list_3c_type) :: nl_3c
291 : TYPE(pw_c1d_gs_type) :: pot_g, rho_g
292 : TYPE(pw_env_type), POINTER :: pw_env_sub
293 : TYPE(pw_poisson_type), POINTER :: poisson_env
294 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
295 : TYPE(pw_r3d_rs_type) :: psi_L, rho_r
296 : TYPE(section_vals_type), POINTER :: qs_section
297 : TYPE(task_list_type), POINTER :: task_list_sub
298 :
299 668 : CALL timeset(routineN, handle)
300 :
301 668 : CALL cite_reference(DelBen2013)
302 :
303 668 : nspins = SIZE(homo)
304 :
305 2004 : ALLOCATE (virtual(nspins))
306 1490 : virtual(:) = nmo - homo(:)
307 668 : gw_corr_lev_total = gw_corr_lev_virt + gw_corr_lev_occ
308 :
309 668 : eri_method = qs_env%mp2_env%eri_method
310 668 : eri_param => qs_env%mp2_env%eri_mme_param
311 668 : do_svd = qs_env%mp2_env%do_svd
312 668 : eps_svd = qs_env%mp2_env%eps_svd
313 668 : potential_type = qs_env%mp2_env%potential_parameter%potential_type
314 668 : ri_metric_type = ri_metric%potential_type
315 668 : omega_pot = qs_env%mp2_env%potential_parameter%omega
316 :
317 : ! whether we need gpw integrals (plus pw stuff)
318 : do_gpw = (eri_method == do_eri_gpw) .OR. &
319 : ((potential_type == do_potential_long .OR. ri_metric_type == do_potential_long) &
320 : .AND. qs_env%mp2_env%eri_method == do_eri_os) &
321 668 : .OR. (ri_metric_type == do_potential_id .AND. qs_env%mp2_env%eri_method == do_eri_mme)
322 :
323 668 : IF (do_svd .AND. calc_forces) THEN
324 0 : CPABORT("SVD not implemented for forces.!")
325 : END IF
326 :
327 668 : do_kpoints_from_Gamma = qs_env%mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma
328 668 : IF (do_kpoints_cubic_RPA .OR. do_kpoints_from_Gamma) THEN
329 : CALL get_qs_env(qs_env=qs_env, &
330 22 : kpoints=kpoints)
331 : END IF
332 22 : IF (do_kpoints_from_Gamma) THEN
333 16 : CALL compute_kpoints(qs_env, kpoints, unit_nr)
334 : END IF
335 :
336 668 : IF (do_bse) THEN
337 42 : IF (.NOT. my_do_gw) THEN
338 0 : CALL cp_abort(__LOCATION__, "BSE calculations require prior GW calculations.")
339 : END IF
340 42 : IF (do_im_time) THEN
341 0 : CALL cp_abort(__LOCATION__, "BSE calculations are not implemented for low-scaling GW.")
342 : END IF
343 : ! GPW integrals have to be implemented later
344 42 : IF (eri_method == do_eri_gpw) THEN
345 : CALL cp_abort(__LOCATION__, &
346 : "BSE calculations are not implemented for GPW integrals. "// &
347 : "This is probably caused by invoking a periodic calculation. "// &
348 0 : "Use PERIODIC NONE for BSE calculations.")
349 : END IF
350 : END IF
351 :
352 668 : ngroup = para_env%num_pe/para_env_sub%num_pe
353 :
354 : ! Preparations for MME method to compute ERIs
355 668 : IF (qs_env%mp2_env%eri_method == do_eri_mme) THEN
356 : ! cell might have changed, so we need to reset parameters
357 126 : CALL cp_eri_mme_set_params(eri_param, cell, qs_kind_set, basis_type_1="ORB", basis_type_2="RI_AUX", para_env=para_env)
358 : END IF
359 :
360 668 : CALL get_cell(cell=cell, periodic=periodic)
361 : ! for minimax Ewald summation, full periodicity is required
362 668 : IF (eri_method == do_eri_mme) THEN
363 126 : CPASSERT(periodic(1) == 1 .AND. periodic(2) == 1 .AND. periodic(3) == 1)
364 : END IF
365 :
366 668 : IF (do_svd .AND. (do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA)) THEN
367 0 : CPABORT("SVD with kpoints not implemented yet!")
368 : END IF
369 :
370 : CALL get_2c_integrals(qs_env, eri_method, eri_param, para_env, para_env_sub, mp2_memory, &
371 : my_Lrows, my_Vrows, fm_matrix_PQ, ngroup, color_sub, dimen_RI, dimen_RI_red, &
372 : kpoints, my_group_L_size, my_group_L_start, my_group_L_end, &
373 : gd_array, calc_PQ_cond_num .AND. .NOT. do_svd, do_svd, eps_svd, &
374 : qs_env%mp2_env%potential_parameter, ri_metric, &
375 : fm_matrix_L_kpoints, fm_matrix_Minv_L_kpoints, fm_matrix_Minv, fm_matrix_Minv_Vtrunc_Minv, &
376 : do_im_time, do_kpoints_from_Gamma .OR. do_kpoints_cubic_RPA, qs_env%mp2_env%mp2_gpw%eps_pgf_orb_S, &
377 1982 : qs_kind_set, sab_orb_sub, calc_forces, unit_nr)
378 :
379 668 : IF (unit_nr > 0) THEN
380 : ASSOCIATE (ri_metric => qs_env%mp2_env%ri_metric)
381 588 : SELECT CASE (ri_metric%potential_type)
382 : CASE (do_potential_coulomb)
383 : WRITE (unit_nr, FMT="(/T3,A,T74,A)") &
384 254 : "RI_INFO| RI metric: ", "COULOMB"
385 : CASE (do_potential_short)
386 : WRITE (unit_nr, FMT="(T3,A,T71,A)") &
387 0 : "RI_INFO| RI metric: ", "SHORTRANGE"
388 : WRITE (unit_nr, '(T3,A,T61,F20.10)') &
389 0 : "RI_INFO| Omega: ", ri_metric%omega
390 0 : rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
391 : WRITE (unit_nr, '(T3,A,T61,F20.10)') &
392 0 : "RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
393 : CASE (do_potential_long)
394 : WRITE (unit_nr, FMT="(T3,A,T72,A)") &
395 8 : "RI_INFO| RI metric: ", "LONGRANGE"
396 : WRITE (unit_nr, '(T3,A,T61,F20.10)') &
397 8 : "RI_INFO| Omega: ", ri_metric%omega
398 : CASE (do_potential_id)
399 : WRITE (unit_nr, FMT="(T3,A,T74,A)") &
400 41 : "RI_INFO| RI metric: ", "OVERLAP"
401 : CASE (do_potential_truncated)
402 : WRITE (unit_nr, FMT="(T3,A,T64,A)") &
403 31 : "RI_INFO| RI metric: ", "TRUNCATED COULOMB"
404 31 : rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
405 : WRITE (unit_nr, '(T3,A,T61,F20.2)') &
406 365 : "RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
407 : END SELECT
408 : END ASSOCIATE
409 : END IF
410 :
411 668 : IF (calc_forces .AND. .NOT. do_im_time) THEN
412 : ! we need (P|Q)^(-1/2) for future use, just save it
413 : ! in a fully (home made) distributed way
414 272 : itmp = get_limit(dimen_RI, para_env_sub%num_pe, para_env_sub%mepos)
415 272 : lll = itmp(2) - itmp(1) + 1
416 1088 : ALLOCATE (qs_env%mp2_env%ri_grad%PQ_half(lll, my_group_L_size))
417 972108 : qs_env%mp2_env%ri_grad%PQ_half(:, :) = my_Lrows(itmp(1):itmp(2), 1:my_group_L_size)
418 272 : IF (.NOT. compare_potential_types(qs_env%mp2_env%ri_metric, qs_env%mp2_env%potential_parameter)) THEN
419 36 : ALLOCATE (qs_env%mp2_env%ri_grad%operator_half(lll, my_group_L_size))
420 41844 : qs_env%mp2_env%ri_grad%operator_half(:, :) = my_Vrows(itmp(1):itmp(2), 1:my_group_L_size)
421 12 : DEALLOCATE (my_Vrows)
422 : END IF
423 : END IF
424 :
425 668 : IF (unit_nr > 0) THEN
426 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
427 334 : "RI_INFO| Number of auxiliary basis functions:", dimen_RI, &
428 334 : "GENERAL_INFO| Number of basis functions:", nmo, &
429 334 : "GENERAL_INFO| Number of occupied orbitals:", homo(1), &
430 668 : "GENERAL_INFO| Number of virtual orbitals:", virtual(1)
431 334 : IF (do_svd) THEN
432 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
433 22 : "RI_INFO| Reduced auxiliary basis set size:", dimen_RI_red
434 : END IF
435 :
436 745 : mem_for_iaK = dimen_RI*REAL(SUM(homo*virtual), KIND=dp)*8.0_dp/(1024_dp**2)
437 745 : mem_for_ijK = dimen_RI*REAL(SUM(homo(1:nspins)**2), KIND=dp)*8.0_dp/(1024_dp**2)
438 745 : mem_for_abK = dimen_RI*REAL(SUM(bse_lev_virt(1:nspins)**2), KIND=dp)*8.0_dp/(1024_dp**2)
439 :
440 334 : IF (.NOT. do_im_time) THEN
441 266 : WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ia|K) integrals:', &
442 532 : mem_for_iaK, ' MiB'
443 266 : IF (my_do_gw .AND. .NOT. do_im_time) THEN
444 35 : mem_for_iaK = dimen_RI*REAL(nmo, KIND=dp)*gw_corr_lev_total*8.0_dp/(1024_dp**2)
445 :
446 35 : WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for G0W0-(nm|K) integrals:', &
447 70 : mem_for_iaK, ' MiB'
448 : END IF
449 : END IF
450 334 : IF (do_bse) THEN
451 21 : WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ij|K) integrals:', &
452 42 : mem_for_ijK, ' MiB'
453 21 : WRITE (unit_nr, '(T3,A,T66,F11.2,A4)') 'RI_INFO| Total memory for (ab|K) integrals:', &
454 42 : mem_for_abK, ' MiB'
455 : END IF
456 334 : CALL m_flush(unit_nr)
457 : END IF
458 :
459 668 : CALL para_env%sync() ! sync to see memory output
460 :
461 : ! in case we do imaginary time, we need the overlap tensor (alpha beta P) or trunc. Coulomb tensor
462 668 : IF (.NOT. do_im_time) THEN
463 :
464 3976 : ALLOCATE (gd_B_virtual(nspins), intermed_mat(nspins))
465 2128 : ALLOCATE (my_B_virtual_start(nspins), my_B_virtual_end(nspins), my_B_size(nspins))
466 1190 : DO ispin = 1, nspins
467 :
468 : CALL create_intermediate_matrices(intermed_mat(ispin), mo_coeff_o(ispin)%matrix, virtual(ispin), homo(ispin), &
469 658 : TRIM(ADJUSTL(cp_to_string(ispin))), blacs_env_sub, para_env_sub)
470 :
471 658 : CALL create_group_dist(gd_B_virtual(ispin), para_env_sub%num_pe, virtual(ispin))
472 : CALL get_group_dist(gd_B_virtual(ispin), para_env_sub%mepos, my_B_virtual_start(ispin), my_B_virtual_end(ispin), &
473 1190 : my_B_size(ispin))
474 :
475 : END DO
476 :
477 : ! in the case of G0W0, we need (K|nm), n,m may be occ or virt (m restricted to corrected levels)
478 532 : IF (my_do_gw) THEN
479 :
480 222 : ALLOCATE (intermed_mat_gw(nspins))
481 152 : DO ispin = 1, nspins
482 : CALL create_intermediate_matrices(intermed_mat_gw(ispin), mo_coeff_gw(ispin)%matrix, &
483 : nmo, gw_corr_lev_total, &
484 : "gw_"//TRIM(ADJUSTL(cp_to_string(ispin))), &
485 152 : blacs_env_sub, para_env_sub)
486 :
487 : END DO
488 :
489 70 : CALL create_group_dist(gd_B_all, para_env_sub%num_pe, nmo)
490 70 : CALL get_group_dist(gd_B_all, para_env_sub%mepos, my_B_all_start, my_B_all_end, my_B_all_size)
491 :
492 70 : IF (do_bse) THEN
493 : ! virt x virt slab size bse_lev_virt(ispin) is per-spin, so gd_B_virt_bse is an array;
494 : ! the occupied count homo(ispin) is per-spin, so gd_B_occ_bse and the bse intermediates are arrays
495 410 : ALLOCATE (intermed_mat_bse_ab(nspins), intermed_mat_bse_ij(nspins), gd_B_occ_bse(nspins), gd_B_virt_bse(nspins))
496 168 : ALLOCATE (my_B_occ_bse_start(nspins), my_B_occ_bse_end(nspins), my_B_occ_bse_size(nspins))
497 168 : ALLOCATE (my_B_virt_bse_start(nspins), my_B_virt_bse_end(nspins), my_B_virt_bse_size(nspins))
498 92 : DO ispin = 1, nspins
499 50 : CALL create_group_dist(gd_B_virt_bse(ispin), para_env_sub%num_pe, bse_lev_virt(ispin))
500 : CALL get_group_dist(gd_B_virt_bse(ispin), para_env_sub%mepos, my_B_virt_bse_start(ispin), &
501 50 : my_B_virt_bse_end(ispin), my_B_virt_bse_size(ispin))
502 : ! virt x virt matrices
503 : CALL create_intermediate_matrices(intermed_mat_bse_ab(ispin), mo_coeff_v_bse(ispin)%matrix, &
504 : bse_lev_virt(ispin), bse_lev_virt(ispin), &
505 50 : "bse_ab_"//TRIM(ADJUSTL(cp_to_string(ispin))), blacs_env_sub, para_env_sub)
506 :
507 : ! occ x occ matrices
508 : ! We do not implement bse_lev_occ here, because the small number of occupied levels
509 : ! does not critically influence the memory
510 : CALL create_intermediate_matrices(intermed_mat_bse_ij(ispin), mo_coeff_o_bse(ispin)%matrix, &
511 : homo(ispin), homo(ispin), &
512 50 : "bse_ij_"//TRIM(ADJUSTL(cp_to_string(ispin))), blacs_env_sub, para_env_sub)
513 :
514 50 : CALL create_group_dist(gd_B_occ_bse(ispin), para_env_sub%num_pe, homo(ispin))
515 : CALL get_group_dist(gd_B_occ_bse(ispin), para_env_sub%mepos, my_B_occ_bse_start(ispin), &
516 92 : my_B_occ_bse_end(ispin), my_B_occ_bse_size(ispin))
517 : END DO
518 :
519 : END IF
520 : END IF
521 :
522 : ! array that will store the (ia|K) integrals
523 2254 : ALLOCATE (BIb_C(nspins))
524 1190 : DO ispin = 1, nspins
525 3284 : ALLOCATE (BIb_C(ispin)%array(my_group_L_size, my_B_size(ispin), homo(ispin)))
526 1763264 : BIb_C(ispin)%array = 0.0_dp
527 : END DO
528 :
529 : ! in the case of GW, we also need (nm|K)
530 532 : IF (my_do_gw) THEN
531 :
532 222 : ALLOCATE (BIb_C_gw(nspins))
533 152 : DO ispin = 1, nspins
534 410 : ALLOCATE (BIb_C_gw(ispin)%array(my_group_L_size, my_B_all_size, gw_corr_lev_total))
535 3358764 : BIb_C_gw(ispin)%array = 0.0_dp
536 : END DO
537 :
538 : END IF
539 :
540 532 : IF (do_bse) THEN
541 :
542 226 : ALLOCATE (BIb_C_bse_ij(nspins), BIb_C_bse_ab(nspins))
543 92 : DO ispin = 1, nspins
544 250 : ALLOCATE (BIb_C_bse_ij(ispin)%array(my_group_L_size, my_B_occ_bse_size(ispin), homo(ispin)))
545 62292 : BIb_C_bse_ij(ispin)%array = 0.0_dp
546 :
547 250 : ALLOCATE (BIb_C_bse_ab(ispin)%array(my_group_L_size, my_B_virt_bse_size(ispin), bse_lev_virt(ispin)))
548 2052906 : BIb_C_bse_ab(ispin)%array = 0.0_dp
549 : END DO
550 :
551 : END IF
552 :
553 532 : CALL timeset(routineN//"_loop", handle2)
554 :
555 : IF (eri_method == do_eri_mme .AND. &
556 532 : (ri_metric%potential_type == do_potential_coulomb .OR. ri_metric%potential_type == do_potential_long) .OR. &
557 : eri_method == do_eri_os .AND. ri_metric%potential_type == do_potential_coulomb) THEN
558 :
559 : ! Add a warning for automatically generated RI_AUX basis sets
560 : ! Tend to be not sufficiently converged
561 182 : IF (qs_env%mp2_env%ri_aux_auto_generated) THEN
562 : CALL cp_warn(__LOCATION__, &
563 : "At least one RI_AUX basis set was not explicitly invoked in &KIND-section. "// &
564 : "Automatically RI-basis sets and ERI_METHOD OS tend to be not converged. "// &
565 0 : "Consider specifying BASIS_SET RI_AUX explicitly with a sufficiently large basis.")
566 : END IF
567 :
568 182 : NULLIFY (mat_munu_local_L)
569 6503 : ALLOCATE (mat_munu_local_L(my_group_L_size))
570 6139 : DO LLL = 1, my_group_L_size
571 5957 : NULLIFY (mat_munu_local_L(LLL)%matrix)
572 5957 : ALLOCATE (mat_munu_local_L(LLL)%matrix)
573 5957 : CALL dbcsr_copy(mat_munu_local_L(LLL)%matrix, mat_munu%matrix)
574 6139 : CALL dbcsr_set(mat_munu_local_L(LLL)%matrix, 0.0_dp)
575 : END DO
576 : CALL mp2_eri_3c_integrate(eri_param, ri_metric, para_env_sub, qs_env, &
577 : first_c=my_group_L_start, last_c=my_group_L_end, &
578 : mat_ab=mat_munu_local_L, &
579 : basis_type_a="ORB", basis_type_b="ORB", &
580 : basis_type_c="RI_AUX", &
581 182 : sab_nl=sab_orb_sub, eri_method=eri_method)
582 :
583 378 : DO ispin = 1, nspins
584 6907 : DO LLL = 1, my_group_L_size
585 : CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat(ispin), &
586 : BIb_C(ispin)%array(LLL, :, :), &
587 : mo_coeff_o(ispin)%matrix, mo_coeff_v(ispin)%matrix, &
588 : eps_filter, &
589 6907 : my_B_virtual_end(ispin), my_B_virtual_start(ispin))
590 : END DO
591 : CALL contract_B_L(BIb_C(ispin)%array, my_Lrows, gd_B_virtual(ispin)%sizes, &
592 : gd_array%sizes, qs_env%mp2_env%eri_blksize, &
593 378 : ngroup, color_sub, para_env, para_env_sub)
594 : END DO
595 :
596 182 : IF (my_do_gw) THEN
597 :
598 152 : DO ispin = 1, nspins
599 3520 : DO LLL = 1, my_group_L_size
600 : CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat_gw(ispin), &
601 : BIb_C_gw(ispin)%array(LLL, :, :), &
602 : mo_coeff_gw(ispin)%matrix, mo_coeff_all(ispin)%matrix, eps_filter, &
603 3520 : my_B_all_end, my_B_all_start)
604 : END DO
605 : CALL contract_B_L(BIb_C_gw(ispin)%array, my_Lrows, gd_B_all%sizes, gd_array%sizes, qs_env%mp2_env%eri_blksize, &
606 152 : ngroup, color_sub, para_env, para_env_sub)
607 : END DO
608 : END IF
609 :
610 182 : IF (do_bse) THEN
611 :
612 92 : DO ispin = 1, nspins
613 : ! B^ab_P matrix elements for BSE
614 2216 : DO LLL = 1, my_group_L_size
615 : CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat_bse_ab(ispin), &
616 : BIb_C_bse_ab(ispin)%array(LLL, :, :), &
617 : mo_coeff_v_bse(ispin)%matrix, mo_coeff_v_bse(ispin)%matrix, eps_filter, &
618 2216 : my_B_all_end, my_B_all_start)
619 : END DO
620 : CALL contract_B_L(BIb_C_bse_ab(ispin)%array, my_Lrows, gd_B_virt_bse(ispin)%sizes, gd_array%sizes, &
621 50 : qs_env%mp2_env%eri_blksize, ngroup, color_sub, para_env, para_env_sub)
622 :
623 : ! B^ij_P matrix elements for BSE
624 2216 : DO LLL = 1, my_group_L_size
625 : CALL ao_to_mo_and_store_B(para_env_sub, mat_munu_local_L(LLL), intermed_mat_bse_ij(ispin), &
626 : BIb_C_bse_ij(ispin)%array(LLL, :, :), &
627 : mo_coeff_o(ispin)%matrix, mo_coeff_o(ispin)%matrix, eps_filter, &
628 2216 : my_B_occ_bse_end(ispin), my_B_occ_bse_start(ispin))
629 : END DO
630 : CALL contract_B_L(BIb_C_bse_ij(ispin)%array, my_Lrows, gd_B_occ_bse(ispin)%sizes, gd_array%sizes, &
631 92 : qs_env%mp2_env%eri_blksize, ngroup, color_sub, para_env, para_env_sub)
632 : END DO
633 :
634 : END IF
635 :
636 6139 : DO LLL = 1, my_group_L_size
637 6139 : CALL dbcsr_release_p(mat_munu_local_L(LLL)%matrix)
638 : END DO
639 182 : DEALLOCATE (mat_munu_local_L)
640 :
641 350 : ELSE IF (do_gpw) THEN
642 :
643 : CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
644 350 : auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
645 :
646 15109 : DO i_counter = 1, my_group_L_size
647 :
648 : CALL mp2_eri_3c_integrate_gpw(psi_L, rho_g, atomic_kind_set, qs_kind_set, cell, dft_control, &
649 : particle_set, pw_env_sub, my_Lrows(:, i_counter), poisson_env, rho_r, pot_g, &
650 14759 : ri_metric, mat_munu, qs_env, task_list_sub)
651 :
652 34272 : DO ispin = 1, nspins
653 : CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, intermed_mat(ispin), &
654 : BIb_C(ispin)%array(i_counter, :, :), &
655 : mo_coeff_o(ispin)%matrix, mo_coeff_v(ispin)%matrix, eps_filter, &
656 34272 : my_B_virtual_end(ispin), my_B_virtual_start(ispin))
657 :
658 : END DO
659 :
660 15109 : IF (my_do_gw) THEN
661 : ! transform (K|mu nu) to (K|nm), n corresponds to corrected GW levels, m is in nmo
662 0 : DO ispin = 1, nspins
663 : CALL ao_to_mo_and_store_B(para_env_sub, mat_munu, intermed_mat_gw(ispin), &
664 : BIb_C_gw(ispin)%array(i_counter, :, :), &
665 : mo_coeff_gw(ispin)%matrix, mo_coeff_all(ispin)%matrix, eps_filter, &
666 0 : my_B_all_end, my_B_all_start)
667 :
668 : END DO
669 : END IF
670 :
671 : END DO
672 :
673 : CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
674 350 : task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
675 : ELSE
676 0 : CPABORT("Integration method not implemented!")
677 : END IF
678 :
679 532 : CALL timestop(handle2)
680 :
681 532 : DEALLOCATE (my_Lrows)
682 :
683 1190 : DO ispin = 1, nspins
684 1190 : CALL release_intermediate_matrices(intermed_mat(ispin))
685 : END DO
686 1190 : DEALLOCATE (intermed_mat)
687 :
688 532 : IF (my_do_gw) THEN
689 152 : DO ispin = 1, nspins
690 152 : CALL release_intermediate_matrices(intermed_mat_gw(ispin))
691 : END DO
692 152 : DEALLOCATE (intermed_mat_gw)
693 : END IF
694 :
695 1064 : IF (do_bse) THEN
696 92 : DO ispin = 1, nspins
697 50 : CALL release_intermediate_matrices(intermed_mat_bse_ab(ispin))
698 92 : CALL release_intermediate_matrices(intermed_mat_bse_ij(ispin))
699 : END DO
700 142 : DEALLOCATE (intermed_mat_bse_ab, intermed_mat_bse_ij)
701 : END IF
702 :
703 : ! imag. time = low-scaling SOS-MP2, RPA, GW
704 : ELSE
705 :
706 136 : memory_info = qs_env%mp2_env%ri_rpa_im_time%memory_info
707 :
708 : ! we need 3 tensors:
709 : ! 1) t_3c_overl_int: 3c overlap integrals, optimized for easy access to integral blocks
710 : ! (atomic blocks)
711 : ! 2) t_3c_O: 3c overlap integrals, optimized for contraction (split blocks)
712 : ! 3) t_3c_M: tensor M, optimized for contraction
713 :
714 136 : CALL get_qs_env(qs_env, natom=natom, nkind=nkind, dft_control=dft_control)
715 :
716 136 : pdims_t3c = 0
717 136 : CALL dbt_pgrid_create(para_env, pdims_t3c, pgrid_t3c_overl)
718 :
719 : ! set up basis
720 544 : ALLOCATE (sizes_RI(natom), sizes_AO(natom))
721 1084 : ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
722 136 : CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
723 136 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_ri_aux)
724 136 : CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
725 136 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao)
726 :
727 : ! make sure we use the QS%EPS_PGF_ORB
728 136 : qs_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS")
729 136 : CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
730 136 : IF (n_rep /= 0) THEN
731 82 : CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=eps_pgf_orb)
732 : ELSE
733 54 : CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=eps_pgf_orb)
734 54 : eps_pgf_orb = SQRT(eps_pgf_orb)
735 : END IF
736 136 : eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
737 :
738 406 : DO ibasis = 1, SIZE(basis_set_ao)
739 270 : orb_basis => basis_set_ao(ibasis)%gto_basis_set
740 270 : CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb)
741 270 : ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
742 406 : CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb)
743 : END DO
744 :
745 136 : cut_memory_int = qs_env%mp2_env%ri_rpa_im_time%cut_memory
746 : CALL create_tensor_batches(sizes_RI, cut_memory_int, starts_array_mc_int, ends_array_mc_int, &
747 136 : starts_array_mc_block_int, ends_array_mc_block_int)
748 :
749 136 : DEALLOCATE (starts_array_mc_int, ends_array_mc_int)
750 :
751 : CALL create_3c_tensor(t_3c_overl_int_template, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, &
752 : sizes_RI, sizes_AO, sizes_AO, map1=[1, 2], map2=[3], &
753 136 : name="O (RI AO | AO)")
754 :
755 136 : CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
756 136 : CALL dbt_mp_environ_pgrid(pgrid_t3c_overl, pdims, pcoord)
757 136 : CALL mp_comm_t3c_2%create(pgrid_t3c_overl%mp_comm_2d, 3, pdims)
758 : CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, &
759 136 : nkind, particle_set, mp_comm_t3c_2, own_comm=.TRUE.)
760 136 : DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
761 :
762 : CALL build_3c_neighbor_lists(nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, &
763 : dist_3d, ri_metric, "RPA_3c_nl", qs_env, &
764 136 : sym_jk=.NOT. do_kpoints_cubic_RPA, own_dist=.TRUE.)
765 :
766 : ! init k points
767 136 : IF (do_kpoints_cubic_RPA) THEN
768 : ! set up new kpoint type with periodic images according to eps_grid from MP2 section
769 : ! instead of eps_pgf_orb from QS section
770 6 : CALL kpoint_init_cell_index(kpoints, nl_3c%jk_list, para_env, dft_control%nimages)
771 6 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") &
772 3 : "3C_OVERLAP_INTEGRALS_INFO| Number of periodic images considered:", dft_control%nimages
773 :
774 6 : nimg = dft_control%nimages
775 : ELSE
776 : nimg = 1
777 : END IF
778 :
779 1800 : ALLOCATE (t_3c_overl_int(nimg, nimg))
780 :
781 296 : DO i = 1, SIZE(t_3c_overl_int, 1)
782 576 : DO j = 1, SIZE(t_3c_overl_int, 2)
783 440 : CALL dbt_create(t_3c_overl_int_template, t_3c_overl_int(i, j))
784 : END DO
785 : END DO
786 :
787 136 : CALL dbt_destroy(t_3c_overl_int_template)
788 :
789 : ! split blocks to improve load balancing for tensor contraction
790 136 : min_bsize = qs_env%mp2_env%ri_rpa_im_time%min_bsize
791 :
792 136 : CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, min_bsize, sizes_AO_split)
793 136 : CALL pgf_block_sizes(atomic_kind_set, basis_set_ri_aux, min_bsize, sizes_RI_split)
794 :
795 136 : pdims_t3c = 0
796 136 : CALL dbt_pgrid_create(para_env, pdims_t3c, pgrid_t3c_M)
797 :
798 : ASSOCIATE (cut_memory => qs_env%mp2_env%ri_rpa_im_time%cut_memory)
799 : CALL create_tensor_batches(sizes_AO_split, cut_memory, starts_array_mc, ends_array_mc, &
800 136 : starts_array_mc_block, ends_array_mc_block)
801 : CALL create_tensor_batches(sizes_RI_split, cut_memory, &
802 : qs_env%mp2_env%ri_rpa_im_time%starts_array_mc_RI, &
803 : qs_env%mp2_env%ri_rpa_im_time%ends_array_mc_RI, &
804 : qs_env%mp2_env%ri_rpa_im_time%starts_array_mc_block_RI, &
805 272 : qs_env%mp2_env%ri_rpa_im_time%ends_array_mc_block_RI)
806 :
807 : END ASSOCIATE
808 136 : cut_memory = qs_env%mp2_env%ri_rpa_im_time%cut_memory
809 :
810 : CALL create_3c_tensor(t_3c_M, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_M, &
811 : sizes_RI_split, sizes_AO_split, sizes_AO_split, &
812 : map1=[1], map2=[2, 3], &
813 136 : name="M (RI | AO AO)")
814 136 : DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
815 136 : CALL dbt_pgrid_destroy(pgrid_t3c_M)
816 :
817 1800 : ALLOCATE (t_3c_O(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2)))
818 288946 : ALLOCATE (t_3c_O_compressed(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2), cut_memory))
819 1714 : ALLOCATE (t_3c_O_ind(SIZE(t_3c_overl_int, 1), SIZE(t_3c_overl_int, 2), cut_memory))
820 : CALL create_3c_tensor(t_3c_O(1, 1), dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, &
821 : sizes_RI_split, sizes_AO_split, sizes_AO_split, &
822 : map1=[1, 2], map2=[3], &
823 136 : name="O (RI AO | AO)")
824 136 : DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
825 136 : CALL dbt_pgrid_destroy(pgrid_t3c_overl)
826 :
827 296 : DO i = 1, SIZE(t_3c_O, 1)
828 576 : DO j = 1, SIZE(t_3c_O, 2)
829 440 : IF (i > 1 .OR. j > 1) CALL dbt_create(t_3c_O(1, 1), t_3c_O(i, j))
830 : END DO
831 : END DO
832 :
833 : ! build integrals in batches and copy to optimized format
834 : ! note: integrals are stored in terms of atomic blocks. To avoid a memory bottleneck,
835 : ! integrals are calculated in batches and copied to optimized format with subatomic blocks
836 :
837 374 : DO cm = 1, cut_memory_int
838 : CALL build_3c_integrals(t_3c_overl_int, &
839 : qs_env%mp2_env%ri_rpa_im_time%eps_filter/2, &
840 : qs_env, &
841 : nl_3c, &
842 : int_eps=qs_env%mp2_env%ri_rpa_im_time%eps_filter/2, &
843 : basis_i=basis_set_ri_aux, &
844 : basis_j=basis_set_ao, basis_k=basis_set_ao, &
845 : potential_parameter=ri_metric, &
846 : do_kpoints=do_kpoints_cubic_RPA, &
847 714 : bounds_i=[starts_array_mc_block_int(cm), ends_array_mc_block_int(cm)], desymmetrize=.FALSE.)
848 238 : CALL timeset(routineN//"_copy_3c", handle4)
849 : ! copy integral tensor t_3c_overl_int to t_3c_O tensor optimized for contraction
850 508 : DO i = 1, SIZE(t_3c_overl_int, 1)
851 938 : DO j = 1, SIZE(t_3c_overl_int, 2)
852 :
853 : CALL dbt_copy(t_3c_overl_int(i, j), t_3c_O(i, j), order=[1, 3, 2], &
854 430 : summation=.TRUE., move_data=.TRUE.)
855 430 : CALL dbt_clear(t_3c_overl_int(i, j))
856 430 : CALL dbt_filter(t_3c_O(i, j), qs_env%mp2_env%ri_rpa_im_time%eps_filter/2)
857 : ! rescaling, probably because of neighbor list
858 700 : IF (do_kpoints_cubic_RPA .AND. cm == cut_memory_int) THEN
859 150 : CALL dbt_scale(t_3c_O(i, j), 0.5_dp)
860 : END IF
861 : END DO
862 : END DO
863 612 : CALL timestop(handle4)
864 : END DO
865 :
866 296 : DO i = 1, SIZE(t_3c_overl_int, 1)
867 576 : DO j = 1, SIZE(t_3c_overl_int, 2)
868 440 : CALL dbt_destroy(t_3c_overl_int(i, j))
869 : END DO
870 : END DO
871 416 : DEALLOCATE (t_3c_overl_int)
872 :
873 136 : CALL timeset(routineN//"_copy_3c", handle4)
874 : ! desymmetrize
875 136 : CALL dbt_create(t_3c_O(1, 1), t_3c_tmp)
876 296 : DO jcell = 1, nimg
877 516 : DO kcell = 1, jcell
878 220 : CALL dbt_copy(t_3c_O(jcell, kcell), t_3c_tmp)
879 220 : CALL dbt_copy(t_3c_tmp, t_3c_O(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.)
880 380 : CALL dbt_filter(t_3c_O(kcell, jcell), qs_env%mp2_env%ri_rpa_im_time%eps_filter)
881 : END DO
882 : END DO
883 296 : DO jcell = 1, nimg
884 356 : DO kcell = jcell + 1, nimg
885 60 : CALL dbt_copy(t_3c_O(jcell, kcell), t_3c_tmp)
886 60 : CALL dbt_copy(t_3c_tmp, t_3c_O(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.)
887 220 : CALL dbt_filter(t_3c_O(kcell, jcell), qs_env%mp2_env%ri_rpa_im_time%eps_filter)
888 : END DO
889 : END DO
890 :
891 136 : CALL dbt_get_info(t_3c_O(1, 1), nfull_total=bounds_3c)
892 136 : CALL get_tensor_occupancy(t_3c_O(1, 1), nze, occ)
893 136 : memory_3c = 0.0_dp
894 :
895 408 : bounds(:, 1) = [1, bounds_3c(1)]
896 408 : bounds(:, 3) = [1, bounds_3c(3)]
897 296 : DO i = 1, SIZE(t_3c_O, 1)
898 576 : DO j = 1, SIZE(t_3c_O, 2)
899 742 : DO i_mem = 1, cut_memory
900 1386 : bounds(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
901 462 : CALL dbt_copy(t_3c_O(i, j), t_3c_tmp, bounds=bounds)
902 :
903 462 : CALL alloc_containers(t_3c_O_compressed(i, j, i_mem), 1)
904 : CALL compress_tensor(t_3c_tmp, t_3c_O_ind(i, j, i_mem)%ind, &
905 : t_3c_O_compressed(i, j, i_mem), &
906 742 : qs_env%mp2_env%ri_rpa_im_time%eps_compress, memory_3c)
907 : END DO
908 440 : CALL dbt_clear(t_3c_O(i, j))
909 : END DO
910 : END DO
911 :
912 136 : CALL para_env%sum(memory_3c)
913 :
914 136 : compression_factor = REAL(nze, dp)*1.0E-06*8.0_dp/memory_3c
915 :
916 136 : IF (unit_nr > 0) THEN
917 : WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") &
918 68 : "MEMORY_INFO| Memory for 3-center integrals (compressed):", memory_3c, ' MiB'
919 :
920 : WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") &
921 68 : "MEMORY_INFO| Compression factor: ", compression_factor
922 : END IF
923 :
924 136 : CALL dbt_destroy(t_3c_tmp)
925 :
926 136 : CALL timestop(handle4)
927 :
928 406 : DO ibasis = 1, SIZE(basis_set_ao)
929 270 : orb_basis => basis_set_ao(ibasis)%gto_basis_set
930 270 : CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb_old)
931 270 : ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
932 406 : CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb_old)
933 : END DO
934 :
935 136 : DEALLOCATE (basis_set_ri_aux, basis_set_ao)
936 :
937 544 : CALL neighbor_list_3c_destroy(nl_3c)
938 :
939 : END IF
940 :
941 668 : CALL timestop(handle)
942 :
943 2672 : END SUBROUTINE mp2_ri_gpw_compute_in
944 :
945 : ! **************************************************************************************************
946 : !> \brief Contract (P|ai) = (R|P) x (R|ai)
947 : !> \param BIb_C (R|ai)
948 : !> \param my_Lrows (R|P)
949 : !> \param sizes_B number of a (virtual) indices per subgroup process
950 : !> \param sizes_L number of P / R (auxiliary) indices per subgroup
951 : !> \param blk_size ...
952 : !> \param ngroup how many subgroups (NG)
953 : !> \param igroup subgroup color
954 : !> \param mp_comm communicator
955 : !> \param para_env_sub ...
956 : ! **************************************************************************************************
957 378 : SUBROUTINE contract_B_L(BIb_C, my_Lrows, sizes_B, sizes_L, blk_size, ngroup, igroup, mp_comm, para_env_sub)
958 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: BIb_C
959 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: my_Lrows
960 : INTEGER, DIMENSION(:), INTENT(IN) :: sizes_B, sizes_L
961 : INTEGER, DIMENSION(2), INTENT(IN) :: blk_size
962 : INTEGER, INTENT(IN) :: ngroup, igroup
963 :
964 : CLASS(mp_comm_type), INTENT(IN) :: mp_comm
965 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
966 :
967 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_B_L'
968 : LOGICAL, PARAMETER :: debug = .FALSE.
969 :
970 : INTEGER :: check_proc, handle, i, iend, ii, ioff, &
971 : istart, loc_a, loc_P, nblk_per_thread
972 378 : INTEGER, ALLOCATABLE, DIMENSION(:) :: block_ind_L_P, block_ind_L_R
973 : INTEGER, DIMENSION(1) :: dist_B_i, map_B_1, map_L_1, map_L_2, &
974 : sizes_i
975 : INTEGER, DIMENSION(2) :: map_B_2, pdims_L
976 : INTEGER, DIMENSION(3) :: pdims_B
977 : LOGICAL :: found
978 756 : INTEGER, DIMENSION(ngroup) :: dist_L_P, dist_L_R
979 756 : INTEGER, DIMENSION(para_env_sub%num_pe) :: dist_B_a
980 6426 : TYPE(dbt_distribution_type) :: dist_B, dist_L
981 1890 : TYPE(dbt_pgrid_type) :: mp_comm_B, mp_comm_L
982 9450 : TYPE(dbt_type) :: tB_in, tB_in_split, tB_out, &
983 9450 : tB_out_split, tL, tL_split
984 :
985 378 : CALL timeset(routineN, handle)
986 :
987 378 : sizes_i(1) = SIZE(BIb_C, 3)
988 :
989 : ASSOCIATE (nproc => para_env_sub%num_pe, iproc => para_env_sub%mepos, iproc_glob => mp_comm%mepos)
990 :
991 : ! local block index for R/P and a
992 378 : loc_P = igroup + 1; loc_a = iproc + 1
993 :
994 0 : CPASSERT(SIZE(sizes_L) == ngroup)
995 378 : CPASSERT(SIZE(sizes_B) == nproc)
996 378 : CPASSERT(sizes_L(loc_P) == SIZE(BIb_C, 1))
997 378 : CPASSERT(sizes_L(loc_P) == SIZE(my_Lrows, 2))
998 378 : CPASSERT(sizes_B(loc_a) == SIZE(BIb_C, 2))
999 :
1000 : ! Tensor distributions as follows:
1001 : ! Process grid NG x Nw
1002 : ! Each process has coordinates (np, nw)
1003 : ! tB_in: (R|ai): R distributed (np), a distributed (nw)
1004 : ! tB_out: (P|ai): P distributed (np), a distributed (nw)
1005 : ! tL: (R|P): R distributed (nw), P distributed (np)
1006 :
1007 : ! define mappings between tensor index and matrix index:
1008 : ! (R|ai) and (P|ai):
1009 378 : map_B_1 = [1] ! index 1 (R or P) maps to 1st matrix index (np distributed)
1010 378 : map_B_2 = [2, 3] ! indices 2, 3 (a, i) map to 2nd matrix index (nw distributed)
1011 : ! (R|P):
1012 378 : map_L_1 = [2] ! index 2 (P) maps to 1st matrix index (np distributed)
1013 378 : map_L_2 = [1] ! index 1 (R) maps to 2nd matrix index (nw distributed)
1014 :
1015 : ! derive nd process grid that is compatible with distributions and 2d process grid
1016 : ! (R|ai) / (P|ai) on process grid NG x Nw x 1
1017 : ! (R|P) on process grid NG x Nw
1018 1512 : pdims_B = [ngroup, nproc, 1]
1019 1134 : pdims_L = [nproc, ngroup]
1020 :
1021 378 : CALL dbt_pgrid_create(mp_comm, pdims_B, mp_comm_B)
1022 378 : CALL dbt_pgrid_create(mp_comm, pdims_L, mp_comm_L)
1023 :
1024 : ! setup distribution vectors such that distribution matches parallel data layout of BIb_C and my_Lrows
1025 378 : dist_B_i = [0]
1026 1524 : dist_B_a = [(i, i=0, nproc - 1)]
1027 2256 : dist_L_R = [(MODULO(i, nproc), i=0, ngroup - 1)] ! R index is replicated in my_Lrows, we impose a cyclic distribution
1028 2256 : dist_L_P = [(i, i=0, ngroup - 1)]
1029 :
1030 : ! create distributions and tensors
1031 378 : CALL dbt_distribution_new(dist_B, mp_comm_B, dist_L_P, dist_B_a, dist_B_i)
1032 378 : CALL dbt_distribution_new(dist_L, mp_comm_L, dist_L_R, dist_L_P)
1033 :
1034 378 : CALL dbt_create(tB_in, "(R|ai)", dist_B, map_B_1, map_B_2, sizes_L, sizes_B, sizes_i)
1035 378 : CALL dbt_create(tB_out, "(P|ai)", dist_B, map_B_1, map_B_2, sizes_L, sizes_B, sizes_i)
1036 378 : CALL dbt_create(tL, "(R|P)", dist_L, map_L_1, map_L_2, sizes_L, sizes_L)
1037 :
1038 : IF (debug) THEN
1039 : ! check that tensor distribution is correct
1040 : CALL dbt_get_stored_coordinates(tB_in, [loc_P, loc_a, 1], check_proc)
1041 : CPASSERT(check_proc == iproc_glob)
1042 : END IF
1043 :
1044 : ! reserve (R|ai) block
1045 378 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tB_in,loc_P,loc_a)
1046 : CALL dbt_reserve_blocks(tB_in, [loc_P], [loc_a], [1])
1047 : !$OMP END PARALLEL
1048 :
1049 : ! reserve (R|P) blocks
1050 : ! in my_Lrows, R index is replicated. For (R|P), we distribute quadratic blocks cyclically over
1051 : ! the processes in a subgroup.
1052 : ! There are NG blocks, so each process holds at most NG/Nw+1 blocks.
1053 1134 : ALLOCATE (block_ind_L_R(ngroup/nproc + 1))
1054 756 : ALLOCATE (block_ind_L_P(ngroup/nproc + 1))
1055 378 : block_ind_L_R(:) = 0; block_ind_L_P(:) = 0
1056 378 : ii = 0
1057 1128 : DO i = 1, ngroup
1058 2250 : CALL dbt_get_stored_coordinates(tL, [i, loc_P], check_proc)
1059 1128 : IF (check_proc == iproc_glob) THEN
1060 747 : ii = ii + 1
1061 747 : block_ind_L_R(ii) = i
1062 747 : block_ind_L_P(ii) = loc_P
1063 : END IF
1064 : END DO
1065 :
1066 : !TODO: Parallelize creation of block list.
1067 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tL,block_ind_L_R,block_ind_L_P,ii) &
1068 378 : !$OMP PRIVATE(nblk_per_thread,istart,iend)
1069 : nblk_per_thread = ii/omp_get_num_threads() + 1
1070 : istart = omp_get_thread_num()*nblk_per_thread + 1
1071 : iend = MIN(istart + nblk_per_thread, ii)
1072 : CALL dbt_reserve_blocks(tL, block_ind_L_R(istart:iend), block_ind_L_P(istart:iend))
1073 : !$OMP END PARALLEL
1074 :
1075 : ! insert (R|ai) block
1076 2646 : CALL dbt_put_block(tB_in, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C)
1077 :
1078 : ! insert (R|P) blocks
1079 378 : ioff = 0
1080 1506 : DO i = 1, ngroup
1081 750 : istart = ioff + 1; iend = ioff + sizes_L(i)
1082 750 : ioff = ioff + sizes_L(i)
1083 2250 : CALL dbt_get_stored_coordinates(tL, [i, loc_P], check_proc)
1084 1128 : IF (check_proc == iproc_glob) THEN
1085 1258330 : CALL dbt_put_block(tL, [i, loc_P], [sizes_L(i), sizes_L(loc_P)], my_Lrows(istart:iend, :))
1086 : END IF
1087 : END DO
1088 : END ASSOCIATE
1089 :
1090 1512 : CALL dbt_split_blocks(tB_in, tB_in_split, [blk_size(2), blk_size(1), blk_size(1)])
1091 1134 : CALL dbt_split_blocks(tL, tL_split, [blk_size(2), blk_size(2)])
1092 1512 : CALL dbt_split_blocks(tB_out, tB_out_split, [blk_size(2), blk_size(1), blk_size(1)])
1093 :
1094 : ! contract
1095 : CALL dbt_contract(alpha=1.0_dp, tensor_1=tB_in_split, tensor_2=tL_split, &
1096 : beta=0.0_dp, tensor_3=tB_out_split, &
1097 : contract_1=[1], notcontract_1=[2, 3], &
1098 : contract_2=[1], notcontract_2=[2], &
1099 378 : map_1=[2, 3], map_2=[1], optimize_dist=.TRUE.)
1100 :
1101 : ! retrieve local block of contraction result (P|ai)
1102 378 : CALL dbt_copy(tB_out_split, tB_out)
1103 :
1104 2646 : CALL dbt_get_block(tB_out, [loc_P, loc_a, 1], SHAPE(BIb_C), BIb_C, found)
1105 378 : CPASSERT(found)
1106 :
1107 : ! cleanup
1108 378 : CALL dbt_destroy(tB_in)
1109 378 : CALL dbt_destroy(tB_in_split)
1110 378 : CALL dbt_destroy(tB_out)
1111 378 : CALL dbt_destroy(tB_out_split)
1112 378 : CALL dbt_destroy(tL)
1113 378 : CALL dbt_destroy(tL_split)
1114 :
1115 378 : CALL dbt_distribution_destroy(dist_B)
1116 378 : CALL dbt_distribution_destroy(dist_L)
1117 :
1118 378 : CALL dbt_pgrid_destroy(mp_comm_B)
1119 378 : CALL dbt_pgrid_destroy(mp_comm_L)
1120 :
1121 378 : CALL timestop(handle)
1122 :
1123 756 : END SUBROUTINE contract_B_L
1124 :
1125 : ! **************************************************************************************************
1126 : !> \brief Encapsulate building of intermediate matrices matrix_ia_jnu(_beta
1127 : !> matrix_ia_jb(_beta),fm_BIb_jb(_beta),matrix_in_jnu(for G0W0) and
1128 : !> fm_BIb_all(for G0W0)
1129 : !> \param intermed_mat ...
1130 : !> \param mo_coeff_templ ...
1131 : !> \param size_1 ...
1132 : !> \param size_2 ...
1133 : !> \param matrix_name_2 ...
1134 : !> \param blacs_env_sub ...
1135 : !> \param para_env_sub ...
1136 : !> \author Jan Wilhelm
1137 : ! **************************************************************************************************
1138 0 : SUBROUTINE create_intermediate_matrices(intermed_mat, mo_coeff_templ, size_1, size_2, &
1139 : matrix_name_2, blacs_env_sub, para_env_sub)
1140 :
1141 : TYPE(intermediate_matrix_type), INTENT(OUT) :: intermed_mat
1142 : TYPE(dbcsr_type), INTENT(INOUT) :: mo_coeff_templ
1143 : INTEGER, INTENT(IN) :: size_1, size_2
1144 : CHARACTER(LEN=*), INTENT(IN) :: matrix_name_2
1145 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
1146 : TYPE(mp_para_env_type), POINTER :: para_env_sub
1147 :
1148 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_intermediate_matrices'
1149 :
1150 : INTEGER :: handle, ncol_local, nfullcols_total, &
1151 : nfullrows_total, nrow_local
1152 840 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
1153 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1154 :
1155 840 : CALL timeset(routineN, handle)
1156 :
1157 : ! initialize and create the matrix (K|jnu)
1158 840 : CALL dbcsr_create(intermed_mat%matrix_ia_jnu, template=mo_coeff_templ)
1159 :
1160 : ! Allocate Sparse matrices: (K|jb)
1161 : CALL cp_dbcsr_m_by_n_from_template(intermed_mat%matrix_ia_jb, template=mo_coeff_templ, m=size_2, n=size_1, &
1162 840 : sym=dbcsr_type_no_symmetry)
1163 :
1164 : ! set all to zero in such a way that the memory is actually allocated
1165 840 : CALL dbcsr_set(intermed_mat%matrix_ia_jnu, 0.0_dp)
1166 840 : CALL dbcsr_set(intermed_mat%matrix_ia_jb, 0.0_dp)
1167 :
1168 : ! create the analogous of matrix_ia_jb in fm type
1169 840 : NULLIFY (fm_struct)
1170 840 : CALL dbcsr_get_info(intermed_mat%matrix_ia_jb, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
1171 : CALL cp_fm_struct_create(fm_struct, context=blacs_env_sub, nrow_global=nfullrows_total, &
1172 840 : ncol_global=nfullcols_total, para_env=para_env_sub)
1173 840 : CALL cp_fm_create(intermed_mat%fm_BIb_jb, fm_struct, name="fm_BIb_jb_"//matrix_name_2)
1174 :
1175 840 : CALL copy_dbcsr_to_fm(intermed_mat%matrix_ia_jb, intermed_mat%fm_BIb_jb)
1176 840 : CALL cp_fm_struct_release(fm_struct)
1177 :
1178 : CALL cp_fm_get_info(matrix=intermed_mat%fm_BIb_jb, &
1179 : nrow_local=nrow_local, &
1180 : ncol_local=ncol_local, &
1181 : row_indices=row_indices, &
1182 840 : col_indices=col_indices)
1183 :
1184 840 : intermed_mat%max_row_col_local = MAX(nrow_local, ncol_local)
1185 840 : CALL para_env_sub%max(intermed_mat%max_row_col_local)
1186 :
1187 3360 : ALLOCATE (intermed_mat%local_col_row_info(0:intermed_mat%max_row_col_local, 2))
1188 31684 : intermed_mat%local_col_row_info = 0
1189 : ! 0,1 nrows
1190 840 : intermed_mat%local_col_row_info(0, 1) = nrow_local
1191 5693 : intermed_mat%local_col_row_info(1:nrow_local, 1) = row_indices(1:nrow_local)
1192 : ! 0,2 ncols
1193 840 : intermed_mat%local_col_row_info(0, 2) = ncol_local
1194 14474 : intermed_mat%local_col_row_info(1:ncol_local, 2) = col_indices(1:ncol_local)
1195 :
1196 840 : intermed_mat%descr = matrix_name_2
1197 :
1198 840 : CALL timestop(handle)
1199 :
1200 2520 : END SUBROUTINE create_intermediate_matrices
1201 :
1202 : ! **************************************************************************************************
1203 : !> \brief Encapsulate ERI postprocessing: AO to MO transformation and store in B matrix.
1204 : !> \param para_env ...
1205 : !> \param mat_munu ...
1206 : !> \param intermed_mat ...
1207 : !> \param BIb_jb ...
1208 : !> \param mo_coeff_o ...
1209 : !> \param mo_coeff_v ...
1210 : !> \param eps_filter ...
1211 : !> \param my_B_end ...
1212 : !> \param my_B_start ...
1213 : ! **************************************************************************************************
1214 33994 : SUBROUTINE ao_to_mo_and_store_B(para_env, mat_munu, intermed_mat, BIb_jb, &
1215 : mo_coeff_o, mo_coeff_v, eps_filter, &
1216 : my_B_end, my_B_start)
1217 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
1218 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
1219 : TYPE(intermediate_matrix_type), INTENT(INOUT) :: intermed_mat
1220 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: BIb_jb
1221 : TYPE(dbcsr_type), POINTER :: mo_coeff_o, mo_coeff_v
1222 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1223 : INTEGER, INTENT(IN) :: my_B_end, my_B_start
1224 :
1225 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ao_to_mo_and_store_B'
1226 :
1227 : INTEGER :: handle
1228 :
1229 33994 : CALL timeset(routineN//"_mult_"//TRIM(intermed_mat%descr), handle)
1230 :
1231 : CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu%matrix, mo_coeff_o, &
1232 33994 : 0.0_dp, intermed_mat%matrix_ia_jnu, filter_eps=eps_filter)
1233 : CALL dbcsr_multiply("T", "N", 1.0_dp, intermed_mat%matrix_ia_jnu, mo_coeff_v, &
1234 33994 : 0.0_dp, intermed_mat%matrix_ia_jb, filter_eps=eps_filter)
1235 33994 : CALL timestop(handle)
1236 :
1237 33994 : CALL timeset(routineN//"_E_Ex_"//TRIM(intermed_mat%descr), handle)
1238 33994 : CALL copy_dbcsr_to_fm(intermed_mat%matrix_ia_jb, intermed_mat%fm_BIb_jb)
1239 :
1240 : CALL grep_my_integrals(para_env, intermed_mat%fm_BIb_jb, BIb_jb, intermed_mat%max_row_col_local, &
1241 : intermed_mat%local_col_row_info, &
1242 33994 : my_B_end, my_B_start)
1243 :
1244 33994 : CALL timestop(handle)
1245 33994 : END SUBROUTINE ao_to_mo_and_store_B
1246 :
1247 : ! **************************************************************************************************
1248 : !> \brief ...
1249 : !> \param intermed_mat ...
1250 : ! **************************************************************************************************
1251 840 : SUBROUTINE release_intermediate_matrices(intermed_mat)
1252 : TYPE(intermediate_matrix_type), INTENT(INOUT) :: intermed_mat
1253 :
1254 840 : CALL dbcsr_release(intermed_mat%matrix_ia_jnu)
1255 840 : CALL dbcsr_release(intermed_mat%matrix_ia_jb)
1256 840 : CALL cp_fm_release(intermed_mat%fm_BIb_jb)
1257 840 : DEALLOCATE (intermed_mat%local_col_row_info)
1258 :
1259 840 : END SUBROUTINE release_intermediate_matrices
1260 :
1261 : ! **************************************************************************************************
1262 : !> \brief ...
1263 : !> \param qs_env ...
1264 : !> \param kpoints ...
1265 : !> \param unit_nr ...
1266 : ! **************************************************************************************************
1267 16 : SUBROUTINE compute_kpoints(qs_env, kpoints, unit_nr)
1268 :
1269 : TYPE(qs_environment_type), POINTER :: qs_env
1270 : TYPE(kpoint_type), POINTER :: kpoints
1271 : INTEGER :: unit_nr
1272 :
1273 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_kpoints'
1274 :
1275 : INTEGER :: handle, i, i_dim, ix, iy, iz, nkp, &
1276 : nkp_extra, nkp_orig
1277 : INTEGER, DIMENSION(3) :: nkp_grid, nkp_grid_extra, periodic
1278 : LOGICAL :: do_extrapolate_kpoints
1279 : TYPE(cell_type), POINTER :: cell
1280 : TYPE(dft_control_type), POINTER :: dft_control
1281 : TYPE(mp_para_env_type), POINTER :: para_env
1282 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1283 16 : POINTER :: sab_orb
1284 :
1285 16 : CALL timeset(routineN, handle)
1286 :
1287 16 : NULLIFY (cell, dft_control, para_env)
1288 16 : CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env, dft_control=dft_control, sab_orb=sab_orb)
1289 16 : CALL get_cell(cell=cell, periodic=periodic)
1290 :
1291 : ! general because we augment a Monkhorst-Pack mesh by additional points in the BZ
1292 16 : kpoints%kp_scheme = "GENERAL"
1293 16 : kpoints%symmetry = .FALSE.
1294 16 : kpoints%verbose = .FALSE.
1295 16 : kpoints%full_grid = .TRUE.
1296 16 : kpoints%use_real_wfn = .FALSE.
1297 16 : kpoints%eps_geo = 1.e-6_dp
1298 64 : nkp_grid(1:3) = qs_env%mp2_env%ri_rpa_im_time%kp_grid(1:3)
1299 16 : do_extrapolate_kpoints = qs_env%mp2_env%ri_rpa_im_time%do_extrapolate_kpoints
1300 :
1301 64 : DO i_dim = 1, 3
1302 48 : IF (periodic(i_dim) == 1) THEN
1303 32 : CPASSERT(MODULO(nkp_grid(i_dim), 2) == 0)
1304 : END IF
1305 64 : IF (periodic(i_dim) == 0) THEN
1306 16 : CPASSERT(nkp_grid(i_dim) == 1)
1307 : END IF
1308 : END DO
1309 :
1310 16 : nkp_orig = nkp_grid(1)*nkp_grid(2)*nkp_grid(3)/2
1311 :
1312 16 : IF (do_extrapolate_kpoints) THEN
1313 :
1314 16 : CPASSERT(qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method == kp_weights_W_uniform)
1315 :
1316 64 : DO i_dim = 1, 3
1317 48 : IF (periodic(i_dim) == 1) nkp_grid_extra(i_dim) = nkp_grid(i_dim) + 2
1318 64 : IF (periodic(i_dim) == 0) nkp_grid_extra(i_dim) = 1
1319 : END DO
1320 :
1321 64 : qs_env%mp2_env%ri_rpa_im_time%kp_grid_extra(1:3) = nkp_grid_extra(1:3)
1322 :
1323 16 : nkp_extra = nkp_grid_extra(1)*nkp_grid_extra(2)*nkp_grid_extra(3)/2
1324 :
1325 : ELSE
1326 :
1327 0 : nkp_grid_extra(1:3) = 0
1328 0 : nkp_extra = 0
1329 :
1330 : END IF
1331 :
1332 16 : nkp = nkp_orig + nkp_extra
1333 :
1334 16 : qs_env%mp2_env%ri_rpa_im_time%nkp_orig = nkp_orig
1335 16 : qs_env%mp2_env%ri_rpa_im_time%nkp_extra = nkp_extra
1336 :
1337 80 : ALLOCATE (kpoints%xkp(3, nkp), kpoints%wkp(nkp))
1338 :
1339 64 : kpoints%nkp_grid(1:3) = nkp_grid(1:3)
1340 16 : kpoints%nkp = nkp
1341 :
1342 32 : ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%wkp_V(nkp))
1343 16 : IF (do_extrapolate_kpoints) THEN
1344 : kpoints%wkp(1:nkp_orig) = 1.0_dp/REAL(nkp_orig, KIND=dp) &
1345 144 : /(1.0_dp - SQRT(REAL(nkp_extra, KIND=dp)/REAL(nkp_orig, KIND=dp)))
1346 : kpoints%wkp(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp) &
1347 304 : /(1.0_dp - SQRT(REAL(nkp_orig, KIND=dp)/REAL(nkp_extra, KIND=dp)))
1348 144 : qs_env%mp2_env%ri_rpa_im_time%wkp_V(1:nkp_orig) = 0.0_dp
1349 304 : qs_env%mp2_env%ri_rpa_im_time%wkp_V(nkp_orig + 1:nkp) = 1.0_dp/REAL(nkp_extra, KIND=dp)
1350 : ELSE
1351 0 : kpoints%wkp(:) = 1.0_dp/REAL(nkp, KIND=dp)
1352 0 : qs_env%mp2_env%ri_rpa_im_time%wkp_V(:) = kpoints%wkp(:)
1353 : END IF
1354 :
1355 16 : i = 0
1356 44 : DO ix = 1, nkp_grid(1)
1357 108 : DO iy = 1, nkp_grid(2)
1358 348 : DO iz = 1, nkp_grid(3)
1359 :
1360 256 : IF (i == nkp_orig) CYCLE
1361 128 : i = i + 1
1362 :
1363 128 : kpoints%xkp(1, i) = REAL(2*ix - nkp_grid(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(1), KIND=dp))
1364 128 : kpoints%xkp(2, i) = REAL(2*iy - nkp_grid(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(2), KIND=dp))
1365 320 : kpoints%xkp(3, i) = REAL(2*iz - nkp_grid(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid(3), KIND=dp))
1366 :
1367 : END DO
1368 : END DO
1369 : END DO
1370 :
1371 52 : DO ix = 1, nkp_grid_extra(1)
1372 148 : DO iy = 1, nkp_grid_extra(2)
1373 708 : DO iz = 1, nkp_grid_extra(3)
1374 :
1375 576 : i = i + 1
1376 576 : IF (i > nkp) CYCLE
1377 :
1378 288 : kpoints%xkp(1, i) = REAL(2*ix - nkp_grid_extra(1) - 1, KIND=dp)/(2._dp*REAL(nkp_grid_extra(1), KIND=dp))
1379 288 : kpoints%xkp(2, i) = REAL(2*iy - nkp_grid_extra(2) - 1, KIND=dp)/(2._dp*REAL(nkp_grid_extra(2), KIND=dp))
1380 672 : kpoints%xkp(3, i) = REAL(2*iz - nkp_grid_extra(3) - 1, KIND=dp)/(2._dp*REAL(nkp_grid_extra(3), KIND=dp))
1381 :
1382 : END DO
1383 : END DO
1384 : END DO
1385 :
1386 16 : CALL kpoint_init_cell_index(kpoints, sab_orb, para_env, dft_control%nimages)
1387 :
1388 16 : CALL set_qs_env(qs_env, kpoints=kpoints)
1389 :
1390 16 : IF (unit_nr > 0) THEN
1391 :
1392 8 : IF (do_extrapolate_kpoints) THEN
1393 8 : WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh for V (leading to Sigma^x):", nkp_grid(1:3)
1394 8 : WRITE (UNIT=unit_nr, FMT="(T3,A,T69)") "KPOINT_INFO| K-point extrapolation for W^c is used (W^c leads to Sigma^c):"
1395 8 : WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh 1 for W^c:", nkp_grid(1:3)
1396 8 : WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh 2 for W^c:", nkp_grid_extra(1:3)
1397 : ELSE
1398 0 : WRITE (UNIT=unit_nr, FMT="(T3,A,T69,3I4)") "KPOINT_INFO| K-point mesh for V and W:", nkp_grid(1:3)
1399 0 : WRITE (UNIT=unit_nr, FMT="(T3,A,T75,I6)") "KPOINT_INFO| Number of kpoints for V and W:", nkp
1400 : END IF
1401 :
1402 8 : SELECT CASE (qs_env%mp2_env%ri_rpa_im_time%kpoint_weights_W_method)
1403 : CASE (kp_weights_W_tailored)
1404 : WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") &
1405 0 : "KPOINT_INFO| K-point weights for W: TAILORED"
1406 : CASE (kp_weights_W_auto)
1407 : WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") &
1408 0 : "KPOINT_INFO| K-point weights for W: AUTO"
1409 : CASE (kp_weights_W_uniform)
1410 : WRITE (UNIT=unit_nr, FMT="(T3,A,T81)") &
1411 8 : "KPOINT_INFO| K-point weights for W: UNIFORM"
1412 : END SELECT
1413 :
1414 : END IF
1415 :
1416 16 : CALL timestop(handle)
1417 :
1418 16 : END SUBROUTINE compute_kpoints
1419 :
1420 : ! **************************************************************************************************
1421 : !> \brief ...
1422 : !> \param para_env_sub ...
1423 : !> \param fm_BIb_jb ...
1424 : !> \param BIb_jb ...
1425 : !> \param max_row_col_local ...
1426 : !> \param local_col_row_info ...
1427 : !> \param my_B_virtual_end ...
1428 : !> \param my_B_virtual_start ...
1429 : ! **************************************************************************************************
1430 33994 : SUBROUTINE grep_my_integrals(para_env_sub, fm_BIb_jb, BIb_jb, max_row_col_local, &
1431 : local_col_row_info, &
1432 : my_B_virtual_end, my_B_virtual_start)
1433 : TYPE(mp_para_env_type), INTENT(IN) :: para_env_sub
1434 : TYPE(cp_fm_type), INTENT(IN) :: fm_BIb_jb
1435 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: BIb_jb
1436 : INTEGER, INTENT(IN) :: max_row_col_local
1437 : INTEGER, ALLOCATABLE, DIMENSION(:, :), INTENT(IN) :: local_col_row_info
1438 : INTEGER, INTENT(IN) :: my_B_virtual_end, my_B_virtual_start
1439 :
1440 : INTEGER :: i_global, iiB, j_global, jjB, ncol_rec, &
1441 : nrow_rec, proc_receive, proc_send, &
1442 : proc_shift
1443 33994 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: rec_col_row_info
1444 33994 : INTEGER, DIMENSION(:), POINTER :: col_indices_rec, row_indices_rec
1445 33994 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_BI, rec_BI
1446 :
1447 135976 : ALLOCATE (rec_col_row_info(0:max_row_col_local, 2))
1448 :
1449 1413728 : rec_col_row_info(:, :) = local_col_row_info
1450 :
1451 33994 : nrow_rec = rec_col_row_info(0, 1)
1452 33994 : ncol_rec = rec_col_row_info(0, 2)
1453 :
1454 101940 : ALLOCATE (row_indices_rec(nrow_rec))
1455 269233 : row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1456 :
1457 101982 : ALLOCATE (col_indices_rec(ncol_rec))
1458 651279 : col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1459 :
1460 : ! accumulate data on BIb_jb buffer starting from myself
1461 651279 : DO jjB = 1, ncol_rec
1462 617285 : j_global = col_indices_rec(jjB)
1463 651279 : IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN
1464 7648657 : DO iiB = 1, nrow_rec
1465 7055972 : i_global = row_indices_rec(iiB)
1466 7648657 : BIb_jb(j_global - my_B_virtual_start + 1, i_global) = fm_BIb_jb%local_data(iiB, jjB)
1467 : END DO
1468 : END IF
1469 : END DO
1470 :
1471 33994 : DEALLOCATE (row_indices_rec)
1472 33994 : DEALLOCATE (col_indices_rec)
1473 :
1474 33994 : IF (para_env_sub%num_pe > 1) THEN
1475 9816 : ALLOCATE (local_BI(nrow_rec, ncol_rec))
1476 156209 : local_BI(1:nrow_rec, 1:ncol_rec) = fm_BIb_jb%local_data(1:nrow_rec, 1:ncol_rec)
1477 :
1478 4908 : DO proc_shift = 1, para_env_sub%num_pe - 1
1479 2454 : proc_send = MODULO(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1480 2454 : proc_receive = MODULO(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1481 :
1482 : ! first exchange information on the local data
1483 2454 : rec_col_row_info = 0
1484 2454 : CALL para_env_sub%sendrecv(local_col_row_info, proc_send, rec_col_row_info, proc_receive)
1485 2454 : nrow_rec = rec_col_row_info(0, 1)
1486 2454 : ncol_rec = rec_col_row_info(0, 2)
1487 :
1488 7362 : ALLOCATE (row_indices_rec(nrow_rec))
1489 7705 : row_indices_rec = rec_col_row_info(1:nrow_rec, 1)
1490 :
1491 7362 : ALLOCATE (col_indices_rec(ncol_rec))
1492 51654 : col_indices_rec = rec_col_row_info(1:ncol_rec, 2)
1493 :
1494 9816 : ALLOCATE (rec_BI(nrow_rec, ncol_rec))
1495 156209 : rec_BI = 0.0_dp
1496 :
1497 : ! then send and receive the real data
1498 309964 : CALL para_env_sub%sendrecv(local_BI, proc_send, rec_BI, proc_receive)
1499 :
1500 : ! accumulate the received data on BIb_jb buffer
1501 51654 : DO jjB = 1, ncol_rec
1502 49200 : j_global = col_indices_rec(jjB)
1503 51654 : IF (j_global >= my_B_virtual_start .AND. j_global <= my_B_virtual_end) THEN
1504 76719 : DO iiB = 1, nrow_rec
1505 52119 : i_global = row_indices_rec(iiB)
1506 76719 : BIb_jb(j_global - my_B_virtual_start + 1, i_global) = rec_BI(iiB, jjB)
1507 : END DO
1508 : END IF
1509 : END DO
1510 :
1511 2454 : DEALLOCATE (col_indices_rec)
1512 2454 : DEALLOCATE (row_indices_rec)
1513 4908 : DEALLOCATE (rec_BI)
1514 : END DO
1515 :
1516 2454 : DEALLOCATE (local_BI)
1517 : END IF
1518 :
1519 33994 : DEALLOCATE (rec_col_row_info)
1520 :
1521 33994 : END SUBROUTINE grep_my_integrals
1522 :
1523 0 : END MODULE mp2_integrals
|
