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