(git:47b47fb)
Loading...
Searching...
No Matches
hfx_energy_potential.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
120CONTAINS
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 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 INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: tmp_task_list_cost
192 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, last_sgf_global, nimages, &
193 tmp_index, first_set_of_atom
194 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, lc_min, ld_max, &
195 ld_min, npgfa, npgfb, npgfc, npgfd, nsgfa, nsgfb, nsgfc, nsgfd, shm_block_offset
196 INTEGER, DIMENSION(:, :), POINTER :: first_sgfb, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, &
197 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 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
200 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 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 REAL(dp), ALLOCATABLE, DIMENSION(:) :: ee_buffer1, ee_buffer2, ee_primitives_tmp, ee_work, &
213 ee_work2, kac_buf, kad_buf, kbc_buf, kbd_buf, pac_buf, pad_buf, pbc_buf, pbd_buf, &
214 primitive_integrals
215 REAL(dp), DIMENSION(:), POINTER :: p_work
216 REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta, full_ks_alpha, &
217 full_ks_beta, max_contraction, ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, shm_pmax_atom, &
218 shm_pmax_block, sphi_b, zeta, zetb, zetc, zetd, &
219 full_ks_alpha_from_gpu, full_ks_beta_from_gpu
220 REAL(dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, &
221 sphi_c_ext_set, sphi_d_ext_set
222 REAL(dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, &
223 sphi_d_ext
224 REAL(dp), DIMENSION(:, :, :), POINTER :: gcc
225 REAL(kind=dp) :: coeffs_kind_max0
226 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 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 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
234 TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches, integral_caches_disk
235 TYPE(hfx_cache_type), POINTER :: maxval_cache, maxval_cache_disk
236 TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers, &
237 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 TYPE(hfx_p_kind), DIMENSION(:), POINTER :: shm_initial_p
244 TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl
245 TYPE(hfx_pgf_product_list), ALLOCATABLE, &
246 DIMENSION(:) :: pgf_product_list
247 TYPE(hfx_potential_type) :: potential_parameter
248 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
249 POINTER :: screen_coeffs_kind, tmp_r_1, tmp_r_2, &
250 tmp_screen_pgf1, tmp_screen_pgf2
251 TYPE(hfx_screen_coeff_type), &
252 DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
253 TYPE(hfx_screen_coeff_type), &
254 DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf
255 TYPE(hfx_screening_type) :: screening_parameter
256 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 DIMENSION(:) :: set_list_ij, set_list_kl
262 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 NULLIFY (dft_control, matrix_ks_aux_fit_hfx)
273
274 CALL timeset(routinen, handle)
275
276 CALL cite_reference(guidon2008)
277 CALL cite_reference(guidon2009)
278
279 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 my_geo_change = geometry_did_change
283 IF (ispin == 2) my_geo_change = .false.
284
285 logger => cp_get_default_logger()
286
287 is_anti_symmetric = dbcsr_get_matrix_type(rho_ao(1, 1)%matrix) == dbcsr_type_antisymmetric
288
289 iw = -1
290
291 IF (my_geo_change) THEN
292 CALL m_memory(memsize_before)
293 CALL para_env%max(memsize_before)
294 iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
295 extension=".scfLog")
296 IF (iw > 0) THEN
297 WRITE (unit=iw, fmt="(/,(T3,A,T60,I21))") &
298 "HFX_MEM_INFO| Est. max. program size before HFX [MiB]:", memsize_before/(1024*1024)
299 CALL m_flush(iw)
300 END IF
301 CALL cp_print_key_finished_output(iw, logger, hfx_section, &
302 "HF_INFO")
303 END IF
304
305 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell)
306
307 NULLIFY (cell_to_index)
308 CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
309 IF (do_kpoints) THEN
310 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
311 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
312 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 nkind = SIZE(atomic_kind_set, 1)
317 l_max = 0
318 DO ikind = 1, nkind
319 l_max = max(l_max, maxval(x_data(1, 1)%basis_parameter(ikind)%lmax))
320 END DO
321 l_max = 4*l_max
322 CALL init_md_ftable(l_max)
323
324 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 IF (l_max > init_t_c_g0_lmax) THEN
327 IF (para_env%is_source()) THEN
328 CALL open_file(unit_number=unit_id, file_name=x_data(1, 1)%potential_parameter%filename)
329 END IF
330 CALL init(l_max, unit_id, para_env%mepos, para_env)
331 IF (para_env%is_source()) THEN
332 CALL close_file(unit_id)
333 END IF
334 init_t_c_g0_lmax = l_max
335 END IF
336 END IF
337
338 n_threads = 1
339!$ n_threads = omp_get_max_threads()
340
341 ! This initialization is needed to prevent a segmentation fault. The correct assigment is done below
342 my_nspins = 0
343 IF (PRESENT(nspins)) my_nspins = nspins
344
345 shm_neris_total = 0
346 shm_nprim_ints = 0
347 shm_neris_onthefly = 0
348 shm_storage_counter_integrals = 0
349 shm_stor_count_int_disk = 0
350 shm_neris_incore = 0
351 shm_neris_disk = 0
352 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!$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 CALL timestop(handle)
2336 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 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 zeta, zetb, zetc, zetd, &
2426 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 nsgfl_a, nsgfl_b, nsgfl_c, &
2434 nsgfl_d, &
2435 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 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 neighbor_cells, nimages, do_periodic)
2486
2487 cart_estimate = 0.0_dp
2488 neris_tmp = 0
2489 primitive_integrals = 0.0_dp
2490 max_l = la_max + lb_max + lc_max + ld_max
2491
2492 DO list_ij = 1, nelements_ij
2493 zetainv = pgf_list_ij(list_ij)%ZetaInv
2494 ipgf = pgf_list_ij(list_ij)%ipgf
2495 jpgf = pgf_list_ij(list_ij)%jpgf
2496
2497 DO list_kl = 1, nelements_kl
2498 etainv = pgf_list_kl(list_kl)%ZetaInv
2499 kpgf = pgf_list_kl(list_kl)%ipgf
2500 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 potential_parameter, max_l, do_periodic)
2505
2506 s_offset_a = 0
2507 DO la = la_min, la_max
2508 s_offset_b = 0
2509 ncoa = nco(la)
2510 nsgfla = nsgfl_a(la)
2511 nsoa = nso(la)
2512
2513 DO lb = lb_min, lb_max
2514 s_offset_c = 0
2515 ncob = nco(lb)
2516 nsgflb = nsgfl_b(lb)
2517 nsob = nso(lb)
2518
2519 DO lc = lc_min, lc_max
2520 s_offset_d = 0
2521 ncoc = nco(lc)
2522 nsgflc = nsgfl_c(lc)
2523 nsoc = nso(lc)
2524
2525 DO ld = ld_min, ld_max
2526 ncod = nco(ld)
2527 nsgfld = nsgfl_d(ld)
2528 nsod = nso(ld)
2529
2530 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 p_work)
2546 cart_estimate = max(tmp_max, cart_estimate)
2547 s_offset_d = s_offset_d + nsod*nsgfld
2548 END DO !ld
2549 s_offset_c = s_offset_c + nsoc*nsgflc
2550 END DO !lc
2551 s_offset_b = s_offset_b + nsob*nsgflb
2552 END DO !lb
2553 s_offset_a = s_offset_a + nsoa*nsgfla
2554 END DO !la
2555 END DO
2556 END DO
2557
2558 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 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 min_ij = min(i, j)
2581 get_1d_idx = min_ij*n + max(i, j) - (min_ij - 1)*min_ij/2 - n
2582
2583 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 SUBROUTINE update_fock_matrix(ma_max, mb_max, mc_max, md_max, &
2628 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 IF (jatom >= latom) THEN
2654 i = 1
2655 offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
2656 j = offset_bd
2657 DO md = 1, md_max
2658 DO mb = 1, mb_max
2659 pbd(i) = density(j)
2660 i = i + 1
2661 j = j + 1
2662 END DO
2663 END DO
2664 ELSE
2665 i = 1
2666 offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
2667 DO md = 1, md_max
2668 j = offset_bd + md - 1
2669 DO mb = 1, mb_max
2670 pbd(i) = density(j)
2671 i = i + 1
2672 j = j + md_max
2673 END DO
2674 END DO
2675 END IF
2676 IF (jatom >= katom) THEN
2677 i = 1
2678 offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
2679 j = offset_bc
2680 DO mc = 1, mc_max
2681 DO mb = 1, mb_max
2682 pbc(i) = density(j)
2683 i = i + 1
2684 j = j + 1
2685 END DO
2686 END DO
2687 ELSE
2688 i = 1
2689 offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
2690 DO mc = 1, mc_max
2691 j = offset_bc + mc - 1
2692 DO mb = 1, mb_max
2693 pbc(i) = density(j)
2694 i = i + 1
2695 j = j + mc_max
2696 END DO
2697 END DO
2698 END IF
2699 IF (iatom >= latom) THEN
2700 i = 1
2701 offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
2702 j = offset_ad
2703 DO md = 1, md_max
2704 DO ma = 1, ma_max
2705 pad(i) = density(j)
2706 i = i + 1
2707 j = j + 1
2708 END DO
2709 END DO
2710 ELSE
2711 i = 1
2712 offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
2713 DO md = 1, md_max
2714 j = offset_ad + md - 1
2715 DO ma = 1, ma_max
2716 pad(i) = density(j)
2717 i = i + 1
2718 j = j + md_max
2719 END DO
2720 END DO
2721 END IF
2722 IF (iatom >= katom) THEN
2723 i = 1
2724 offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
2725 j = offset_ac
2726 DO mc = 1, mc_max
2727 DO ma = 1, ma_max
2728 pac(i) = density(j)
2729 i = i + 1
2730 j = j + 1
2731 END DO
2732 END DO
2733 ELSE
2734 i = 1
2735 offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
2736 DO mc = 1, mc_max
2737 j = offset_ac + mc - 1
2738 DO ma = 1, ma_max
2739 pac(i) = density(j)
2740 i = i + 1
2741 j = j + mc_max
2742 END DO
2743 END DO
2744 END IF
2745
2746 CALL contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, fac*symm_fac)
2747 IF (jatom >= latom) THEN
2748 i = 1
2749 j = offset_bd
2750 DO md = 1, md_max
2751 DO mb = 1, mb_max
2752 ki = kbd(i)
2753!$OMP ATOMIC
2754 ks(j) = ks(j) + ki
2755 i = i + 1
2756 j = j + 1
2757 END DO
2758 END DO
2759 ELSE
2760 i = 1
2761 DO md = 1, md_max
2762 j = offset_bd + md - 1
2763 DO mb = 1, mb_max
2764 ki = kbd(i)
2765!$OMP ATOMIC
2766 ks(j) = ks(j) + ki
2767 i = i + 1
2768 j = j + md_max
2769 END DO
2770 END DO
2771 END IF
2772 IF (jatom >= katom) THEN
2773 i = 1
2774 j = offset_bc
2775 DO mc = 1, mc_max
2776 DO mb = 1, mb_max
2777 ki = kbc(i)
2778!$OMP ATOMIC
2779 ks(j) = ks(j) + ki
2780 i = i + 1
2781 j = j + 1
2782 END DO
2783 END DO
2784 ELSE
2785 i = 1
2786 DO mc = 1, mc_max
2787 j = offset_bc + mc - 1
2788 DO mb = 1, mb_max
2789 ki = kbc(i)
2790!$OMP ATOMIC
2791 ks(j) = ks(j) + ki
2792 i = i + 1
2793 j = j + mc_max
2794 END DO
2795 END DO
2796 END IF
2797 IF (iatom >= latom) THEN
2798 i = 1
2799 j = offset_ad
2800 DO md = 1, md_max
2801 DO ma = 1, ma_max
2802 ki = kad(i)
2803!$OMP ATOMIC
2804 ks(j) = ks(j) + ki
2805 i = i + 1
2806 j = j + 1
2807 END DO
2808 END DO
2809 ELSE
2810 i = 1
2811 DO md = 1, md_max
2812 j = offset_ad + md - 1
2813 DO ma = 1, ma_max
2814 ki = kad(i)
2815!$OMP ATOMIC
2816 ks(j) = ks(j) + ki
2817 i = i + 1
2818 j = j + md_max
2819 END DO
2820 END DO
2821 END IF
2822 IF (iatom >= katom) THEN
2823 i = 1
2824 j = offset_ac
2825 DO mc = 1, mc_max
2826 DO ma = 1, ma_max
2827 ki = kac(i)
2828!$OMP ATOMIC
2829 ks(j) = ks(j) + ki
2830 i = i + 1
2831 j = j + 1
2832 END DO
2833 END DO
2834 ELSE
2835 i = 1
2836 DO mc = 1, mc_max
2837 j = offset_ac + mc - 1
2838 DO ma = 1, ma_max
2839 ki = kac(i)
2840!$OMP ATOMIC
2841 ks(j) = ks(j) + ki
2842 i = i + 1
2843 j = j + mc_max
2844 END DO
2845 END DO
2846 END IF
2847 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 SUBROUTINE update_fock_matrix_as(ma_max, mb_max, mc_max, md_max, &
2886 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 IF (jatom >= latom) THEN
2911 i = 1
2912 offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
2913 j = offset_bd
2914 DO md = 1, md_max
2915 DO mb = 1, mb_max
2916 pbd(i) = +density(j)
2917 i = i + 1
2918 j = j + 1
2919 END DO
2920 END DO
2921 ELSE
2922 i = 1
2923 offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
2924 DO md = 1, md_max
2925 j = offset_bd + md - 1
2926 DO mb = 1, mb_max
2927 pbd(i) = -density(j)
2928 i = i + 1
2929 j = j + md_max
2930 END DO
2931 END DO
2932 END IF
2933 IF (jatom >= katom) THEN
2934 i = 1
2935 offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
2936 j = offset_bc
2937 DO mc = 1, mc_max
2938 DO mb = 1, mb_max
2939 pbc(i) = -density(j)
2940 i = i + 1
2941 j = j + 1
2942 END DO
2943 END DO
2944 ELSE
2945 i = 1
2946 offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
2947 DO mc = 1, mc_max
2948 j = offset_bc + mc - 1
2949 DO mb = 1, mb_max
2950 pbc(i) = density(j)
2951 i = i + 1
2952 j = j + mc_max
2953 END DO
2954 END DO
2955 END IF
2956 IF (iatom >= latom) THEN
2957 i = 1
2958 offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
2959 j = offset_ad
2960 DO md = 1, md_max
2961 DO ma = 1, ma_max
2962 pad(i) = -density(j)
2963 i = i + 1
2964 j = j + 1
2965 END DO
2966 END DO
2967 ELSE
2968 i = 1
2969 offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
2970 DO md = 1, md_max
2971 j = offset_ad + md - 1
2972 DO ma = 1, ma_max
2973 pad(i) = density(j)
2974 i = i + 1
2975 j = j + md_max
2976 END DO
2977 END DO
2978 END IF
2979 IF (iatom >= katom) THEN
2980 i = 1
2981 offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
2982 j = offset_ac
2983 DO mc = 1, mc_max
2984 DO ma = 1, ma_max
2985 pac(i) = +density(j)
2986 i = i + 1
2987 j = j + 1
2988 END DO
2989 END DO
2990 ELSE
2991 i = 1
2992 offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
2993 DO mc = 1, mc_max
2994 j = offset_ac + mc - 1
2995 DO ma = 1, ma_max
2996 pac(i) = -density(j)
2997 i = i + 1
2998 j = j + mc_max
2999 END DO
3000 END DO
3001 END IF
3002
3003 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 IF (jatom >= latom) THEN
3006 i = 1
3007 j = offset_bd
3008 DO md = 1, md_max
3009 DO mb = 1, mb_max
3010!$OMP ATOMIC
3011 ks(j) = ks(j) + kbd(i)
3012 i = i + 1
3013 j = j + 1
3014 END DO
3015 END DO
3016 ELSE
3017 i = 1
3018 DO md = 1, md_max
3019 j = offset_bd + md - 1
3020 DO mb = 1, mb_max
3021!$OMP ATOMIC
3022 ks(j) = ks(j) - kbd(i)
3023 i = i + 1
3024 j = j + md_max
3025 END DO
3026 END DO
3027 END IF
3028 IF (jatom >= katom) THEN
3029 i = 1
3030 j = offset_bc
3031 DO mc = 1, mc_max
3032 DO mb = 1, mb_max
3033!$OMP ATOMIC
3034 ks(j) = ks(j) - kbc(i)
3035 i = i + 1
3036 j = j + 1
3037 END DO
3038 END DO
3039 ELSE
3040 i = 1
3041 DO mc = 1, mc_max
3042 j = offset_bc + mc - 1
3043 DO mb = 1, mb_max
3044!$OMP ATOMIC
3045 ks(j) = ks(j) + kbc(i)
3046 i = i + 1
3047 j = j + mc_max
3048 END DO
3049 END DO
3050 END IF
3051 IF (iatom >= latom) THEN
3052 i = 1
3053 j = offset_ad
3054 DO md = 1, md_max
3055 DO ma = 1, ma_max
3056!$OMP ATOMIC
3057 ks(j) = ks(j) - kad(i)
3058 i = i + 1
3059 j = j + 1
3060 END DO
3061 END DO
3062 ELSE
3063 i = 1
3064 DO md = 1, md_max
3065 j = offset_ad + md - 1
3066 DO ma = 1, ma_max
3067!$OMP ATOMIC
3068 ks(j) = ks(j) + kad(i)
3069 i = i + 1
3070 j = j + md_max
3071 END DO
3072 END DO
3073 END IF
3074! XXXXXXXXXXXXXXXXXXXXXXXX
3075 IF (iatom >= katom) THEN
3076 i = 1
3077 j = offset_ac
3078 DO mc = 1, mc_max
3079 DO ma = 1, ma_max
3080!$OMP ATOMIC
3081 ks(j) = ks(j) + kac(i)
3082 i = i + 1
3083 j = j + 1
3084 END DO
3085 END DO
3086 ELSE
3087 i = 1
3088 DO mc = 1, mc_max
3089 j = offset_ac + mc - 1
3090 DO ma = 1, ma_max
3091!$OMP ATOMIC
3092 ks(j) = ks(j) - kac(i)
3093 i = i + 1
3094 j = j + mc_max
3095 END DO
3096 END DO
3097 END IF
3098
3099 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!> \param output_unit ...
3119! **************************************************************************************************
3120 SUBROUTINE print_integrals(i, j, k, l, set_offsets, atom_offsets, iset, jset, kset, lset, &
3121 ma_max, mb_max, mc_max, md_max, prim, output_unit)
3122 INTEGER :: i, j, k, l
3123 INTEGER, DIMENSION(:, :, :, :), POINTER :: set_offsets
3124 INTEGER, DIMENSION(:, :), POINTER :: atom_offsets
3125 INTEGER :: iset, jset, kset, lset, ma_max, mb_max, &
3126 mc_max, md_max
3127 REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
3128 INTENT(IN) :: prim
3129 INTEGER, OPTIONAL :: output_unit
3130
3131 REAL(kind=dp), PARAMETER :: eps_int = 1.0e-13_dp
3132
3133 INTEGER :: iint, ma, mb, mc, md, my_output_unit
3134
3135 my_output_unit = 99
3136 IF (PRESENT(output_unit)) my_output_unit = output_unit
3137
3138 iint = 0
3139 DO md = 1, md_max
3140 DO mc = 1, mc_max
3141 DO mb = 1, mb_max
3142 DO ma = 1, ma_max
3143 iint = iint + 1
3144 IF (abs(prim(iint)) > eps_int) THEN
3145 IF (my_output_unit > 0) THEN
3146 WRITE (my_output_unit, *) atom_offsets(i, 1) + ma + set_offsets(iset, 1, i, 1) - 1, &
3147 atom_offsets(j, 1) + ma + set_offsets(jset, 1, j, 1) - 1, &
3148 atom_offsets(k, 1) + ma + set_offsets(kset, 1, k, 1) - 1, &
3149 atom_offsets(l, 1) + ma + set_offsets(lset, 1, l, 1) - 1, &
3150 prim(iint)
3151 END IF
3152 END IF
3153 END DO
3154 END DO
3155 END DO
3156 END DO
3157
3158 END SUBROUTINE print_integrals
3159
3160
3161! **************************************************************************************************
3162!> \brief This routine calculates the maximum density matrix element, when
3163!> screening on an initial density matrix is applied. Due to symmetry of
3164!> the ERI's, there are always 4 matrix elements to be considered.
3165!> CASE 0-15 belong to an energy calculation (linear screening)
3166!> CASE 16-31 belong to a force calculation (square screening)
3167!> \param ptr_p_1 Pointers to atomic density matrices
3168!> \param ptr_p_2 Pointers to atomic density matrices
3169!> \param ptr_p_3 Pointers to atomic density matrices
3170!> \param ptr_p_4 Pointers to atomic density matrices
3171!> \param iset Current set
3172!> \param jset Current set
3173!> \param kset Current set
3174!> \param lset Current set
3175!> \param pmax_val value to be calculated
3176!> \param swap_id Defines how the matrices are accessed
3177!> \par History
3178!> 06.2009 created [Manuel Guidon]
3179!> \author Manuel Guidon
3180! **************************************************************************************************
3181PURE SUBROUTINE get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, iset, jset, kset, lset, pmax_val, swap_id)
3182
3183 REAL(dp), DIMENSION(:, :), POINTER :: ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4
3184 INTEGER, INTENT(IN) :: iset, jset, kset, lset
3185
3186 REAL(dp), INTENT(OUT) :: pmax_val
3187 INTEGER, INTENT(IN) :: swap_id
3188
3189 REAL(dp) :: pmax_1, pmax_2, pmax_3, pmax_4
3190
3191 SELECT CASE (swap_id)
3192 CASE (0)
3193 pmax_1 = ptr_p_1(kset, iset)
3194 pmax_2 = ptr_p_2(lset, jset)
3195 pmax_3 = ptr_p_3(lset, iset)
3196 pmax_4 = ptr_p_4(kset, jset)
3197 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3198 CASE (1)
3199 pmax_1 = ptr_p_1(iset, kset)
3200 pmax_2 = ptr_p_2(lset, jset)
3201 pmax_3 = ptr_p_3(lset, iset)
3202 pmax_4 = ptr_p_4(kset, jset)
3203 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3204 CASE (2)
3205 pmax_1 = ptr_p_1(kset, iset)
3206 pmax_2 = ptr_p_2(jset, lset)
3207 pmax_3 = ptr_p_3(lset, iset)
3208 pmax_4 = ptr_p_4(kset, jset)
3209 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3210 CASE (3)
3211 pmax_1 = ptr_p_1(iset, kset)
3212 pmax_2 = ptr_p_2(jset, lset)
3213 pmax_3 = ptr_p_3(lset, iset)
3214 pmax_4 = ptr_p_4(kset, jset)
3215 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3216 CASE (4)
3217 pmax_1 = ptr_p_1(kset, iset)
3218 pmax_2 = ptr_p_2(lset, jset)
3219 pmax_3 = ptr_p_3(iset, lset)
3220 pmax_4 = ptr_p_4(kset, jset)
3221 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3222 CASE (5)
3223 pmax_1 = ptr_p_1(iset, kset)
3224 pmax_2 = ptr_p_2(lset, jset)
3225 pmax_3 = ptr_p_3(iset, lset)
3226 pmax_4 = ptr_p_4(kset, jset)
3227 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3228 CASE (6)
3229 pmax_1 = ptr_p_1(kset, iset)
3230 pmax_2 = ptr_p_2(jset, lset)
3231 pmax_3 = ptr_p_3(iset, lset)
3232 pmax_4 = ptr_p_4(kset, jset)
3233 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3234 CASE (7)
3235 pmax_1 = ptr_p_1(iset, kset)
3236 pmax_2 = ptr_p_2(jset, lset)
3237 pmax_3 = ptr_p_3(iset, lset)
3238 pmax_4 = ptr_p_4(kset, jset)
3239 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3240 CASE (8)
3241 pmax_1 = ptr_p_1(kset, iset)
3242 pmax_2 = ptr_p_2(lset, jset)
3243 pmax_3 = ptr_p_3(lset, iset)
3244 pmax_4 = ptr_p_4(jset, kset)
3245 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3246 CASE (9)
3247 pmax_1 = ptr_p_1(iset, kset)
3248 pmax_2 = ptr_p_2(lset, jset)
3249 pmax_3 = ptr_p_3(lset, iset)
3250 pmax_4 = ptr_p_4(jset, kset)
3251 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3252 CASE (10)
3253 pmax_1 = ptr_p_1(kset, iset)
3254 pmax_2 = ptr_p_2(jset, lset)
3255 pmax_3 = ptr_p_3(lset, iset)
3256 pmax_4 = ptr_p_4(jset, kset)
3257 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3258 CASE (11)
3259 pmax_1 = ptr_p_1(iset, kset)
3260 pmax_2 = ptr_p_2(jset, lset)
3261 pmax_3 = ptr_p_3(lset, iset)
3262 pmax_4 = ptr_p_4(jset, kset)
3263 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3264 CASE (12)
3265 pmax_1 = ptr_p_1(kset, iset)
3266 pmax_2 = ptr_p_2(lset, jset)
3267 pmax_3 = ptr_p_3(iset, lset)
3268 pmax_4 = ptr_p_4(jset, kset)
3269 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3270 CASE (13)
3271 pmax_1 = ptr_p_1(iset, kset)
3272 pmax_2 = ptr_p_2(lset, jset)
3273 pmax_3 = ptr_p_3(iset, lset)
3274 pmax_4 = ptr_p_4(jset, kset)
3275 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3276 CASE (14)
3277 pmax_1 = ptr_p_1(kset, iset)
3278 pmax_2 = ptr_p_2(jset, lset)
3279 pmax_3 = ptr_p_3(iset, lset)
3280 pmax_4 = ptr_p_4(jset, kset)
3281 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3282 CASE (15)
3283 pmax_1 = ptr_p_1(iset, kset)
3284 pmax_2 = ptr_p_2(jset, lset)
3285 pmax_3 = ptr_p_3(iset, lset)
3286 pmax_4 = ptr_p_4(jset, kset)
3287 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3288 CASE (16)
3289 pmax_1 = ptr_p_1(kset, iset)
3290 pmax_2 = ptr_p_2(lset, jset)
3291 pmax_3 = ptr_p_3(lset, iset)
3292 pmax_4 = ptr_p_4(kset, jset)
3293 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3294 CASE (17)
3295 pmax_1 = ptr_p_1(iset, kset)
3296 pmax_2 = ptr_p_2(lset, jset)
3297 pmax_3 = ptr_p_3(lset, iset)
3298 pmax_4 = ptr_p_4(kset, jset)
3299 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3300 CASE (18)
3301 pmax_1 = ptr_p_1(kset, iset)
3302 pmax_2 = ptr_p_2(jset, lset)
3303 pmax_3 = ptr_p_3(lset, iset)
3304 pmax_4 = ptr_p_4(kset, jset)
3305 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3306 CASE (19)
3307 pmax_1 = ptr_p_1(iset, kset)
3308 pmax_2 = ptr_p_2(jset, lset)
3309 pmax_3 = ptr_p_3(lset, iset)
3310 pmax_4 = ptr_p_4(kset, jset)
3311 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3312 CASE (20)
3313 pmax_1 = ptr_p_1(kset, iset)
3314 pmax_2 = ptr_p_2(lset, jset)
3315 pmax_3 = ptr_p_3(iset, lset)
3316 pmax_4 = ptr_p_4(kset, jset)
3317 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3318 CASE (21)
3319 pmax_1 = ptr_p_1(iset, kset)
3320 pmax_2 = ptr_p_2(lset, jset)
3321 pmax_3 = ptr_p_3(iset, lset)
3322 pmax_4 = ptr_p_4(kset, jset)
3323 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3324 CASE (22)
3325 pmax_1 = ptr_p_1(kset, iset)
3326 pmax_2 = ptr_p_2(jset, lset)
3327 pmax_3 = ptr_p_3(iset, lset)
3328 pmax_4 = ptr_p_4(kset, jset)
3329 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3330 CASE (23)
3331 pmax_1 = ptr_p_1(iset, kset)
3332 pmax_2 = ptr_p_2(jset, lset)
3333 pmax_3 = ptr_p_3(iset, lset)
3334 pmax_4 = ptr_p_4(kset, jset)
3335 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3336 CASE (24)
3337 pmax_1 = ptr_p_1(kset, iset)
3338 pmax_2 = ptr_p_2(lset, jset)
3339 pmax_3 = ptr_p_3(lset, iset)
3340 pmax_4 = ptr_p_4(jset, kset)
3341 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3342 CASE (25)
3343 pmax_1 = ptr_p_1(iset, kset)
3344 pmax_2 = ptr_p_2(lset, jset)
3345 pmax_3 = ptr_p_3(lset, iset)
3346 pmax_4 = ptr_p_4(jset, kset)
3347 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3348 CASE (26)
3349 pmax_1 = ptr_p_1(kset, iset)
3350 pmax_2 = ptr_p_2(jset, lset)
3351 pmax_3 = ptr_p_3(lset, iset)
3352 pmax_4 = ptr_p_4(jset, kset)
3353 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3354 CASE (27)
3355 pmax_1 = ptr_p_1(iset, kset)
3356 pmax_2 = ptr_p_2(jset, lset)
3357 pmax_3 = ptr_p_3(lset, iset)
3358 pmax_4 = ptr_p_4(jset, kset)
3359 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3360 CASE (28)
3361 pmax_1 = ptr_p_1(kset, iset)
3362 pmax_2 = ptr_p_2(lset, jset)
3363 pmax_3 = ptr_p_3(iset, lset)
3364 pmax_4 = ptr_p_4(jset, kset)
3365 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3366 CASE (29)
3367 pmax_1 = ptr_p_1(iset, kset)
3368 pmax_2 = ptr_p_2(lset, jset)
3369 pmax_3 = ptr_p_3(iset, lset)
3370 pmax_4 = ptr_p_4(jset, kset)
3371 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3372 CASE (30)
3373 pmax_1 = ptr_p_1(kset, iset)
3374 pmax_2 = ptr_p_2(jset, lset)
3375 pmax_3 = ptr_p_3(iset, lset)
3376 pmax_4 = ptr_p_4(jset, kset)
3377 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3378 CASE (31)
3379 pmax_1 = ptr_p_1(iset, kset)
3380 pmax_2 = ptr_p_2(jset, lset)
3381 pmax_3 = ptr_p_3(iset, lset)
3382 pmax_4 = ptr_p_4(jset, kset)
3383 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3384 END SELECT
3385
3386END SUBROUTINE get_pmax_val
3387
3388
3389END MODULE hfx_energy_potential
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cell_types::pbc
Definition cell_types.F:105
libgint_wrapper::libgint_get_fock_matrix
Definition libgint_wrapper.F:61
libgint_wrapper::libgint_set_density
Definition libgint_wrapper.F:56
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:599
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::guidon2008
integer, save, public guidon2008
Definition bibliography.F:36
bibliography::guidon2009
integer, save, public guidon2009
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp2k_info
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:22
cp2k_info::cp2k_flags
character(len=10 *default_string_length) function, public cp2k_flags()
list all compile time options that influence the capabilities of cp2k. All new flags should be added ...
Definition cp2k_info.F:89
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_dot_threadsafe
subroutine, public dbcsr_dot_threadsafe(matrix_a, matrix_b, result)
...
Definition cp_dbcsr_api.F:1325
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:196
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:63
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
gamma::init_md_ftable
subroutine, public init_md_ftable(nmax)
Initialize a table of F_n(t) values in the range 0 <= t <= 12 with a stepsize of 0....
Definition gamma.F:540
hfx_communication
Routines for data exchange between MPI processes.
Definition hfx_communication.F:14
hfx_communication::distribute_ks_matrix
subroutine, public distribute_ks_matrix(para_env, full_ks, ks_matrix, number_of_p_entries, block_offset, kind_of, basis_parameter, off_diag_fac, diag_fac)
Distributes the local full Kohn-Sham matrix to all CPUS
Definition hfx_communication.F:203
hfx_communication::get_atomic_block_maps
subroutine, public get_atomic_block_maps(matrix, basis_parameter, kind_of, is_assoc_atomic_block, number_of_p_entries, para_env, atomic_block_offset, set_offset, block_offset, map_atoms_to_cpus, nkind)
create a several maps array that reflects the ks matrix sparsity
Definition hfx_communication.F:427
hfx_communication::get_full_density
subroutine, public get_full_density(para_env, full_density, rho, number_of_p_entries, block_offset, kind_of, basis_parameter, get_max_vals_spin, rho_beta, antisymmetric)
Collects full density matrix from all CPUs
Definition hfx_communication.F:76
hfx_compression_methods
routines and types for Hartree-Fock-Exchange
Definition hfx_compression_methods.F:14
hfx_compression_methods::hfx_add_single_cache_element
subroutine, public hfx_add_single_cache_element(value, nbits, cache, container, memory_usage, use_disk_storage, max_val_memory)
This routine adds an int_8 value to a cache. If the cache is full a compression routine is invoked an...
Definition hfx_compression_methods.F:72
hfx_compression_methods::hfx_get_mult_cache_elements
subroutine, public hfx_get_mult_cache_elements(values, nints, nbits, cache, container, eps_schwarz, pmax_entry, memory_usage, use_disk_storage)
This routine returns a bunch real values from a cache. If the cache is empty a decompression routine ...
Definition hfx_compression_methods.F:425
hfx_compression_methods::hfx_flush_last_cache
subroutine, public hfx_flush_last_cache(nbits, cache, container, memory_usage, use_disk_storage)
This routine compresses the last probably not yet compressed cache into a container
Definition hfx_compression_methods.F:309
hfx_compression_methods::hfx_get_single_cache_element
subroutine, public hfx_get_single_cache_element(value, nbits, cache, container, memory_usage, use_disk_storage)
This routine returns an int_8 value from a cache. If the cache is empty a decompression routine is in...
Definition hfx_compression_methods.F:168
hfx_compression_methods::hfx_decompress_first_cache
subroutine, public hfx_decompress_first_cache(nbits, cache, container, memory_usage, use_disk_storage)
This routine decompresses the first bunch of data in a container and copies them into a cache
Definition hfx_compression_methods.F:286
hfx_compression_methods::hfx_add_mult_cache_elements
subroutine, public hfx_add_mult_cache_elements(values, nints, nbits, cache, container, eps_schwarz, pmax_entry, memory_usage, use_disk_storage)
This routine adds an a few real values to a cache. If the cache is full a compression routine is invo...
Definition hfx_compression_methods.F:345
hfx_compression_methods::hfx_reset_cache_and_container
subroutine, public hfx_reset_cache_and_container(cache, container, memory_usage, do_disk_storage)
This routine resets the containers list pointer to the first element and moves the element counters o...
Definition hfx_compression_methods.F:255
hfx_contract_block
routines to contract density matrix blocks with the for center integrals to yield the Kohn-Sham matri...
Definition hfx_contract_block.F:15
hfx_contract_block::contract_block
subroutine, public contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, scale)
...
Definition hfx_contract_block.F:42
hfx_energy_potential
Routines to calculate HFX energy and potential.
Definition hfx_energy_potential.F:15
hfx_energy_potential::coulomb4
subroutine, public coulomb4(lib, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, zeta, zetb, zetc, zetd, primitive_integrals, potential_parameter, neighbor_cells, screen1, screen2, eps_schwarz, max_contraction_val, cart_estimate, cell, neris_tmp, log10_pmax, log10_eps_schwarz, r1_pgf, r2_pgf, pgf1, pgf2, pgf_list_ij, pgf_list_kl, pgf_product_list, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, nimages, do_periodic, p_work)
calculates two-electron integrals of a quartet/shell using the library lib_int in the periodic case
Definition hfx_energy_potential.F:2438
hfx_energy_potential::integrate_four_center
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin, nspins)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Definition hfx_energy_potential.F:150
hfx_libint_interface
Interface to the Libint-Library.
Definition hfx_libint_interface.F:15
hfx_libint_interface::evaluate_eri
subroutine, public evaluate_eri(libint, nproducts, pgf_product_list, n_a, n_b, n_c, n_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd, primitives, max_contraction, tmp_max, eps_schwarz, neris, zetainv, etainv, s_offset_a, s_offset_b, s_offset_c, s_offset_d, nl_a, nl_b, nl_c, nl_d, nsoa, nsob, nsoc, nsod, sphi_a, sphi_b, sphi_c, sphi_d, work, work2, buffer1, buffer2, primitives_tmp, p_work)
Evaluate electron repulsion integrals for a primitive quartet.
Definition hfx_libint_interface.F:1464
hfx_load_balance_methods
Routines for optimizing load balance between processes in HFX calculations.
Definition hfx_load_balance_methods.F:14
hfx_load_balance_methods::collect_load_balance_info
subroutine, public collect_load_balance_info(para_env, x_data, iw, n_threads, i_thread, eval_type)
...
Definition hfx_load_balance_methods.F:2331
hfx_load_balance_methods::hfx_load_balance
subroutine, public hfx_load_balance(x_data, eps_schwarz, particle_set, max_set, para_env, coeffs_set, coeffs_kind, is_assoc_atomic_block_global, do_periodic, load_balance_parameter, kind_of, basis_parameter, pmax_set, pmax_atom, i_thread, n_threads, cell, do_p_screening, map_atom_to_kind_atom, nkind, eval_type, pmax_block, use_virial)
Distributes the computation of eri's to all available processes.
Definition hfx_load_balance_methods.F:130
hfx_load_balance_methods::hfx_update_load_balance
subroutine, public hfx_update_load_balance(x_data, para_env, load_balance_parameter, i_thread, n_threads, eval_type)
Cheap way of redistributing the eri's.
Definition hfx_load_balance_methods.F:1483
hfx_pair_list_methods
Routines for optimizing load balance between processes in HFX calculations.
Definition hfx_pair_list_methods.F:15
hfx_pair_list_methods::build_pgf_product_list
subroutine, public build_pgf_product_list(list1, list2, product_list, nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, potential_parameter, m_max, do_periodic)
...
Definition hfx_pair_list_methods.F:196
hfx_pair_list_methods::build_atomic_pair_list
subroutine, public build_atomic_pair_list(natom, atomic_pair_list, kind_of, basis_parameter, particle_set, do_periodic, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, blocks)
...
Definition hfx_pair_list_methods.F:686
hfx_pair_list_methods::build_pair_list
subroutine, public build_pair_list(natom, list, set_list, i_start, i_end, j_start, j_end, kind_of, basis_parameter, particle_set, do_periodic, coeffs_set, coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, pmax_blocks, atomic_pair_list)
...
Definition hfx_pair_list_methods.F:591
hfx_pair_list_methods::build_pair_list_pgf
subroutine, public build_pair_list_pgf(npgfa, npgfb, list, zeta, zetb, screen1, screen2, pgf, r_pgf, log10_pmax, log10_eps_schwarz, ra, rb, nelements, neighbor_cells, nimages, do_periodic)
...
Definition hfx_pair_list_methods.F:472
hfx_pair_list_methods::pgf_product_list_size
integer, save, public pgf_product_list_size
Definition hfx_pair_list_methods.F:48
hfx_screening_methods
Several screening methods used in HFX calcualtions.
Definition hfx_screening_methods.F:14
hfx_screening_methods::calc_screening_functions
subroutine, public calc_screening_functions(qs_env, basis_parameter, lib, potential_parameter, coeffs_set, coeffs_kind, coeffs_pgf, radii_pgf, max_set, max_pgf, n_threads, i_thread, p_work)
calculates screening functions for schwarz screening
Definition hfx_screening_methods.F:304
hfx_screening_methods::update_pmax_mat
subroutine, public update_pmax_mat(pmax_set, map_atom_to_kind_atom, set_offset, atomic_block_offset, pmax_atom, full_density_alpha, full_density_beta, natom, kind_of, basis_parameter, nkind, is_assoc_atomic_block)
updates the maximum of the density matrix in compressed form for screening purposes
Definition hfx_screening_methods.F:138
hfx_screening_methods::calc_pair_dist_radii
subroutine, public calc_pair_dist_radii(qs_env, basis_parameter, radii_pgf, max_set, max_pgf, eps_schwarz, n_threads, i_thread)
calculates radius functions for longrange screening
Definition hfx_screening_methods.F:660
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:16
hfx_types::hfx_init_container
subroutine, public hfx_init_container(container, memory_usage, do_disk_storage)
This routine deletes all list entries in a container in order to deallocate the memory.
Definition hfx_types.F:2598
hfx_types::init_t_c_g0_lmax
integer, save, public init_t_c_g0_lmax
Definition hfx_types.F:136
hfx_types::log_zero
real(dp), parameter, public log_zero
Definition hfx_types.F:121
hfx_types::alloc_containers
subroutine, public alloc_containers(data, bin_size)
...
Definition hfx_types.F:2989
hfx_types::hfx_create_neighbor_cells
subroutine, public hfx_create_neighbor_cells(x_data, pbc_shells, cell, i_thread, nkp_grid)
This routine computes the neighbor cells that are taken into account in periodic runs
Definition hfx_types.F:2120
hfx_types::dealloc_containers
subroutine, public dealloc_containers(data, memory_usage)
...
Definition hfx_types.F:2957
hfx_types::hfx_reset_memory_usage_counter
subroutine, public hfx_reset_memory_usage_counter(memory_parameter, subtr_size_mb)
resets the maximum memory usage for a HFX calculation subtracting all relevant buffers from the input...
Definition hfx_types.F:2684
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::hfx_library_is_libint
integer, parameter, public hfx_library_is_libint
Definition input_constants.F:885
input_constants::hfx_library_is_both
integer, parameter, public hfx_library_is_both
Definition input_constants.F:885
input_constants::hfx_do_eval_energy
integer, parameter, public hfx_do_eval_energy
Definition input_constants.F:877
input_constants::hfx_library_is_libgint
integer, parameter, public hfx_library_is_libgint
Definition input_constants.F:885
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:855
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:855
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:855
input_constants::do_potential_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
Definition input_constants.F:855
input_constants::do_potential_tshpsc
integer, parameter, public do_potential_tshpsc
Definition input_constants.F:855
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:587
libgint_wrapper
Interface to the LibGint-Library.
Definition libgint_wrapper.F:14
libgint_wrapper::libgint_update_env
subroutine, public libgint_update_env(fac, memory_parameter, do_periodic, cell, actual_x_data, nneighbors, max_pgf, natom, kind_of, particle_set, basis_parameter)
Initialize and update the libGint computational environment. Must be called at least once after geo_c...
Definition libgint_wrapper.F:117
libgint_wrapper::cp_libgint_init
subroutine, public cp_libgint_init(actual_x_data)
Sets libGint internal enviroment, must be called at least once before libGint can be used.
Definition libgint_wrapper.F:72
libgint_wrapper::libgint_update_fock_matrix
subroutine, public libgint_update_fock_matrix(symm_fac, iatom, jatom, katom, latom, iset, jset, kset, lset, atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, nsgfa, nsgfb, nsgfc, nsgfd, la_min, la_max, lb_min, lb_max, lc_min, lc_max, ld_min, ld_max, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d)
The previous coulomb_4 function assigned an integral between primitive gaussian. Now we assign the ns...
Definition libgint_wrapper.F:509
libgint_wrapper::libgint_coulomb4
subroutine, public libgint_coulomb4(iatom, jatom, katom, latom, iset, jset, kset, lset, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, potential_parameter, screen1, screen2, log10_pmax, log10_eps_schwarz, pgf1, pgf2, neighbor_cells, cell, do_periodic, screened)
Assign two-electron integrals of a quartet/shell to libGint.
Definition libgint_wrapper.F:319
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:442
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:126
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:143
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
offload_api
Fortran API for the offload package, which is written in C.
Definition offload_api.F:12
offload_api::offload_set_chosen_device
subroutine, public offload_set_chosen_device(device_id)
Selects the chosen device to be used.
Definition offload_api.F:132
offload_api::offload_activate_chosen_device
subroutine, public offload_activate_chosen_device()
Activates the device selected via offload_set_chosen_device().
Definition offload_api.F:175
offload_api::offload_get_chosen_device
integer function, public offload_get_chosen_device()
Returns the chosen device.
Definition offload_api.F:152
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:58
t_c_g0::init
subroutine, public init(nder, iunit, mepos, group)
...
Definition t_c_g0.F:1361
util
All kind of helpful little routines.
Definition util.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:787
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
hfx_types::hfx_basis_info_type
Definition hfx_types.F:289
hfx_types::hfx_basis_type
Definition hfx_types.F:267
hfx_types::hfx_cache_type
Definition hfx_types.F:245
hfx_types::hfx_cell_type
Definition hfx_types.F:208
hfx_types::hfx_container_type
Definition hfx_types.F:257
hfx_types::hfx_distribution
Definition hfx_types.F:214
hfx_types::hfx_general_type
Definition hfx_types.F:201
hfx_types::hfx_load_balance_type
Definition hfx_types.F:191
hfx_types::hfx_memory_type
Definition hfx_types.F:160
hfx_types::hfx_p_kind
Definition hfx_types.F:301
hfx_types::hfx_pgf_list
Definition hfx_types.F:323
hfx_types::hfx_pgf_product_list
Definition hfx_types.F:334
hfx_types::hfx_potential_type
Definition hfx_types.F:141
hfx_types::hfx_screen_coeff_type
Definition hfx_types.F:296
hfx_types::hfx_screening_type
Definition hfx_types.F:152
hfx_types::hfx_task_list_type
Definition hfx_types.F:350
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:514
hfx_types::pair_list_type
Definition hfx_types.F:239
hfx_types::pair_set_list_type
Definition hfx_types.F:234
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:198
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:71
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137