Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate HFX energy and potential
10 : !> \par History
11 : !> 11.2006 created [Manuel Guidon]
12 : !> 10.2025 added libGint interface [M Puligheddu]
13 : !> \author Manuel Guidon
14 : ! **************************************************************************************************
15 : MODULE hfx_energy_potential
16 : USE admm_types, ONLY: get_admm_env
17 : USE atomic_kind_types, ONLY: atomic_kind_type, &
18 : get_atomic_kind_set
19 : USE bibliography, ONLY: cite_reference, &
20 : guidon2008, &
21 : guidon2009
22 : USE cell_types, ONLY: cell_type, &
23 : pbc
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_files, ONLY: close_file, &
26 : open_file
27 : USE cp_log_handling, ONLY: cp_get_default_logger, &
28 : cp_logger_type
29 : USE cp_output_handling, ONLY: cp_p_file, &
30 : cp_print_key_finished_output, &
31 : cp_print_key_should_output, &
32 : cp_print_key_unit_nr
33 : USE message_passing, ONLY: mp_para_env_type
34 : USE cp_dbcsr_api, ONLY: dbcsr_copy, &
35 : dbcsr_get_matrix_type, &
36 : dbcsr_p_type, &
37 : dbcsr_type_antisymmetric, &
38 : dbcsr_dot_threadsafe
39 : USE gamma, ONLY: init_md_ftable
40 : USE hfx_communication, ONLY: distribute_ks_matrix, &
41 : get_atomic_block_maps, &
42 : get_full_density
43 : USE hfx_compression_methods, ONLY: hfx_add_mult_cache_elements, &
44 : hfx_add_single_cache_element, &
45 : hfx_decompress_first_cache, &
46 : hfx_flush_last_cache, &
47 : hfx_get_mult_cache_elements, &
48 : hfx_get_single_cache_element, &
49 : hfx_reset_cache_and_container
50 : USE hfx_contract_block, ONLY: contract_block
51 : USE hfx_libint_interface, ONLY: evaluate_eri
52 : USE hfx_load_balance_methods, ONLY: collect_load_balance_info, &
53 : hfx_load_balance, &
54 : hfx_update_load_balance
55 : USE hfx_pair_list_methods, ONLY: build_atomic_pair_list, &
56 : build_pair_list, &
57 : build_pair_list_pgf, &
58 : build_pgf_product_list, &
59 : pgf_product_list_size
60 : USE hfx_screening_methods, ONLY: calc_pair_dist_radii, &
61 : calc_screening_functions, &
62 : update_pmax_mat
63 : USE hfx_types, ONLY: &
64 : alloc_containers, dealloc_containers, hfx_basis_info_type, hfx_basis_type, hfx_cache_type, &
65 : hfx_cell_type, hfx_container_type, hfx_create_neighbor_cells, hfx_distribution, &
66 : hfx_general_type, hfx_init_container, hfx_load_balance_type, hfx_memory_type, hfx_p_kind, &
67 : hfx_pgf_list, hfx_pgf_product_list, hfx_potential_type, hfx_reset_memory_usage_counter, &
68 : hfx_screen_coeff_type, hfx_screening_type, hfx_task_list_type, hfx_type, init_t_c_g0_lmax, &
69 : log_zero, pair_list_type, pair_set_list_type
70 :
71 : USE input_constants, ONLY: do_potential_mix_cl_trunc, &
72 : do_potential_truncated, &
73 : do_potential_TShPSC, &
74 : do_potential_short, &
75 : do_potential_coulomb, &
76 : hfx_do_eval_energy, &
77 : hfx_library_is_libint, hfx_library_is_libGint, hfx_library_is_both
78 :
79 : USE input_section_types, ONLY: section_vals_type
80 : USE kinds, ONLY: default_string_length, &
81 : dp, &
82 : int_8
83 : USE kpoint_types, ONLY: get_kpoint_info, &
84 : kpoint_type
85 : USE libint_wrapper, ONLY: cp_libint_t
86 : USE machine, ONLY: m_flush, &
87 : m_memory, &
88 : m_walltime
89 : USE mathconstants, ONLY: fac
90 : USE orbital_pointers, ONLY: nco, &
91 : ncoset, &
92 : nso
93 : USE particle_types, ONLY: particle_type
94 : USE qs_environment_types, ONLY: get_qs_env, &
95 : qs_environment_type
96 : USE qs_ks_types, ONLY: get_ks_env, &
97 : qs_ks_env_type
98 : USE t_c_g0, ONLY: init
99 : USE util, ONLY: sort
100 : use, intrinsic :: iso_c_binding, only: c_ptr
101 :
102 : use cp2k_info, only: cp2k_flags
103 : use offload_api, only: offload_get_chosen_device, offload_set_chosen_device, offload_activate_chosen_device
104 : use libGint_wrapper, only: cp_libGint_init, libGint_update_env, libGint_set_density, libGint_coulomb4, &
105 : libGint_update_fock_matrix, libGint_get_fock_matrix
106 :
107 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
108 :
109 : #include "./base/base_uses.f90"
110 :
111 : IMPLICIT NONE
112 : PRIVATE
113 :
114 : PUBLIC :: integrate_four_center, coulomb4
115 :
116 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_energy_potential'
117 :
118 : !***
119 :
120 : CONTAINS
121 :
122 : ! **************************************************************************************************
123 : !> \brief computes four center integrals for a full basis set and updates the
124 : !> Kohn-Sham-Matrix and energy. Uses all 8 eri symmetries
125 : !> \param qs_env ...
126 : !> \param x_data ...
127 : !> \param ks_matrix ...
128 : !> \param ehfx energy calculated with the updated HFX matrix
129 : !> \param rho_ao density matrix in ao basis
130 : !> \param hfx_section input_section HFX
131 : !> \param para_env ...
132 : !> \param geometry_did_change flag that indicates we have to recalc integrals
133 : !> \param irep Index for the HFX replica
134 : !> \param distribute_fock_matrix Flag that indicates whether to communicate the
135 : !> new fock matrix or not
136 : !> \param ispin ...
137 : !> \par History
138 : !> 06.2007 created [Manuel Guidon]
139 : !> 08.2007 optimized load balance [Manuel Guidon]
140 : !> 09.2007 new parallelization [Manuel Guidon]
141 : !> 02.2009 completely rewritten screening part [Manuel Guidon]
142 : !> 12.2017 major bug fix. removed wrong cycle that was caussing segfault.
143 : !> see https://groups.google.com/forum/#!topic/cp2k/pc6B14XOALY
144 : !> [Tobias Binninger + Valery Weber]
145 : !> \author Manuel Guidon
146 : ! **************************************************************************************************
147 42353 : SUBROUTINE integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, &
148 : geometry_did_change, irep, distribute_fock_matrix, &
149 : ispin, nspins)
150 :
151 : TYPE(qs_environment_type), POINTER :: qs_env
152 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
153 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
154 : REAL(KIND=dp), INTENT(OUT) :: ehfx
155 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
156 : TYPE(section_vals_type), POINTER :: hfx_section
157 : TYPE(mp_para_env_type), POINTER :: para_env
158 : LOGICAL :: geometry_did_change
159 : INTEGER :: irep
160 : LOGICAL, INTENT(IN) :: distribute_fock_matrix
161 : INTEGER, INTENT(IN) :: ispin
162 : INTEGER, INTENT(IN), OPTIONAL :: nspins
163 :
164 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_four_center'
165 :
166 : CHARACTER(len=default_string_length) :: eps_scaling_str, eps_schwarz_min_str
167 : INTEGER :: act_atomic_block_offset, act_set_offset, atomic_offset_ac, atomic_offset_ad, &
168 : atomic_offset_bc, atomic_offset_bd, bin, bits_max_val, buffer_left, buffer_size, &
169 : buffer_start, cache_size, current_counter, handle, handle_bin, handle_dist_ks, &
170 : handle_getP, handle_load, handle_main, i, inla, nla, i_list_ij, i_list_kl, i_set_list_ij, &
171 : i_set_list_ij_start, i_set_list_ij_stop, i_set_list_kl, i_set_list_kl_start, &
172 : i_set_list_kl_stop, i_thread, iatom, iatom_block, iatom_end, iatom_start, ikind, img, &
173 : iset, iw, j, jatom, jatom_block, jatom_end, jatom_start, jkind, jset, katom, katom_block, &
174 : katom_end
175 : INTEGER :: katom_start, kind_kind_idx, kkind, kset, l_max, latom, latom_block, latom_end, &
176 : latom_start, lkind, lset, ma, max_am, max_pgf, max_set, mb, my_bin_id, my_bin_size, &
177 : my_thread_id, n_threads, natom, nbits, ncob, ncos_max, nints, nkimages, nkind, &
178 : nneighbors, nseta, nsetb, nsgf_max, my_nspins, pa, sgfb, shm_task_counter, shm_total_bins, &
179 : sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, &
180 : sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, swap_id, tmp_i4, unit_id
181 : INTEGER(int_8) :: atom_block, counter, estimate_to_store_int, max_val_memory, &
182 : mb_size_buffers, mb_size_f, mb_size_p, mem_compression_counter, &
183 : mem_compression_counter_disk, mem_eris, mem_eris_disk, mem_max_val, memsize_after, &
184 : memsize_before, my_current_counter, my_istart, n_processes, nblocks, ncpu, neris_disk, &
185 : neris_incore, neris_onthefly, neris_tmp, neris_total, nprim_ints, &
186 : shm_mem_compression_counter, shm_neris_disk, shm_neris_incore, shm_neris_onthefly, &
187 : shm_neris_total, shm_nprim_ints, shm_stor_count_int_disk, shm_stor_count_max_val, &
188 : shm_storage_counter_integrals, stor_count_int_disk
189 : INTEGER(int_8) :: stor_count_max_val, storage_counter_integrals, subtr_size_mb, tmp_block, &
190 : tmp_i8(8)
191 42353 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: tmp_task_list_cost
192 84706 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, last_sgf_global, nimages, &
193 84706 : tmp_index, first_set_of_atom
194 42353 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, lc_min, ld_max, &
195 84706 : ld_min, npgfa, npgfb, npgfc, npgfd, nsgfa, nsgfb, nsgfc, nsgfd, shm_block_offset
196 84706 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfb, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, &
197 42353 : offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, shm_atomic_block_offset
198 : INTEGER, DIMENSION(:, :), POINTER, SAVE :: shm_is_assoc_atomic_block
199 42353 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
200 42353 : INTEGER, DIMENSION(:, :, :, :), POINTER :: shm_set_offset
201 : INTEGER, SAVE :: shm_number_of_p_entries
202 : LOGICAL :: bins_left, buffer_overflow, do_disk_storage, do_dynamic_load_balancing, do_it, &
203 : do_kpoints, do_p_screening, do_periodic, do_print_load_balance_info, is_anti_symmetric, &
204 : ks_fully_occ, my_geo_change, treat_lsd_in_core, use_disk_storage
205 42353 : LOGICAL, DIMENSION(:, :), POINTER :: shm_atomic_pair_list
206 : REAL(dp) :: afac, bintime_start, bintime_stop, cartesian_estimate, compression_factor, &
207 : compression_factor_disk, ene_x_aa, ene_x_aa_diag, ene_x_bb, ene_x_bb_diag, eps_schwarz, &
208 : eps_storage, etmp, fac, hf_fraction, ln_10, log10_eps_schwarz, log10_pmax, &
209 : max_contraction_val, max_val1, max_val2, max_val2_set, pmax_atom, pmax_blocks, &
210 : pmax_entry, ra(3), rab2, rb(3), rc(3), rcd2, rd(3), screen_kind_ij, screen_kind_kl, &
211 : spherical_estimate, symm_fac
212 42353 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: ee_buffer1, ee_buffer2, ee_primitives_tmp, ee_work, &
213 84706 : ee_work2, kac_buf, kad_buf, kbc_buf, kbd_buf, pac_buf, pad_buf, pbc_buf, pbd_buf, &
214 42353 : primitive_integrals
215 42353 : REAL(dp), DIMENSION(:), POINTER :: p_work
216 42353 : REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta, full_ks_alpha, &
217 84706 : full_ks_beta, max_contraction, ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, shm_pmax_atom, &
218 42353 : shm_pmax_block, sphi_b, zeta, zetb, zetc, zetd, &
219 42353 : full_ks_alpha_from_gpu, full_ks_beta_from_gpu
220 42353 : REAL(dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, &
221 42353 : sphi_c_ext_set, sphi_d_ext_set
222 42353 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, &
223 42353 : sphi_d_ext
224 42353 : REAL(dp), DIMENSION(:, :, :), POINTER :: gcc
225 : REAL(KIND=dp) :: coeffs_kind_max0
226 42353 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
227 : TYPE(cell_type), POINTER :: cell
228 : TYPE(cp_libint_t) :: private_lib
229 : TYPE(cp_logger_type), POINTER :: logger
230 42353 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks_aux_fit_hfx
231 : TYPE(dft_control_type), POINTER :: dft_control
232 : TYPE(hfx_basis_info_type), POINTER :: basis_info
233 42353 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
234 42353 : TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches, integral_caches_disk
235 : TYPE(hfx_cache_type), POINTER :: maxval_cache, maxval_cache_disk
236 42353 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers, &
237 42353 : integral_containers_disk
238 : TYPE(hfx_container_type), POINTER :: maxval_container, maxval_container_disk
239 : TYPE(hfx_distribution), POINTER :: distribution_energy
240 : TYPE(hfx_general_type) :: general_parameter
241 : TYPE(hfx_load_balance_type), POINTER :: load_balance_parameter
242 : TYPE(hfx_memory_type), POINTER :: memory_parameter
243 42353 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: shm_initial_p
244 42353 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl
245 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
246 42353 : DIMENSION(:) :: pgf_product_list
247 : TYPE(hfx_potential_type) :: potential_parameter
248 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
249 42353 : POINTER :: screen_coeffs_kind, tmp_R_1, tmp_R_2, &
250 42353 : tmp_screen_pgf1, tmp_screen_pgf2
251 : TYPE(hfx_screen_coeff_type), &
252 42353 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
253 : TYPE(hfx_screen_coeff_type), &
254 42353 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf
255 : TYPE(hfx_screening_type) :: screening_parameter
256 42353 : TYPE(hfx_task_list_type), DIMENSION(:), POINTER :: shm_task_list, tmp_task_list
257 : TYPE(hfx_type), POINTER :: actual_x_data, shm_master_x_data
258 : TYPE(kpoint_type), POINTER :: kpoints
259 : TYPE(pair_list_type) :: list_ij, list_kl
260 : TYPE(pair_set_list_type), ALLOCATABLE, &
261 42353 : DIMENSION(:) :: set_list_ij, set_list_kl
262 42353 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
263 : TYPE(qs_ks_env_type), POINTER :: ks_env
264 :
265 : logical :: screened
266 : real(dp) :: max_abs_delta_KS
267 : integer :: iatom_set, jatom_set, katom_set, latom_set, shared_dev
268 : logical :: use_libint, use_libgint, use_only_alpha_spin
269 : logical :: cp2k_was_compiled_with_cuda_offload, cp2k_was_compiled_with_libint
270 : logical :: cp2k_was_compiled_with_libGint
271 :
272 42353 : NULLIFY (dft_control, matrix_ks_aux_fit_hfx)
273 :
274 42353 : CALL timeset(routineN, handle)
275 :
276 42353 : CALL cite_reference(Guidon2008)
277 42353 : CALL cite_reference(Guidon2009)
278 :
279 42353 : ehfx = 0.0_dp
280 :
281 : ! This is not very clean, but effective. ispin can only be 2 if we do the beta spin part in core
282 42353 : my_geo_change = geometry_did_change
283 42353 : IF (ispin == 2) my_geo_change = .FALSE.
284 :
285 42353 : logger => cp_get_default_logger()
286 :
287 42353 : is_anti_symmetric = dbcsr_get_matrix_type(rho_ao(1, 1)%matrix) == dbcsr_type_antisymmetric
288 :
289 42353 : iw = -1
290 :
291 42353 : IF (my_geo_change) THEN
292 2536 : CALL m_memory(memsize_before)
293 2536 : CALL para_env%max(memsize_before)
294 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
295 2536 : extension=".scfLog")
296 2536 : IF (iw > 0) THEN
297 : WRITE (UNIT=iw, FMT="(/,(T3,A,T60,I21))") &
298 476 : "HFX_MEM_INFO| Est. max. program size before HFX [MiB]:", memsize_before/(1024*1024)
299 476 : CALL m_flush(iw)
300 : END IF
301 : CALL cp_print_key_finished_output(iw, logger, hfx_section, &
302 2536 : "HF_INFO")
303 : END IF
304 :
305 42353 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell)
306 :
307 42353 : NULLIFY (cell_to_index)
308 42353 : CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
309 42353 : IF (do_kpoints) THEN
310 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
311 0 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
312 0 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
313 : END IF
314 :
315 : !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe
316 42353 : nkind = SIZE(atomic_kind_set, 1)
317 42353 : l_max = 0
318 120251 : DO ikind = 1, nkind
319 352507 : l_max = MAX(l_max, MAXVAL(x_data(1, 1)%basis_parameter(ikind)%lmax))
320 : END DO
321 42353 : l_max = 4*l_max
322 42353 : CALL init_md_ftable(l_max)
323 :
324 42353 : IF (x_data(1, 1)%potential_parameter%potential_type == do_potential_truncated .OR. &
325 : x_data(1, 1)%potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
326 10026 : IF (l_max > init_t_c_g0_lmax) THEN
327 320 : IF (para_env%is_source()) THEN
328 160 : CALL open_file(unit_number=unit_id, file_name=x_data(1, 1)%potential_parameter%filename)
329 : END IF
330 320 : CALL init(l_max, unit_id, para_env%mepos, para_env)
331 320 : IF (para_env%is_source()) THEN
332 160 : CALL close_file(unit_id)
333 : END IF
334 320 : init_t_c_g0_lmax = l_max
335 : END IF
336 : END IF
337 :
338 42353 : n_threads = 1
339 42353 : !$ n_threads = omp_get_max_threads()
340 :
341 : ! This initialization is needed to prevent a segmentation fault. The correct assigment is done below
342 42353 : my_nspins = 0
343 42353 : IF (PRESENT(nspins)) my_nspins = nspins
344 :
345 42353 : shm_neris_total = 0
346 42353 : shm_nprim_ints = 0
347 42353 : shm_neris_onthefly = 0
348 42353 : shm_storage_counter_integrals = 0
349 42353 : shm_stor_count_int_disk = 0
350 42353 : shm_neris_incore = 0
351 42353 : shm_neris_disk = 0
352 42353 : shm_stor_count_max_val = 0
353 :
354 : !$OMP PARALLEL DEFAULT(OMP_DEFAULT_NONE_WITH_OOP) &
355 : !$OMP SHARED(qs_env,&
356 : !$OMP x_data,&
357 : !$OMP ks_matrix,&
358 : !$OMP ehfx,&
359 : !$OMP rho_ao,&
360 : !$OMP matrix_ks_aux_fit_hfx,&
361 : !$OMP hfx_section,&
362 : !$OMP para_env,&
363 : !$OMP my_geo_change,&
364 : !$OMP irep,&
365 : !$OMP distribute_fock_matrix,&
366 : !$OMP cell_to_index,&
367 : !$OMP ncoset,&
368 : !$OMP nso,&
369 : !$OMP nco,&
370 : !$OMP full_ks_alpha,&
371 : !$OMP full_ks_alpha_from_gpu,&
372 : !$OMP full_ks_beta_from_gpu,&
373 : !$OMP full_ks_beta,&
374 : !$OMP n_threads,&
375 : !$OMP full_density_alpha,&
376 : !$OMP full_density_beta,&
377 : !$OMP shm_initial_p,&
378 : !$OMP shm_is_assoc_atomic_block,&
379 : !$OMP shm_number_of_p_entries,&
380 : !$OMP shm_neris_total,&
381 : !$OMP shm_neris_onthefly,&
382 : !$OMP shm_storage_counter_integrals,&
383 : !$OMP shm_stor_count_int_disk,&
384 : !$OMP shm_neris_incore,&
385 : !$OMP shm_neris_disk,&
386 : !$OMP shm_nprim_ints,&
387 : !$OMP shm_stor_count_max_val,&
388 : !$OMP cell,&
389 : !$OMP screen_coeffs_set,&
390 : !$OMP screen_coeffs_kind,&
391 : !$OMP screen_coeffs_pgf,&
392 : !$OMP pgf_product_list_size,&
393 : !$OMP radii_pgf,&
394 : !$OMP nkind,&
395 : !$OMP ispin,&
396 : !$OMP is_anti_symmetric,&
397 : !$OMP shm_atomic_block_offset,&
398 : !$OMP shm_set_offset,&
399 : !$OMP shm_block_offset,&
400 : !$OMP shm_task_counter,&
401 : !$OMP shm_task_list,&
402 : !$OMP shm_total_bins,&
403 : !$OMP shm_master_x_data,&
404 : !$OMP shm_pmax_atom,&
405 : !$OMP shm_pmax_block,&
406 : !$OMP shm_atomic_pair_list,&
407 : !$OMP shm_mem_compression_counter,&
408 : !$OMP shared_dev, &
409 : !$OMP do_print_load_balance_info,use_libgint,use_libint, &
410 : !$OMP cp2k_was_compiled_with_cuda_offload, &
411 : !$OMP cp2k_was_compiled_with_libint, &
412 : !$OMP cp2k_was_compiled_with_libGint, &
413 : !$OMP my_nspins) &
414 : !$OMP PRIVATE(ln_10,i_thread,actual_x_data,do_periodic,screening_parameter,potential_parameter,&
415 : !$OMP general_parameter,load_balance_parameter,memory_parameter,cache_size,bits_max_val,&
416 : !$OMP basis_parameter,basis_info,treat_lsd_in_core,ncpu,n_processes,neris_total,neris_incore,&
417 : !$OMP neris_disk,neris_onthefly,mem_eris,mem_eris_disk,mem_max_val,compression_factor,&
418 : !$OMP compression_factor_disk,nprim_ints,neris_tmp,max_val_memory,max_am,do_p_screening,&
419 : !$OMP max_set,particle_set,atomic_kind_set,natom,kind_of,ncos_max,nsgf_max,ikind,&
420 : !$OMP nseta,npgfa,la_max,nsgfa,primitive_integrals,pbd_buf,pbc_buf,pad_buf,pac_buf,kbd_buf,kbc_buf,&
421 : !$OMP kad_buf,kac_buf,ee_work,ee_work2,ee_buffer1,ee_buffer2,ee_primitives_tmp,max_contraction,&
422 : !$OMP max_pgf,jkind,lb_max,nsetb,npgfb,first_sgfb,sphi_b,nsgfb,ncob,sgfb,nneighbors,pgf_list_ij,pgf_list_kl,&
423 : !$OMP pgf_product_list,nimages,ks_fully_occ,subtr_size_mb,use_disk_storage,counter,do_disk_storage,&
424 : !$OMP maxval_container_disk,maxval_cache_disk,integral_containers_disk,integral_caches_disk,eps_schwarz,&
425 : !$OMP log10_eps_schwarz,eps_storage,hf_fraction,buffer_overflow,logger,private_lib,last_sgf_global,handle_getp,&
426 : !$OMP p_work,fac,handle_load,do_dynamic_load_balancing,my_bin_size,maxval_container,integral_containers,maxval_cache,&
427 : !$OMP integral_caches,tmp_task_list,tmp_task_list_cost,tmp_index,handle_main,coeffs_kind_max0,set_list_ij,&
428 : !$OMP set_list_kl,iatom_start,iatom_end,jatom_start,jatom_end,nblocks,bins_left,do_it,distribution_energy,&
429 : !$OMP my_thread_id,my_bin_id,handle_bin,bintime_start,my_istart,my_current_counter,latom_block,tmp_block,&
430 : !$OMP katom_block,katom_start,katom_end,latom_start,latom_end,pmax_blocks,list_ij,list_kl,i_set_list_ij_start,&
431 : !$OMP i_set_list_ij_stop,ra,rb,rab2,la_min,zeta,sphi_a_ext,nsgfl_a,sphi_a_u1,sphi_a_u2,sphi_a_u3,&
432 : !$OMP lb_min,zetb,sphi_b_ext,nsgfl_b,sphi_b_u1,sphi_b_u2,sphi_b_u3,katom,latom,i_set_list_kl_start,i_set_list_kl_stop,&
433 : !$OMP kkind,lkind,rc,rd,rcd2,pmax_atom,screen_kind_ij,screen_kind_kl,symm_fac,lc_max,lc_min,npgfc,zetc,nsgfc,sphi_c_ext,&
434 : !$OMP nsgfl_c,sphi_c_u1,sphi_c_u2,sphi_c_u3,ld_max,ld_min,npgfd,zetd,nsgfd,sphi_d_ext,nsgfl_d,sphi_d_u1,sphi_d_u2,&
435 : !$OMP sphi_d_u3,atomic_offset_bd,atomic_offset_bc,atomic_offset_ad,atomic_offset_ac,offset_bd_set,offset_bc_set,&
436 : !$OMP offset_ad_set,offset_ac_set,swap_id,kind_kind_idx,ptr_p_1,ptr_p_2,ptr_p_3,ptr_p_4,mem_compression_counter,&
437 : !$OMP mem_compression_counter_disk,max_val1,sphi_a_ext_set,sphi_b_ext_set,kset,lset,max_val2_set,max_val2,&
438 : !$OMP sphi_c_ext_set,sphi_d_ext_set,pmax_entry,log10_pmax,current_counter,nints,estimate_to_store_int,&
439 : !$OMP spherical_estimate,nbits,buffer_left,buffer_start,buffer_size,max_contraction_val,tmp_r_1,tmp_r_2,&
440 : !$OMP tmp_screen_pgf1,tmp_screen_pgf2,cartesian_estimate,bintime_stop,iw,memsize_after,storage_counter_integrals,&
441 : !$OMP stor_count_int_disk,stor_count_max_val,ene_x_aa,ene_x_bb,mb_size_p,mb_size_f,mb_size_buffers,afac,ene_x_aa_diag,&
442 : !$OMP ene_x_bb_diag,act_atomic_block_offset,act_set_offset,j,handle_dist_ks,tmp_i8,tmp_i4,dft_control,&
443 : !$OMP etmp,nkimages,img,bin,eps_scaling_str,eps_schwarz_min_str, &
444 : !$OMP i,inla,nla,gcc,first_set_of_atom,max_abs_delta_KS, screened, &
445 42353 : !$OMP iatom_set,jatom_set,katom_set,latom_set,use_only_alpha_spin)
446 :
447 : ln_10 = LOG(10.0_dp)
448 : i_thread = 0
449 : !$ i_thread = omp_get_thread_num()
450 :
451 : actual_x_data => x_data(irep, i_thread + 1)
452 : !$OMP MASTER
453 : shm_master_x_data => x_data(irep, 1)
454 : !$OMP END MASTER
455 : !$OMP BARRIER
456 :
457 : do_periodic = actual_x_data%periodic_parameter%do_periodic
458 :
459 : IF (do_periodic) THEN
460 : ! ** Rebuild neighbor lists in case the cell has changed (i.e. NPT MD)
461 : actual_x_data%periodic_parameter%number_of_shells = actual_x_data%periodic_parameter%mode
462 : CALL hfx_create_neighbor_cells(actual_x_data, actual_x_data%periodic_parameter%number_of_shells_from_input, &
463 : cell, i_thread)
464 : END IF
465 :
466 : screening_parameter = actual_x_data%screening_parameter
467 : potential_parameter = actual_x_data%potential_parameter
468 :
469 : general_parameter = actual_x_data%general_parameter
470 : load_balance_parameter => actual_x_data%load_balance_parameter
471 : memory_parameter => actual_x_data%memory_parameter
472 :
473 : cache_size = memory_parameter%cache_size
474 : bits_max_val = memory_parameter%bits_max_val
475 :
476 : basis_parameter => actual_x_data%basis_parameter
477 : basis_info => actual_x_data%basis_info
478 :
479 : treat_lsd_in_core = general_parameter%treat_lsd_in_core
480 :
481 : ncpu = para_env%num_pe
482 : n_processes = ncpu*n_threads
483 :
484 : !! initialize some counters
485 : neris_total = 0_int_8
486 : neris_incore = 0_int_8
487 : neris_disk = 0_int_8
488 : neris_onthefly = 0_int_8
489 : mem_eris = 0_int_8
490 : mem_eris_disk = 0_int_8
491 : mem_max_val = 0_int_8
492 : compression_factor = 0.0_dp
493 : compression_factor_disk = 0.0_dp
494 : nprim_ints = 0_int_8
495 : neris_tmp = 0_int_8
496 : max_val_memory = 1_int_8
497 :
498 : ! Set the integral engine
499 : IF (general_parameter%hfx_library == hfx_library_is_libint) THEN
500 : use_libint = .true.
501 : use_libGint = .false.
502 : END IF
503 : IF (general_parameter%hfx_library == hfx_library_is_libGint) THEN
504 : use_libint = .false.
505 : use_libGint = .true.
506 : END IF
507 : IF (general_parameter%hfx_library == hfx_library_is_both) THEN
508 : use_libint = .true.
509 : use_libGint = .true.
510 : END IF
511 :
512 : max_am = basis_info%max_am
513 :
514 : CALL get_qs_env(qs_env=qs_env, &
515 : atomic_kind_set=atomic_kind_set, &
516 : particle_set=particle_set, &
517 : dft_control=dft_control)
518 : IF (dft_control%do_admm) CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_hfx=matrix_ks_aux_fit_hfx)
519 :
520 : do_p_screening = screening_parameter%do_initial_p_screening
521 : ! Special treatment for MP2 with initial density screening
522 : IF (do_p_screening) THEN
523 : IF (ASSOCIATED(qs_env%mp2_env)) THEN
524 : IF ((qs_env%mp2_env%ri_grad%free_hfx_buffer)) THEN
525 : do_p_screening = ((qs_env%mp2_env%p_screen) .AND. (qs_env%mp2_env%not_last_hfx))
526 : ELSE
527 : do_p_screening = .FALSE.
528 : END IF
529 : END IF
530 : END IF
531 : max_set = basis_info%max_set
532 : natom = SIZE(particle_set, 1)
533 :
534 : ! Number of image matrices in k-point calculations (nkimages==1 -> no kpoints)
535 : nkimages = dft_control%nimages
536 : CPASSERT(nkimages == 1)
537 :
538 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
539 :
540 : !! precompute maximum nco and allocate scratch
541 : ncos_max = 0
542 : nsgf_max = 0
543 : DO iatom = 1, natom
544 : ikind = kind_of(iatom)
545 : nseta = basis_parameter(ikind)%nset
546 : npgfa => basis_parameter(ikind)%npgf
547 : la_max => basis_parameter(ikind)%lmax
548 : nsgfa => basis_parameter(ikind)%nsgf
549 : DO iset = 1, nseta
550 : ncos_max = MAX(ncos_max, ncoset(la_max(iset)))
551 : nsgf_max = MAX(nsgf_max, nsgfa(iset))
552 : END DO
553 : END DO
554 : !! Allocate the arrays for the integrals.
555 : ALLOCATE (primitive_integrals(nsgf_max**4))
556 : primitive_integrals = 0.0_dp
557 :
558 : ALLOCATE (pbd_buf(nsgf_max**2))
559 : ALLOCATE (pbc_buf(nsgf_max**2))
560 : ALLOCATE (pad_buf(nsgf_max**2))
561 : ALLOCATE (pac_buf(nsgf_max**2))
562 : ALLOCATE (kbd_buf(nsgf_max**2))
563 : ALLOCATE (kbc_buf(nsgf_max**2))
564 : ALLOCATE (kad_buf(nsgf_max**2))
565 : ALLOCATE (kac_buf(nsgf_max**2))
566 : ALLOCATE (ee_work(ncos_max**4))
567 : ALLOCATE (ee_work2(ncos_max**4))
568 : ALLOCATE (ee_buffer1(ncos_max**4))
569 : ALLOCATE (ee_buffer2(ncos_max**4))
570 : ALLOCATE (ee_primitives_tmp(nsgf_max**4))
571 :
572 : IF (my_nspins == 0) my_nspins = dft_control%nspins
573 :
574 : ALLOCATE (max_contraction(max_set, natom))
575 :
576 : max_contraction = 0.0_dp
577 : max_pgf = 0
578 : DO jatom = 1, natom
579 : jkind = kind_of(jatom)
580 : lb_max => basis_parameter(jkind)%lmax
581 : nsetb = basis_parameter(jkind)%nset
582 : npgfb => basis_parameter(jkind)%npgf
583 : first_sgfb => basis_parameter(jkind)%first_sgf
584 : sphi_b => basis_parameter(jkind)%sphi
585 : nsgfb => basis_parameter(jkind)%nsgf
586 : DO jset = 1, nsetb
587 : ! takes the primitive to contracted transformation into account
588 : ncob = npgfb(jset)*ncoset(lb_max(jset))
589 : sgfb = first_sgfb(1, jset)
590 : ! if the primitives are assumed to be all of max_val2, max_val2*p2s_b becomes
591 : ! the maximum value after multiplication with sphi_b
592 : max_contraction(jset, jatom) = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
593 : max_pgf = MAX(max_pgf, npgfb(jset))
594 : END DO
595 : END DO
596 :
597 : ! ** Allocate buffers for pgf_lists
598 : nneighbors = SIZE(actual_x_data%neighbor_cells)
599 : ALLOCATE (pgf_list_ij(max_pgf**2))
600 : ALLOCATE (pgf_list_kl(max_pgf**2))
601 : ! the size of pgf_product_list is allocated and resized as needed. The initial guess grows as needed
602 : !$OMP ATOMIC READ
603 : tmp_i4 = pgf_product_list_size
604 : ALLOCATE (pgf_product_list(tmp_i4))
605 : ALLOCATE (nimages(max_pgf**2))
606 :
607 : DO i = 1, max_pgf**2
608 : ALLOCATE (pgf_list_ij(i)%image_list(nneighbors))
609 : ALLOCATE (pgf_list_kl(i)%image_list(nneighbors))
610 : END DO
611 : !$OMP BARRIER
612 : !$OMP MASTER
613 : !! Calculate helper array that stores if a certain atomic pair is associated in the KS matrix
614 : IF (my_geo_change) THEN
615 : CALL get_atomic_block_maps(ks_matrix(1, 1)%matrix, basis_parameter, kind_of, &
616 : shm_master_x_data%is_assoc_atomic_block, &
617 : shm_master_x_data%number_of_p_entries, &
618 : para_env, &
619 : shm_master_x_data%atomic_block_offset, &
620 : shm_master_x_data%set_offset, &
621 : shm_master_x_data%block_offset, &
622 : shm_master_x_data%map_atoms_to_cpus, &
623 : nkind)
624 :
625 : shm_is_assoc_atomic_block => shm_master_x_data%is_assoc_atomic_block
626 :
627 : !! Get occupation of KS-matrix
628 : ks_fully_occ = .TRUE.
629 : outer: DO iatom = 1, natom
630 : DO jatom = iatom, natom
631 : IF (shm_is_assoc_atomic_block(jatom, iatom) == 0) THEN
632 : ks_fully_occ = .FALSE.
633 : EXIT outer
634 : END IF
635 : END DO
636 : END DO outer
637 :
638 : IF (.NOT. ks_fully_occ) THEN
639 : CALL cp_warn(__LOCATION__, &
640 : "The Kohn Sham matrix is not 100% occupied. This "// &
641 : "may result in incorrect Hartree-Fock results. Setting "// &
642 : "MIN_PAIR_LIST_RADIUS to -1 in the QS section ensures a "// &
643 : "fully occupied KS matrix. For more information "// &
644 : "see FAQ: https://www.cp2k.org/faq:hfx_eps_warning")
645 : END IF
646 : END IF
647 :
648 : ! ** Set pointers
649 : shm_number_of_p_entries = shm_master_x_data%number_of_p_entries
650 : shm_is_assoc_atomic_block => shm_master_x_data%is_assoc_atomic_block
651 : shm_atomic_block_offset => shm_master_x_data%atomic_block_offset
652 : shm_set_offset => shm_master_x_data%set_offset
653 : shm_block_offset => shm_master_x_data%block_offset
654 : !$OMP END MASTER
655 : !$OMP BARRIER
656 :
657 : ! ** Reset storage counter given by MAX_MEMORY by subtracting all buffers
658 : ! ** Fock and density Matrices (shared)
659 : subtr_size_mb = 2_int_8*shm_block_offset(ncpu + 1)
660 : ! ** if non restricted ==> alpha/beta spin
661 : IF (.NOT. treat_lsd_in_core) THEN
662 : IF (my_nspins == 2) subtr_size_mb = subtr_size_mb*2_int_8
663 : END IF
664 : ! ** Initial P only MAX(alpha,beta) (shared)
665 : IF (do_p_screening .OR. screening_parameter%do_p_screening_forces) THEN
666 : subtr_size_mb = subtr_size_mb + memory_parameter%size_p_screen
667 : END IF
668 : ! ** In core forces require their own initial P
669 : IF (screening_parameter%do_p_screening_forces) THEN
670 : IF (memory_parameter%treat_forces_in_core) THEN
671 : subtr_size_mb = subtr_size_mb + memory_parameter%size_p_screen
672 : END IF
673 : END IF
674 : ! ** primitive buffer (not shared by the threads)
675 : subtr_size_mb = subtr_size_mb + nsgf_max**4*n_threads
676 : ! ** density + fock buffers
677 : subtr_size_mb = subtr_size_mb + 8_int_8*nsgf_max**2*n_threads
678 : ! ** screening functions (shared)
679 : ! ** coeffs_pgf
680 : subtr_size_mb = subtr_size_mb + max_pgf**2*max_set**2*nkind**2
681 : ! ** coeffs_set
682 : subtr_size_mb = subtr_size_mb + max_set**2*nkind**2
683 : ! ** coeffs_kind
684 : subtr_size_mb = subtr_size_mb + nkind**2
685 : ! ** radii_pgf
686 : subtr_size_mb = subtr_size_mb + max_pgf**2*max_set**2*nkind**2
687 :
688 : ! ** is_assoc (shared)
689 : subtr_size_mb = subtr_size_mb + natom**2
690 :
691 : ! ** pmax_atom (shared)
692 : IF (do_p_screening) THEN
693 : subtr_size_mb = subtr_size_mb + natom**2
694 : END IF
695 : IF (screening_parameter%do_p_screening_forces) THEN
696 : IF (memory_parameter%treat_forces_in_core) THEN
697 : subtr_size_mb = subtr_size_mb + natom**2
698 : END IF
699 : END IF
700 :
701 : ! ** Convert into MiB's
702 : subtr_size_mb = subtr_size_mb*8_int_8/1024_int_8/1024_int_8
703 :
704 : ! ** Subtracting all these buffers from MAX_MEMORY yields the amount
705 : ! ** of RAM that is left for the compressed integrals. When using threads
706 : ! ** all the available memory is shared among all n_threads. i.e. the faster
707 : ! ** ones can steal from the slower ones
708 :
709 : CALL hfx_reset_memory_usage_counter(memory_parameter, subtr_size_mb)
710 :
711 : use_disk_storage = .FALSE.
712 : counter = 0_int_8
713 : do_disk_storage = memory_parameter%do_disk_storage
714 : IF (do_disk_storage) THEN
715 : maxval_container_disk => actual_x_data%store_ints%maxval_container_disk
716 : maxval_cache_disk => actual_x_data%store_ints%maxval_cache_disk
717 :
718 : integral_containers_disk => actual_x_data%store_ints%integral_containers_disk
719 : integral_caches_disk => actual_x_data%store_ints%integral_caches_disk
720 : END IF
721 :
722 : IF (my_geo_change) THEN
723 : memory_parameter%ram_counter = HUGE(memory_parameter%ram_counter)
724 : END IF
725 :
726 : IF (my_geo_change) THEN
727 : memory_parameter%recalc_forces = .TRUE.
728 : ELSE
729 : IF (.NOT. memory_parameter%treat_forces_in_core) memory_parameter%recalc_forces = .TRUE.
730 : END IF
731 :
732 : !! Get screening parameter
733 : eps_schwarz = screening_parameter%eps_schwarz
734 : IF (eps_schwarz <= 0.0_dp) THEN
735 : log10_eps_schwarz = log_zero
736 : ELSE
737 : log10_eps_schwarz = LOG10(eps_schwarz)
738 : END IF
739 : !! get storage epsilon
740 : eps_storage = eps_schwarz*memory_parameter%eps_storage_scaling
741 :
742 : !! If we have a hybrid functional, we may need only a fraction of exact exchange
743 : hf_fraction = general_parameter%fraction
744 :
745 : !! The number of integrals that fit into the given MAX_MEMORY
746 :
747 : !! Parameters related to the potential 1/r, erf(wr)/r, erfc(wr/r)
748 : potential_parameter = actual_x_data%potential_parameter
749 :
750 : !! Variable to check if we calculate the integrals in-core or on the fly
751 : !! If TRUE -> on the fly
752 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
753 : buffer_overflow = .FALSE.
754 : ELSE
755 : buffer_overflow = .TRUE.
756 : END IF
757 : logger => cp_get_default_logger()
758 :
759 : private_lib = actual_x_data%lib
760 :
761 : !! Helper array to map local basis function indices to global ones
762 : ALLOCATE (last_sgf_global(0:natom))
763 : last_sgf_global(0) = 0
764 : DO iatom = 1, natom
765 : ikind = kind_of(iatom)
766 : last_sgf_global(iatom) = last_sgf_global(iatom - 1) + basis_parameter(ikind)%nsgf_total
767 : END DO
768 : !$OMP BARRIER
769 : !$OMP MASTER
770 : !! Let master thread get the density (avoid problems with MPI)
771 : !! Get the full density from all the processors
772 : NULLIFY (full_density_alpha, full_density_beta)
773 : ALLOCATE (full_density_alpha(shm_block_offset(ncpu + 1), nkimages))
774 : IF (.NOT. treat_lsd_in_core .OR. my_nspins == 1) THEN
775 : CALL timeset(routineN//"_getP", handle_getP)
776 : DO img = 1, nkimages
777 : CALL get_full_density(para_env, full_density_alpha(:, img), rho_ao(ispin, img)%matrix, shm_number_of_p_entries, &
778 : shm_master_x_data%block_offset, &
779 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., antisymmetric=is_anti_symmetric)
780 : END DO
781 :
782 : IF (my_nspins == 2) THEN
783 : ALLOCATE (full_density_beta(shm_block_offset(ncpu + 1), nkimages))
784 : DO img = 1, nkimages
785 : CALL get_full_density(para_env, full_density_beta(:, img), rho_ao(2, img)%matrix, shm_number_of_p_entries, &
786 : shm_master_x_data%block_offset, &
787 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., antisymmetric=is_anti_symmetric)
788 : END DO
789 : END IF
790 : CALL timestop(handle_getP)
791 :
792 : !! Calculate the max values of the density matrix actual_pmax stores the data from the actual density matrix
793 : !! and x_data%initial_p stores the same for the initial guess. The initial guess is updated only in the case of
794 : !! a changed geometry
795 : NULLIFY (shm_initial_p)
796 : IF (do_p_screening) THEN
797 : shm_initial_p => shm_master_x_data%initial_p
798 : shm_pmax_atom => shm_master_x_data%pmax_atom
799 : IF (my_geo_change) THEN
800 : CALL update_pmax_mat(shm_master_x_data%initial_p, &
801 : shm_master_x_data%map_atom_to_kind_atom, &
802 : shm_master_x_data%set_offset, &
803 : shm_master_x_data%atomic_block_offset, &
804 : shm_pmax_atom, &
805 : full_density_alpha, full_density_beta, &
806 : natom, kind_of, basis_parameter, &
807 : nkind, shm_master_x_data%is_assoc_atomic_block)
808 : END IF
809 : END IF
810 : ELSE
811 : IF (do_p_screening) THEN
812 : CALL timeset(routineN//"_getP", handle_getP)
813 : DO img = 1, nkimages
814 : CALL get_full_density(para_env, full_density_alpha(:, img), rho_ao(1, img)%matrix, shm_number_of_p_entries, &
815 : shm_master_x_data%block_offset, &
816 : kind_of, basis_parameter, get_max_vals_spin=.TRUE., &
817 : rho_beta=rho_ao(2, img)%matrix, antisymmetric=is_anti_symmetric)
818 : END DO
819 : CALL timestop(handle_getP)
820 :
821 : !! Calculate the max values of the density matrix actual_pmax stores the data from the actual density matrix
822 : !! and x_data%initial_p stores the same for the initial guess. The initial guess is updated only in the case of
823 : !! a changed geometry
824 : NULLIFY (shm_initial_p)
825 : shm_initial_p => actual_x_data%initial_p
826 : shm_pmax_atom => shm_master_x_data%pmax_atom
827 : IF (my_geo_change) THEN
828 : CALL update_pmax_mat(shm_master_x_data%initial_p, &
829 : shm_master_x_data%map_atom_to_kind_atom, &
830 : shm_master_x_data%set_offset, &
831 : shm_master_x_data%atomic_block_offset, &
832 : shm_pmax_atom, &
833 : full_density_alpha, full_density_beta, &
834 : natom, kind_of, basis_parameter, &
835 : nkind, shm_master_x_data%is_assoc_atomic_block)
836 : END IF
837 : END IF
838 : ! ** Now get the density(ispin)
839 : DO img = 1, nkimages
840 : CALL get_full_density(para_env, full_density_alpha(:, img), rho_ao(ispin, img)%matrix, shm_number_of_p_entries, &
841 : shm_master_x_data%block_offset, &
842 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
843 : antisymmetric=is_anti_symmetric)
844 : END DO
845 : END IF
846 :
847 : NULLIFY (full_ks_alpha, full_ks_beta)
848 : ALLOCATE (shm_master_x_data%full_ks_alpha(shm_block_offset(ncpu + 1), nkimages))
849 : full_ks_alpha => shm_master_x_data%full_ks_alpha
850 : full_ks_alpha = 0.0_dp
851 :
852 : IF (.NOT. treat_lsd_in_core) THEN
853 : IF (my_nspins == 2) THEN
854 : ALLOCATE (shm_master_x_data%full_ks_beta(shm_block_offset(ncpu + 1), nkimages))
855 : full_ks_beta => shm_master_x_data%full_ks_beta
856 : full_ks_beta = 0.0_dp
857 : END IF
858 : END IF
859 :
860 : !$OMP END MASTER
861 : !$OMP BARRIER
862 :
863 : !! Initialize schwarz screening matrices for near field estimates and boxing screening matrices
864 : !! for far field estimates. The update is only performed if the geomtry of the system changed.
865 : !! If the system is periodic, then the corresponding routines are called and some variables
866 : !! are initialized
867 :
868 : IF (.NOT. shm_master_x_data%screen_funct_is_initialized) THEN
869 : CALL calc_pair_dist_radii(qs_env, basis_parameter, &
870 : shm_master_x_data%pair_dist_radii_pgf, max_set, max_pgf, eps_schwarz, &
871 : n_threads, i_thread)
872 : !$OMP BARRIER
873 : CALL calc_screening_functions(qs_env, basis_parameter, private_lib, shm_master_x_data%potential_parameter, &
874 : shm_master_x_data%screen_funct_coeffs_set, &
875 : shm_master_x_data%screen_funct_coeffs_kind, &
876 : shm_master_x_data%screen_funct_coeffs_pgf, &
877 : shm_master_x_data%pair_dist_radii_pgf, &
878 : max_set, max_pgf, n_threads, i_thread, p_work)
879 :
880 : !$OMP MASTER
881 : shm_master_x_data%screen_funct_is_initialized = .TRUE.
882 : !$OMP END MASTER
883 : END IF
884 : !$OMP BARRIER
885 :
886 : !$OMP MASTER
887 : screen_coeffs_set => shm_master_x_data%screen_funct_coeffs_set
888 : screen_coeffs_kind => shm_master_x_data%screen_funct_coeffs_kind
889 : screen_coeffs_pgf => shm_master_x_data%screen_funct_coeffs_pgf
890 : radii_pgf => shm_master_x_data%pair_dist_radii_pgf
891 : !$OMP END MASTER
892 : !$OMP BARRIER
893 :
894 : !! Initialize a prefactor depending on the fraction and the number of spins
895 : IF (my_nspins == 1) THEN
896 : fac = 0.5_dp*hf_fraction
897 : ELSE
898 : fac = 1.0_dp*hf_fraction
899 : END IF
900 :
901 : !! Call routines that distribute the load on all processes. If we want to screen on a initial density matrix, there is
902 : !! an optional parameter. Of course, this is only done if the geometry did change
903 : !$OMP BARRIER
904 : !$OMP MASTER
905 : CALL timeset(routineN//"_load", handle_load)
906 : !$OMP END MASTER
907 : !$OMP BARRIER
908 : IF (my_geo_change) THEN
909 : IF (actual_x_data%b_first_load_balance_energy) THEN
910 : CALL hfx_load_balance(actual_x_data, eps_schwarz, particle_set, max_set, para_env, &
911 : screen_coeffs_set, screen_coeffs_kind, &
912 : shm_is_assoc_atomic_block, do_periodic, load_balance_parameter, &
913 : kind_of, basis_parameter, shm_initial_p, shm_pmax_atom, i_thread, n_threads, &
914 : cell, do_p_screening, actual_x_data%map_atom_to_kind_atom, &
915 : nkind, hfx_do_eval_energy, shm_pmax_block, use_virial=.FALSE.)
916 : actual_x_data%b_first_load_balance_energy = .FALSE.
917 : ELSE
918 : CALL hfx_update_load_balance(actual_x_data, para_env, &
919 : load_balance_parameter, &
920 : i_thread, n_threads, hfx_do_eval_energy)
921 : END IF
922 : END IF
923 : !$OMP BARRIER
924 : !$OMP MASTER
925 : CALL timestop(handle_load)
926 : !$OMP END MASTER
927 : !$OMP BARRIER
928 :
929 : use_only_alpha_spin = treat_lsd_in_core .or. .not. my_nspins == 2
930 :
931 : ! Set up the gpu enviroment for libGint
932 : if (use_libGint) then
933 :
934 : call offload_activate_chosen_device()
935 :
936 : ! At least once, and before any other call from all OMP threads
937 : call cp_libGint_init(actual_x_data)
938 :
939 : ! At least at geo change
940 : call libGint_update_env(fac, memory_parameter, do_periodic, cell, actual_x_data, &
941 : nneighbors, max_pgf, natom, kind_of, particle_set, basis_parameter)
942 :
943 : ! Every cycle, comunicate the updated density to libGint
944 : if (use_only_alpha_spin) then
945 : call libGint_set_density(full_density_alpha)
946 : else
947 : call libGint_set_density(full_density_alpha, full_density_beta)
948 : end if
949 : end if
950 : !!!
951 :
952 : !! Start calculating integrals of the form (ab|cd) or (ij|kl)
953 : !! In order to do so, there is a main four-loop structure that takes into account the two symmetries
954 : !!
955 : !! (ab|cd) = (ba|cd) = (ab|dc) = (ba|dc)
956 : !!
957 : !! by iterating in the following way
958 : !!
959 : !! DO iatom=1,natom and DO katom=1,natom
960 : !! DO jatom=iatom,natom DO latom=katom,natom
961 : !!
962 : !! The third symmetry
963 : !!
964 : !! (ab|cd) = (cd|ab)
965 : !!
966 : !! is taken into account by the following criterion:
967 : !!
968 : !! IF(katom+latom<=iatom+jatom) THEN
969 : !! IF( ((iatom+jatom)==(katom+latom) ) .AND.(katom<iatom)) CYCLE
970 : !!
971 : !! Depending on the degeneracy of an integral the exchange contribution is multiplied by a corresponding
972 : !! factor ( symm_fac ).
973 : !!
974 : !! If one quartet does not pass the screening we CYCLE on the outer most possible loop. Thats why we use
975 : !! different hierarchies of short range screening matrices.
976 : !!
977 : !! If we do a parallel run, each process owns a unique array of workloads. Here, a workload is
978 : !! defined as :
979 : !!
980 : !! istart, jstart, kstart, lstart, number_of_atom_quartets, initial_cpu_id
981 : !!
982 : !! This tells the process where to start the main loops and how many bunches of integrals it has to
983 : !! calculate. The original parallelization is a simple modulo distribution that is binned and
984 : !! optimized in the load_balance routines. Since the Monte Carlo routines can swap processors,
985 : !! we need to know which was the initial cpu_id.
986 : !! Furthermore, the indices iatom, jatom, katom, latom have to be set to istart, jstart, kstart and
987 : !! lstart only the first time the loop is executed. All subsequent loops have to start with one or
988 : !! iatom and katom respectively. Therefore, we use flags like first_j_loop etc.
989 :
990 : !! LibGint follows the same patter, distributing integrals over mpi and openmp with a shared-over-mpi
991 : !! density and fock matrix. The main differences are that the fock matrix is not in a valid state
992 : !! until after the libGint_get_fock_matrix has returned -in approximately a thousand lines-,
993 : !! and the integrals are not saved in memory, but recomputed each cycle
994 :
995 : do_dynamic_load_balancing = .TRUE.
996 :
997 : IF (n_threads == 1 .OR. do_disk_storage) do_dynamic_load_balancing = .FALSE.
998 :
999 : IF (do_dynamic_load_balancing) THEN
1000 : my_bin_size = SIZE(actual_x_data%distribution_energy)
1001 : ELSE
1002 : my_bin_size = 1
1003 : END IF
1004 : !! We do not need the containers if MAX_MEM = 0
1005 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
1006 : !! IF new md step -> reinitialize containers
1007 : IF (my_geo_change) THEN
1008 : CALL dealloc_containers(actual_x_data%store_ints, memory_parameter%actual_memory_usage)
1009 : CALL alloc_containers(actual_x_data%store_ints, my_bin_size)
1010 :
1011 : DO bin = 1, my_bin_size
1012 : maxval_container => actual_x_data%store_ints%maxval_container(bin)
1013 : integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
1014 : CALL hfx_init_container(maxval_container, memory_parameter%actual_memory_usage, .FALSE.)
1015 : DO i = 1, 64
1016 : CALL hfx_init_container(integral_containers(i), memory_parameter%actual_memory_usage, .FALSE.)
1017 : END DO
1018 : END DO
1019 : END IF
1020 :
1021 : !! Decompress the first cache for maxvals and integrals
1022 : IF (.NOT. my_geo_change) THEN
1023 : DO bin = 1, my_bin_size
1024 : maxval_cache => actual_x_data%store_ints%maxval_cache(bin)
1025 : maxval_container => actual_x_data%store_ints%maxval_container(bin)
1026 : integral_caches => actual_x_data%store_ints%integral_caches(:, bin)
1027 : integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
1028 : CALL hfx_decompress_first_cache(bits_max_val, maxval_cache, maxval_container, &
1029 : memory_parameter%actual_memory_usage, .FALSE.)
1030 : DO i = 1, 64
1031 : CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), &
1032 : memory_parameter%actual_memory_usage, .FALSE.)
1033 : END DO
1034 : END DO
1035 : END IF
1036 : END IF
1037 :
1038 : !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here
1039 : !$OMP CRITICAL(hfxenergy_io_critical)
1040 : !! If we do disk storage, we have to initialize the containers/caches
1041 : IF (do_disk_storage) THEN
1042 : !! IF new md step -> reinitialize containers
1043 : IF (my_geo_change) THEN
1044 : CALL hfx_init_container(maxval_container_disk, memory_parameter%actual_memory_usage_disk, do_disk_storage)
1045 : DO i = 1, 64
1046 : CALL hfx_init_container(integral_containers_disk(i), memory_parameter%actual_memory_usage_disk, do_disk_storage)
1047 : END DO
1048 : END IF
1049 : !! Decompress the first cache for maxvals and integrals
1050 : IF (.NOT. my_geo_change) THEN
1051 : CALL hfx_decompress_first_cache(bits_max_val, maxval_cache_disk, maxval_container_disk, &
1052 : memory_parameter%actual_memory_usage_disk, .TRUE.)
1053 : DO i = 1, 64
1054 : CALL hfx_decompress_first_cache(i, integral_caches_disk(i), integral_containers_disk(i), &
1055 : memory_parameter%actual_memory_usage_disk, .TRUE.)
1056 : END DO
1057 : END IF
1058 : END IF
1059 : !$OMP END CRITICAL(hfxenergy_io_critical)
1060 :
1061 : !$OMP BARRIER
1062 : !$OMP MASTER
1063 :
1064 : IF (do_dynamic_load_balancing) THEN
1065 : ! ** Lets construct the task list
1066 : shm_total_bins = 0
1067 : DO i = 1, n_threads
1068 : shm_total_bins = shm_total_bins + SIZE(x_data(irep, i)%distribution_energy)
1069 : END DO
1070 : ALLOCATE (tmp_task_list(shm_total_bins))
1071 : shm_task_counter = 0
1072 : DO i = 1, n_threads
1073 : DO bin = 1, SIZE(x_data(irep, i)%distribution_energy)
1074 : shm_task_counter = shm_task_counter + 1
1075 : tmp_task_list(shm_task_counter)%thread_id = i
1076 : tmp_task_list(shm_task_counter)%bin_id = bin
1077 : tmp_task_list(shm_task_counter)%cost = x_data(irep, i)%distribution_energy(bin)%cost
1078 : END DO
1079 : END DO
1080 :
1081 : ! ** Now sort the task list
1082 : ALLOCATE (tmp_task_list_cost(shm_total_bins))
1083 : ALLOCATE (tmp_index(shm_total_bins))
1084 :
1085 : DO i = 1, shm_total_bins
1086 : tmp_task_list_cost(i) = tmp_task_list(i)%cost
1087 : END DO
1088 :
1089 : CALL sort(tmp_task_list_cost, shm_total_bins, tmp_index)
1090 :
1091 : ALLOCATE (shm_master_x_data%task_list(shm_total_bins))
1092 :
1093 : DO i = 1, shm_total_bins
1094 : shm_master_x_data%task_list(i) = tmp_task_list(tmp_index(shm_total_bins - i + 1))
1095 : END DO
1096 :
1097 : shm_task_list => shm_master_x_data%task_list
1098 : shm_task_counter = 0
1099 :
1100 : DEALLOCATE (tmp_task_list_cost, tmp_index, tmp_task_list)
1101 : END IF
1102 : !$OMP END MASTER
1103 : !$OMP BARRIER
1104 :
1105 : IF (my_bin_size > 0) THEN
1106 : maxval_container => actual_x_data%store_ints%maxval_container(1)
1107 : maxval_cache => actual_x_data%store_ints%maxval_cache(1)
1108 :
1109 : integral_containers => actual_x_data%store_ints%integral_containers(:, 1)
1110 : integral_caches => actual_x_data%store_ints%integral_caches(:, 1)
1111 : END IF
1112 :
1113 : !$OMP BARRIER
1114 : !$OMP MASTER
1115 : CALL timeset(routineN//"_main", handle_main)
1116 : !$OMP END MASTER
1117 : !$OMP BARRIER
1118 :
1119 : coeffs_kind_max0 = MAXVAL(screen_coeffs_kind(:, :)%x(2))
1120 : ALLOCATE (set_list_ij((max_set*load_balance_parameter%block_size)**2))
1121 : ALLOCATE (set_list_kl((max_set*load_balance_parameter%block_size)**2))
1122 :
1123 : !$OMP BARRIER
1124 : !$OMP MASTER
1125 :
1126 : !! precalculate maximum density matrix elements in blocks
1127 : actual_x_data%pmax_block = 0.0_dp
1128 : shm_pmax_block => actual_x_data%pmax_block
1129 : IF (do_p_screening) THEN
1130 : DO iatom_block = 1, SIZE(actual_x_data%blocks)
1131 : iatom_start = actual_x_data%blocks(iatom_block)%istart
1132 : iatom_end = actual_x_data%blocks(iatom_block)%iend
1133 : DO jatom_block = 1, SIZE(actual_x_data%blocks)
1134 : jatom_start = actual_x_data%blocks(jatom_block)%istart
1135 : jatom_end = actual_x_data%blocks(jatom_block)%iend
1136 : shm_pmax_block(iatom_block, jatom_block) = MAXVAL(shm_pmax_atom(iatom_start:iatom_end, jatom_start:jatom_end))
1137 : END DO
1138 : END DO
1139 : END IF
1140 : shm_atomic_pair_list => actual_x_data%atomic_pair_list
1141 : IF (my_geo_change) THEN
1142 : CALL build_atomic_pair_list(natom, shm_atomic_pair_list, kind_of, basis_parameter, particle_set, &
1143 : do_periodic, screen_coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
1144 : actual_x_data%blocks)
1145 : END IF
1146 :
1147 : my_bin_size = SIZE(actual_x_data%distribution_energy)
1148 : ! reset timings for the new SCF round
1149 : IF (my_geo_change) THEN
1150 : DO bin = 1, my_bin_size
1151 : actual_x_data%distribution_energy(bin)%time_first_scf = 0.0_dp
1152 : actual_x_data%distribution_energy(bin)%time_other_scf = 0.0_dp
1153 : actual_x_data%distribution_energy(bin)%time_forces = 0.0_dp
1154 : END DO
1155 : END IF
1156 : !$OMP END MASTER
1157 : !$OMP BARRIER
1158 :
1159 : my_bin_size = SIZE(actual_x_data%distribution_energy)
1160 : nblocks = load_balance_parameter%nblocks
1161 :
1162 : bins_left = .TRUE.
1163 : do_it = .TRUE.
1164 : bin = 0
1165 : ! Main cycle starts here
1166 : DO WHILE (bins_left)
1167 : ! Check if this thread can run this bin
1168 : IF (.NOT. do_dynamic_load_balancing) THEN
1169 : bin = bin + 1
1170 : IF (bin > my_bin_size) THEN
1171 : do_it = .FALSE.
1172 : bins_left = .FALSE.
1173 : ELSE
1174 : do_it = .TRUE.
1175 : bins_left = .TRUE.
1176 : distribution_energy => actual_x_data%distribution_energy(bin)
1177 : END IF
1178 : ELSE
1179 : !$OMP CRITICAL(hfxenergy_critical)
1180 : shm_task_counter = shm_task_counter + 1
1181 : IF (shm_task_counter <= shm_total_bins) THEN
1182 : my_thread_id = shm_task_list(shm_task_counter)%thread_id
1183 : my_bin_id = shm_task_list(shm_task_counter)%bin_id
1184 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
1185 : maxval_cache => x_data(irep, my_thread_id)%store_ints%maxval_cache(my_bin_id)
1186 : maxval_container => x_data(irep, my_thread_id)%store_ints%maxval_container(my_bin_id)
1187 : integral_caches => x_data(irep, my_thread_id)%store_ints%integral_caches(:, my_bin_id)
1188 : integral_containers => x_data(irep, my_thread_id)%store_ints%integral_containers(:, my_bin_id)
1189 : END IF
1190 : distribution_energy => x_data(irep, my_thread_id)%distribution_energy(my_bin_id)
1191 : do_it = .TRUE.
1192 : bins_left = .TRUE.
1193 : IF (my_geo_change) THEN
1194 : distribution_energy%ram_counter = HUGE(distribution_energy%ram_counter)
1195 : END IF
1196 : counter = 0_Int_8
1197 : ELSE
1198 : do_it = .FALSE.
1199 : bins_left = .FALSE.
1200 : END IF
1201 : !$OMP END CRITICAL(hfxenergy_critical)
1202 : END IF
1203 : IF (.NOT. do_it) CYCLE
1204 :
1205 : !$OMP MASTER
1206 : CALL timeset(routineN//"_bin", handle_bin)
1207 : !$OMP END MASTER
1208 :
1209 : bintime_start = m_walltime()
1210 : my_istart = distribution_energy%istart
1211 : my_current_counter = 0
1212 : IF (distribution_energy%number_of_atom_quartets == 0 .OR. &
1213 : my_istart == -1_int_8) my_istart = nblocks**4
1214 : atomic_blocks: DO atom_block = my_istart, nblocks**4 - 1, n_processes
1215 : latom_block = INT(MODULO(atom_block, nblocks)) + 1
1216 : tmp_block = atom_block/nblocks
1217 : katom_block = INT(MODULO(tmp_block, nblocks)) + 1
1218 : IF (latom_block < katom_block) CYCLE atomic_blocks
1219 : tmp_block = tmp_block/nblocks
1220 : jatom_block = INT(MODULO(tmp_block, nblocks)) + 1
1221 : tmp_block = tmp_block/nblocks
1222 : iatom_block = INT(MODULO(tmp_block, nblocks)) + 1
1223 : IF (jatom_block < iatom_block) CYCLE atomic_blocks
1224 : my_current_counter = my_current_counter + 1
1225 : IF (my_current_counter > distribution_energy%number_of_atom_quartets) EXIT atomic_blocks
1226 :
1227 : iatom_start = actual_x_data%blocks(iatom_block)%istart
1228 : iatom_end = actual_x_data%blocks(iatom_block)%iend
1229 : jatom_start = actual_x_data%blocks(jatom_block)%istart
1230 : jatom_end = actual_x_data%blocks(jatom_block)%iend
1231 : katom_start = actual_x_data%blocks(katom_block)%istart
1232 : katom_end = actual_x_data%blocks(katom_block)%iend
1233 : latom_start = actual_x_data%blocks(latom_block)%istart
1234 : latom_end = actual_x_data%blocks(latom_block)%iend
1235 :
1236 : pmax_blocks = MAX(shm_pmax_block(katom_block, iatom_block), &
1237 : shm_pmax_block(latom_block, jatom_block), &
1238 : shm_pmax_block(latom_block, iatom_block), &
1239 : shm_pmax_block(katom_block, jatom_block))
1240 :
1241 : IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE atomic_blocks
1242 :
1243 : CALL build_pair_list(natom, list_ij, set_list_ij, iatom_start, iatom_end, &
1244 : jatom_start, jatom_end, &
1245 : kind_of, basis_parameter, particle_set, &
1246 : do_periodic, screen_coeffs_set, screen_coeffs_kind, &
1247 : coeffs_kind_max0, log10_eps_schwarz, cell, pmax_blocks, &
1248 : shm_atomic_pair_list)
1249 :
1250 : CALL build_pair_list(natom, list_kl, set_list_kl, katom_start, katom_end, &
1251 : latom_start, latom_end, &
1252 : kind_of, basis_parameter, particle_set, &
1253 : do_periodic, screen_coeffs_set, screen_coeffs_kind, &
1254 : coeffs_kind_max0, log10_eps_schwarz, cell, pmax_blocks, &
1255 : shm_atomic_pair_list)
1256 :
1257 : DO i_list_ij = 1, list_ij%n_element
1258 :
1259 : iatom = list_ij%elements(i_list_ij)%pair(1)
1260 : jatom = list_ij%elements(i_list_ij)%pair(2)
1261 : i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
1262 : i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
1263 : ikind = list_ij%elements(i_list_ij)%kind_pair(1)
1264 : jkind = list_ij%elements(i_list_ij)%kind_pair(2)
1265 : ra = list_ij%elements(i_list_ij)%r1
1266 : rb = list_ij%elements(i_list_ij)%r2
1267 : rab2 = list_ij%elements(i_list_ij)%dist2
1268 :
1269 : la_max => basis_parameter(ikind)%lmax
1270 : la_min => basis_parameter(ikind)%lmin
1271 : npgfa => basis_parameter(ikind)%npgf
1272 : nseta = basis_parameter(ikind)%nset
1273 : zeta => basis_parameter(ikind)%zet
1274 : nsgfa => basis_parameter(ikind)%nsgf
1275 : sphi_a_ext => basis_parameter(ikind)%sphi_ext(:, :, :, :)
1276 : nsgfl_a => basis_parameter(ikind)%nsgfl
1277 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
1278 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
1279 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
1280 :
1281 : lb_max => basis_parameter(jkind)%lmax
1282 : lb_min => basis_parameter(jkind)%lmin
1283 : npgfb => basis_parameter(jkind)%npgf
1284 : nsetb = basis_parameter(jkind)%nset
1285 : zetb => basis_parameter(jkind)%zet
1286 : nsgfb => basis_parameter(jkind)%nsgf
1287 : sphi_b_ext => basis_parameter(jkind)%sphi_ext(:, :, :, :)
1288 : nsgfl_b => basis_parameter(jkind)%nsgfl
1289 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
1290 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
1291 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
1292 :
1293 : DO i_list_kl = 1, list_kl%n_element
1294 : katom = list_kl%elements(i_list_kl)%pair(1)
1295 : latom = list_kl%elements(i_list_kl)%pair(2)
1296 :
1297 : IF (.NOT. (katom + latom <= iatom + jatom)) CYCLE
1298 : IF (((iatom + jatom) == (katom + latom)) .AND. (katom < iatom)) CYCLE
1299 : i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
1300 : i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
1301 : kkind = list_kl%elements(i_list_kl)%kind_pair(1)
1302 : lkind = list_kl%elements(i_list_kl)%kind_pair(2)
1303 : rc = list_kl%elements(i_list_kl)%r1
1304 : rd = list_kl%elements(i_list_kl)%r2
1305 : rcd2 = list_kl%elements(i_list_kl)%dist2
1306 :
1307 : IF (do_p_screening) THEN
1308 : pmax_atom = MAX(shm_pmax_atom(katom, iatom), &
1309 : shm_pmax_atom(latom, jatom), &
1310 : shm_pmax_atom(latom, iatom), &
1311 : shm_pmax_atom(katom, jatom))
1312 : ELSE
1313 : pmax_atom = 0.0_dp
1314 : END IF
1315 :
1316 : screen_kind_ij = screen_coeffs_kind(jkind, ikind)%x(1)*rab2 + &
1317 : screen_coeffs_kind(jkind, ikind)%x(2)
1318 : screen_kind_kl = screen_coeffs_kind(lkind, kkind)%x(1)*rcd2 + &
1319 : screen_coeffs_kind(lkind, kkind)%x(2)
1320 :
1321 : IF (screen_kind_ij + screen_kind_kl + pmax_atom < log10_eps_schwarz) CYCLE
1322 :
1323 : !! we want to be consistent with the KS matrix. If none of the atomic indices
1324 : !! is associated cycle
1325 : IF (.NOT. (shm_is_assoc_atomic_block(latom, iatom) >= 1 .AND. &
1326 : shm_is_assoc_atomic_block(katom, iatom) >= 1 .AND. &
1327 : shm_is_assoc_atomic_block(katom, jatom) >= 1 .AND. &
1328 : shm_is_assoc_atomic_block(latom, jatom) >= 1)) CYCLE
1329 :
1330 : !! calculate symmetry_factor according to degeneracy of atomic quartet
1331 : symm_fac = 0.5_dp
1332 : IF (iatom == jatom) symm_fac = symm_fac*2.0_dp
1333 : IF (katom == latom) symm_fac = symm_fac*2.0_dp
1334 : IF (iatom == katom .AND. jatom == latom .AND. iatom /= jatom .AND. katom /= latom) symm_fac = symm_fac*2.0_dp
1335 : IF (iatom == katom .AND. iatom == jatom .AND. katom == latom) symm_fac = symm_fac*2.0_dp
1336 : symm_fac = 1.0_dp/symm_fac
1337 :
1338 : lc_max => basis_parameter(kkind)%lmax
1339 : lc_min => basis_parameter(kkind)%lmin
1340 : npgfc => basis_parameter(kkind)%npgf
1341 : zetc => basis_parameter(kkind)%zet
1342 : nsgfc => basis_parameter(kkind)%nsgf
1343 : sphi_c_ext => basis_parameter(kkind)%sphi_ext(:, :, :, :)
1344 : nsgfl_c => basis_parameter(kkind)%nsgfl
1345 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
1346 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
1347 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
1348 :
1349 : ld_max => basis_parameter(lkind)%lmax
1350 : ld_min => basis_parameter(lkind)%lmin
1351 : npgfd => basis_parameter(lkind)%npgf
1352 : zetd => basis_parameter(lkind)%zet
1353 : nsgfd => basis_parameter(lkind)%nsgf
1354 : sphi_d_ext => basis_parameter(lkind)%sphi_ext(:, :, :, :)
1355 : nsgfl_d => basis_parameter(lkind)%nsgfl
1356 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
1357 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
1358 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
1359 :
1360 : ! Note: shm_atomic_block_offset is already symmetric
1361 : ! see get_atomic_block_maps in hfx_communication.F
1362 : atomic_offset_bd = shm_atomic_block_offset(jatom, latom)
1363 : atomic_offset_bc = shm_atomic_block_offset(jatom, katom)
1364 : atomic_offset_ad = shm_atomic_block_offset(iatom, latom)
1365 : atomic_offset_ac = shm_atomic_block_offset(iatom, katom)
1366 :
1367 : IF (jatom < latom) THEN
1368 : offset_bd_set => shm_set_offset(:, :, lkind, jkind)
1369 : ELSE
1370 : offset_bd_set => shm_set_offset(:, :, jkind, lkind)
1371 : END IF
1372 : IF (jatom < katom) THEN
1373 : offset_bc_set => shm_set_offset(:, :, kkind, jkind)
1374 : ELSE
1375 : offset_bc_set => shm_set_offset(:, :, jkind, kkind)
1376 : END IF
1377 : IF (iatom < latom) THEN
1378 : offset_ad_set => shm_set_offset(:, :, lkind, ikind)
1379 : ELSE
1380 : offset_ad_set => shm_set_offset(:, :, ikind, lkind)
1381 : END IF
1382 : IF (iatom < katom) THEN
1383 : offset_ac_set => shm_set_offset(:, :, kkind, ikind)
1384 : ELSE
1385 : offset_ac_set => shm_set_offset(:, :, ikind, kkind)
1386 : END IF
1387 :
1388 : IF (do_p_screening) THEN
1389 : swap_id = 0
1390 : kind_kind_idx = INT(get_1D_idx(kkind, ikind, INT(nkind, int_8)))
1391 : IF (ikind >= kkind) THEN
1392 : ptr_p_1 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1393 : actual_x_data%map_atom_to_kind_atom(katom), &
1394 : actual_x_data%map_atom_to_kind_atom(iatom))
1395 : ELSE
1396 : ptr_p_1 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1397 : actual_x_data%map_atom_to_kind_atom(iatom), &
1398 : actual_x_data%map_atom_to_kind_atom(katom))
1399 : swap_id = swap_id + 1
1400 : END IF
1401 : kind_kind_idx = INT(get_1D_idx(lkind, jkind, INT(nkind, int_8)))
1402 : IF (jkind >= lkind) THEN
1403 : ptr_p_2 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1404 : actual_x_data%map_atom_to_kind_atom(latom), &
1405 : actual_x_data%map_atom_to_kind_atom(jatom))
1406 : ELSE
1407 : ptr_p_2 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1408 : actual_x_data%map_atom_to_kind_atom(jatom), &
1409 : actual_x_data%map_atom_to_kind_atom(latom))
1410 : swap_id = swap_id + 2
1411 : END IF
1412 : kind_kind_idx = INT(get_1D_idx(lkind, ikind, INT(nkind, int_8)))
1413 : IF (ikind >= lkind) THEN
1414 : ptr_p_3 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1415 : actual_x_data%map_atom_to_kind_atom(latom), &
1416 : actual_x_data%map_atom_to_kind_atom(iatom))
1417 : ELSE
1418 : ptr_p_3 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1419 : actual_x_data%map_atom_to_kind_atom(iatom), &
1420 : actual_x_data%map_atom_to_kind_atom(latom))
1421 : swap_id = swap_id + 4
1422 : END IF
1423 : kind_kind_idx = INT(get_1D_idx(kkind, jkind, INT(nkind, int_8)))
1424 : IF (jkind >= kkind) THEN
1425 : ptr_p_4 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1426 : actual_x_data%map_atom_to_kind_atom(katom), &
1427 : actual_x_data%map_atom_to_kind_atom(jatom))
1428 : ELSE
1429 : ptr_p_4 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1430 : actual_x_data%map_atom_to_kind_atom(jatom), &
1431 : actual_x_data%map_atom_to_kind_atom(katom))
1432 : swap_id = swap_id + 8
1433 : END IF
1434 : END IF
1435 :
1436 : !! At this stage, check for memory used in compression
1437 : IF (my_geo_change) THEN
1438 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
1439 : ! ** We know the maximum amount of integrals that we can store per MPI-process
1440 : ! ** Now we need to sum the current memory usage among all openMP threads
1441 : ! ** We can just read what is currently stored on the corresponding x_data type
1442 : ! ** If this is thread i and it tries to read the data from thread j, that is
1443 : ! ** currently writing that data, we just dont care, because the possible error
1444 : ! ** is of the order of CACHE_SIZE
1445 : mem_compression_counter = 0
1446 : DO i = 1, n_threads
1447 : !$OMP ATOMIC READ
1448 : tmp_i4 = x_data(irep, i)%memory_parameter%actual_memory_usage
1449 : mem_compression_counter = mem_compression_counter + &
1450 : tmp_i4*memory_parameter%cache_size
1451 : END DO
1452 : IF (mem_compression_counter > memory_parameter%max_compression_counter) THEN
1453 : buffer_overflow = .TRUE.
1454 : IF (do_dynamic_load_balancing) THEN
1455 : distribution_energy%ram_counter = counter
1456 : ELSE
1457 : memory_parameter%ram_counter = counter
1458 : END IF
1459 : ELSE
1460 : counter = counter + 1
1461 : buffer_overflow = .FALSE.
1462 : END IF
1463 : END IF
1464 : ELSE
1465 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
1466 : IF (do_dynamic_load_balancing) THEN
1467 : IF (distribution_energy%ram_counter == counter) THEN
1468 : buffer_overflow = .TRUE.
1469 : ELSE
1470 : counter = counter + 1
1471 : buffer_overflow = .FALSE.
1472 : END IF
1473 :
1474 : ELSE
1475 : IF (memory_parameter%ram_counter == counter) THEN
1476 : buffer_overflow = .TRUE.
1477 : ELSE
1478 : counter = counter + 1
1479 : buffer_overflow = .FALSE.
1480 : END IF
1481 : END IF
1482 : END IF
1483 : END IF
1484 :
1485 : IF (buffer_overflow .AND. do_disk_storage) THEN
1486 : use_disk_storage = .TRUE.
1487 : buffer_overflow = .FALSE.
1488 : END IF
1489 :
1490 : IF (use_disk_storage) THEN
1491 : !$OMP ATOMIC READ
1492 : tmp_i4 = memory_parameter%actual_memory_usage_disk
1493 : mem_compression_counter_disk = tmp_i4*memory_parameter%cache_size
1494 : IF (mem_compression_counter_disk > memory_parameter%max_compression_counter_disk) THEN
1495 : buffer_overflow = .TRUE.
1496 : use_disk_storage = .FALSE.
1497 : END IF
1498 : END IF
1499 :
1500 : DO i_set_list_ij = i_set_list_ij_start, i_set_list_ij_stop
1501 : iset = set_list_ij(i_set_list_ij)%pair(1)
1502 : jset = set_list_ij(i_set_list_ij)%pair(2)
1503 :
1504 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1505 : max_val1 = screen_coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + &
1506 : screen_coeffs_set(jset, iset, jkind, ikind)%x(2)
1507 :
1508 : IF (max_val1 + screen_kind_kl + pmax_atom < log10_eps_schwarz) CYCLE
1509 :
1510 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
1511 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
1512 : DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
1513 : kset = set_list_kl(i_set_list_kl)%pair(1)
1514 : lset = set_list_kl(i_set_list_kl)%pair(2)
1515 :
1516 : max_val2_set = (screen_coeffs_set(lset, kset, lkind, kkind)%x(1)*rcd2 + &
1517 : screen_coeffs_set(lset, kset, lkind, kkind)%x(2))
1518 : max_val2 = max_val1 + max_val2_set
1519 :
1520 : !! Near field screening
1521 : IF (max_val2 + pmax_atom < log10_eps_schwarz) CYCLE
1522 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
1523 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
1524 : !! get max_vals if we screen on initial density
1525 : IF (do_p_screening) THEN
1526 : CALL get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, &
1527 : iset, jset, kset, lset, &
1528 : pmax_entry, swap_id)
1529 : ELSE
1530 : pmax_entry = 0.0_dp
1531 : END IF
1532 : log10_pmax = pmax_entry
1533 : max_val2 = max_val2 + log10_pmax
1534 : IF (max_val2 < log10_eps_schwarz) CYCLE
1535 : pmax_entry = EXP(log10_pmax*ln_10)
1536 :
1537 : ! libGint can only run in direct mode at the moment
1538 : ! so there is no need to check the buffers
1539 : if (use_libGint) then
1540 : tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind)
1541 : tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind)
1542 : ! TODO a more modern approach in which the call is more like
1543 : ! call coulom_4_gpu( HFX_data, iset, jset, kset, lset )
1544 : call libGint_coulomb4( &
1545 : iatom, jatom, katom, latom, iset, jset, kset, lset, &
1546 : ra, rb, rc, rd, &
1547 : npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
1548 : potential_parameter, &
1549 : screen_coeffs_set(jset, iset, jkind, ikind)%x, &
1550 : screen_coeffs_set(lset, kset, lkind, kkind)%x, &
1551 : log10_pmax, log10_eps_schwarz, &
1552 : tmp_screen_pgf1, tmp_screen_pgf2, &
1553 : actual_x_data%neighbor_cells, cell, do_periodic, screened)
1554 :
1555 : if (.not. screened) then
1556 :
1557 : ! TODO a more modern approach in which the call is more like
1558 : ! call update_fock_matrix_gpu( HFX_data, iset, jset, kset, lset )
1559 : call libGint_update_fock_matrix( &
1560 : symm_fac, &
1561 : iatom, jatom, katom, latom, &
1562 : iset, jset, kset, lset, &
1563 : atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, &
1564 : offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, &
1565 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1566 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
1567 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
1568 : nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d)
1569 :
1570 : end if
1571 : end if
1572 : ! Run using libint as integral engine
1573 : if_use_libint: if (use_libint) then
1574 : !! store current number of integrals, update total number and number of integrals in buffer
1575 : current_counter = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)
1576 : IF (buffer_overflow) THEN
1577 : neris_onthefly = neris_onthefly + current_counter
1578 : END IF
1579 :
1580 : !! Get integrals from buffer and update Kohn-Sham matrix
1581 : IF (.NOT. buffer_overflow .AND. .NOT. my_geo_change) THEN
1582 : nints = current_counter
1583 : IF (.NOT. use_disk_storage) THEN
1584 : CALL hfx_get_single_cache_element( &
1585 : estimate_to_store_int, 6, &
1586 : maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
1587 : use_disk_storage)
1588 : ELSE
1589 : CALL hfx_get_single_cache_element( &
1590 : estimate_to_store_int, 6, &
1591 : maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, &
1592 : use_disk_storage)
1593 : END IF
1594 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
1595 : IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE
1596 : nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1
1597 : buffer_left = nints
1598 : buffer_start = 1
1599 : IF (.NOT. use_disk_storage) THEN
1600 : neris_incore = neris_incore + INT(nints, int_8)
1601 : ELSE
1602 : neris_disk = neris_disk + INT(nints, int_8)
1603 : END IF
1604 : DO WHILE (buffer_left > 0)
1605 : buffer_size = MIN(buffer_left, cache_size)
1606 : IF (.NOT. use_disk_storage) THEN
1607 : CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), &
1608 : buffer_size, nbits, &
1609 : integral_caches(nbits), &
1610 : integral_containers(nbits), &
1611 : eps_storage, pmax_entry, &
1612 : memory_parameter%actual_memory_usage, &
1613 : use_disk_storage)
1614 : ELSE
1615 : CALL hfx_get_mult_cache_elements(primitive_integrals(buffer_start), &
1616 : buffer_size, nbits, &
1617 : integral_caches_disk(nbits), &
1618 : integral_containers_disk(nbits), &
1619 : eps_storage, pmax_entry, &
1620 : memory_parameter%actual_memory_usage_disk, &
1621 : use_disk_storage)
1622 : END IF
1623 : buffer_left = buffer_left - buffer_size
1624 : buffer_start = buffer_start + buffer_size
1625 : END DO
1626 : END IF
1627 : !! Calculate integrals if we run out of buffer or the geometry did change
1628 : IF (my_geo_change .OR. buffer_overflow) THEN
1629 :
1630 : max_contraction_val = max_contraction(iset, iatom)* &
1631 : max_contraction(jset, jatom)* &
1632 : max_contraction(kset, katom)* &
1633 : max_contraction(lset, latom)*pmax_entry
1634 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
1635 : tmp_R_2 => radii_pgf(:, :, lset, kset, lkind, kkind)
1636 : tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind)
1637 : tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind)
1638 :
1639 : CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), &
1640 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
1641 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
1642 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1643 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
1644 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
1645 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
1646 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
1647 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
1648 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
1649 : primitive_integrals, &
1650 : potential_parameter, &
1651 : actual_x_data%neighbor_cells, screen_coeffs_set(jset, iset, jkind, ikind)%x, &
1652 : screen_coeffs_set(lset, kset, lkind, kkind)%x, eps_schwarz, &
1653 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
1654 : log10_pmax, log10_eps_schwarz, &
1655 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
1656 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
1657 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
1658 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
1659 : sphi_a_ext_set, &
1660 : sphi_b_ext_set, &
1661 : sphi_c_ext_set, &
1662 : sphi_d_ext_set, &
1663 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
1664 : nimages, do_periodic, p_work)
1665 :
1666 : nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)
1667 : neris_total = neris_total + nints
1668 : nprim_ints = nprim_ints + neris_tmp
1669 :
1670 : !! Compress the array for storage
1671 : spherical_estimate = 0.0_dp
1672 : DO i = 1, nints
1673 : spherical_estimate = MAX(spherical_estimate, ABS(primitive_integrals(i)))
1674 : END DO
1675 :
1676 : IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate)
1677 : estimate_to_store_int = EXPONENT(spherical_estimate)
1678 : estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
1679 :
1680 : IF (.NOT. buffer_overflow .AND. my_geo_change) THEN
1681 : IF (.NOT. use_disk_storage) THEN
1682 : CALL hfx_add_single_cache_element( &
1683 : estimate_to_store_int, 6, &
1684 : maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
1685 : use_disk_storage, max_val_memory)
1686 : ELSE
1687 : CALL hfx_add_single_cache_element( &
1688 : estimate_to_store_int, 6, &
1689 : maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, &
1690 : use_disk_storage)
1691 : END IF
1692 : END IF
1693 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
1694 : IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE
1695 : IF (.NOT. buffer_overflow) THEN
1696 : nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1
1697 :
1698 : ! In case of a tight eps_storage threshold the number of significant
1699 : ! bits in the integer number NINT(value*pmax_entry/eps_storage) may
1700 : ! exceed the width of the storage element. As the compression algorithm
1701 : ! is designed for IEEE 754 double precision numbers, a 64-bit signed
1702 : ! integer variable which is used to store the result of this float-to-
1703 : ! integer conversion (we have no wish to use more memory for storing
1704 : ! compressed ERIs than it is needed for uncompressed ERIs) may overflow.
1705 : ! Abort with a meaningful message when it happens.
1706 : !
1707 : ! The magic number 63 stands for the number of magnitude bits
1708 : ! (64 bits minus one sign bit).
1709 : IF (nbits > 63) THEN
1710 : WRITE (eps_schwarz_min_str, '(ES10.3E2)') &
1711 : spherical_estimate*pmax_entry/ &
1712 : (SET_EXPONENT(1.0_dp, 63)*memory_parameter%eps_storage_scaling)
1713 :
1714 : WRITE (eps_scaling_str, '(ES10.3E2)') &
1715 : spherical_estimate*pmax_entry/(SET_EXPONENT(1.0_dp, 63)*eps_schwarz)
1716 :
1717 : CALL cp_abort(__LOCATION__, &
1718 : "Overflow during ERI's compression. Please use a larger "// &
1719 : "EPS_SCHWARZ threshold (above "//TRIM(ADJUSTL(eps_schwarz_min_str))// &
1720 : ") or increase the EPS_STORAGE_SCALING factor above "// &
1721 : TRIM(ADJUSTL(eps_scaling_str))//".")
1722 : END IF
1723 :
1724 : buffer_left = nints
1725 : buffer_start = 1
1726 : IF (.NOT. use_disk_storage) THEN
1727 : neris_incore = neris_incore + INT(nints, int_8)
1728 : ELSE
1729 : neris_disk = neris_disk + INT(nints, int_8)
1730 : END IF
1731 : DO WHILE (buffer_left > 0)
1732 : buffer_size = MIN(buffer_left, CACHE_SIZE)
1733 : IF (.NOT. use_disk_storage) THEN
1734 : CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start), &
1735 : buffer_size, nbits, &
1736 : integral_caches(nbits), &
1737 : integral_containers(nbits), &
1738 : eps_storage, pmax_entry, &
1739 : memory_parameter%actual_memory_usage, &
1740 : use_disk_storage)
1741 : ELSE
1742 : CALL hfx_add_mult_cache_elements(primitive_integrals(buffer_start), &
1743 : buffer_size, nbits, &
1744 : integral_caches_disk(nbits), &
1745 : integral_containers_disk(nbits), &
1746 : eps_storage, pmax_entry, &
1747 : memory_parameter%actual_memory_usage_disk, &
1748 : use_disk_storage)
1749 : END IF
1750 : buffer_left = buffer_left - buffer_size
1751 : buffer_start = buffer_start + buffer_size
1752 : END DO
1753 : ELSE
1754 : !! In order to be consistent with in-core part, round all the eris wrt. eps_schwarz
1755 : DO i = 1, nints
1756 : primitive_integrals(i) = primitive_integrals(i)*pmax_entry
1757 : IF (ABS(primitive_integrals(i)) > eps_storage) THEN
1758 : primitive_integrals(i) = ANINT(primitive_integrals(i)/eps_storage, dp)*eps_storage/pmax_entry
1759 : ELSE
1760 : primitive_integrals(i) = 0.0_dp
1761 : END IF
1762 : END DO
1763 : END IF
1764 : END IF
1765 : !!! DEBUG, print out primitive integrals and indices. Only works serial no OMP !!!
1766 : IF (.FALSE.) THEN
1767 : CALL print_integrals( &
1768 : iatom, jatom, katom, latom, shm_set_offset, shm_atomic_block_offset, &
1769 : iset, jset, kset, lset, nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), primitive_integrals)
1770 : END IF
1771 : IF (.NOT. is_anti_symmetric) THEN
1772 : !! Update Kohn-Sham matrix
1773 : CALL update_fock_matrix( &
1774 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1775 : fac, symm_fac, full_density_alpha(:, 1), full_ks_alpha(:, 1), &
1776 : primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, &
1777 : kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, &
1778 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, &
1779 : atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac)
1780 : IF (.NOT. treat_lsd_in_core) THEN
1781 : IF (my_nspins == 2) THEN
1782 : CALL update_fock_matrix( &
1783 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1784 : fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), &
1785 : primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, &
1786 : kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, &
1787 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, &
1788 : atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac)
1789 : END IF
1790 : END IF
1791 : ELSE
1792 : !! Update Kohn-Sham matrix
1793 : CALL update_fock_matrix_as( &
1794 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1795 : fac, symm_fac, full_density_alpha(:, 1), full_ks_alpha(:, 1), &
1796 : primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, &
1797 : kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, &
1798 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, &
1799 : atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac)
1800 : IF (.NOT. treat_lsd_in_core) THEN
1801 : IF (my_nspins == 2) THEN
1802 : CALL update_fock_matrix_as( &
1803 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1804 : fac, symm_fac, full_density_beta(:, 1), full_ks_beta(:, 1), &
1805 : primitive_integrals, pbd_buf, pbc_buf, pad_buf, pac_buf, kbd_buf, &
1806 : kbc_buf, kad_buf, kac_buf, iatom, jatom, katom, latom, &
1807 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, offset_ac_set, &
1808 : atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, atomic_offset_ac)
1809 : END IF
1810 : END IF
1811 : END IF
1812 : END IF if_use_libint
1813 : END DO ! i_set_list_kl
1814 : END DO ! i_set_list_ij
1815 : IF (do_disk_storage) THEN
1816 : buffer_overflow = .TRUE.
1817 : END IF
1818 : END DO !i_list_ij
1819 : END DO !i_list_kl
1820 : END DO atomic_blocks
1821 : bintime_stop = m_walltime()
1822 : IF (my_geo_change) THEN
1823 : distribution_energy%time_first_scf = bintime_stop - bintime_start
1824 : ELSE
1825 : distribution_energy%time_other_scf = &
1826 : distribution_energy%time_other_scf + bintime_stop - bintime_start
1827 : END IF
1828 : !$OMP MASTER
1829 : CALL timestop(handle_bin)
1830 : !$OMP END MASTER
1831 : END DO !bin
1832 :
1833 : !!!
1834 : ! IF running both libint and libGint, uses the values from libint to check libGint values
1835 : IF (use_libGint .and. use_libint) THEN
1836 : ! barrier needed after allocate full_ks_alpha(beta)_from_gpu before libGint get K
1837 : ! to feed libgint a valid position for full_ks_alpha(beta)_from_gpu
1838 : ALLOCATE (full_ks_alpha_from_gpu, Mold=full_ks_alpha)
1839 : !$OMP BARRIER
1840 : IF (use_only_alpha_spin) THEN
1841 : call libGint_get_fock_matrix(full_ks_alpha_from_gpu)
1842 : ELSE
1843 : ALLOCATE (full_ks_beta_from_gpu, Mold=full_ks_beta)
1844 : !$OMP BARRIER
1845 : CALL libGint_get_fock_matrix(full_ks_alpha_from_gpu, full_ks_beta_from_gpu)
1846 : END IF
1847 :
1848 : !$OMP BARRIER
1849 : !$omp single
1850 : max_abs_delta_KS = 0.0_dp
1851 : DO i = 1, SIZE(full_ks_alpha)
1852 : max_abs_delta_KS = MAX(max_abs_delta_ks, ABS(full_ks_alpha(i, 1) - full_ks_alpha_from_gpu(i, 1)))
1853 : END DO
1854 : IF (iw > 0) THEN
1855 : WRITE (UNIT=iw, FMT=*) " ------------------------------------ "
1856 : WRITE (UNIT=iw, FMT=*) " | rank ", para_env%mepos
1857 : WRITE (UNIT=iw, FMT=*) " | max abs dif: ", max_abs_delta_KS
1858 : WRITE (UNIT=iw, FMT=*) " ------------------------------------ "
1859 : END IF
1860 :
1861 : ! Print all F if delta is too large or (suspiciously) too small
1862 : IF (iw > 0) THEN
1863 : if (max_abs_delta_KS > 100.0_dp*eps_schwarz .or. max_abs_delta_KS < 1.0E-20_dp) THEN
1864 : WRITE (UNIT=iw, FMT=*) " K from CPU | K from GPU | Delta | Ratio "
1865 : DO i = 1, SIZE(full_ks_alpha)
1866 : write (UNIT=iw, FMT=*) i, full_ks_alpha(i, 1), full_ks_alpha_from_gpu(i, 1), &
1867 : full_ks_alpha(i, 1) - full_ks_alpha_from_gpu(i, 1), &
1868 : full_ks_alpha(i, 1)/full_ks_alpha_from_gpu(i, 1)
1869 : END DO
1870 : END IF
1871 : END IF
1872 :
1873 : IF (.NOT. use_only_alpha_spin) THEN
1874 : max_abs_delta_KS = 0.0_dp
1875 : DO i = 1, size(full_ks_beta)
1876 : max_abs_delta_KS = max(max_abs_delta_ks, abs(full_ks_beta(i, 1) - full_ks_beta_from_gpu(i, 1)))
1877 : END DO
1878 : IF (iw > 0) THEN
1879 : WRITE (UNIT=iw, FMT=*) " ------------------------------------ "
1880 : WRITE (UNIT=iw, FMT=*) " | rank ", para_env%mepos
1881 : WRITE (UNIT=iw, FMT=*) " | max abs dif B: ", max_abs_delta_KS
1882 : WRITE (UNIT=iw, FMT=*) " ------------------------------------ "
1883 : END IF
1884 :
1885 : ! Print all F if delta is too large or (suspiciously) too small
1886 : IF (max_abs_delta_KS > 100.0_dp*eps_schwarz .or. max_abs_delta_KS < 1.0E-18_dp) THEN
1887 : IF (iw > 0) THEN
1888 : WRITE (UNIT=iw, FMT=*) " BETA K from CPU | K from GPU | Delta | Ratio "
1889 : DO i = 1, size(full_ks_beta)
1890 : WRITE (UNIT=iw, FMT=*) i, full_ks_beta(i, 1), full_ks_beta_from_gpu(i, 1), &
1891 : full_ks_beta(i, 1) - full_ks_beta_from_gpu(i, 1), &
1892 : full_ks_beta(i, 1)/full_ks_beta_from_gpu(i, 1)
1893 : END DO
1894 : END IF
1895 : END IF
1896 : END IF
1897 : !$omp end single
1898 : END IF
1899 :
1900 : ! Copy the KS matrix from device to host.
1901 : IF (use_libGint .AND. (.NOT. use_libint)) THEN
1902 : IF (use_only_alpha_spin) THEN
1903 : CALL libGint_get_fock_matrix(full_ks_alpha)
1904 : ELSE
1905 : CALL libGint_get_fock_matrix(full_ks_alpha, full_ks_beta)
1906 : END IF
1907 : END IF
1908 :
1909 : !$OMP MASTER
1910 : logger => cp_get_default_logger()
1911 : do_print_load_balance_info = .FALSE.
1912 : do_print_load_balance_info = BTEST(cp_print_key_should_output(logger%iter_info, hfx_section, &
1913 : "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO"), cp_p_file)
1914 : !$OMP END MASTER
1915 : !$OMP BARRIER
1916 : IF (do_print_load_balance_info) THEN
1917 : iw = -1
1918 : !$OMP MASTER
1919 : iw = cp_print_key_unit_nr(logger, hfx_section, "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO", &
1920 : extension=".scfLog")
1921 : !$OMP END MASTER
1922 :
1923 : CALL collect_load_balance_info(para_env, actual_x_data, iw, n_threads, i_thread, &
1924 : hfx_do_eval_energy)
1925 :
1926 : !$OMP MASTER
1927 : CALL cp_print_key_finished_output(iw, logger, hfx_section, &
1928 : "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO")
1929 : !$OMP END MASTER
1930 : END IF
1931 :
1932 : !$OMP BARRIER
1933 : !$OMP MASTER
1934 : CALL m_memory(memsize_after)
1935 : !$OMP END MASTER
1936 : !$OMP BARRIER
1937 :
1938 : DEALLOCATE (primitive_integrals)
1939 : !$OMP BARRIER
1940 : !! Get some number about ERIS
1941 : !$OMP ATOMIC
1942 : shm_neris_total = shm_neris_total + neris_total
1943 : !$OMP ATOMIC
1944 : shm_neris_onthefly = shm_neris_onthefly + neris_onthefly
1945 : !$OMP ATOMIC
1946 : shm_nprim_ints = shm_nprim_ints + nprim_ints
1947 :
1948 : storage_counter_integrals = memory_parameter%actual_memory_usage* &
1949 : memory_parameter%cache_size
1950 : stor_count_int_disk = memory_parameter%actual_memory_usage_disk* &
1951 : memory_parameter%cache_size
1952 : stor_count_max_val = max_val_memory*memory_parameter%cache_size
1953 : !$OMP ATOMIC
1954 : shm_storage_counter_integrals = shm_storage_counter_integrals + storage_counter_integrals
1955 : !$OMP ATOMIC
1956 : shm_stor_count_int_disk = shm_stor_count_int_disk + stor_count_int_disk
1957 : !$OMP ATOMIC
1958 : shm_neris_incore = shm_neris_incore + neris_incore
1959 : !$OMP ATOMIC
1960 : shm_neris_disk = shm_neris_disk + neris_disk
1961 : !$OMP ATOMIC
1962 : shm_stor_count_max_val = shm_stor_count_max_val + stor_count_max_val
1963 : !$OMP BARRIER
1964 :
1965 : ! ** Calculate how much memory has already been used (might be needed for in-core forces
1966 : !$OMP MASTER
1967 : shm_mem_compression_counter = 0
1968 : DO i = 1, n_threads
1969 : !$OMP ATOMIC READ
1970 : tmp_i4 = x_data(irep, i)%memory_parameter%actual_memory_usage
1971 : shm_mem_compression_counter = shm_mem_compression_counter + &
1972 : tmp_i4*memory_parameter%cache_size
1973 : END DO
1974 : !$OMP END MASTER
1975 : !$OMP BARRIER
1976 : actual_x_data%memory_parameter%final_comp_counter_energy = shm_mem_compression_counter
1977 :
1978 : !$OMP MASTER
1979 : !! Calculate the exchange energies from the Kohn-Sham matrix. Before we can go on, we have to symmetrize.
1980 : ene_x_aa = 0.0_dp
1981 : ene_x_bb = 0.0_dp
1982 :
1983 : mb_size_p = shm_block_offset(ncpu + 1)/1024/128
1984 : mb_size_f = shm_block_offset(ncpu + 1)/1024/128
1985 : IF (.NOT. treat_lsd_in_core) THEN
1986 : IF (my_nspins == 2) THEN
1987 : mb_size_f = mb_size_f*2
1988 : mb_size_p = mb_size_p*2
1989 : END IF
1990 : END IF
1991 : !! size of primitive_integrals(not shared)
1992 : mb_size_buffers = INT(nsgf_max, int_8)**4*n_threads
1993 : !! fock density buffers (not shared)
1994 : mb_size_buffers = mb_size_buffers + INT(nsgf_max, int_8)**2*n_threads
1995 : subtr_size_mb = subtr_size_mb + 8_int_8*nsgf_max**2*n_threads
1996 : !! size of screening functions (shared)
1997 : mb_size_buffers = mb_size_buffers + max_pgf**2*max_set**2*nkind**2 &
1998 : + max_set**2*nkind**2 &
1999 : + nkind**2 &
2000 : + max_pgf**2*max_set**2*nkind**2
2001 : !! is_assoc (shared)
2002 : mb_size_buffers = mb_size_buffers + natom**2
2003 : ! ** pmax_atom (shared)
2004 : IF (do_p_screening) THEN
2005 : mb_size_buffers = mb_size_buffers + natom**2
2006 : END IF
2007 : IF (screening_parameter%do_p_screening_forces) THEN
2008 : IF (memory_parameter%treat_forces_in_core) THEN
2009 : mb_size_buffers = mb_size_buffers + natom**2
2010 : END IF
2011 : END IF
2012 : ! ** Initial P only MAX(alpha,beta) (shared)
2013 : IF (do_p_screening .OR. screening_parameter%do_p_screening_forces) THEN
2014 : mb_size_buffers = mb_size_buffers + memory_parameter%size_p_screen
2015 : END IF
2016 : ! ** In core forces require their own initial P
2017 : IF (screening_parameter%do_p_screening_forces) THEN
2018 : IF (memory_parameter%treat_forces_in_core) THEN
2019 : mb_size_buffers = mb_size_buffers + memory_parameter%size_p_screen
2020 : END IF
2021 : END IF
2022 :
2023 : !! mb
2024 : mb_size_buffers = mb_size_buffers/1024/128
2025 :
2026 : afac = 1.0_dp
2027 : IF (is_anti_symmetric) afac = -1.0_dp
2028 : CALL timestop(handle_main)
2029 : ene_x_aa_diag = 0.0_dp
2030 : ene_x_bb_diag = 0.0_dp
2031 : DO iatom = 1, natom
2032 : ikind = kind_of(iatom)
2033 : nseta = basis_parameter(ikind)%nset
2034 : nsgfa => basis_parameter(ikind)%nsgf
2035 : jatom = iatom
2036 : jkind = kind_of(jatom)
2037 : nsetb = basis_parameter(jkind)%nset
2038 : nsgfb => basis_parameter(jkind)%nsgf
2039 : act_atomic_block_offset = shm_atomic_block_offset(jatom, iatom)
2040 : DO img = 1, nkimages
2041 : DO iset = 1, nseta
2042 : DO jset = 1, nsetb
2043 : act_set_offset = shm_set_offset(jset, iset, jkind, ikind)
2044 : i = act_set_offset + act_atomic_block_offset - 1
2045 : DO ma = 1, nsgfa(iset)
2046 : j = shm_set_offset(iset, jset, jkind, ikind) + act_atomic_block_offset - 1 + ma - 1
2047 : DO mb = 1, nsgfb(jset)
2048 : IF (i > j) THEN
2049 : full_ks_alpha(i, img) = (full_ks_alpha(i, img) + full_ks_alpha(j, img)*afac)
2050 : full_ks_alpha(j, img) = full_ks_alpha(i, img)*afac
2051 : IF (.NOT. treat_lsd_in_core .AND. my_nspins == 2) THEN
2052 : full_ks_beta(i, img) = (full_ks_beta(i, img) + full_ks_beta(j, img)*afac)
2053 : full_ks_beta(j, img) = full_ks_beta(i, img)*afac
2054 : END IF
2055 : END IF
2056 : ! ** For adiabatically rescaled functionals we need the energy coming from the diagonal elements
2057 : IF (i == j) THEN
2058 : ene_x_aa_diag = ene_x_aa_diag + full_ks_alpha(i, img)*full_density_alpha(i, img)
2059 : IF (.NOT. treat_lsd_in_core .AND. my_nspins == 2) THEN
2060 : ene_x_bb_diag = ene_x_bb_diag + full_ks_beta(i, img)*full_density_beta(i, img)
2061 : END IF
2062 : END IF
2063 : i = i + 1
2064 : j = j + nsgfa(iset)
2065 : END DO
2066 : END DO
2067 : END DO
2068 : END DO
2069 : END DO
2070 : END DO
2071 :
2072 : CALL para_env%sync()
2073 : afac = 1.0_dp
2074 : IF (is_anti_symmetric) afac = 0._dp
2075 : IF (distribute_fock_matrix) THEN
2076 : !! Distribute the current KS-matrix to all the processes
2077 : CALL timeset(routineN//"_dist_KS", handle_dist_ks)
2078 : DO img = 1, nkimages
2079 : CALL distribute_ks_matrix(para_env, full_ks_alpha(:, img), ks_matrix(ispin, img)%matrix, shm_number_of_p_entries, &
2080 : shm_block_offset, kind_of, basis_parameter, &
2081 : off_diag_fac=0.5_dp, diag_fac=afac)
2082 : END DO
2083 :
2084 : NULLIFY (full_ks_alpha)
2085 : DEALLOCATE (shm_master_x_data%full_ks_alpha)
2086 : IF (.NOT. treat_lsd_in_core) THEN
2087 : IF (my_nspins == 2) THEN
2088 : DO img = 1, nkimages
2089 : CALL distribute_ks_matrix(para_env, full_ks_beta(:, img), ks_matrix(2, img)%matrix, shm_number_of_p_entries, &
2090 : shm_block_offset, kind_of, basis_parameter, &
2091 : off_diag_fac=0.5_dp, diag_fac=afac)
2092 : END DO
2093 : NULLIFY (full_ks_beta)
2094 : DEALLOCATE (shm_master_x_data%full_ks_beta)
2095 : END IF
2096 : END IF
2097 : CALL timestop(handle_dist_ks)
2098 : END IF
2099 :
2100 : IF (distribute_fock_matrix) THEN
2101 : !! ** Calculate the exchange energy
2102 : ene_x_aa = 0.0_dp
2103 : DO img = 1, nkimages
2104 : CALL dbcsr_dot_threadsafe(ks_matrix(ispin, img)%matrix, rho_ao(ispin, img)%matrix, etmp)
2105 : ene_x_aa = ene_x_aa + etmp
2106 : END DO
2107 : !for ADMMS, we need the exchange matrix k(d) for both spins
2108 : IF (dft_control%do_admm) THEN
2109 : CPASSERT(nkimages == 1)
2110 : CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, ks_matrix(ispin, 1)%matrix, &
2111 : name="HF exch. part of matrix_ks_aux_fit for ADMMS")
2112 : END IF
2113 :
2114 : ene_x_bb = 0.0_dp
2115 : IF (my_nspins == 2 .AND. .NOT. treat_lsd_in_core) THEN
2116 : DO img = 1, nkimages
2117 : CALL dbcsr_dot_threadsafe(ks_matrix(2, img)%matrix, rho_ao(2, img)%matrix, etmp)
2118 : ene_x_bb = ene_x_bb + etmp
2119 : END DO
2120 : !for ADMMS, we need the exchange matrix k(d) for both spins
2121 : IF (dft_control%do_admm) THEN
2122 : CPASSERT(nkimages == 1)
2123 : CALL dbcsr_copy(matrix_ks_aux_fit_hfx(2)%matrix, ks_matrix(2, 1)%matrix, &
2124 : name="HF exch. part of matrix_ks_aux_fit for ADMMS")
2125 : END IF
2126 : END IF
2127 :
2128 : !! Update energy type
2129 : ehfx = 0.5_dp*(ene_x_aa + ene_x_bb)
2130 : ELSE
2131 : ! ** It is easier to correct the following expression by the diagonal energy contribution,
2132 : ! ** than explicitly going throuhg the diagonal elements
2133 : DO img = 1, nkimages
2134 : DO pa = 1, SIZE(full_ks_alpha, 1)
2135 : ene_x_aa = ene_x_aa + full_ks_alpha(pa, img)*full_density_alpha(pa, img)
2136 : END DO
2137 : END DO
2138 : ! ** Now correct
2139 : ene_x_aa = (ene_x_aa + ene_x_aa_diag)*0.5_dp
2140 : IF (my_nspins == 2) THEN
2141 : DO img = 1, nkimages
2142 : DO pa = 1, SIZE(full_ks_beta, 1)
2143 : ene_x_bb = ene_x_bb + full_ks_beta(pa, img)*full_density_beta(pa, img)
2144 : END DO
2145 : END DO
2146 : ! ** Now correct
2147 : ene_x_bb = (ene_x_bb + ene_x_bb_diag)*0.5_dp
2148 : END IF
2149 : CALL para_env%sum(ene_x_aa)
2150 : IF (my_nspins == 2) CALL para_env%sum(ene_x_bb)
2151 : ehfx = 0.5_dp*(ene_x_aa + ene_x_bb)
2152 : END IF
2153 :
2154 : !! Print some memeory information if this is the first step
2155 : IF (my_geo_change) THEN
2156 : tmp_i8(1:8) = [shm_storage_counter_integrals, shm_neris_onthefly, shm_neris_incore, shm_neris_disk, &
2157 : shm_neris_total, shm_stor_count_int_disk, shm_nprim_ints, shm_stor_count_max_val]
2158 : CALL para_env%sum(tmp_i8)
2159 : shm_storage_counter_integrals = tmp_i8(1)
2160 : shm_neris_onthefly = tmp_i8(2)
2161 : shm_neris_incore = tmp_i8(3)
2162 : shm_neris_disk = tmp_i8(4)
2163 : shm_neris_total = tmp_i8(5)
2164 : shm_stor_count_int_disk = tmp_i8(6)
2165 : shm_nprim_ints = tmp_i8(7)
2166 : shm_stor_count_max_val = tmp_i8(8)
2167 : CALL para_env%max(memsize_after)
2168 : mem_eris = (shm_storage_counter_integrals + 128*1024 - 1)/1024/128
2169 : compression_factor = REAL(shm_neris_incore, dp)/REAL(shm_storage_counter_integrals, dp)
2170 : mem_eris_disk = (shm_stor_count_int_disk + 128*1024 - 1)/1024/128
2171 : compression_factor_disk = REAL(shm_neris_disk, dp)/REAL(shm_stor_count_int_disk, dp)
2172 : mem_max_val = (shm_stor_count_max_val + 128*1024 - 1)/1024/128
2173 :
2174 : IF (shm_neris_incore == 0) THEN
2175 : mem_eris = 0
2176 : compression_factor = 0.0_dp
2177 : END IF
2178 : IF (shm_neris_disk == 0) THEN
2179 : mem_eris_disk = 0
2180 : compression_factor_disk = 0.0_dp
2181 : END IF
2182 :
2183 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2184 : extension=".scfLog")
2185 : IF (iw > 0) THEN
2186 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2187 : "HFX_MEM_INFO| Number of cart. primitive ERI's calculated: ", shm_nprim_ints
2188 :
2189 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2190 : "HFX_MEM_INFO| Number of sph. ERI's calculated: ", shm_neris_total
2191 :
2192 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2193 : "HFX_MEM_INFO| Number of sph. ERI's stored in-core: ", shm_neris_incore
2194 :
2195 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2196 : "HFX_MEM_INFO| Number of sph. ERI's stored on disk: ", shm_neris_disk
2197 :
2198 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2199 : "HFX_MEM_INFO| Number of sph. ERI's calculated on the fly: ", shm_neris_onthefly
2200 :
2201 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2202 : "HFX_MEM_INFO| Total memory consumption ERI's RAM [MiB]: ", mem_eris
2203 :
2204 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2205 : "HFX_MEM_INFO| Whereof max-vals [MiB]: ", mem_max_val
2206 :
2207 : WRITE (UNIT=iw, FMT="((T3,A,T60,F21.2))") &
2208 : "HFX_MEM_INFO| Total compression factor ERI's RAM: ", compression_factor
2209 :
2210 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2211 : "HFX_MEM_INFO| Total memory consumption ERI's disk [MiB]: ", mem_eris_disk
2212 :
2213 : WRITE (UNIT=iw, FMT="((T3,A,T60,F21.2))") &
2214 : "HFX_MEM_INFO| Total compression factor ERI's disk: ", compression_factor_disk
2215 :
2216 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2217 : "HFX_MEM_INFO| Size of density/Fock matrix [MiB]: ", 2_int_8*mb_size_p
2218 :
2219 : IF (do_periodic) THEN
2220 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2221 : "HFX_MEM_INFO| Size of buffers [MiB]: ", mb_size_buffers
2222 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2223 : "HFX_MEM_INFO| Number of periodic image cells considered: ", SIZE(shm_master_x_data%neighbor_cells)
2224 : ELSE
2225 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
2226 : "HFX_MEM_INFO| Size of buffers [MiB]: ", mb_size_buffers
2227 : END IF
2228 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21),/)") &
2229 : "HFX_MEM_INFO| Est. max. program size after HFX [MiB]:", memsize_after/(1024*1024)
2230 : CALL m_flush(iw)
2231 : END IF
2232 :
2233 : CALL cp_print_key_finished_output(iw, logger, hfx_section, &
2234 : "HF_INFO")
2235 : END IF
2236 : !$OMP END MASTER
2237 :
2238 : !! flush caches if the geometry changed
2239 : IF (do_dynamic_load_balancing) THEN
2240 : my_bin_size = SIZE(actual_x_data%distribution_energy)
2241 : ELSE
2242 : my_bin_size = 1
2243 : END IF
2244 :
2245 : IF (my_geo_change) THEN
2246 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
2247 : DO bin = 1, my_bin_size
2248 : maxval_cache => actual_x_data%store_ints%maxval_cache(bin)
2249 : maxval_container => actual_x_data%store_ints%maxval_container(bin)
2250 : integral_caches => actual_x_data%store_ints%integral_caches(:, bin)
2251 : integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
2252 : CALL hfx_flush_last_cache(bits_max_val, maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
2253 : .FALSE.)
2254 : DO i = 1, 64
2255 : CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
2256 : memory_parameter%actual_memory_usage, .FALSE.)
2257 : END DO
2258 : END DO
2259 : END IF
2260 : END IF
2261 : !! reset all caches except we calculate all on the fly
2262 : IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
2263 : DO bin = 1, my_bin_size
2264 : maxval_cache => actual_x_data%store_ints%maxval_cache(bin)
2265 : maxval_container => actual_x_data%store_ints%maxval_container(bin)
2266 : integral_caches => actual_x_data%store_ints%integral_caches(:, bin)
2267 : integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
2268 :
2269 : CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_parameter%actual_memory_usage, .FALSE.)
2270 : DO i = 1, 64
2271 : CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
2272 : memory_parameter%actual_memory_usage, &
2273 : .FALSE.)
2274 : END DO
2275 : END DO
2276 : END IF
2277 :
2278 : !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here
2279 : !$OMP CRITICAL(hfxenergy_out_critical)
2280 : IF (do_disk_storage) THEN
2281 : !! flush caches if the geometry changed
2282 : IF (my_geo_change) THEN
2283 : CALL hfx_flush_last_cache(bits_max_val, maxval_cache_disk, maxval_container_disk, &
2284 : memory_parameter%actual_memory_usage_disk, .TRUE.)
2285 : DO i = 1, 64
2286 : CALL hfx_flush_last_cache(i, integral_caches_disk(i), integral_containers_disk(i), &
2287 : memory_parameter%actual_memory_usage_disk, .TRUE.)
2288 : END DO
2289 : END IF
2290 : !! reset all caches except we calculate all on the fly
2291 : CALL hfx_reset_cache_and_container(maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, &
2292 : do_disk_storage)
2293 : DO i = 1, 64
2294 : CALL hfx_reset_cache_and_container(integral_caches_disk(i), integral_containers_disk(i), &
2295 : memory_parameter%actual_memory_usage_disk, do_disk_storage)
2296 : END DO
2297 : END IF
2298 : !$OMP END CRITICAL(hfxenergy_out_critical)
2299 : !$OMP BARRIER
2300 : !! Clean up
2301 : DEALLOCATE (last_sgf_global)
2302 : !$OMP MASTER
2303 : DEALLOCATE (full_density_alpha)
2304 : IF (.NOT. treat_lsd_in_core) THEN
2305 : IF (my_nspins == 2) THEN
2306 : DEALLOCATE (full_density_beta)
2307 : END IF
2308 : END IF
2309 : IF (do_dynamic_load_balancing) THEN
2310 : DEALLOCATE (shm_master_x_data%task_list)
2311 : END IF
2312 : !$OMP END MASTER
2313 : DEALLOCATE (pbd_buf, pbc_buf, pad_buf, pac_buf)
2314 : DEALLOCATE (kbd_buf, kbc_buf, kad_buf, kac_buf)
2315 : DEALLOCATE (set_list_ij, set_list_kl)
2316 :
2317 : DO i = 1, max_pgf**2
2318 : DEALLOCATE (pgf_list_ij(i)%image_list)
2319 : DEALLOCATE (pgf_list_kl(i)%image_list)
2320 : END DO
2321 :
2322 : DEALLOCATE (pgf_list_ij)
2323 : DEALLOCATE (pgf_list_kl)
2324 : DEALLOCATE (pgf_product_list)
2325 :
2326 : DEALLOCATE (max_contraction, kind_of)
2327 :
2328 : DEALLOCATE (ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
2329 :
2330 : DEALLOCATE (nimages)
2331 :
2332 : !$OMP BARRIER
2333 : !$OMP END PARALLEL
2334 :
2335 42353 : CALL timestop(handle)
2336 87840122 : END SUBROUTINE integrate_four_center
2337 :
2338 : ! **************************************************************************************************
2339 : !> \brief calculates two-electron integrals of a quartet/shell using the library
2340 : !> lib_int in the periodic case
2341 : !> \param lib ...
2342 : !> \param ra ...
2343 : !> \param rb ...
2344 : !> \param rc ...
2345 : !> \param rd ...
2346 : !> \param npgfa ...
2347 : !> \param npgfb ...
2348 : !> \param npgfc ...
2349 : !> \param npgfd ...
2350 : !> \param la_min ...
2351 : !> \param la_max ...
2352 : !> \param lb_min ...
2353 : !> \param lb_max ...
2354 : !> \param lc_min ...
2355 : !> \param lc_max ...
2356 : !> \param ld_min ...
2357 : !> \param ld_max ...
2358 : !> \param nsgfa ...
2359 : !> \param nsgfb ...
2360 : !> \param nsgfc ...
2361 : !> \param nsgfd ...
2362 : !> \param sphi_a_u1 ...
2363 : !> \param sphi_a_u2 ...
2364 : !> \param sphi_a_u3 ...
2365 : !> \param sphi_b_u1 ...
2366 : !> \param sphi_b_u2 ...
2367 : !> \param sphi_b_u3 ...
2368 : !> \param sphi_c_u1 ...
2369 : !> \param sphi_c_u2 ...
2370 : !> \param sphi_c_u3 ...
2371 : !> \param sphi_d_u1 ...
2372 : !> \param sphi_d_u2 ...
2373 : !> \param sphi_d_u3 ...
2374 : !> \param zeta ...
2375 : !> \param zetb ...
2376 : !> \param zetc ...
2377 : !> \param zetd ...
2378 : !> \param primitive_integrals array of primitive_integrals
2379 : !> \param potential_parameter contains info for libint
2380 : !> \param neighbor_cells Periodic images
2381 : !> \param screen1 set based coefficients for near field screening
2382 : !> \param screen2 set based coefficients for near field screening
2383 : !> \param eps_schwarz threshold
2384 : !> \param max_contraction_val maximum multiplication factor for cart -> sph
2385 : !> \param cart_estimate maximum calculated integral value
2386 : !> \param cell cell
2387 : !> \param neris_tmp counter for calculated cart integrals
2388 : !> \param log10_pmax logarithm of initial p matrix max element
2389 : !> \param log10_eps_schwarz log of threshold
2390 : !> \param R1_pgf coefficients for radii of product distribution function
2391 : !> \param R2_pgf coefficients for radii of product distribution function
2392 : !> \param pgf1 schwarz coefficients pgf basid
2393 : !> \param pgf2 schwarz coefficients pgf basid
2394 : !> \param pgf_list_ij ...
2395 : !> \param pgf_list_kl ...
2396 : !> \param pgf_product_list ...
2397 : !> \param nsgfl_a ...
2398 : !> \param nsgfl_b ...
2399 : !> \param nsgfl_c ...
2400 : !> \param nsgfl_d ...
2401 : !> \param sphi_a_ext ...
2402 : !> \param sphi_b_ext ...
2403 : !> \param sphi_c_ext ...
2404 : !> \param sphi_d_ext ...
2405 : !> \param ee_work ...
2406 : !> \param ee_work2 ...
2407 : !> \param ee_buffer1 ...
2408 : !> \param ee_buffer2 ...
2409 : !> \param ee_primitives_tmp ...
2410 : !> \param nimages ...
2411 : !> \param do_periodic ...
2412 : !> \param p_work ...
2413 : !> \par History
2414 : !> 11.2006 created [Manuel Guidon]
2415 : !> 02.2009 completely rewritten screening part [Manuel Guidon]
2416 : !> \author Manuel Guidon
2417 : ! **************************************************************************************************
2418 8586686 : SUBROUTINE coulomb4(lib, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, &
2419 : la_min, la_max, lb_min, lb_max, &
2420 : lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, &
2421 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
2422 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
2423 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
2424 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
2425 8586686 : zeta, zetb, zetc, zetd, &
2426 8586686 : primitive_integrals, &
2427 : potential_parameter, neighbor_cells, &
2428 : screen1, screen2, eps_schwarz, max_contraction_val, &
2429 : cart_estimate, cell, neris_tmp, log10_pmax, &
2430 : log10_eps_schwarz, R1_pgf, R2_pgf, pgf1, pgf2, &
2431 : pgf_list_ij, pgf_list_kl, &
2432 : pgf_product_list, &
2433 8586686 : nsgfl_a, nsgfl_b, nsgfl_c, &
2434 8586686 : nsgfl_d, &
2435 8586686 : sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, &
2436 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
2437 : nimages, do_periodic, p_work)
2438 :
2439 : TYPE(cp_libint_t) :: lib
2440 : REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3)
2441 : INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, &
2442 : lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, &
2443 : sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, &
2444 : sphi_d_u3
2445 : REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
2446 : REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
2447 : REAL(dp), DIMENSION(1:npgfc), INTENT(IN) :: zetc
2448 : REAL(dp), DIMENSION(1:npgfd), INTENT(IN) :: zetd
2449 : REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitive_integrals
2450 : TYPE(hfx_potential_type) :: potential_parameter
2451 : TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
2452 : REAL(dp), INTENT(IN) :: screen1(2), screen2(2), eps_schwarz, &
2453 : max_contraction_val
2454 : REAL(dp) :: cart_estimate
2455 : TYPE(cell_type), POINTER :: cell
2456 : INTEGER(int_8) :: neris_tmp
2457 : REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
2458 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
2459 : POINTER :: R1_pgf, R2_pgf, pgf1, pgf2
2460 : TYPE(hfx_pgf_list), DIMENSION(*) :: pgf_list_ij, pgf_list_kl
2461 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
2462 : DIMENSION(:), INTENT(INOUT) :: pgf_product_list
2463 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
2464 : REAL(dp), INTENT(IN) :: sphi_a_ext(sphi_a_u1, sphi_a_u2, sphi_a_u3), &
2465 : sphi_b_ext(sphi_b_u1, sphi_b_u2, sphi_b_u3), sphi_c_ext(sphi_c_u1, sphi_c_u2, sphi_c_u3), &
2466 : sphi_d_ext(sphi_d_u1, sphi_d_u2, sphi_d_u3)
2467 : REAL(dp), DIMENSION(*) :: ee_work, ee_work2, ee_buffer1, &
2468 : ee_buffer2, ee_primitives_tmp
2469 : INTEGER, DIMENSION(*) :: nimages
2470 : LOGICAL, INTENT(IN) :: do_periodic
2471 : REAL(dp), DIMENSION(:), POINTER :: p_work
2472 :
2473 : INTEGER :: ipgf, jpgf, kpgf, la, lb, lc, ld, list_ij, list_kl, lpgf, max_l, ncoa, ncob, &
2474 : ncoc, ncod, nelements_ij, nelements_kl, nproducts, nsgfla, nsgflb, nsgflc, nsgfld, nsoa, &
2475 : nsob, nsoc, nsod, s_offset_a, s_offset_b, s_offset_c, s_offset_d
2476 : REAL(dp) :: EtaInv, tmp_max, ZetaInv
2477 :
2478 : CALL build_pair_list_pgf(npgfa, npgfb, pgf_list_ij, zeta, zetb, screen1, screen2, &
2479 : pgf1, R1_pgf, log10_pmax, log10_eps_schwarz, ra, rb, &
2480 : nelements_ij, &
2481 8586686 : neighbor_cells, nimages, do_periodic)
2482 : CALL build_pair_list_pgf(npgfc, npgfd, pgf_list_kl, zetc, zetd, screen2, screen1, &
2483 : pgf2, R2_pgf, log10_pmax, log10_eps_schwarz, rc, rd, &
2484 : nelements_kl, &
2485 8586686 : neighbor_cells, nimages, do_periodic)
2486 :
2487 8586686 : cart_estimate = 0.0_dp
2488 8586686 : neris_tmp = 0
2489 569844009 : primitive_integrals = 0.0_dp
2490 8586686 : max_l = la_max + lb_max + lc_max + ld_max
2491 :
2492 29539852 : DO list_ij = 1, nelements_ij
2493 20953166 : ZetaInv = pgf_list_ij(list_ij)%ZetaInv
2494 20953166 : ipgf = pgf_list_ij(list_ij)%ipgf
2495 20953166 : jpgf = pgf_list_ij(list_ij)%jpgf
2496 :
2497 106155984 : DO list_kl = 1, nelements_kl
2498 76616132 : EtaInv = pgf_list_kl(list_kl)%ZetaInv
2499 76616132 : kpgf = pgf_list_kl(list_kl)%ipgf
2500 76616132 : lpgf = pgf_list_kl(list_kl)%jpgf
2501 :
2502 : CALL build_pgf_product_list(pgf_list_ij(list_ij), pgf_list_kl(list_kl), pgf_product_list, &
2503 : nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, &
2504 76616132 : potential_parameter, max_l, do_periodic)
2505 :
2506 76616132 : s_offset_a = 0
2507 188964830 : DO la = la_min, la_max
2508 91395532 : s_offset_b = 0
2509 91395532 : ncoa = nco(la)
2510 91395532 : nsgfla = nsgfl_a(la)
2511 91395532 : nsoa = nso(la)
2512 :
2513 195344826 : DO lb = lb_min, lb_max
2514 103949294 : s_offset_c = 0
2515 103949294 : ncob = nco(lb)
2516 103949294 : nsgflb = nsgfl_b(lb)
2517 103949294 : nsob = nso(lb)
2518 :
2519 244587535 : DO lc = lc_min, lc_max
2520 140638241 : s_offset_d = 0
2521 140638241 : ncoc = nco(lc)
2522 140638241 : nsgflc = nsgfl_c(lc)
2523 140638241 : nsoc = nso(lc)
2524 :
2525 327053899 : DO ld = ld_min, ld_max
2526 186415658 : ncod = nco(ld)
2527 186415658 : nsgfld = nsgfl_d(ld)
2528 186415658 : nsod = nso(ld)
2529 :
2530 186415658 : tmp_max = 0.0_dp
2531 : CALL evaluate_eri(lib, nproducts, pgf_product_list, &
2532 : la, lb, lc, ld, &
2533 : ncoa, ncob, ncoc, ncod, &
2534 : nsgfa, nsgfb, nsgfc, nsgfd, &
2535 : primitive_integrals, &
2536 : max_contraction_val, tmp_max, eps_schwarz, &
2537 : neris_tmp, ZetaInv, EtaInv, &
2538 : s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
2539 : nsgfla, nsgflb, nsgflc, nsgfld, nsoa, nsob, nsoc, nsod, &
2540 : sphi_a_ext(1, la + 1, ipgf), &
2541 : sphi_b_ext(1, lb + 1, jpgf), &
2542 : sphi_c_ext(1, lc + 1, kpgf), &
2543 : sphi_d_ext(1, ld + 1, lpgf), &
2544 : ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
2545 186415658 : p_work)
2546 186415658 : cart_estimate = MAX(tmp_max, cart_estimate)
2547 327053899 : s_offset_d = s_offset_d + nsod*nsgfld
2548 : END DO !ld
2549 244587535 : s_offset_c = s_offset_c + nsoc*nsgflc
2550 : END DO !lc
2551 195344826 : s_offset_b = s_offset_b + nsob*nsgflb
2552 : END DO !lb
2553 168011664 : s_offset_a = s_offset_a + nsoa*nsgfla
2554 : END DO !la
2555 : END DO
2556 : END DO
2557 :
2558 8586686 : END SUBROUTINE coulomb4
2559 :
2560 : ! **************************************************************************************************
2561 : !> \brief Given a 2d index pair, this function returns a 1d index pair for
2562 : !> a symmetric upper triangle NxN matrix
2563 : !> The compiler should inline this function, therefore it appears in
2564 : !> several modules
2565 : !> \param i 2d index
2566 : !> \param j 2d index
2567 : !> \param N matrix size
2568 : !> \return ...
2569 : !> \par History
2570 : !> 03.2009 created [Manuel Guidon]
2571 : !> \author Manuel Guidon
2572 : ! **************************************************************************************************
2573 28264 : PURE FUNCTION get_1D_idx(i, j, N)
2574 : INTEGER, INTENT(IN) :: i, j
2575 : INTEGER(int_8), INTENT(IN) :: N
2576 : INTEGER(int_8) :: get_1D_idx
2577 :
2578 : INTEGER(int_8) :: min_ij
2579 :
2580 28264 : min_ij = MIN(i, j)
2581 28264 : get_1D_idx = min_ij*N + MAX(i, j) - (min_ij - 1)*min_ij/2 - N
2582 :
2583 28264 : END FUNCTION get_1D_idx
2584 :
2585 : ! **************************************************************************************************
2586 : !> \brief This routine prefetches density/fock matrix elements and stores them
2587 : !> in cache friendly arrays. These buffers are then used to update the
2588 : !> fock matrix
2589 : !> \param ma_max Size of matrix blocks
2590 : !> \param mb_max Size of matrix blocks
2591 : !> \param mc_max Size of matrix blocks
2592 : !> \param md_max Size of matrix blocks
2593 : !> \param fac multiplication factor (spin)
2594 : !> \param symm_fac multiplication factor (symmetry)
2595 : !> \param density upper triangular density matrix
2596 : !> \param ks upper triangular fock matrix
2597 : !> \param prim primitive integrals
2598 : !> \param pbd buffer that will contain P(b,d)
2599 : !> \param pbc buffer that will contain P(b,c)
2600 : !> \param pad buffer that will contain P(a,d)
2601 : !> \param pac buffer that will contain P(a,c)
2602 : !> \param kbd buffer for KS(b,d)
2603 : !> \param kbc buffer for KS(b,c)
2604 : !> \param kad buffer for KS(a,d)
2605 : !> \param kac buffer for KS(a,c)
2606 : !> \param iatom ...
2607 : !> \param jatom ...
2608 : !> \param katom ...
2609 : !> \param latom ...
2610 : !> \param iset ...
2611 : !> \param jset ...
2612 : !> \param kset ...
2613 : !> \param lset ...
2614 : !> \param offset_bd_set ...
2615 : !> \param offset_bc_set ...
2616 : !> \param offset_ad_set ...
2617 : !> \param offset_ac_set ...
2618 : !> \param atomic_offset_bd ...
2619 : !> \param atomic_offset_bc ...
2620 : !> \param atomic_offset_ad ...
2621 : !> \param atomic_offset_ac ...
2622 : !> \par History
2623 : !> 03.2009 created [Manuel Guidon]
2624 : !> \author Manuel Guidon
2625 : ! **************************************************************************************************
2626 :
2627 80800944 : SUBROUTINE update_fock_matrix(ma_max, mb_max, mc_max, md_max, &
2628 80800944 : fac, symm_fac, density, ks, prim, &
2629 : pbd, pbc, pad, pac, kbd, kbc, kad, kac, &
2630 : iatom, jatom, katom, latom, &
2631 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, &
2632 : offset_ac_set, atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, &
2633 : atomic_offset_ac)
2634 :
2635 : INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
2636 : REAL(dp), INTENT(IN) :: fac, symm_fac
2637 : REAL(dp), DIMENSION(:), INTENT(IN) :: density
2638 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: ks
2639 : REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
2640 : INTENT(IN) :: prim
2641 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: pbd, pbc, pad, pac, kbd, kbc, kad, kac
2642 : INTEGER, INTENT(IN) :: iatom, jatom, katom, latom, iset, jset, &
2643 : kset, lset
2644 : INTEGER, DIMENSION(:, :), POINTER, INTENT(IN) :: offset_bd_set, offset_bc_set, &
2645 : offset_ad_set, offset_ac_set
2646 : INTEGER, INTENT(IN) :: atomic_offset_bd, atomic_offset_bc, &
2647 : atomic_offset_ad, atomic_offset_ac
2648 :
2649 : INTEGER :: i, j, ma, mb, mc, md, offset_ac, &
2650 : offset_ad, offset_bc, offset_bd
2651 : REAL(dp) :: ki
2652 :
2653 80800944 : IF (jatom >= latom) THEN
2654 79591598 : i = 1
2655 79591598 : offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
2656 79591598 : j = offset_bd
2657 238897667 : DO md = 1, md_max
2658 533439716 : DO mb = 1, mb_max
2659 294542049 : pbd(i) = density(j)
2660 294542049 : i = i + 1
2661 453848118 : j = j + 1
2662 : END DO
2663 : END DO
2664 : ELSE
2665 1209346 : i = 1
2666 1209346 : offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
2667 3552704 : DO md = 1, md_max
2668 2343358 : j = offset_bd + md - 1
2669 8698035 : DO mb = 1, mb_max
2670 5145331 : pbd(i) = density(j)
2671 5145331 : i = i + 1
2672 7488689 : j = j + md_max
2673 : END DO
2674 : END DO
2675 : END IF
2676 80800944 : IF (jatom >= katom) THEN
2677 80800944 : i = 1
2678 80800944 : offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
2679 80800944 : j = offset_bc
2680 267469659 : DO mc = 1, mc_max
2681 602635807 : DO mb = 1, mb_max
2682 335166148 : pbc(i) = density(j)
2683 335166148 : i = i + 1
2684 521834863 : j = j + 1
2685 : END DO
2686 : END DO
2687 : ELSE
2688 0 : i = 1
2689 0 : offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
2690 0 : DO mc = 1, mc_max
2691 0 : j = offset_bc + mc - 1
2692 0 : DO mb = 1, mb_max
2693 0 : pbc(i) = density(j)
2694 0 : i = i + 1
2695 0 : j = j + mc_max
2696 : END DO
2697 : END DO
2698 : END IF
2699 80800944 : IF (iatom >= latom) THEN
2700 58153030 : i = 1
2701 58153030 : offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
2702 58153030 : j = offset_ad
2703 182420283 : DO md = 1, md_max
2704 450517005 : DO ma = 1, ma_max
2705 268096722 : pad(i) = density(j)
2706 268096722 : i = i + 1
2707 392363975 : j = j + 1
2708 : END DO
2709 : END DO
2710 : ELSE
2711 22647914 : i = 1
2712 22647914 : offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
2713 60030088 : DO md = 1, md_max
2714 37382174 : j = offset_ad + md - 1
2715 149516166 : DO ma = 1, ma_max
2716 89486078 : pad(i) = density(j)
2717 89486078 : i = i + 1
2718 126868252 : j = j + md_max
2719 : END DO
2720 : END DO
2721 : END IF
2722 80800944 : IF (iatom >= katom) THEN
2723 77070995 : i = 1
2724 77070995 : offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
2725 77070995 : j = offset_ac
2726 256969660 : DO mc = 1, mc_max
2727 650439078 : DO ma = 1, ma_max
2728 393469418 : pac(i) = density(j)
2729 393469418 : i = i + 1
2730 573368083 : j = j + 1
2731 : END DO
2732 : END DO
2733 : ELSE
2734 3729949 : i = 1
2735 3729949 : offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
2736 10499999 : DO mc = 1, mc_max
2737 6770050 : j = offset_ac + mc - 1
2738 28342867 : DO ma = 1, ma_max
2739 17842868 : pac(i) = density(j)
2740 17842868 : i = i + 1
2741 24612918 : j = j + mc_max
2742 : END DO
2743 : END DO
2744 : END IF
2745 :
2746 80800944 : CALL contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, fac*symm_fac)
2747 80800944 : IF (jatom >= latom) THEN
2748 : i = 1
2749 : j = offset_bd
2750 238897667 : DO md = 1, md_max
2751 533439716 : DO mb = 1, mb_max
2752 294542049 : ki = kbd(i)
2753 294542049 : !$OMP ATOMIC
2754 : ks(j) = ks(j) + ki
2755 294542049 : i = i + 1
2756 453848118 : j = j + 1
2757 : END DO
2758 : END DO
2759 : ELSE
2760 : i = 1
2761 3552704 : DO md = 1, md_max
2762 2343358 : j = offset_bd + md - 1
2763 8698035 : DO mb = 1, mb_max
2764 5145331 : ki = kbd(i)
2765 5145331 : !$OMP ATOMIC
2766 : ks(j) = ks(j) + ki
2767 5145331 : i = i + 1
2768 7488689 : j = j + md_max
2769 : END DO
2770 : END DO
2771 : END IF
2772 80800944 : IF (jatom >= katom) THEN
2773 : i = 1
2774 : j = offset_bc
2775 267469659 : DO mc = 1, mc_max
2776 602635807 : DO mb = 1, mb_max
2777 335166148 : ki = kbc(i)
2778 335166148 : !$OMP ATOMIC
2779 : ks(j) = ks(j) + ki
2780 335166148 : i = i + 1
2781 521834863 : j = j + 1
2782 : END DO
2783 : END DO
2784 : ELSE
2785 : i = 1
2786 0 : DO mc = 1, mc_max
2787 0 : j = offset_bc + mc - 1
2788 0 : DO mb = 1, mb_max
2789 0 : ki = kbc(i)
2790 0 : !$OMP ATOMIC
2791 : ks(j) = ks(j) + ki
2792 0 : i = i + 1
2793 0 : j = j + mc_max
2794 : END DO
2795 : END DO
2796 : END IF
2797 80800944 : IF (iatom >= latom) THEN
2798 : i = 1
2799 : j = offset_ad
2800 182420283 : DO md = 1, md_max
2801 450517005 : DO ma = 1, ma_max
2802 268096722 : ki = kad(i)
2803 268096722 : !$OMP ATOMIC
2804 : ks(j) = ks(j) + ki
2805 268096722 : i = i + 1
2806 392363975 : j = j + 1
2807 : END DO
2808 : END DO
2809 : ELSE
2810 : i = 1
2811 60030088 : DO md = 1, md_max
2812 37382174 : j = offset_ad + md - 1
2813 149516166 : DO ma = 1, ma_max
2814 89486078 : ki = kad(i)
2815 89486078 : !$OMP ATOMIC
2816 : ks(j) = ks(j) + ki
2817 89486078 : i = i + 1
2818 126868252 : j = j + md_max
2819 : END DO
2820 : END DO
2821 : END IF
2822 80800944 : IF (iatom >= katom) THEN
2823 : i = 1
2824 : j = offset_ac
2825 256969660 : DO mc = 1, mc_max
2826 650439078 : DO ma = 1, ma_max
2827 393469418 : ki = kac(i)
2828 393469418 : !$OMP ATOMIC
2829 : ks(j) = ks(j) + ki
2830 393469418 : i = i + 1
2831 573368083 : j = j + 1
2832 : END DO
2833 : END DO
2834 : ELSE
2835 : i = 1
2836 10499999 : DO mc = 1, mc_max
2837 6770050 : j = offset_ac + mc - 1
2838 28342867 : DO ma = 1, ma_max
2839 17842868 : ki = kac(i)
2840 17842868 : !$OMP ATOMIC
2841 : ks(j) = ks(j) + ki
2842 17842868 : i = i + 1
2843 24612918 : j = j + mc_max
2844 : END DO
2845 : END DO
2846 : END IF
2847 80800944 : END SUBROUTINE update_fock_matrix
2848 :
2849 : ! **************************************************************************************************
2850 : !> \brief ...
2851 : !> \param ma_max ...
2852 : !> \param mb_max ...
2853 : !> \param mc_max ...
2854 : !> \param md_max ...
2855 : !> \param fac ...
2856 : !> \param symm_fac ...
2857 : !> \param density ...
2858 : !> \param ks ...
2859 : !> \param prim ...
2860 : !> \param pbd ...
2861 : !> \param pbc ...
2862 : !> \param pad ...
2863 : !> \param pac ...
2864 : !> \param kbd ...
2865 : !> \param kbc ...
2866 : !> \param kad ...
2867 : !> \param kac ...
2868 : !> \param iatom ...
2869 : !> \param jatom ...
2870 : !> \param katom ...
2871 : !> \param latom ...
2872 : !> \param iset ...
2873 : !> \param jset ...
2874 : !> \param kset ...
2875 : !> \param lset ...
2876 : !> \param offset_bd_set ...
2877 : !> \param offset_bc_set ...
2878 : !> \param offset_ad_set ...
2879 : !> \param offset_ac_set ...
2880 : !> \param atomic_offset_bd ...
2881 : !> \param atomic_offset_bc ...
2882 : !> \param atomic_offset_ad ...
2883 : !> \param atomic_offset_ac ...
2884 : ! **************************************************************************************************
2885 5695767 : SUBROUTINE update_fock_matrix_as(ma_max, mb_max, mc_max, md_max, &
2886 5695767 : fac, symm_fac, density, ks, prim, &
2887 : pbd, pbc, pad, pac, kbd, kbc, kad, kac, &
2888 : iatom, jatom, katom, latom, &
2889 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, &
2890 : offset_ac_set, atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, &
2891 : atomic_offset_ac)
2892 :
2893 : INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
2894 : REAL(dp), INTENT(IN) :: fac, symm_fac
2895 : REAL(dp), DIMENSION(:), INTENT(IN) :: density
2896 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: ks
2897 : REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
2898 : INTENT(IN) :: prim
2899 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: pbd, pbc, pad, pac, kbd, kbc, kad, kac
2900 : INTEGER, INTENT(IN) :: iatom, jatom, katom, latom, iset, jset, &
2901 : kset, lset
2902 : INTEGER, DIMENSION(:, :), POINTER :: offset_bd_set, offset_bc_set, &
2903 : offset_ad_set, offset_ac_set
2904 : INTEGER, INTENT(IN) :: atomic_offset_bd, atomic_offset_bc, &
2905 : atomic_offset_ad, atomic_offset_ac
2906 :
2907 : INTEGER :: i, j, ma, mb, mc, md, offset_ac, &
2908 : offset_ad, offset_bc, offset_bd
2909 :
2910 5695767 : IF (jatom >= latom) THEN
2911 5682255 : i = 1
2912 5682255 : offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
2913 5682255 : j = offset_bd
2914 15565048 : DO md = 1, md_max
2915 30051681 : DO mb = 1, mb_max
2916 14486633 : pbd(i) = +density(j)
2917 14486633 : i = i + 1
2918 24369426 : j = j + 1
2919 : END DO
2920 : END DO
2921 : ELSE
2922 13512 : i = 1
2923 13512 : offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
2924 28504 : DO md = 1, md_max
2925 14992 : j = offset_bd + md - 1
2926 47164 : DO mb = 1, mb_max
2927 18660 : pbd(i) = -density(j)
2928 18660 : i = i + 1
2929 33652 : j = j + md_max
2930 : END DO
2931 : END DO
2932 : END IF
2933 5695767 : IF (jatom >= katom) THEN
2934 5695767 : i = 1
2935 5695767 : offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
2936 5695767 : j = offset_bc
2937 17488404 : DO mc = 1, mc_max
2938 34264985 : DO mb = 1, mb_max
2939 16776581 : pbc(i) = -density(j)
2940 16776581 : i = i + 1
2941 28569218 : j = j + 1
2942 : END DO
2943 : END DO
2944 : ELSE
2945 0 : i = 1
2946 0 : offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
2947 0 : DO mc = 1, mc_max
2948 0 : j = offset_bc + mc - 1
2949 0 : DO mb = 1, mb_max
2950 0 : pbc(i) = density(j)
2951 0 : i = i + 1
2952 0 : j = j + mc_max
2953 : END DO
2954 : END DO
2955 : END IF
2956 5695767 : IF (iatom >= latom) THEN
2957 4264733 : i = 1
2958 4264733 : offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
2959 4264733 : j = offset_ad
2960 12447084 : DO md = 1, md_max
2961 27986983 : DO ma = 1, ma_max
2962 15539899 : pad(i) = -density(j)
2963 15539899 : i = i + 1
2964 23722250 : j = j + 1
2965 : END DO
2966 : END DO
2967 : ELSE
2968 1431034 : i = 1
2969 1431034 : offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
2970 3146468 : DO md = 1, md_max
2971 1715434 : j = offset_ad + md - 1
2972 6978906 : DO ma = 1, ma_max
2973 3832438 : pad(i) = density(j)
2974 3832438 : i = i + 1
2975 5547872 : j = j + md_max
2976 : END DO
2977 : END DO
2978 : END IF
2979 5695767 : IF (iatom >= katom) THEN
2980 5487993 : i = 1
2981 5487993 : offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
2982 5487993 : j = offset_ac
2983 16983984 : DO mc = 1, mc_max
2984 39281281 : DO ma = 1, ma_max
2985 22297297 : pac(i) = +density(j)
2986 22297297 : i = i + 1
2987 33793288 : j = j + 1
2988 : END DO
2989 : END DO
2990 : ELSE
2991 207774 : i = 1
2992 207774 : offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
2993 504420 : DO mc = 1, mc_max
2994 296646 : j = offset_ac + mc - 1
2995 1265292 : DO ma = 1, ma_max
2996 760872 : pac(i) = -density(j)
2997 760872 : i = i + 1
2998 1057518 : j = j + mc_max
2999 : END DO
3000 : END DO
3001 : END IF
3002 :
3003 5695767 : CALL contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, fac*symm_fac)
3004 :
3005 5695767 : IF (jatom >= latom) THEN
3006 : i = 1
3007 : j = offset_bd
3008 15565048 : DO md = 1, md_max
3009 30051681 : DO mb = 1, mb_max
3010 14486633 : !$OMP ATOMIC
3011 : ks(j) = ks(j) + kbd(i)
3012 14486633 : i = i + 1
3013 24369426 : j = j + 1
3014 : END DO
3015 : END DO
3016 : ELSE
3017 : i = 1
3018 28504 : DO md = 1, md_max
3019 14992 : j = offset_bd + md - 1
3020 47164 : DO mb = 1, mb_max
3021 18660 : !$OMP ATOMIC
3022 : ks(j) = ks(j) - kbd(i)
3023 18660 : i = i + 1
3024 33652 : j = j + md_max
3025 : END DO
3026 : END DO
3027 : END IF
3028 5695767 : IF (jatom >= katom) THEN
3029 : i = 1
3030 : j = offset_bc
3031 17488404 : DO mc = 1, mc_max
3032 34264985 : DO mb = 1, mb_max
3033 16776581 : !$OMP ATOMIC
3034 : ks(j) = ks(j) - kbc(i)
3035 16776581 : i = i + 1
3036 28569218 : j = j + 1
3037 : END DO
3038 : END DO
3039 : ELSE
3040 : i = 1
3041 0 : DO mc = 1, mc_max
3042 0 : j = offset_bc + mc - 1
3043 0 : DO mb = 1, mb_max
3044 0 : !$OMP ATOMIC
3045 : ks(j) = ks(j) + kbc(i)
3046 0 : i = i + 1
3047 0 : j = j + mc_max
3048 : END DO
3049 : END DO
3050 : END IF
3051 5695767 : IF (iatom >= latom) THEN
3052 : i = 1
3053 : j = offset_ad
3054 12447084 : DO md = 1, md_max
3055 27986983 : DO ma = 1, ma_max
3056 15539899 : !$OMP ATOMIC
3057 : ks(j) = ks(j) - kad(i)
3058 15539899 : i = i + 1
3059 23722250 : j = j + 1
3060 : END DO
3061 : END DO
3062 : ELSE
3063 : i = 1
3064 3146468 : DO md = 1, md_max
3065 1715434 : j = offset_ad + md - 1
3066 6978906 : DO ma = 1, ma_max
3067 3832438 : !$OMP ATOMIC
3068 : ks(j) = ks(j) + kad(i)
3069 3832438 : i = i + 1
3070 5547872 : j = j + md_max
3071 : END DO
3072 : END DO
3073 : END IF
3074 : ! XXXXXXXXXXXXXXXXXXXXXXXX
3075 5695767 : IF (iatom >= katom) THEN
3076 : i = 1
3077 : j = offset_ac
3078 16983984 : DO mc = 1, mc_max
3079 39281281 : DO ma = 1, ma_max
3080 22297297 : !$OMP ATOMIC
3081 : ks(j) = ks(j) + kac(i)
3082 22297297 : i = i + 1
3083 33793288 : j = j + 1
3084 : END DO
3085 : END DO
3086 : ELSE
3087 : i = 1
3088 504420 : DO mc = 1, mc_max
3089 296646 : j = offset_ac + mc - 1
3090 1265292 : DO ma = 1, ma_max
3091 760872 : !$OMP ATOMIC
3092 : ks(j) = ks(j) - kac(i)
3093 760872 : i = i + 1
3094 1057518 : j = j + mc_max
3095 : END DO
3096 : END DO
3097 : END IF
3098 :
3099 5695767 : END SUBROUTINE update_fock_matrix_as
3100 :
3101 : ! **************************************************************************************************
3102 : !> \brief ...
3103 : !> \param i ...
3104 : !> \param j ...
3105 : !> \param k ...
3106 : !> \param l ...
3107 : !> \param set_offsets ...
3108 : !> \param atom_offsets ...
3109 : !> \param iset ...
3110 : !> \param jset ...
3111 : !> \param kset ...
3112 : !> \param lset ...
3113 : !> \param ma_max ...
3114 : !> \param mb_max ...
3115 : !> \param mc_max ...
3116 : !> \param md_max ...
3117 : !> \param prim ...
3118 : ! **************************************************************************************************
3119 0 : SUBROUTINE print_integrals(i, j, k, l, set_offsets, atom_offsets, iset, jset, kset, lset, ma_max, mb_max, mc_max, md_max, prim)
3120 : INTEGER :: i, j, k, l
3121 : INTEGER, DIMENSION(:, :, :, :), POINTER :: set_offsets
3122 : INTEGER, DIMENSION(:, :), POINTER :: atom_offsets
3123 : INTEGER :: iset, jset, kset, lset, ma_max, mb_max, &
3124 : mc_max, md_max
3125 : REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
3126 : INTENT(IN) :: prim
3127 :
3128 : INTEGER :: iint, ma, mb, mc, md
3129 :
3130 0 : iint = 0
3131 0 : DO md = 1, md_max
3132 0 : DO mc = 1, mc_max
3133 0 : DO mb = 1, mb_max
3134 0 : DO ma = 1, ma_max
3135 0 : iint = iint + 1
3136 0 : IF (ABS(prim(iint)) > 0.0000000000001) THEN
3137 0 : WRITE (99, *) atom_offsets(i, 1) + ma + set_offsets(iset, 1, i, 1) - 1, &
3138 0 : atom_offsets(j, 1) + ma + set_offsets(jset, 1, j, 1) - 1, &
3139 0 : atom_offsets(k, 1) + ma + set_offsets(kset, 1, k, 1) - 1, &
3140 0 : atom_offsets(l, 1) + ma + set_offsets(lset, 1, l, 1) - 1, &
3141 0 : prim(iint)
3142 : END IF
3143 : END DO
3144 : END DO
3145 : END DO
3146 : END DO
3147 :
3148 0 : END SUBROUTINE print_integrals
3149 :
3150 : #:include "hfx_get_pmax_val.fypp"
3151 :
3152 : END MODULE hfx_energy_potential
|
