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