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