(git:ed6f26b)
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-2025 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 cp_dbcsr_api, ONLY: dbcsr_copy, &
34 dbcsr_get_matrix_type, &
35 dbcsr_p_type, &
36 dbcsr_type_antisymmetric, &
37 dbcsr_dot_threadsafe
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_threadsafe(ks_matrix(ispin, img)%matrix, rho_ao(ispin, img)%matrix, etmp)
1929 ene_x_aa = ene_x_aa + etmp
1930 END DO
1931 !for ADMMS, we need the exchange matrix k(d) for both spins
1932 IF (dft_control%do_admm) THEN
1933 cpassert(nkimages == 1)
1934 CALL dbcsr_copy(matrix_ks_aux_fit_hfx(ispin)%matrix, ks_matrix(ispin, 1)%matrix, &
1935 name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1936 END IF
1937
1938 ene_x_bb = 0.0_dp
1939 IF (nspins == 2 .AND. .NOT. treat_lsd_in_core) THEN
1940 DO img = 1, nkimages
1941 CALL dbcsr_dot_threadsafe(ks_matrix(2, img)%matrix, rho_ao(2, img)%matrix, etmp)
1942 ene_x_bb = ene_x_bb + etmp
1943 END DO
1944 !for ADMMS, we need the exchange matrix k(d) for both spins
1945 IF (dft_control%do_admm) THEN
1946 cpassert(nkimages == 1)
1947 CALL dbcsr_copy(matrix_ks_aux_fit_hfx(2)%matrix, ks_matrix(2, 1)%matrix, &
1948 name="HF exch. part of matrix_ks_aux_fit for ADMMS")
1949 END IF
1950 END IF
1951
1952 !! Update energy type
1953 ehfx = 0.5_dp*(ene_x_aa + ene_x_bb)
1954 ELSE
1955 ! ** It is easier to correct the following expression by the diagonal energy contribution,
1956 ! ** than explicitly going throuhg the diagonal elements
1957 DO img = 1, nkimages
1958 DO pa = 1, SIZE(full_ks_alpha, 1)
1959 ene_x_aa = ene_x_aa + full_ks_alpha(pa, img)*full_density_alpha(pa, img)
1960 END DO
1961 END DO
1962 ! ** Now correct
1963 ene_x_aa = (ene_x_aa + ene_x_aa_diag)*0.5_dp
1964 IF (nspins == 2) THEN
1965 DO img = 1, nkimages
1966 DO pa = 1, SIZE(full_ks_beta, 1)
1967 ene_x_bb = ene_x_bb + full_ks_beta(pa, img)*full_density_beta(pa, img)
1968 END DO
1969 END DO
1970 ! ** Now correct
1971 ene_x_bb = (ene_x_bb + ene_x_bb_diag)*0.5_dp
1972 END IF
1973 CALL para_env%sum(ene_x_aa)
1974 IF (nspins == 2) CALL para_env%sum(ene_x_bb)
1975 ehfx = 0.5_dp*(ene_x_aa + ene_x_bb)
1976 END IF
1977
1978 !! Print some memeory information if this is the first step
1979 IF (my_geo_change) THEN
1980 tmp_i8(1:8) = (/shm_storage_counter_integrals, shm_neris_onthefly, shm_neris_incore, shm_neris_disk, &
1981 shm_neris_total, shm_stor_count_int_disk, shm_nprim_ints, shm_stor_count_max_val/)
1982 CALL para_env%sum(tmp_i8)
1983 shm_storage_counter_integrals = tmp_i8(1)
1984 shm_neris_onthefly = tmp_i8(2)
1985 shm_neris_incore = tmp_i8(3)
1986 shm_neris_disk = tmp_i8(4)
1987 shm_neris_total = tmp_i8(5)
1988 shm_stor_count_int_disk = tmp_i8(6)
1989 shm_nprim_ints = tmp_i8(7)
1990 shm_stor_count_max_val = tmp_i8(8)
1991 CALL para_env%max(memsize_after)
1992 mem_eris = (shm_storage_counter_integrals + 128*1024 - 1)/1024/128
1993 compression_factor = real(shm_neris_incore, dp)/real(shm_storage_counter_integrals, dp)
1994 mem_eris_disk = (shm_stor_count_int_disk + 128*1024 - 1)/1024/128
1995 compression_factor_disk = real(shm_neris_disk, dp)/real(shm_stor_count_int_disk, dp)
1996 mem_max_val = (shm_stor_count_max_val + 128*1024 - 1)/1024/128
1997
1998 IF (shm_neris_incore == 0) THEN
1999 mem_eris = 0
2000 compression_factor = 0.0_dp
2001 END IF
2002 IF (shm_neris_disk == 0) THEN
2003 mem_eris_disk = 0
2004 compression_factor_disk = 0.0_dp
2005 END IF
2006
2007 iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2008 extension=".scfLog")
2009 IF (iw > 0) THEN
2010 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2011 "HFX_MEM_INFO| Number of cart. primitive ERI's calculated: ", shm_nprim_ints
2012
2013 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2014 "HFX_MEM_INFO| Number of sph. ERI's calculated: ", shm_neris_total
2015
2016 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2017 "HFX_MEM_INFO| Number of sph. ERI's stored in-core: ", shm_neris_incore
2018
2019 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2020 "HFX_MEM_INFO| Number of sph. ERI's stored on disk: ", shm_neris_disk
2021
2022 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2023 "HFX_MEM_INFO| Number of sph. ERI's calculated on the fly: ", shm_neris_onthefly
2024
2025 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2026 "HFX_MEM_INFO| Total memory consumption ERI's RAM [MiB]: ", mem_eris
2027
2028 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2029 "HFX_MEM_INFO| Whereof max-vals [MiB]: ", mem_max_val
2030
2031 WRITE (unit=iw, fmt="((T3,A,T60,F21.2))") &
2032 "HFX_MEM_INFO| Total compression factor ERI's RAM: ", compression_factor
2033
2034 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2035 "HFX_MEM_INFO| Total memory consumption ERI's disk [MiB]: ", mem_eris_disk
2036
2037 WRITE (unit=iw, fmt="((T3,A,T60,F21.2))") &
2038 "HFX_MEM_INFO| Total compression factor ERI's disk: ", compression_factor_disk
2039
2040 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2041 "HFX_MEM_INFO| Size of density/Fock matrix [MiB]: ", 2_int_8*mb_size_p
2042
2043 IF (do_periodic) THEN
2044 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2045 "HFX_MEM_INFO| Size of buffers [MiB]: ", mb_size_buffers
2046 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2047 "HFX_MEM_INFO| Number of periodic image cells considered: ", SIZE(shm_master_x_data%neighbor_cells)
2048 ELSE
2049 WRITE (unit=iw, fmt="((T3,A,T60,I21))") &
2050 "HFX_MEM_INFO| Size of buffers [MiB]: ", mb_size_buffers
2051 END IF
2052 WRITE (unit=iw, fmt="((T3,A,T60,I21),/)") &
2053 "HFX_MEM_INFO| Est. max. program size after HFX [MiB]:", memsize_after/(1024*1024)
2054 CALL m_flush(iw)
2055 END IF
2056
2057 CALL cp_print_key_finished_output(iw, logger, hfx_section, &
2058 "HF_INFO")
2059 END IF
2060!$OMP END MASTER
2061
2062 !! flush caches if the geometry changed
2063 IF (do_dynamic_load_balancing) THEN
2064 my_bin_size = SIZE(actual_x_data%distribution_energy)
2065 ELSE
2066 my_bin_size = 1
2067 END IF
2068
2069 IF (my_geo_change) THEN
2070 IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
2071 DO bin = 1, my_bin_size
2072 maxval_cache => actual_x_data%store_ints%maxval_cache(bin)
2073 maxval_container => actual_x_data%store_ints%maxval_container(bin)
2074 integral_caches => actual_x_data%store_ints%integral_caches(:, bin)
2075 integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
2076 CALL hfx_flush_last_cache(bits_max_val, maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
2077 .false.)
2078 DO i = 1, 64
2079 CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
2080 memory_parameter%actual_memory_usage, .false.)
2081 END DO
2082 END DO
2083 END IF
2084 END IF
2085 !! reset all caches except we calculate all on the fly
2086 IF (.NOT. memory_parameter%do_all_on_the_fly) THEN
2087 DO bin = 1, my_bin_size
2088 maxval_cache => actual_x_data%store_ints%maxval_cache(bin)
2089 maxval_container => actual_x_data%store_ints%maxval_container(bin)
2090 integral_caches => actual_x_data%store_ints%integral_caches(:, bin)
2091 integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
2092
2093 CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_parameter%actual_memory_usage, .false.)
2094 DO i = 1, 64
2095 CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
2096 memory_parameter%actual_memory_usage, &
2097 .false.)
2098 END DO
2099 END DO
2100 END IF
2101
2102 !! Since the I/O routines are no thread-safe, i.e. the procedure to get the unit number, put a lock here
2103!$OMP CRITICAL(hfxenergy_out_critical)
2104 IF (do_disk_storage) THEN
2105 !! flush caches if the geometry changed
2106 IF (my_geo_change) THEN
2107 CALL hfx_flush_last_cache(bits_max_val, maxval_cache_disk, maxval_container_disk, &
2108 memory_parameter%actual_memory_usage_disk, .true.)
2109 DO i = 1, 64
2110 CALL hfx_flush_last_cache(i, integral_caches_disk(i), integral_containers_disk(i), &
2111 memory_parameter%actual_memory_usage_disk, .true.)
2112 END DO
2113 END IF
2114 !! reset all caches except we calculate all on the fly
2115 CALL hfx_reset_cache_and_container(maxval_cache_disk, maxval_container_disk, memory_parameter%actual_memory_usage_disk, &
2116 do_disk_storage)
2117 DO i = 1, 64
2118 CALL hfx_reset_cache_and_container(integral_caches_disk(i), integral_containers_disk(i), &
2119 memory_parameter%actual_memory_usage_disk, do_disk_storage)
2120 END DO
2121 END IF
2122!$OMP END CRITICAL(hfxenergy_out_critical)
2123!$OMP BARRIER
2124 !! Clean up
2125 DEALLOCATE (last_sgf_global)
2126!$OMP MASTER
2127 DEALLOCATE (full_density_alpha)
2128 IF (.NOT. treat_lsd_in_core) THEN
2129 IF (nspins == 2) THEN
2130 DEALLOCATE (full_density_beta)
2131 END IF
2132 END IF
2133 IF (do_dynamic_load_balancing) THEN
2134 DEALLOCATE (shm_master_x_data%task_list)
2135 END IF
2136!$OMP END MASTER
2137 DEALLOCATE (pbd_buf, pbc_buf, pad_buf, pac_buf)
2138 DEALLOCATE (kbd_buf, kbc_buf, kad_buf, kac_buf)
2139 DEALLOCATE (set_list_ij, set_list_kl)
2140
2141 DO i = 1, max_pgf**2
2142 DEALLOCATE (pgf_list_ij(i)%image_list)
2143 DEALLOCATE (pgf_list_kl(i)%image_list)
2144 END DO
2145
2146 DEALLOCATE (pgf_list_ij)
2147 DEALLOCATE (pgf_list_kl)
2148 DEALLOCATE (pgf_product_list)
2149
2150 DEALLOCATE (max_contraction, kind_of)
2151
2152 DEALLOCATE (ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp)
2153
2154 DEALLOCATE (nimages)
2155
2156!$OMP BARRIER
2157!$OMP END PARALLEL
2158
2159 CALL timestop(handle)
2160 END SUBROUTINE integrate_four_center
2161
2162! **************************************************************************************************
2163!> \brief calculates two-electron integrals of a quartet/shell using the library
2164!> lib_int in the periodic case
2165!> \param lib ...
2166!> \param ra ...
2167!> \param rb ...
2168!> \param rc ...
2169!> \param rd ...
2170!> \param npgfa ...
2171!> \param npgfb ...
2172!> \param npgfc ...
2173!> \param npgfd ...
2174!> \param la_min ...
2175!> \param la_max ...
2176!> \param lb_min ...
2177!> \param lb_max ...
2178!> \param lc_min ...
2179!> \param lc_max ...
2180!> \param ld_min ...
2181!> \param ld_max ...
2182!> \param nsgfa ...
2183!> \param nsgfb ...
2184!> \param nsgfc ...
2185!> \param nsgfd ...
2186!> \param sphi_a_u1 ...
2187!> \param sphi_a_u2 ...
2188!> \param sphi_a_u3 ...
2189!> \param sphi_b_u1 ...
2190!> \param sphi_b_u2 ...
2191!> \param sphi_b_u3 ...
2192!> \param sphi_c_u1 ...
2193!> \param sphi_c_u2 ...
2194!> \param sphi_c_u3 ...
2195!> \param sphi_d_u1 ...
2196!> \param sphi_d_u2 ...
2197!> \param sphi_d_u3 ...
2198!> \param zeta ...
2199!> \param zetb ...
2200!> \param zetc ...
2201!> \param zetd ...
2202!> \param primitive_integrals array of primitive_integrals
2203!> \param potential_parameter contains info for libint
2204!> \param neighbor_cells Periodic images
2205!> \param screen1 set based coefficients for near field screening
2206!> \param screen2 set based coefficients for near field screening
2207!> \param eps_schwarz threshold
2208!> \param max_contraction_val maximum multiplication factor for cart -> sph
2209!> \param cart_estimate maximum calculated integral value
2210!> \param cell cell
2211!> \param neris_tmp counter for calculated cart integrals
2212!> \param log10_pmax logarithm of initial p matrix max element
2213!> \param log10_eps_schwarz log of threshold
2214!> \param R1_pgf coefficients for radii of product distribution function
2215!> \param R2_pgf coefficients for radii of product distribution function
2216!> \param pgf1 schwarz coefficients pgf basid
2217!> \param pgf2 schwarz coefficients pgf basid
2218!> \param pgf_list_ij ...
2219!> \param pgf_list_kl ...
2220!> \param pgf_product_list ...
2221!> \param nsgfl_a ...
2222!> \param nsgfl_b ...
2223!> \param nsgfl_c ...
2224!> \param nsgfl_d ...
2225!> \param sphi_a_ext ...
2226!> \param sphi_b_ext ...
2227!> \param sphi_c_ext ...
2228!> \param sphi_d_ext ...
2229!> \param ee_work ...
2230!> \param ee_work2 ...
2231!> \param ee_buffer1 ...
2232!> \param ee_buffer2 ...
2233!> \param ee_primitives_tmp ...
2234!> \param nimages ...
2235!> \param do_periodic ...
2236!> \param p_work ...
2237!> \par History
2238!> 11.2006 created [Manuel Guidon]
2239!> 02.2009 completely rewritten screening part [Manuel Guidon]
2240!> \author Manuel Guidon
2241! **************************************************************************************************
2242 SUBROUTINE coulomb4(lib, ra, rb, rc, rd, npgfa, npgfb, npgfc, npgfd, &
2243 la_min, la_max, lb_min, lb_max, &
2244 lc_min, lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, &
2245 sphi_a_u1, sphi_a_u2, sphi_a_u3, &
2246 sphi_b_u1, sphi_b_u2, sphi_b_u3, &
2247 sphi_c_u1, sphi_c_u2, sphi_c_u3, &
2248 sphi_d_u1, sphi_d_u2, sphi_d_u3, &
2249 zeta, zetb, zetc, zetd, &
2250 primitive_integrals, &
2251 potential_parameter, neighbor_cells, &
2252 screen1, screen2, eps_schwarz, max_contraction_val, &
2253 cart_estimate, cell, neris_tmp, log10_pmax, &
2254 log10_eps_schwarz, R1_pgf, R2_pgf, pgf1, pgf2, &
2255 pgf_list_ij, pgf_list_kl, &
2256 pgf_product_list, &
2257 nsgfl_a, nsgfl_b, nsgfl_c, &
2258 nsgfl_d, &
2259 sphi_a_ext, sphi_b_ext, sphi_c_ext, sphi_d_ext, &
2260 ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
2261 nimages, do_periodic, p_work)
2262
2263 TYPE(cp_libint_t) :: lib
2264 REAL(dp), INTENT(IN) :: ra(3), rb(3), rc(3), rd(3)
2265 INTEGER, INTENT(IN) :: npgfa, npgfb, npgfc, npgfd, la_min, la_max, lb_min, lb_max, lc_min, &
2266 lc_max, ld_min, ld_max, nsgfa, nsgfb, nsgfc, nsgfd, sphi_a_u1, sphi_a_u2, sphi_a_u3, &
2267 sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, &
2268 sphi_d_u3
2269 REAL(dp), DIMENSION(1:npgfa), INTENT(IN) :: zeta
2270 REAL(dp), DIMENSION(1:npgfb), INTENT(IN) :: zetb
2271 REAL(dp), DIMENSION(1:npgfc), INTENT(IN) :: zetc
2272 REAL(dp), DIMENSION(1:npgfd), INTENT(IN) :: zetd
2273 REAL(dp), DIMENSION(nsgfa, nsgfb, nsgfc, nsgfd) :: primitive_integrals
2274 TYPE(hfx_potential_type) :: potential_parameter
2275 TYPE(hfx_cell_type), DIMENSION(:), POINTER :: neighbor_cells
2276 REAL(dp), INTENT(IN) :: screen1(2), screen2(2), eps_schwarz, &
2277 max_contraction_val
2278 REAL(dp) :: cart_estimate
2279 TYPE(cell_type), POINTER :: cell
2280 INTEGER(int_8) :: neris_tmp
2281 REAL(dp), INTENT(IN) :: log10_pmax, log10_eps_schwarz
2282 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
2283 POINTER :: r1_pgf, r2_pgf, pgf1, pgf2
2284 TYPE(hfx_pgf_list), DIMENSION(*) :: pgf_list_ij, pgf_list_kl
2285 TYPE(hfx_pgf_product_list), ALLOCATABLE, &
2286 DIMENSION(:), INTENT(INOUT) :: pgf_product_list
2287 INTEGER, DIMENSION(0:), INTENT(IN) :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d
2288 REAL(dp), INTENT(IN) :: sphi_a_ext(sphi_a_u1, sphi_a_u2, sphi_a_u3), &
2289 sphi_b_ext(sphi_b_u1, sphi_b_u2, sphi_b_u3), sphi_c_ext(sphi_c_u1, sphi_c_u2, sphi_c_u3), &
2290 sphi_d_ext(sphi_d_u1, sphi_d_u2, sphi_d_u3)
2291 REAL(dp), DIMENSION(*) :: ee_work, ee_work2, ee_buffer1, &
2292 ee_buffer2, ee_primitives_tmp
2293 INTEGER, DIMENSION(*) :: nimages
2294 LOGICAL, INTENT(IN) :: do_periodic
2295 REAL(dp), DIMENSION(:), POINTER :: p_work
2296
2297 INTEGER :: ipgf, jpgf, kpgf, la, lb, lc, ld, list_ij, list_kl, lpgf, max_l, ncoa, ncob, &
2298 ncoc, ncod, nelements_ij, nelements_kl, nproducts, nsgfla, nsgflb, nsgflc, nsgfld, nsoa, &
2299 nsob, nsoc, nsod, s_offset_a, s_offset_b, s_offset_c, s_offset_d
2300 REAL(dp) :: etainv, tmp_max, zetainv
2301
2302 CALL build_pair_list_pgf(npgfa, npgfb, pgf_list_ij, zeta, zetb, screen1, screen2, &
2303 pgf1, r1_pgf, log10_pmax, log10_eps_schwarz, ra, rb, &
2304 nelements_ij, &
2305 neighbor_cells, nimages, do_periodic)
2306 CALL build_pair_list_pgf(npgfc, npgfd, pgf_list_kl, zetc, zetd, screen2, screen1, &
2307 pgf2, r2_pgf, log10_pmax, log10_eps_schwarz, rc, rd, &
2308 nelements_kl, &
2309 neighbor_cells, nimages, do_periodic)
2310
2311 cart_estimate = 0.0_dp
2312 neris_tmp = 0
2313 primitive_integrals = 0.0_dp
2314 max_l = la_max + lb_max + lc_max + ld_max
2315 DO list_ij = 1, nelements_ij
2316 zetainv = pgf_list_ij(list_ij)%ZetaInv
2317 ipgf = pgf_list_ij(list_ij)%ipgf
2318 jpgf = pgf_list_ij(list_ij)%jpgf
2319
2320 DO list_kl = 1, nelements_kl
2321 etainv = pgf_list_kl(list_kl)%ZetaInv
2322 kpgf = pgf_list_kl(list_kl)%ipgf
2323 lpgf = pgf_list_kl(list_kl)%jpgf
2324
2325 CALL build_pgf_product_list(pgf_list_ij(list_ij), pgf_list_kl(list_kl), pgf_product_list, &
2326 nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, &
2327 potential_parameter, max_l, do_periodic)
2328
2329 s_offset_a = 0
2330 DO la = la_min, la_max
2331 s_offset_b = 0
2332 ncoa = nco(la)
2333 nsgfla = nsgfl_a(la)
2334 nsoa = nso(la)
2335
2336 DO lb = lb_min, lb_max
2337 s_offset_c = 0
2338 ncob = nco(lb)
2339 nsgflb = nsgfl_b(lb)
2340 nsob = nso(lb)
2341
2342 DO lc = lc_min, lc_max
2343 s_offset_d = 0
2344 ncoc = nco(lc)
2345 nsgflc = nsgfl_c(lc)
2346 nsoc = nso(lc)
2347
2348 DO ld = ld_min, ld_max
2349 ncod = nco(ld)
2350 nsgfld = nsgfl_d(ld)
2351 nsod = nso(ld)
2352
2353 tmp_max = 0.0_dp
2354 CALL evaluate_eri(lib, nproducts, pgf_product_list, &
2355 la, lb, lc, ld, &
2356 ncoa, ncob, ncoc, ncod, &
2357 nsgfa, nsgfb, nsgfc, nsgfd, &
2358 primitive_integrals, &
2359 max_contraction_val, tmp_max, eps_schwarz, &
2360 neris_tmp, zetainv, etainv, &
2361 s_offset_a, s_offset_b, s_offset_c, s_offset_d, &
2362 nsgfla, nsgflb, nsgflc, nsgfld, nsoa, nsob, nsoc, nsod, &
2363 sphi_a_ext(1, la + 1, ipgf), &
2364 sphi_b_ext(1, lb + 1, jpgf), &
2365 sphi_c_ext(1, lc + 1, kpgf), &
2366 sphi_d_ext(1, ld + 1, lpgf), &
2367 ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, &
2368 p_work)
2369 cart_estimate = max(tmp_max, cart_estimate)
2370 s_offset_d = s_offset_d + nsod*nsgfld
2371 END DO !ld
2372 s_offset_c = s_offset_c + nsoc*nsgflc
2373 END DO !lc
2374 s_offset_b = s_offset_b + nsob*nsgflb
2375 END DO !lb
2376 s_offset_a = s_offset_a + nsoa*nsgfla
2377 END DO !la
2378 END DO
2379 END DO
2380
2381 END SUBROUTINE coulomb4
2382
2383! **************************************************************************************************
2384!> \brief Given a 2d index pair, this function returns a 1d index pair for
2385!> a symmetric upper triangle NxN matrix
2386!> The compiler should inline this function, therefore it appears in
2387!> several modules
2388!> \param i 2d index
2389!> \param j 2d index
2390!> \param N matrix size
2391!> \return ...
2392!> \par History
2393!> 03.2009 created [Manuel Guidon]
2394!> \author Manuel Guidon
2395! **************************************************************************************************
2396 PURE FUNCTION get_1d_idx(i, j, N)
2397 INTEGER, INTENT(IN) :: i, j
2398 INTEGER(int_8), INTENT(IN) :: n
2399 INTEGER(int_8) :: get_1d_idx
2400
2401 INTEGER(int_8) :: min_ij
2402
2403 min_ij = min(i, j)
2404 get_1d_idx = min_ij*n + max(i, j) - (min_ij - 1)*min_ij/2 - n
2405
2406 END FUNCTION get_1d_idx
2407
2408! **************************************************************************************************
2409!> \brief This routine prefetches density/fock matrix elements and stores them
2410!> in cache friendly arrays. These buffers are then used to update the
2411!> fock matrix
2412!> \param ma_max Size of matrix blocks
2413!> \param mb_max Size of matrix blocks
2414!> \param mc_max Size of matrix blocks
2415!> \param md_max Size of matrix blocks
2416!> \param fac multiplication factor (spin)
2417!> \param symm_fac multiplication factor (symmetry)
2418!> \param density upper triangular density matrix
2419!> \param ks upper triangular fock matrix
2420!> \param prim primitive integrals
2421!> \param pbd buffer that will contain P(b,d)
2422!> \param pbc buffer that will contain P(b,c)
2423!> \param pad buffer that will contain P(a,d)
2424!> \param pac buffer that will contain P(a,c)
2425!> \param kbd buffer for KS(b,d)
2426!> \param kbc buffer for KS(b,c)
2427!> \param kad buffer for KS(a,d)
2428!> \param kac buffer for KS(a,c)
2429!> \param iatom ...
2430!> \param jatom ...
2431!> \param katom ...
2432!> \param latom ...
2433!> \param iset ...
2434!> \param jset ...
2435!> \param kset ...
2436!> \param lset ...
2437!> \param offset_bd_set ...
2438!> \param offset_bc_set ...
2439!> \param offset_ad_set ...
2440!> \param offset_ac_set ...
2441!> \param atomic_offset_bd ...
2442!> \param atomic_offset_bc ...
2443!> \param atomic_offset_ad ...
2444!> \param atomic_offset_ac ...
2445!> \par History
2446!> 03.2009 created [Manuel Guidon]
2447!> \author Manuel Guidon
2448! **************************************************************************************************
2449
2450 SUBROUTINE update_fock_matrix(ma_max, mb_max, mc_max, md_max, &
2451 fac, symm_fac, density, ks, prim, &
2452 pbd, pbc, pad, pac, kbd, kbc, kad, kac, &
2453 iatom, jatom, katom, latom, &
2454 iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, &
2455 offset_ac_set, atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, &
2456 atomic_offset_ac)
2457
2458 INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
2459 REAL(dp), INTENT(IN) :: fac, symm_fac
2460 REAL(dp), DIMENSION(:), INTENT(IN) :: density
2461 REAL(dp), DIMENSION(:), INTENT(INOUT) :: ks
2462 REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
2463 INTENT(IN) :: prim
2464 REAL(dp), DIMENSION(*), INTENT(INOUT) :: pbd, pbc, pad, pac, kbd, kbc, kad, kac
2465 INTEGER, INTENT(IN) :: iatom, jatom, katom, latom, iset, jset, &
2466 kset, lset
2467 INTEGER, DIMENSION(:, :), POINTER, INTENT(IN) :: offset_bd_set, offset_bc_set, &
2468 offset_ad_set, offset_ac_set
2469 INTEGER, INTENT(IN) :: atomic_offset_bd, atomic_offset_bc, &
2470 atomic_offset_ad, atomic_offset_ac
2471
2472 INTEGER :: i, j, ma, mb, mc, md, offset_ac, &
2473 offset_ad, offset_bc, offset_bd
2474 REAL(dp) :: ki
2475
2476 IF (jatom >= latom) THEN
2477 i = 1
2478 offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
2479 j = offset_bd
2480 DO md = 1, md_max
2481 DO mb = 1, mb_max
2482 pbd(i) = density(j)
2483 i = i + 1
2484 j = j + 1
2485 END DO
2486 END DO
2487 ELSE
2488 i = 1
2489 offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
2490 DO md = 1, md_max
2491 j = offset_bd + md - 1
2492 DO mb = 1, mb_max
2493 pbd(i) = density(j)
2494 i = i + 1
2495 j = j + md_max
2496 END DO
2497 END DO
2498 END IF
2499 IF (jatom >= katom) THEN
2500 i = 1
2501 offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
2502 j = offset_bc
2503 DO mc = 1, mc_max
2504 DO mb = 1, mb_max
2505 pbc(i) = density(j)
2506 i = i + 1
2507 j = j + 1
2508 END DO
2509 END DO
2510 ELSE
2511 i = 1
2512 offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
2513 DO mc = 1, mc_max
2514 j = offset_bc + mc - 1
2515 DO mb = 1, mb_max
2516 pbc(i) = density(j)
2517 i = i + 1
2518 j = j + mc_max
2519 END DO
2520 END DO
2521 END IF
2522 IF (iatom >= latom) THEN
2523 i = 1
2524 offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
2525 j = offset_ad
2526 DO md = 1, md_max
2527 DO ma = 1, ma_max
2528 pad(i) = density(j)
2529 i = i + 1
2530 j = j + 1
2531 END DO
2532 END DO
2533 ELSE
2534 i = 1
2535 offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
2536 DO md = 1, md_max
2537 j = offset_ad + md - 1
2538 DO ma = 1, ma_max
2539 pad(i) = density(j)
2540 i = i + 1
2541 j = j + md_max
2542 END DO
2543 END DO
2544 END IF
2545 IF (iatom >= katom) THEN
2546 i = 1
2547 offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
2548 j = offset_ac
2549 DO mc = 1, mc_max
2550 DO ma = 1, ma_max
2551 pac(i) = density(j)
2552 i = i + 1
2553 j = j + 1
2554 END DO
2555 END DO
2556 ELSE
2557 i = 1
2558 offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
2559 DO mc = 1, mc_max
2560 j = offset_ac + mc - 1
2561 DO ma = 1, ma_max
2562 pac(i) = density(j)
2563 i = i + 1
2564 j = j + mc_max
2565 END DO
2566 END DO
2567 END IF
2568
2569 CALL contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, fac*symm_fac)
2570 IF (jatom >= latom) THEN
2571 i = 1
2572 j = offset_bd
2573 DO md = 1, md_max
2574 DO mb = 1, mb_max
2575 ki = kbd(i)
2576!$OMP ATOMIC
2577 ks(j) = ks(j) + ki
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 ki = kbd(i)
2588!$OMP ATOMIC
2589 ks(j) = ks(j) + ki
2590 i = i + 1
2591 j = j + md_max
2592 END DO
2593 END DO
2594 END IF
2595 IF (jatom >= katom) THEN
2596 i = 1
2597 j = offset_bc
2598 DO mc = 1, mc_max
2599 DO mb = 1, mb_max
2600 ki = kbc(i)
2601!$OMP ATOMIC
2602 ks(j) = ks(j) + ki
2603 i = i + 1
2604 j = j + 1
2605 END DO
2606 END DO
2607 ELSE
2608 i = 1
2609 DO mc = 1, mc_max
2610 j = offset_bc + mc - 1
2611 DO mb = 1, mb_max
2612 ki = kbc(i)
2613!$OMP ATOMIC
2614 ks(j) = ks(j) + ki
2615 i = i + 1
2616 j = j + mc_max
2617 END DO
2618 END DO
2619 END IF
2620 IF (iatom >= latom) THEN
2621 i = 1
2622 j = offset_ad
2623 DO md = 1, md_max
2624 DO ma = 1, ma_max
2625 ki = kad(i)
2626!$OMP ATOMIC
2627 ks(j) = ks(j) + ki
2628 i = i + 1
2629 j = j + 1
2630 END DO
2631 END DO
2632 ELSE
2633 i = 1
2634 DO md = 1, md_max
2635 j = offset_ad + md - 1
2636 DO ma = 1, ma_max
2637 ki = kad(i)
2638!$OMP ATOMIC
2639 ks(j) = ks(j) + ki
2640 i = i + 1
2641 j = j + md_max
2642 END DO
2643 END DO
2644 END IF
2645 IF (iatom >= katom) THEN
2646 i = 1
2647 j = offset_ac
2648 DO mc = 1, mc_max
2649 DO ma = 1, ma_max
2650 ki = kac(i)
2651!$OMP ATOMIC
2652 ks(j) = ks(j) + ki
2653 i = i + 1
2654 j = j + 1
2655 END DO
2656 END DO
2657 ELSE
2658 i = 1
2659 DO mc = 1, mc_max
2660 j = offset_ac + mc - 1
2661 DO ma = 1, ma_max
2662 ki = kac(i)
2663!$OMP ATOMIC
2664 ks(j) = ks(j) + ki
2665 i = i + 1
2666 j = j + mc_max
2667 END DO
2668 END DO
2669 END IF
2670 END SUBROUTINE update_fock_matrix
2671
2672! **************************************************************************************************
2673!> \brief ...
2674!> \param ma_max ...
2675!> \param mb_max ...
2676!> \param mc_max ...
2677!> \param md_max ...
2678!> \param fac ...
2679!> \param symm_fac ...
2680!> \param density ...
2681!> \param ks ...
2682!> \param prim ...
2683!> \param pbd ...
2684!> \param pbc ...
2685!> \param pad ...
2686!> \param pac ...
2687!> \param kbd ...
2688!> \param kbc ...
2689!> \param kad ...
2690!> \param kac ...
2691!> \param iatom ...
2692!> \param jatom ...
2693!> \param katom ...
2694!> \param latom ...
2695!> \param iset ...
2696!> \param jset ...
2697!> \param kset ...
2698!> \param lset ...
2699!> \param offset_bd_set ...
2700!> \param offset_bc_set ...
2701!> \param offset_ad_set ...
2702!> \param offset_ac_set ...
2703!> \param atomic_offset_bd ...
2704!> \param atomic_offset_bc ...
2705!> \param atomic_offset_ad ...
2706!> \param atomic_offset_ac ...
2707! **************************************************************************************************
2708 SUBROUTINE update_fock_matrix_as(ma_max, mb_max, mc_max, md_max, &
2709 fac, symm_fac, density, ks, prim, &
2710 pbd, pbc, pad, pac, kbd, kbc, kad, kac, &
2711 iatom, jatom, katom, latom, &
2712 iset, jset, kset, lset, offset_bd_set, offset_bc_set, offset_ad_set, &
2713 offset_ac_set, atomic_offset_bd, atomic_offset_bc, atomic_offset_ad, &
2714 atomic_offset_ac)
2715
2716 INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
2717 REAL(dp), INTENT(IN) :: fac, symm_fac
2718 REAL(dp), DIMENSION(:), INTENT(IN) :: density
2719 REAL(dp), DIMENSION(:), INTENT(INOUT) :: ks
2720 REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
2721 INTENT(IN) :: prim
2722 REAL(dp), DIMENSION(*), INTENT(INOUT) :: pbd, pbc, pad, pac, kbd, kbc, kad, kac
2723 INTEGER, INTENT(IN) :: iatom, jatom, katom, latom, iset, jset, &
2724 kset, lset
2725 INTEGER, DIMENSION(:, :), POINTER :: offset_bd_set, offset_bc_set, &
2726 offset_ad_set, offset_ac_set
2727 INTEGER, INTENT(IN) :: atomic_offset_bd, atomic_offset_bc, &
2728 atomic_offset_ad, atomic_offset_ac
2729
2730 INTEGER :: i, j, ma, mb, mc, md, offset_ac, &
2731 offset_ad, offset_bc, offset_bd
2732
2733 IF (jatom >= latom) THEN
2734 i = 1
2735 offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
2736 j = offset_bd
2737 DO md = 1, md_max
2738 DO mb = 1, mb_max
2739 pbd(i) = +density(j)
2740 i = i + 1
2741 j = j + 1
2742 END DO
2743 END DO
2744 ELSE
2745 i = 1
2746 offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
2747 DO md = 1, md_max
2748 j = offset_bd + md - 1
2749 DO mb = 1, mb_max
2750 pbd(i) = -density(j)
2751 i = i + 1
2752 j = j + md_max
2753 END DO
2754 END DO
2755 END IF
2756 IF (jatom >= katom) THEN
2757 i = 1
2758 offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
2759 j = offset_bc
2760 DO mc = 1, mc_max
2761 DO mb = 1, mb_max
2762 pbc(i) = -density(j)
2763 i = i + 1
2764 j = j + 1
2765 END DO
2766 END DO
2767 ELSE
2768 i = 1
2769 offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
2770 DO mc = 1, mc_max
2771 j = offset_bc + mc - 1
2772 DO mb = 1, mb_max
2773 pbc(i) = density(j)
2774 i = i + 1
2775 j = j + mc_max
2776 END DO
2777 END DO
2778 END IF
2779 IF (iatom >= latom) THEN
2780 i = 1
2781 offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
2782 j = offset_ad
2783 DO md = 1, md_max
2784 DO ma = 1, ma_max
2785 pad(i) = -density(j)
2786 i = i + 1
2787 j = j + 1
2788 END DO
2789 END DO
2790 ELSE
2791 i = 1
2792 offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
2793 DO md = 1, md_max
2794 j = offset_ad + md - 1
2795 DO ma = 1, ma_max
2796 pad(i) = density(j)
2797 i = i + 1
2798 j = j + md_max
2799 END DO
2800 END DO
2801 END IF
2802 IF (iatom >= katom) THEN
2803 i = 1
2804 offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
2805 j = offset_ac
2806 DO mc = 1, mc_max
2807 DO ma = 1, ma_max
2808 pac(i) = +density(j)
2809 i = i + 1
2810 j = j + 1
2811 END DO
2812 END DO
2813 ELSE
2814 i = 1
2815 offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
2816 DO mc = 1, mc_max
2817 j = offset_ac + mc - 1
2818 DO ma = 1, ma_max
2819 pac(i) = -density(j)
2820 i = i + 1
2821 j = j + mc_max
2822 END DO
2823 END DO
2824 END IF
2825
2826 CALL contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, fac*symm_fac)
2827
2828 IF (jatom >= latom) THEN
2829 i = 1
2830 j = offset_bd
2831 DO md = 1, md_max
2832 DO mb = 1, mb_max
2833!$OMP ATOMIC
2834 ks(j) = ks(j) + kbd(i)
2835 i = i + 1
2836 j = j + 1
2837 END DO
2838 END DO
2839 ELSE
2840 i = 1
2841 DO md = 1, md_max
2842 j = offset_bd + md - 1
2843 DO mb = 1, mb_max
2844!$OMP ATOMIC
2845 ks(j) = ks(j) - kbd(i)
2846 i = i + 1
2847 j = j + md_max
2848 END DO
2849 END DO
2850 END IF
2851 IF (jatom >= katom) THEN
2852 i = 1
2853 j = offset_bc
2854 DO mc = 1, mc_max
2855 DO mb = 1, mb_max
2856!$OMP ATOMIC
2857 ks(j) = ks(j) - kbc(i)
2858 i = i + 1
2859 j = j + 1
2860 END DO
2861 END DO
2862 ELSE
2863 i = 1
2864 DO mc = 1, mc_max
2865 j = offset_bc + mc - 1
2866 DO mb = 1, mb_max
2867!$OMP ATOMIC
2868 ks(j) = ks(j) + kbc(i)
2869 i = i + 1
2870 j = j + mc_max
2871 END DO
2872 END DO
2873 END IF
2874 IF (iatom >= latom) THEN
2875 i = 1
2876 j = offset_ad
2877 DO md = 1, md_max
2878 DO ma = 1, ma_max
2879!$OMP ATOMIC
2880 ks(j) = ks(j) - kad(i)
2881 i = i + 1
2882 j = j + 1
2883 END DO
2884 END DO
2885 ELSE
2886 i = 1
2887 DO md = 1, md_max
2888 j = offset_ad + md - 1
2889 DO ma = 1, ma_max
2890!$OMP ATOMIC
2891 ks(j) = ks(j) + kad(i)
2892 i = i + 1
2893 j = j + md_max
2894 END DO
2895 END DO
2896 END IF
2897! XXXXXXXXXXXXXXXXXXXXXXXX
2898 IF (iatom >= katom) THEN
2899 i = 1
2900 j = offset_ac
2901 DO mc = 1, mc_max
2902 DO ma = 1, ma_max
2903!$OMP ATOMIC
2904 ks(j) = ks(j) + kac(i)
2905 i = i + 1
2906 j = j + 1
2907 END DO
2908 END DO
2909 ELSE
2910 i = 1
2911 DO mc = 1, mc_max
2912 j = offset_ac + mc - 1
2913 DO ma = 1, ma_max
2914!$OMP ATOMIC
2915 ks(j) = ks(j) - kac(i)
2916 i = i + 1
2917 j = j + mc_max
2918 END DO
2919 END DO
2920 END IF
2921
2922 END SUBROUTINE update_fock_matrix_as
2923
2924! **************************************************************************************************
2925!> \brief ...
2926!> \param i ...
2927!> \param j ...
2928!> \param k ...
2929!> \param l ...
2930!> \param set_offsets ...
2931!> \param atom_offsets ...
2932!> \param iset ...
2933!> \param jset ...
2934!> \param kset ...
2935!> \param lset ...
2936!> \param ma_max ...
2937!> \param mb_max ...
2938!> \param mc_max ...
2939!> \param md_max ...
2940!> \param prim ...
2941! **************************************************************************************************
2942 SUBROUTINE print_integrals(i, j, k, l, set_offsets, atom_offsets, iset, jset, kset, lset, ma_max, mb_max, mc_max, md_max, prim)
2943 INTEGER :: i, j, k, l
2944 INTEGER, DIMENSION(:, :, :, :), POINTER :: set_offsets
2945 INTEGER, DIMENSION(:, :), POINTER :: atom_offsets
2946 INTEGER :: iset, jset, kset, lset, ma_max, mb_max, &
2947 mc_max, md_max
2948 REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
2949 INTENT(IN) :: prim
2950
2951 INTEGER :: iint, ma, mb, mc, md
2952
2953 iint = 0
2954 DO md = 1, md_max
2955 DO mc = 1, mc_max
2956 DO mb = 1, mb_max
2957 DO ma = 1, ma_max
2958 iint = iint + 1
2959 IF (abs(prim(iint)) .GT. 0.0000000000001) &
2960 WRITE (99, *) atom_offsets(i, 1) + ma + set_offsets(iset, 1, i, 1) - 1, &
2961 atom_offsets(j, 1) + ma + set_offsets(jset, 1, j, 1) - 1, &
2962 atom_offsets(k, 1) + ma + set_offsets(kset, 1, k, 1) - 1, &
2963 atom_offsets(l, 1) + ma + set_offsets(lset, 1, l, 1) - 1, &
2964 prim(iint)
2965 END DO
2966 END DO
2967 END DO
2968 END DO
2969
2970 END SUBROUTINE print_integrals
2971
2972
2973! **************************************************************************************************
2974!> \brief This routine calculates the maximum density matrix element, when
2975!> screening on an initial density matrix is applied. Due to symmetry of
2976!> the ERI's, there are always 4 matrix elements to be considered.
2977!> CASE 0-15 belong to an energy calculation (linear screening)
2978!> CASE 16-31 belong to a force calculation (square screening)
2979!> \param ptr_p_1 Pointers to atomic density matrices
2980!> \param ptr_p_2 Pointers to atomic density matrices
2981!> \param ptr_p_3 Pointers to atomic density matrices
2982!> \param ptr_p_4 Pointers to atomic density matrices
2983!> \param iset Current set
2984!> \param jset Current set
2985!> \param kset Current set
2986!> \param lset Current set
2987!> \param pmax_val value to be calculated
2988!> \param swap_id Defines how the matrices are accessed
2989!> \par History
2990!> 06.2009 created [Manuel Guidon]
2991!> \author Manuel Guidon
2992! **************************************************************************************************
2993PURE SUBROUTINE get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, iset, jset, kset, lset, pmax_val, swap_id)
2994
2995 REAL(dp), DIMENSION(:, :), POINTER :: ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4
2996 INTEGER, INTENT(IN) :: iset, jset, kset, lset
2997
2998 REAL(dp), INTENT(OUT) :: pmax_val
2999 INTEGER, INTENT(IN) :: swap_id
3000
3001 REAL(dp) :: pmax_1, pmax_2, pmax_3, pmax_4
3002
3003 SELECT CASE (swap_id)
3004 CASE (0)
3005 pmax_1 = ptr_p_1(kset, iset)
3006 pmax_2 = ptr_p_2(lset, jset)
3007 pmax_3 = ptr_p_3(lset, iset)
3008 pmax_4 = ptr_p_4(kset, jset)
3009 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3010 CASE (1)
3011 pmax_1 = ptr_p_1(iset, kset)
3012 pmax_2 = ptr_p_2(lset, jset)
3013 pmax_3 = ptr_p_3(lset, iset)
3014 pmax_4 = ptr_p_4(kset, jset)
3015 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3016 CASE (2)
3017 pmax_1 = ptr_p_1(kset, iset)
3018 pmax_2 = ptr_p_2(jset, lset)
3019 pmax_3 = ptr_p_3(lset, iset)
3020 pmax_4 = ptr_p_4(kset, jset)
3021 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3022 CASE (3)
3023 pmax_1 = ptr_p_1(iset, kset)
3024 pmax_2 = ptr_p_2(jset, lset)
3025 pmax_3 = ptr_p_3(lset, iset)
3026 pmax_4 = ptr_p_4(kset, jset)
3027 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3028 CASE (4)
3029 pmax_1 = ptr_p_1(kset, iset)
3030 pmax_2 = ptr_p_2(lset, jset)
3031 pmax_3 = ptr_p_3(iset, lset)
3032 pmax_4 = ptr_p_4(kset, jset)
3033 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3034 CASE (5)
3035 pmax_1 = ptr_p_1(iset, kset)
3036 pmax_2 = ptr_p_2(lset, jset)
3037 pmax_3 = ptr_p_3(iset, lset)
3038 pmax_4 = ptr_p_4(kset, jset)
3039 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3040 CASE (6)
3041 pmax_1 = ptr_p_1(kset, iset)
3042 pmax_2 = ptr_p_2(jset, lset)
3043 pmax_3 = ptr_p_3(iset, lset)
3044 pmax_4 = ptr_p_4(kset, jset)
3045 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3046 CASE (7)
3047 pmax_1 = ptr_p_1(iset, kset)
3048 pmax_2 = ptr_p_2(jset, lset)
3049 pmax_3 = ptr_p_3(iset, lset)
3050 pmax_4 = ptr_p_4(kset, jset)
3051 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3052 CASE (8)
3053 pmax_1 = ptr_p_1(kset, iset)
3054 pmax_2 = ptr_p_2(lset, jset)
3055 pmax_3 = ptr_p_3(lset, iset)
3056 pmax_4 = ptr_p_4(jset, kset)
3057 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3058 CASE (9)
3059 pmax_1 = ptr_p_1(iset, kset)
3060 pmax_2 = ptr_p_2(lset, jset)
3061 pmax_3 = ptr_p_3(lset, iset)
3062 pmax_4 = ptr_p_4(jset, kset)
3063 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3064 CASE (10)
3065 pmax_1 = ptr_p_1(kset, iset)
3066 pmax_2 = ptr_p_2(jset, lset)
3067 pmax_3 = ptr_p_3(lset, iset)
3068 pmax_4 = ptr_p_4(jset, kset)
3069 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3070 CASE (11)
3071 pmax_1 = ptr_p_1(iset, kset)
3072 pmax_2 = ptr_p_2(jset, lset)
3073 pmax_3 = ptr_p_3(lset, iset)
3074 pmax_4 = ptr_p_4(jset, kset)
3075 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3076 CASE (12)
3077 pmax_1 = ptr_p_1(kset, iset)
3078 pmax_2 = ptr_p_2(lset, jset)
3079 pmax_3 = ptr_p_3(iset, lset)
3080 pmax_4 = ptr_p_4(jset, kset)
3081 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3082 CASE (13)
3083 pmax_1 = ptr_p_1(iset, kset)
3084 pmax_2 = ptr_p_2(lset, jset)
3085 pmax_3 = ptr_p_3(iset, lset)
3086 pmax_4 = ptr_p_4(jset, kset)
3087 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3088 CASE (14)
3089 pmax_1 = ptr_p_1(kset, iset)
3090 pmax_2 = ptr_p_2(jset, lset)
3091 pmax_3 = ptr_p_3(iset, lset)
3092 pmax_4 = ptr_p_4(jset, kset)
3093 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3094 CASE (15)
3095 pmax_1 = ptr_p_1(iset, kset)
3096 pmax_2 = ptr_p_2(jset, lset)
3097 pmax_3 = ptr_p_3(iset, lset)
3098 pmax_4 = ptr_p_4(jset, kset)
3099 pmax_val = max(pmax_1, pmax_2, pmax_3, pmax_4)
3100 CASE (16)
3101 pmax_1 = ptr_p_1(kset, iset)
3102 pmax_2 = ptr_p_2(lset, jset)
3103 pmax_3 = ptr_p_3(lset, iset)
3104 pmax_4 = ptr_p_4(kset, jset)
3105 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3106 CASE (17)
3107 pmax_1 = ptr_p_1(iset, kset)
3108 pmax_2 = ptr_p_2(lset, jset)
3109 pmax_3 = ptr_p_3(lset, iset)
3110 pmax_4 = ptr_p_4(kset, jset)
3111 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3112 CASE (18)
3113 pmax_1 = ptr_p_1(kset, iset)
3114 pmax_2 = ptr_p_2(jset, lset)
3115 pmax_3 = ptr_p_3(lset, iset)
3116 pmax_4 = ptr_p_4(kset, jset)
3117 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3118 CASE (19)
3119 pmax_1 = ptr_p_1(iset, kset)
3120 pmax_2 = ptr_p_2(jset, lset)
3121 pmax_3 = ptr_p_3(lset, iset)
3122 pmax_4 = ptr_p_4(kset, jset)
3123 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3124 CASE (20)
3125 pmax_1 = ptr_p_1(kset, iset)
3126 pmax_2 = ptr_p_2(lset, jset)
3127 pmax_3 = ptr_p_3(iset, lset)
3128 pmax_4 = ptr_p_4(kset, jset)
3129 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3130 CASE (21)
3131 pmax_1 = ptr_p_1(iset, kset)
3132 pmax_2 = ptr_p_2(lset, jset)
3133 pmax_3 = ptr_p_3(iset, lset)
3134 pmax_4 = ptr_p_4(kset, jset)
3135 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3136 CASE (22)
3137 pmax_1 = ptr_p_1(kset, iset)
3138 pmax_2 = ptr_p_2(jset, lset)
3139 pmax_3 = ptr_p_3(iset, lset)
3140 pmax_4 = ptr_p_4(kset, jset)
3141 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3142 CASE (23)
3143 pmax_1 = ptr_p_1(iset, kset)
3144 pmax_2 = ptr_p_2(jset, lset)
3145 pmax_3 = ptr_p_3(iset, lset)
3146 pmax_4 = ptr_p_4(kset, jset)
3147 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3148 CASE (24)
3149 pmax_1 = ptr_p_1(kset, iset)
3150 pmax_2 = ptr_p_2(lset, jset)
3151 pmax_3 = ptr_p_3(lset, iset)
3152 pmax_4 = ptr_p_4(jset, kset)
3153 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3154 CASE (25)
3155 pmax_1 = ptr_p_1(iset, kset)
3156 pmax_2 = ptr_p_2(lset, jset)
3157 pmax_3 = ptr_p_3(lset, iset)
3158 pmax_4 = ptr_p_4(jset, kset)
3159 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3160 CASE (26)
3161 pmax_1 = ptr_p_1(kset, iset)
3162 pmax_2 = ptr_p_2(jset, lset)
3163 pmax_3 = ptr_p_3(lset, iset)
3164 pmax_4 = ptr_p_4(jset, kset)
3165 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3166 CASE (27)
3167 pmax_1 = ptr_p_1(iset, kset)
3168 pmax_2 = ptr_p_2(jset, lset)
3169 pmax_3 = ptr_p_3(lset, iset)
3170 pmax_4 = ptr_p_4(jset, kset)
3171 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3172 CASE (28)
3173 pmax_1 = ptr_p_1(kset, iset)
3174 pmax_2 = ptr_p_2(lset, jset)
3175 pmax_3 = ptr_p_3(iset, lset)
3176 pmax_4 = ptr_p_4(jset, kset)
3177 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3178 CASE (29)
3179 pmax_1 = ptr_p_1(iset, kset)
3180 pmax_2 = ptr_p_2(lset, jset)
3181 pmax_3 = ptr_p_3(iset, lset)
3182 pmax_4 = ptr_p_4(jset, kset)
3183 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3184 CASE (30)
3185 pmax_1 = ptr_p_1(kset, iset)
3186 pmax_2 = ptr_p_2(jset, lset)
3187 pmax_3 = ptr_p_3(iset, lset)
3188 pmax_4 = ptr_p_4(jset, kset)
3189 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3190 CASE (31)
3191 pmax_1 = ptr_p_1(iset, kset)
3192 pmax_2 = ptr_p_2(jset, lset)
3193 pmax_3 = ptr_p_3(iset, lset)
3194 pmax_4 = ptr_p_4(jset, kset)
3195 pmax_val = max(pmax_1 + pmax_2, pmax_3 + pmax_4)
3196 END SELECT
3197
3198END SUBROUTINE get_pmax_val
3199
3200
3201END 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:120
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:21
bibliography::guidon2008
integer, save, public guidon2008
Definition bibliography.F:36
bibliography::guidon2009
integer, save, public guidon2009
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:847
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_dot_threadsafe
subroutine, public dbcsr_dot_threadsafe(matrix_a, matrix_b, result)
...
Definition cp_dbcsr_api.F:1302
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:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
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:2262
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:2331
hfx_load_balance_methods::hfx_load_balance
subroutine, public hfx_load_balance(x_data, eps_schwarz, particle_set, max_set, para_env, coeffs_set, coeffs_kind, is_assoc_atomic_block_global, do_periodic, load_balance_parameter, kind_of, basis_parameter, pmax_set, pmax_atom, i_thread, n_threads, cell, do_p_screening, map_atom_to_kind_atom, nkind, eval_type, pmax_block, use_virial)
Distributes the computation of eri's to all available processes.
Definition hfx_load_balance_methods.F:130
hfx_load_balance_methods::hfx_update_load_balance
subroutine, public hfx_update_load_balance(x_data, para_env, load_balance_parameter, i_thread, n_threads, eval_type)
Cheap way of redistributing the eri's.
Definition hfx_load_balance_methods.F:1483
hfx_pair_list_methods
Routines for optimizing load balance between processes in HFX calculations.
Definition hfx_pair_list_methods.F:15
hfx_pair_list_methods::build_pgf_product_list
subroutine, public build_pgf_product_list(list1, list2, product_list, nproducts, log10_pmax, log10_eps_schwarz, neighbor_cells, cell, potential_parameter, m_max, do_periodic)
...
Definition hfx_pair_list_methods.F: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:2522
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:2905
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:2873
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:2608
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:829
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:807
input_constants::do_potential_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
Definition input_constants.F:807
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:365
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:440
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
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:147
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_pp, 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, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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_pp, 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:333
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:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
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:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
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:217
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:136