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