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