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