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