Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate derivatives with respect to basis function origin
10 : !> \par History
11 : !> 04.2008 created [Manuel Guidon]
12 : !> \author Manuel Guidon
13 : ! **************************************************************************************************
14 : MODULE hfx_derivatives
15 : USE admm_types, ONLY: admm_type
16 : USE atomic_kind_types, ONLY: atomic_kind_type, &
17 : get_atomic_kind_set
18 : USE cell_types, ONLY: cell_type, &
19 : pbc
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_files, ONLY: close_file, &
22 : open_file
23 : USE cp_log_handling, ONLY: cp_get_default_logger, &
24 : cp_logger_type
25 : USE cp_output_handling, ONLY: cp_p_file, &
26 : cp_print_key_finished_output, &
27 : cp_print_key_should_output, &
28 : cp_print_key_unit_nr
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE cp_dbcsr_api, ONLY: dbcsr_get_matrix_type, &
31 : dbcsr_p_type, &
32 : dbcsr_type_antisymmetric
33 : USE four_center_integral_interface, ONLY: four_center_integral_context_type
34 : USE four_center_integral_methods, ONLY: forces4, &
35 : pgf_product_list_size
36 : USE gamma, ONLY: init_md_ftable
37 : USE hfx_communication, ONLY: get_full_density
38 : USE hfx_compression_methods, ONLY: hfx_add_mult_cache_elements, &
39 : hfx_add_single_cache_element, &
40 : hfx_decompress_first_cache, &
41 : hfx_flush_last_cache, &
42 : hfx_get_mult_cache_elements, &
43 : hfx_get_single_cache_element, &
44 : hfx_reset_cache_and_container
45 : USE hfx_load_balance_methods, ONLY: collect_load_balance_info, &
46 : hfx_load_balance, &
47 : hfx_update_load_balance
48 : USE hfx_pair_list_methods, ONLY: build_atomic_pair_list, &
49 : build_pair_list
50 : USE hfx_screening_methods, ONLY: update_pmax_mat
51 : USE hfx_types, ONLY: &
52 : alloc_containers, dealloc_containers, hfx_basis_info_type, hfx_basis_type, hfx_cache_type, &
53 : hfx_cell_type, hfx_container_type, hfx_distribution, hfx_general_type, hfx_init_container, &
54 : hfx_load_balance_type, hfx_memory_type, hfx_p_kind, hfx_pgf_list, hfx_pgf_product_list, &
55 : hfx_potential_type, hfx_screen_coeff_type, hfx_screening_type, hfx_task_list_type, &
56 : hfx_type, init_t_c_g0_lmax, log_zero, pair_list_type, pair_set_list_type
57 : USE input_constants, ONLY: do_potential_mix_cl_trunc, &
58 : do_potential_truncated, &
59 : hfx_do_eval_forces
60 : USE input_section_types, ONLY: section_vals_type
61 : USE kinds, ONLY: dp, &
62 : int_8
63 : USE machine, ONLY: m_walltime
64 : USE mathconstants, ONLY: fac
65 : USE orbital_pointers, ONLY: nco, &
66 : ncoset, &
67 : nso
68 : USE particle_types, ONLY: particle_type
69 : USE qs_environment_types, ONLY: get_qs_env, &
70 : qs_environment_type
71 : USE qs_force_types, ONLY: qs_force_type
72 : USE t_c_g0, ONLY: init
73 : USE util, ONLY: sort
74 : USE virial_types, ONLY: virial_type
75 :
76 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
77 :
78 : #include "./base/base_uses.f90"
79 :
80 : IMPLICIT NONE
81 : PRIVATE
82 :
83 : PUBLIC :: derivatives_four_center
84 :
85 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_derivatives'
86 :
87 : !***
88 :
89 : CONTAINS
90 :
91 : ! **************************************************************************************************
92 : !> \brief computes four center derivatives for a full basis set and updates the
93 : !> forces%fock_4c arrays. Uses all 8 eri symmetries
94 : !> \param qs_env ...
95 : !> \param rho_ao density matrix
96 : !> \param rho_ao_resp relaxed density matrix from response
97 : !> \param hfx_section HFX input section
98 : !> \param para_env para_env
99 : !> \param irep ID of HFX replica
100 : !> \param use_virial ...
101 : !> \param adiabatic_rescale_factor parameter used for MCY3 hybrid
102 : !> \param resp_only Calculate only forces from response density
103 : !> \par History
104 : !> 06.2007 created [Manuel Guidon]
105 : !> 08.2007 optimized load balance [Manuel Guidon]
106 : !> 02.2009 completely rewritten screening part [Manuel Guidon]
107 : !> 12.2017 major bug fix. removed wrong cycle that was causing segfault.
108 : !> see https://groups.google.com/forum/#!topic/cp2k/pc6B14XOALY
109 : !> [Tobias Binninger + Valery Weber]
110 : !> \author Manuel Guidon
111 : ! **************************************************************************************************
112 2012 : SUBROUTINE derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, &
113 : irep, use_virial, adiabatic_rescale_factor, resp_only, &
114 : external_x_data, nspins)
115 :
116 : TYPE(qs_environment_type), POINTER :: qs_env
117 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
118 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_resp
119 : TYPE(section_vals_type), POINTER :: hfx_section
120 : TYPE(mp_para_env_type), POINTER :: para_env
121 : INTEGER, INTENT(IN) :: irep
122 : LOGICAL, INTENT(IN) :: use_virial
123 : REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
124 : LOGICAL, INTENT(IN), OPTIONAL :: resp_only
125 : TYPE(hfx_type), DIMENSION(:, :), POINTER, OPTIONAL :: external_x_data
126 : INTEGER, INTENT(IN), OPTIONAL :: nspins
127 :
128 : CHARACTER(LEN=*), PARAMETER :: routineN = 'derivatives_four_center'
129 :
130 : INTEGER :: atomic_offset_ac, atomic_offset_ad, atomic_offset_bc, atomic_offset_bd, bin, &
131 : bits_max_val, buffer_left, buffer_size, buffer_start, cache_size, coord, current_counter, &
132 : forces_map(4, 2), handle, handle_bin, handle_getP, handle_load, handle_main, i, i_atom, &
133 : i_list_ij, i_list_kl, i_set_list_ij, i_set_list_ij_start, i_set_list_ij_stop, &
134 : i_set_list_kl, i_set_list_kl_start, i_set_list_kl_stop, i_thread, iatom, iatom_block, &
135 : iatom_end, iatom_start, ikind, iset, iw, j, j_atom, jatom, jatom_block, jatom_end, &
136 : jatom_start, jkind, jset, k, k_atom, katom, katom_block, katom_end, katom_start
137 : INTEGER :: kind_kind_idx, kkind, kset, l_atom, l_max, latom, latom_block, latom_end, &
138 : latom_start, lkind, lset, max_am, max_pgf, max_set, my_bin_id, my_bin_size, my_thread_id, &
139 : n_threads, natom, nbits, ncob, ncos_max, nints, nkimages, nkind, nneighbors, nseta, &
140 : nsetb, nsgf_max, my_nspins, sgfb, shm_task_counter, shm_total_bins, sphi_a_u1, sphi_a_u2, &
141 : sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3, sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, &
142 : sphi_d_u2, sphi_d_u3, swap_id, tmp_i4, unit_id
143 : INTEGER(int_8) :: atom_block, counter, estimate_to_store_int, max_val_memory, &
144 : mem_compression_counter, mem_eris, mem_max_val, my_current_counter, my_istart, &
145 : n_processes, nblocks, ncpu, neris_incore, neris_onthefly, neris_tmp, neris_total, &
146 : nprim_ints, shm_neris_incore, shm_neris_onthefly, shm_neris_total, shm_nprim_ints, &
147 : shm_stor_count_max_val, shm_storage_counter_integrals, stor_count_max_val, &
148 : storage_counter_integrals, tmp_block, tmp_i8(6)
149 2012 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:) :: tmp_task_list_cost
150 2012 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, last_sgf_global, &
151 2012 : nimages, tmp_index
152 2012 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, lc_min, ld_max, &
153 4024 : ld_min, npgfa, npgfb, npgfc, npgfd, nsgfa, nsgfb, nsgfc, nsgfd, shm_block_offset
154 4024 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfb, nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d, &
155 2012 : offset_ac_set, offset_ad_set, offset_bc_set, offset_bd_set, shm_atomic_block_offset
156 : INTEGER, DIMENSION(:, :), POINTER, SAVE :: shm_is_assoc_atomic_block
157 2012 : INTEGER, DIMENSION(:, :, :, :), POINTER :: shm_set_offset
158 : INTEGER, SAVE :: shm_number_of_p_entries
159 : LOGICAL :: bins_left, buffer_overflow, do_dynamic_load_balancing, do_it, do_periodic, &
160 : do_print_load_balance_info, is_anti_symmetric, my_resp_only, screen_pmat_forces, &
161 : treat_forces_in_core, use_disk_storage, with_resp_density
162 2012 : LOGICAL, DIMENSION(:, :), POINTER :: shm_atomic_pair_list
163 : REAL(dp) :: bintime_start, bintime_stop, cartesian_estimate, cartesian_estimate_virial, &
164 : compression_factor, eps_schwarz, eps_storage, fac, hf_fraction, ln_10, log10_eps_schwarz, &
165 : log10_pmax, log_2, max_contraction_val, max_val1, max_val2, max_val2_set, &
166 : my_adiabatic_rescale_factor, pmax_atom, pmax_blocks, pmax_entry, pmax_tmp, ra(3), rab2, &
167 : rb(3), rc(3), rcd2, rd(3), spherical_estimate, spherical_estimate_virial, symm_fac, &
168 : tmp_virial(3, 3)
169 2012 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: ede_buffer1, ede_buffer2, ede_primitives_tmp, &
170 2012 : ede_primitives_tmp_virial, ede_work, ede_work2, ede_work2_virial, ede_work_forces, &
171 4024 : ede_work_virial, pac_buf, pac_buf_beta, pac_buf_resp, pac_buf_resp_beta, pad_buf, &
172 2012 : pad_buf_beta, pad_buf_resp, pad_buf_resp_beta, pbc_buf, pbc_buf_beta, pbc_buf_resp, &
173 2012 : pbc_buf_resp_beta, pbd_buf, pbd_buf_beta, pbd_buf_resp, pbd_buf_resp_beta
174 2012 : REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: primitive_forces, primitive_forces_virial
175 2012 : REAL(dp), DIMENSION(:), POINTER :: full_density_resp, &
176 4024 : full_density_resp_beta, T2
177 2012 : REAL(dp), DIMENSION(:, :), POINTER :: full_density_alpha, full_density_beta, &
178 2012 : max_contraction, ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, shm_pmax_atom, shm_pmax_block, &
179 2012 : sphi_b, zeta, zetb, zetc, zetd
180 2012 : REAL(dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, &
181 2012 : sphi_c_ext_set, sphi_d_ext_set
182 2012 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, &
183 2012 : sphi_d_ext
184 : REAL(KIND=dp) :: coeffs_kind_max0
185 : TYPE(admm_type), POINTER :: admm_env
186 2012 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
187 : TYPE(cell_type), POINTER :: cell
188 : TYPE(four_center_integral_context_type) :: private_context
189 : TYPE(cp_logger_type), POINTER :: logger
190 : TYPE(dft_control_type), POINTER :: dft_control
191 : TYPE(hfx_basis_info_type), POINTER :: basis_info
192 2012 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
193 2012 : TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
194 : TYPE(hfx_cache_type), POINTER :: maxval_cache
195 2012 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
196 : TYPE(hfx_container_type), POINTER :: maxval_container
197 : TYPE(hfx_distribution), POINTER :: distribution_forces
198 : TYPE(hfx_general_type) :: general_parameter
199 : TYPE(hfx_load_balance_type), POINTER :: load_balance_parameter
200 : TYPE(hfx_memory_type), POINTER :: memory_parameter
201 2012 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: shm_initial_p
202 2012 : TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl
203 : TYPE(hfx_pgf_product_list), ALLOCATABLE, &
204 2012 : DIMENSION(:) :: pgf_product_list
205 : TYPE(hfx_potential_type) :: potential_parameter
206 : TYPE(hfx_screen_coeff_type), DIMENSION(:, :), &
207 2012 : POINTER :: screen_coeffs_kind, tmp_R_1, tmp_R_2, &
208 2012 : tmp_screen_pgf1, tmp_screen_pgf2
209 : TYPE(hfx_screen_coeff_type), &
210 2012 : DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set
211 : TYPE(hfx_screen_coeff_type), &
212 2012 : DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf
213 : TYPE(hfx_screening_type) :: screening_parameter
214 2012 : TYPE(hfx_task_list_type), DIMENSION(:), POINTER :: shm_task_list, tmp_task_list
215 : TYPE(hfx_type), POINTER :: actual_x_data
216 2012 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: my_x_data
217 : TYPE(pair_list_type) :: list_ij, list_kl
218 : TYPE(pair_set_list_type), ALLOCATABLE, &
219 2012 : DIMENSION(:) :: set_list_ij, set_list_kl
220 2012 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
221 2012 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
222 : TYPE(virial_type), POINTER :: virial
223 :
224 2012 : NULLIFY (dft_control, admm_env)
225 :
226 2012 : with_resp_density = .FALSE.
227 684 : IF (ASSOCIATED(rho_ao_resp)) with_resp_density = .TRUE.
228 :
229 2012 : my_resp_only = .FALSE.
230 2012 : IF (PRESENT(resp_only)) my_resp_only = resp_only
231 :
232 2012 : is_anti_symmetric = dbcsr_get_matrix_type(rho_ao(1, 1)%matrix) == dbcsr_type_antisymmetric
233 :
234 : CALL get_qs_env(qs_env, &
235 : admm_env=admm_env, &
236 : particle_set=particle_set, &
237 : atomic_kind_set=atomic_kind_set, &
238 : cell=cell, &
239 2012 : dft_control=dft_control)
240 :
241 2012 : IF (PRESENT(nspins)) THEN
242 252 : my_nspins = nspins
243 : ELSE
244 1760 : my_nspins = dft_control%nspins
245 : END IF
246 2012 : nkimages = dft_control%nimages
247 2012 : CPASSERT(nkimages == 1)
248 :
249 2012 : my_x_data => qs_env%x_data
250 2012 : IF (PRESENT(external_x_data)) my_x_data => external_x_data
251 :
252 : !! One atom systems have no contribution to forces
253 2012 : IF (SIZE(particle_set, 1) == 1) THEN
254 4 : IF (.NOT. use_virial) RETURN
255 : END IF
256 :
257 2008 : CALL timeset(routineN, handle)
258 : !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe
259 2008 : nkind = SIZE(atomic_kind_set, 1)
260 2008 : l_max = 0
261 5786 : DO ikind = 1, nkind
262 16122 : l_max = MAX(l_max, MAXVAL(my_x_data(1, 1)%basis_parameter(ikind)%lmax))
263 : END DO
264 2008 : l_max = 4*l_max + 1
265 2008 : CALL init_md_ftable(l_max)
266 :
267 2008 : IF (my_x_data(1, 1)%potential_parameter%potential_type == do_potential_truncated .OR. &
268 : my_x_data(1, 1)%potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
269 752 : IF (l_max > init_t_c_g0_lmax) THEN
270 208 : CALL open_file(unit_number=unit_id, file_name=my_x_data(1, 1)%potential_parameter%filename)
271 208 : CALL init(l_max, unit_id, para_env%mepos, para_env)
272 208 : CALL close_file(unit_id)
273 208 : init_t_c_g0_lmax = l_max
274 : END IF
275 : END IF
276 :
277 2008 : n_threads = 1
278 2008 : !$ n_threads = omp_get_max_threads()
279 :
280 2008 : shm_neris_total = 0
281 2008 : shm_nprim_ints = 0
282 2008 : shm_neris_onthefly = 0
283 2008 : shm_storage_counter_integrals = 0
284 2008 : shm_neris_incore = 0
285 2008 : shm_stor_count_max_val = 0
286 :
287 : !! get force array
288 2008 : CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
289 :
290 2008 : my_adiabatic_rescale_factor = 1.0_dp
291 2008 : IF (PRESENT(adiabatic_rescale_factor)) THEN
292 260 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
293 : END IF
294 :
295 : !$OMP PARALLEL DEFAULT(NONE) &
296 : !$OMP SHARED(qs_env,&
297 : !$OMP rho_ao,&
298 : !$OMP rho_ao_resp,&
299 : !$OMP with_resp_density,&
300 : !$OMP hfx_section,&
301 : !$OMP para_env,&
302 : !$OMP irep,&
303 : !$OMP my_resp_only,&
304 : !$OMP ncoset,&
305 : !$OMP nco,&
306 : !$OMP nso,&
307 : !$OMP n_threads,&
308 : !$OMP full_density_alpha,&
309 : !$OMP full_density_resp,&
310 : !$OMP full_density_resp_beta,&
311 : !$OMP full_density_beta,&
312 : !$OMP shm_initial_p,&
313 : !$OMP shm_is_assoc_atomic_block,&
314 : !$OMP shm_number_of_p_entries,&
315 : !$OMP force,&
316 : !$OMP virial,&
317 : !$OMP my_adiabatic_rescale_factor,&
318 : !$OMP shm_neris_total,&
319 : !$OMP shm_neris_onthefly,&
320 : !$OMP shm_storage_counter_integrals,&
321 : !$OMP shm_neris_incore,&
322 : !$OMP shm_nprim_ints,&
323 : !$OMP shm_stor_count_max_val,&
324 : !$OMP cell,&
325 : !$OMP screen_coeffs_set,&
326 : !$OMP screen_coeffs_kind,&
327 : !$OMP screen_coeffs_pgf,&
328 : !$OMP pgf_product_list_size,&
329 : !$OMP nkind,&
330 : !$OMP is_anti_symmetric,&
331 : !$OMP radii_pgf,&
332 : !$OMP shm_atomic_block_offset,&
333 : !$OMP shm_set_offset,&
334 : !$OMP shm_block_offset,&
335 : !$OMP shm_task_counter,&
336 : !$OMP shm_task_list,&
337 : !$OMP shm_total_bins,&
338 : !$OMP shm_pmax_block,&
339 : !$OMP shm_atomic_pair_list,&
340 : !$OMP shm_pmax_atom,&
341 : !$OMP use_virial,&
342 : !$OMP do_print_load_balance_info,&
343 : !$OMP nkimages,my_nspins,my_x_data)&
344 : !$OMP PRIVATE(actual_x_data,atomic_kind_set,atomic_offset_ac,atomic_offset_ad,atomic_offset_bc,&
345 : !$OMP atomic_offset_bd,atom_of_kind,basis_info,&
346 : !$OMP basis_parameter,bin,bins_left,bintime_start,bintime_stop,bits_max_val,buffer_left,buffer_overflow,buffer_size,&
347 : !$OMP buffer_start,cache_size,cartesian_estimate,cartesian_estimate_virial,&
348 : !$OMP coeffs_kind_max0,compression_factor,counter,current_counter,&
349 : !$OMP distribution_forces,do_dynamic_load_balancing,do_it,do_periodic,ede_buffer1,ede_buffer2,&
350 : !$OMP ede_primitives_tmp,ede_primitives_tmp_virial,&
351 : !$OMP ede_work,ede_work2,ede_work2_virial,ede_work_forces,ede_work_virial,eps_schwarz,eps_storage,estimate_to_store_int,&
352 : !$OMP fac,first_sgfb,forces_map,general_parameter,handle_bin,handle_getp,handle_load,&
353 : !$OMP handle_main,hf_fraction,i_atom,iatom_end,iatom_start,ikind,integral_caches,integral_containers,&
354 : !$OMP i_set_list_ij_start,i_set_list_ij_stop,i_set_list_kl_start,i_set_list_kl_stop,i_thread,iw,j_atom,jatom_end,&
355 : !$OMP jatom_start,jkind,katom,k_atom,katom_block,katom_end,katom_start,kind_kind_idx,&
356 : !$OMP kind_of,kkind,kset,la_max,la_min,last_sgf_global,latom,l_atom,&
357 : !$OMP latom_block,latom_end,latom_start,lb_max,lb_min,lc_max,lc_min,ld_max,&
358 : !$OMP ld_min,list_ij,list_kl,lkind,ln_10,load_balance_parameter,log10_eps_schwarz,log10_pmax,&
359 : !$OMP log_2,logger,lset,max_am,max_contraction,max_contraction_val,max_pgf,max_set,&
360 : !$OMP max_val1,max_val2,max_val2_set,maxval_cache,maxval_container,max_val_memory,mem_compression_counter,mem_eris,&
361 : !$OMP mem_max_val,memory_parameter,my_bin_id,my_bin_size,my_current_counter,my_istart,my_thread_id,natom,&
362 : !$OMP nbits,nblocks,ncob,ncos_max,ncpu,neris_incore,neris_onthefly,neris_tmp,&
363 : !$OMP neris_total,nimages,nints,nneighbors,npgfa,npgfb,npgfc,npgfd,&
364 : !$OMP nprim_ints,n_processes,nseta,nsetb,nsgfa,nsgfb,nsgfc,nsgfd,&
365 : !$OMP nsgfl_a,nsgfl_b,nsgfl_c,nsgfl_d,nsgf_max,offset_ac_set,offset_ad_set,&
366 : !$OMP offset_bc_set,offset_bd_set,pac_buf,pac_buf_beta,pac_buf_resp,pac_buf_resp_beta,pad_buf,pad_buf_beta,pad_buf_resp,&
367 : !$OMP pad_buf_resp_beta,particle_set,pbc_buf,pbc_buf_beta,pbc_buf_resp,pbc_buf_resp_beta,pbd_buf,pbd_buf_beta, &
368 : !$OMP pbd_buf_resp, pbd_buf_resp_beta, pgf_list_ij,&
369 : !$OMP pgf_list_kl,pgf_product_list,pmax_atom,pmax_blocks,pmax_entry,pmax_tmp,potential_parameter,primitive_forces,&
370 : !$OMP primitive_forces_virial,private_context,ptr_p_1,ptr_p_2,ptr_p_3,ptr_p_4,ra,rab2,&
371 : !$OMP rb,rc,rcd2,rd,screening_parameter,screen_pmat_forces,set_list_ij,set_list_kl,&
372 : !$OMP sgfb,spherical_estimate,spherical_estimate_virial,sphi_a_ext,sphi_a_ext_set,sphi_a_u1,sphi_a_u2,sphi_a_u3,&
373 : !$OMP sphi_b,sphi_b_ext,sphi_b_ext_set,sphi_b_u1,sphi_b_u2,sphi_b_u3,sphi_c_ext,sphi_c_ext_set,&
374 : !$OMP sphi_c_u1,sphi_c_u2,sphi_c_u3,sphi_d_ext,sphi_d_ext_set,sphi_d_u1,sphi_d_u2,sphi_d_u3,&
375 : !$OMP storage_counter_integrals,stor_count_max_val,swap_id,symm_fac,t2,tmp_block,tmp_i8, tmp_i4,&
376 : !$OMP tmp_index,tmp_r_1,tmp_r_2,tmp_screen_pgf1,tmp_screen_pgf2,tmp_task_list,tmp_task_list_cost,tmp_virial,&
377 2008 : !$OMP treat_forces_in_core,use_disk_storage,zeta,zetb,zetc,zetd)
378 :
379 : i_thread = 0
380 : !$ i_thread = omp_get_thread_num()
381 :
382 : ln_10 = LOG(10.0_dp)
383 : log_2 = LOG10(2.0_dp)
384 : actual_x_data => my_x_data(irep, i_thread + 1)
385 : do_periodic = actual_x_data%periodic_parameter%do_periodic
386 :
387 : screening_parameter = actual_x_data%screening_parameter
388 : general_parameter = actual_x_data%general_parameter
389 : potential_parameter = actual_x_data%potential_parameter
390 : basis_info => actual_x_data%basis_info
391 :
392 : load_balance_parameter => actual_x_data%load_balance_parameter
393 : basis_parameter => actual_x_data%basis_parameter
394 :
395 : ! IF(with_resp_density) THEN
396 : ! ! here we also do a copy of the original load_balance_parameter
397 : ! ! since if MP2 forces are required after the calculation of the HFX
398 : ! ! derivatives we need to calculate again the SCF energy.
399 : ! ALLOCATE(load_balance_parameter_energy)
400 : ! load_balance_parameter_energy = actual_x_data%load_balance_parameter
401 : ! END IF
402 :
403 : memory_parameter => actual_x_data%memory_parameter
404 :
405 : cache_size = memory_parameter%cache_size
406 : bits_max_val = memory_parameter%bits_max_val
407 :
408 : use_disk_storage = .FALSE.
409 :
410 : CALL get_qs_env(qs_env=qs_env, &
411 : atomic_kind_set=atomic_kind_set, &
412 : particle_set=particle_set)
413 :
414 : max_set = basis_info%max_set
415 : max_am = basis_info%max_am
416 : natom = SIZE(particle_set, 1)
417 :
418 : treat_forces_in_core = memory_parameter%treat_forces_in_core
419 :
420 : hf_fraction = general_parameter%fraction
421 : hf_fraction = hf_fraction*my_adiabatic_rescale_factor
422 : eps_schwarz = screening_parameter%eps_schwarz_forces
423 : IF (eps_schwarz < 0.0_dp) THEN
424 : log10_eps_schwarz = log_zero
425 : ELSE
426 : !! ** make eps_schwarz tighter in case of a stress calculation
427 : IF (use_virial) eps_schwarz = eps_schwarz/actual_x_data%periodic_parameter%R_max_stress
428 : log10_eps_schwarz = LOG10(eps_schwarz)
429 : END IF
430 : screen_pmat_forces = screening_parameter%do_p_screening_forces
431 : ! don't do density screening in case of response forces
432 : IF (with_resp_density) screen_pmat_forces = .FALSE.
433 :
434 : eps_storage = eps_schwarz*memory_parameter%eps_storage_scaling
435 :
436 : counter = 0_int_8
437 :
438 : !! Copy the libint_guy
439 : private_context = actual_x_data%eri4_context
440 :
441 : !! Get screening parameter
442 :
443 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
444 : atom_of_kind=atom_of_kind, &
445 : kind_of=kind_of)
446 :
447 : !! Create helper arrray for mapping local basis functions to global ones
448 : ALLOCATE (last_sgf_global(0:natom))
449 : last_sgf_global(0) = 0
450 : DO iatom = 1, natom
451 : ikind = kind_of(iatom)
452 : last_sgf_global(iatom) = last_sgf_global(iatom - 1) + basis_parameter(ikind)%nsgf_total
453 : END DO
454 :
455 : ALLOCATE (max_contraction(max_set, natom))
456 :
457 : max_contraction = 0.0_dp
458 : max_pgf = 0
459 : DO jatom = 1, natom
460 : jkind = kind_of(jatom)
461 : lb_max => basis_parameter(jkind)%lmax
462 : npgfb => basis_parameter(jkind)%npgf
463 : nsetb = basis_parameter(jkind)%nset
464 : first_sgfb => basis_parameter(jkind)%first_sgf
465 : sphi_b => basis_parameter(jkind)%sphi
466 : nsgfb => basis_parameter(jkind)%nsgf
467 : DO jset = 1, nsetb
468 : ! takes the primitive to contracted transformation into account
469 : ncob = npgfb(jset)*ncoset(lb_max(jset))
470 : sgfb = first_sgfb(1, jset)
471 : ! if the primitives are assumed to be all of max_val2, max_val2*p2s_b becomes
472 : ! the maximum value after multiplication with sphi_b
473 : max_contraction(jset, jatom) = MAXVAL([(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)])
474 : max_pgf = MAX(max_pgf, npgfb(jset))
475 : END DO
476 : END DO
477 :
478 : ncpu = para_env%num_pe
479 : n_processes = ncpu*n_threads
480 :
481 : !! initialize some counters
482 : neris_total = 0_int_8
483 : neris_incore = 0_int_8
484 : neris_onthefly = 0_int_8
485 : mem_eris = 0_int_8
486 : mem_max_val = 0_int_8
487 : compression_factor = 0.0_dp
488 : nprim_ints = 0_int_8
489 : neris_tmp = 0_int_8
490 : max_val_memory = 1_int_8
491 :
492 : !$OMP MASTER
493 : !! Set pointer for is_assoc helper
494 : shm_is_assoc_atomic_block => actual_x_data%is_assoc_atomic_block
495 : shm_number_of_p_entries = actual_x_data%number_of_p_entries
496 : shm_atomic_block_offset => actual_x_data%atomic_block_offset
497 : shm_set_offset => actual_x_data%set_offset
498 : shm_block_offset => actual_x_data%block_offset
499 :
500 : NULLIFY (full_density_alpha)
501 : NULLIFY (full_density_beta)
502 : ALLOCATE (full_density_alpha(shm_block_offset(ncpu + 1), nkimages))
503 :
504 : !! Get the full density from all the processors
505 : CALL timeset(routineN//"_getP", handle_getP)
506 : CALL get_full_density(para_env, full_density_alpha(:, 1), rho_ao(1, 1)%matrix, shm_number_of_p_entries, &
507 : shm_block_offset, &
508 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
509 : antisymmetric=is_anti_symmetric)
510 : ! for now only closed shell case
511 : IF (with_resp_density) THEN
512 : NULLIFY (full_density_resp)
513 : ALLOCATE (full_density_resp(shm_block_offset(ncpu + 1)))
514 : CALL get_full_density(para_env, full_density_resp, rho_ao_resp(1)%matrix, shm_number_of_p_entries, &
515 : shm_block_offset, &
516 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
517 : antisymmetric=is_anti_symmetric)
518 : END IF
519 :
520 : IF (my_nspins == 2) THEN
521 : ALLOCATE (full_density_beta(shm_block_offset(ncpu + 1), 1))
522 : CALL get_full_density(para_env, full_density_beta(:, 1), rho_ao(2, 1)%matrix, shm_number_of_p_entries, &
523 : shm_block_offset, &
524 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
525 : antisymmetric=is_anti_symmetric)
526 : ! With resp density
527 : IF (with_resp_density) THEN
528 : NULLIFY (full_density_resp_beta)
529 : ALLOCATE (full_density_resp_beta(shm_block_offset(ncpu + 1)))
530 : CALL get_full_density(para_env, full_density_resp_beta, rho_ao_resp(2)%matrix, shm_number_of_p_entries, &
531 : shm_block_offset, &
532 : kind_of, basis_parameter, get_max_vals_spin=.FALSE., &
533 : antisymmetric=is_anti_symmetric)
534 : END IF
535 :
536 : END IF
537 : CALL timestop(handle_getP)
538 :
539 : !! Calculate max entries for screening on actual density. If screen_p_mat_forces = FALSE, the
540 : !! matrix is initialized to 1.0
541 : IF (screen_pmat_forces) THEN
542 : NULLIFY (shm_initial_p)
543 : shm_initial_p => actual_x_data%initial_p_forces
544 : shm_pmax_atom => actual_x_data%pmax_atom_forces
545 : IF (memory_parameter%recalc_forces) THEN
546 : CALL update_pmax_mat(shm_initial_p, &
547 : actual_x_data%map_atom_to_kind_atom, &
548 : actual_x_data%set_offset, &
549 : actual_x_data%atomic_block_offset, &
550 : shm_pmax_atom, &
551 : full_density_alpha, full_density_beta, &
552 : natom, kind_of, basis_parameter, &
553 : nkind, actual_x_data%is_assoc_atomic_block, my_nspins)
554 : END IF
555 : END IF
556 :
557 : ! restore as full density the HF density
558 : ! maybe in the future
559 : IF (with_resp_density .AND. .NOT. my_resp_only) THEN
560 : full_density_alpha(:, 1) = full_density_alpha(:, 1) - full_density_resp
561 : IF (my_nspins == 2) THEN
562 : full_density_beta(:, 1) = &
563 : full_density_beta(:, 1) - full_density_resp_beta
564 : END IF
565 : ! full_density_resp=full_density+full_density_resp
566 : END IF
567 :
568 : screen_coeffs_set => actual_x_data%screen_funct_coeffs_set
569 : screen_coeffs_kind => actual_x_data%screen_funct_coeffs_kind
570 : screen_coeffs_pgf => actual_x_data%screen_funct_coeffs_pgf
571 : radii_pgf => actual_x_data%pair_dist_radii_pgf
572 : shm_pmax_block => actual_x_data%pmax_block
573 :
574 : !$OMP END MASTER
575 : !$OMP BARRIER
576 :
577 : !$OMP MASTER
578 : CALL timeset(routineN//"_load", handle_load)
579 : !$OMP END MASTER
580 : !! Load balance the work
581 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
582 : IF (actual_x_data%b_first_load_balance_forces) THEN
583 : CALL hfx_load_balance(actual_x_data, eps_schwarz, particle_set, max_set, para_env, &
584 : screen_coeffs_set, screen_coeffs_kind, &
585 : shm_is_assoc_atomic_block, do_periodic, &
586 : load_balance_parameter, kind_of, basis_parameter, shm_initial_p, &
587 : shm_pmax_atom, i_thread, n_threads, &
588 : cell, screen_pmat_forces, actual_x_data%map_atom_to_kind_atom, &
589 : nkind, hfx_do_eval_forces, shm_pmax_block, use_virial)
590 : actual_x_data%b_first_load_balance_forces = .FALSE.
591 : ELSE
592 : CALL hfx_update_load_balance(actual_x_data, para_env, &
593 : load_balance_parameter, &
594 : i_thread, n_threads, hfx_do_eval_forces)
595 : END IF
596 : END IF
597 : !$OMP MASTER
598 : CALL timestop(handle_load)
599 : !$OMP END MASTER
600 :
601 : !! precompute maximum nco and allocate scratch
602 : ncos_max = 0
603 : nsgf_max = 0
604 : DO iatom = 1, natom
605 : ikind = kind_of(iatom)
606 : nseta = basis_parameter(ikind)%nset
607 : npgfa => basis_parameter(ikind)%npgf
608 : la_max => basis_parameter(ikind)%lmax
609 : nsgfa => basis_parameter(ikind)%nsgf
610 : DO iset = 1, nseta
611 : ncos_max = MAX(ncos_max, ncoset(la_max(iset)))
612 : nsgf_max = MAX(nsgf_max, nsgfa(iset))
613 : END DO
614 : END DO
615 :
616 : !! Allocate work arrays
617 : ALLOCATE (primitive_forces(12*nsgf_max**4))
618 : primitive_forces = 0.0_dp
619 :
620 : ! ** Allocate buffers for pgf_lists
621 : nneighbors = SIZE(actual_x_data%neighbor_cells)
622 : ALLOCATE (pgf_list_ij(max_pgf**2))
623 : ALLOCATE (pgf_list_kl(max_pgf**2))
624 : !$OMP ATOMIC READ
625 : tmp_i4 = pgf_product_list_size
626 : ALLOCATE (pgf_product_list(tmp_i4))
627 : ALLOCATE (nimages(max_pgf**2))
628 :
629 : DO i = 1, max_pgf**2
630 : ALLOCATE (pgf_list_ij(i)%image_list(nneighbors))
631 : ALLOCATE (pgf_list_kl(i)%image_list(nneighbors))
632 : END DO
633 :
634 : ALLOCATE (pbd_buf(nsgf_max**2))
635 : ALLOCATE (pbc_buf(nsgf_max**2))
636 : ALLOCATE (pad_buf(nsgf_max**2))
637 : ALLOCATE (pac_buf(nsgf_max**2))
638 : IF (with_resp_density) THEN
639 : ALLOCATE (pbd_buf_resp(nsgf_max**2))
640 : ALLOCATE (pbc_buf_resp(nsgf_max**2))
641 : ALLOCATE (pad_buf_resp(nsgf_max**2))
642 : ALLOCATE (pac_buf_resp(nsgf_max**2))
643 : END IF
644 : IF (my_nspins == 2) THEN
645 : ALLOCATE (pbd_buf_beta(nsgf_max**2))
646 : ALLOCATE (pbc_buf_beta(nsgf_max**2))
647 : ALLOCATE (pad_buf_beta(nsgf_max**2))
648 : ALLOCATE (pac_buf_beta(nsgf_max**2))
649 : IF (with_resp_density) THEN
650 : ALLOCATE (pbd_buf_resp_beta(nsgf_max**2))
651 : ALLOCATE (pbc_buf_resp_beta(nsgf_max**2))
652 : ALLOCATE (pad_buf_resp_beta(nsgf_max**2))
653 : ALLOCATE (pac_buf_resp_beta(nsgf_max**2))
654 : END IF
655 : END IF
656 :
657 : ALLOCATE (ede_work(ncos_max**4*12))
658 : ALLOCATE (ede_work2(ncos_max**4*12))
659 : ALLOCATE (ede_work_forces(ncos_max**4*12))
660 : ALLOCATE (ede_buffer1(ncos_max**4))
661 : ALLOCATE (ede_buffer2(ncos_max**4))
662 : ALLOCATE (ede_primitives_tmp(nsgf_max**4))
663 :
664 : IF (use_virial) THEN
665 : ALLOCATE (primitive_forces_virial(12*nsgf_max**4*3))
666 : primitive_forces_virial = 0.0_dp
667 : ALLOCATE (ede_work_virial(ncos_max**4*12*3))
668 : ALLOCATE (ede_work2_virial(ncos_max**4*12*3))
669 : ALLOCATE (ede_primitives_tmp_virial(nsgf_max**4))
670 : tmp_virial = 0.0_dp
671 : ELSE
672 : ! dummy allocation
673 : ALLOCATE (primitive_forces_virial(1))
674 : primitive_forces_virial = 0.0_dp
675 : ALLOCATE (ede_work_virial(1))
676 : ALLOCATE (ede_work2_virial(1))
677 : ALLOCATE (ede_primitives_tmp_virial(1))
678 : END IF
679 :
680 : !! Start caluclating integrals of the form (ab|cd) or (ij|kl)
681 : !! In order to do so, there is a main four-loop structure that takes into account the two symmetries
682 : !!
683 : !! (ab|cd) = (ba|cd) = (ab|dc) = (ba|dc)
684 : !!
685 : !! by iterating in the following way
686 : !!
687 : !! DO iatom=1,natom and DO katom=1,natom
688 : !! DO jatom=iatom,natom DO latom=katom,natom
689 : !!
690 : !! The third symmetry
691 : !!
692 : !! (ab|cd) = (cd|ab)
693 : !!
694 : !! is taken into account by the following criterion:
695 : !!
696 : !! IF(katom+latom<=iatom+jatom) THEN
697 : !! IF( ((iatom+jatom)==(katom+latom) ) .AND.(katom<iatom)) CYCLE
698 : !!
699 : !! Furthermore, if iatom==jatom==katom==latom we cycle, because the derivatives are zero anyway.
700 : !!
701 : !! Depending on the degeneracy of an integral the exchange contribution is multiplied by a corresponding
702 : !! factor ( symm_fac ).
703 : !!
704 : !! If one quartet does not pass the screening we CYCLE on the outer most possible loop. Thats why we use
705 : !! different hierarchies of short range screening matrices.
706 : !!
707 : !! If we do a parallel run, each process owns a unique array of workloads. Here, a workload is
708 : !! defined as :
709 : !!
710 : !! istart, jstart, kstart, lstart, number_of_atom_quartets, initial_cpu_id
711 : !!
712 : !! This tells the process where to start the main loops and how many bunches of integrals it has to
713 : !! calculate. The original parallelization is a simple modulo distribution that is binned and
714 : !! optimized in the load_balance routines. Since the Monte Carlo routines can swap processors,
715 : !! we need to know which was the initial cpu_id.
716 : !! Furthermore, the indices iatom, jatom, katom, latom have to be set to istart, jstart, kstart and
717 : !! lstart only the first time the loop is executed. All subsequent loops have to start with one or
718 : !! iatom and katom respectively. Therefore, we use flags like first_j_loop etc.
719 :
720 : !$OMP BARRIER
721 : !$OMP MASTER
722 : CALL timeset(routineN//"_main", handle_main)
723 : !$OMP END MASTER
724 : !$OMP BARRIER
725 :
726 : coeffs_kind_max0 = MAXVAL(screen_coeffs_kind(:, :)%x(2))
727 : ALLOCATE (set_list_ij((max_set*load_balance_parameter%block_size)**2))
728 : ALLOCATE (set_list_kl((max_set*load_balance_parameter%block_size)**2))
729 :
730 : !$OMP BARRIER
731 :
732 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
733 : actual_x_data%distribution_forces%ram_counter = HUGE(distribution_forces%ram_counter)
734 : memory_parameter%ram_counter_forces = HUGE(memory_parameter%ram_counter_forces)
735 : END IF
736 :
737 : IF (treat_forces_in_core) THEN
738 : buffer_overflow = .FALSE.
739 : ELSE
740 : buffer_overflow = .TRUE.
741 : END IF
742 :
743 : do_dynamic_load_balancing = .TRUE.
744 : IF (n_threads == 1) do_dynamic_load_balancing = .FALSE.
745 :
746 : IF (do_dynamic_load_balancing) THEN
747 : my_bin_size = SIZE(actual_x_data%distribution_forces)
748 : ELSE
749 : my_bin_size = 1
750 : END IF
751 : !! We do not need the containers if MAX_MEM = 0
752 : IF (treat_forces_in_core) THEN
753 : !! IF new md step -> reinitialize containers
754 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
755 : CALL dealloc_containers(actual_x_data%store_forces, memory_parameter%actual_memory_usage)
756 : CALL alloc_containers(actual_x_data%store_forces, my_bin_size)
757 :
758 : DO bin = 1, my_bin_size
759 : maxval_container => actual_x_data%store_forces%maxval_container(bin)
760 : integral_containers => actual_x_data%store_forces%integral_containers(:, bin)
761 : CALL hfx_init_container(maxval_container, memory_parameter%actual_memory_usage, .FALSE.)
762 : DO i = 1, 64
763 : CALL hfx_init_container(integral_containers(i), memory_parameter%actual_memory_usage, .FALSE.)
764 : END DO
765 : END DO
766 : END IF
767 :
768 : !! Decompress the first cache for maxvals and integrals
769 : IF (.NOT. actual_x_data%memory_parameter%recalc_forces) THEN
770 : DO bin = 1, my_bin_size
771 : maxval_cache => actual_x_data%store_forces%maxval_cache(bin)
772 : maxval_container => actual_x_data%store_forces%maxval_container(bin)
773 : integral_caches => actual_x_data%store_forces%integral_caches(:, bin)
774 : integral_containers => actual_x_data%store_forces%integral_containers(:, bin)
775 : CALL hfx_decompress_first_cache(bits_max_val, maxval_cache, maxval_container, &
776 : memory_parameter%actual_memory_usage, .FALSE.)
777 : DO i = 1, 64
778 : CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), &
779 : memory_parameter%actual_memory_usage, .FALSE.)
780 : END DO
781 : END DO
782 : END IF
783 : END IF
784 :
785 : !$OMP BARRIER
786 : !$OMP MASTER
787 : IF (do_dynamic_load_balancing) THEN
788 : ! ** Lets construct the task list
789 : shm_total_bins = 0
790 : DO i = 1, n_threads
791 : shm_total_bins = shm_total_bins + SIZE(my_x_data(irep, i)%distribution_forces)
792 : END DO
793 : ALLOCATE (tmp_task_list(shm_total_bins))
794 : shm_task_counter = 0
795 : DO i = 1, n_threads
796 : DO bin = 1, SIZE(my_x_data(irep, i)%distribution_forces)
797 : shm_task_counter = shm_task_counter + 1
798 : tmp_task_list(shm_task_counter)%thread_id = i
799 : tmp_task_list(shm_task_counter)%bin_id = bin
800 : tmp_task_list(shm_task_counter)%cost = my_x_data(irep, i)%distribution_forces(bin)%cost
801 : END DO
802 : END DO
803 :
804 : ! ** Now sort the task list
805 : ALLOCATE (tmp_task_list_cost(shm_total_bins))
806 : ALLOCATE (tmp_index(shm_total_bins))
807 :
808 : DO i = 1, shm_total_bins
809 : tmp_task_list_cost(i) = tmp_task_list(i)%cost
810 : END DO
811 :
812 : CALL sort(tmp_task_list_cost, shm_total_bins, tmp_index)
813 :
814 : ALLOCATE (actual_x_data%task_list(shm_total_bins))
815 :
816 : DO i = 1, shm_total_bins
817 : actual_x_data%task_list(i) = tmp_task_list(tmp_index(shm_total_bins - i + 1))
818 : END DO
819 :
820 : shm_task_list => actual_x_data%task_list
821 : shm_task_counter = 0
822 :
823 : DEALLOCATE (tmp_task_list_cost, tmp_index, tmp_task_list)
824 : END IF
825 :
826 : !! precalculate maximum density matrix elements in blocks
827 : shm_pmax_block => actual_x_data%pmax_block
828 : actual_x_data%pmax_block = 0.0_dp
829 : IF (screen_pmat_forces) THEN
830 : DO iatom_block = 1, SIZE(actual_x_data%blocks)
831 : iatom_start = actual_x_data%blocks(iatom_block)%istart
832 : iatom_end = actual_x_data%blocks(iatom_block)%iend
833 : DO jatom_block = 1, SIZE(actual_x_data%blocks)
834 : jatom_start = actual_x_data%blocks(jatom_block)%istart
835 : jatom_end = actual_x_data%blocks(jatom_block)%iend
836 : shm_pmax_block(iatom_block, jatom_block) = MAXVAL(shm_pmax_atom(iatom_start:iatom_end, jatom_start:jatom_end))
837 : END DO
838 : END DO
839 : END IF
840 : shm_atomic_pair_list => actual_x_data%atomic_pair_list_forces
841 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
842 : CALL build_atomic_pair_list(natom, shm_atomic_pair_list, kind_of, basis_parameter, particle_set, &
843 : do_periodic, screen_coeffs_kind, coeffs_kind_max0, log10_eps_schwarz, cell, &
844 : actual_x_data%blocks)
845 : END IF
846 :
847 : my_bin_size = SIZE(actual_x_data%distribution_forces)
848 : DO bin = 1, my_bin_size
849 : actual_x_data%distribution_forces(bin)%time_forces = 0.0_dp
850 : END DO
851 : !$OMP END MASTER
852 : !$OMP BARRIER
853 :
854 : IF (treat_forces_in_core) THEN
855 : IF (my_bin_size > 0) THEN
856 : maxval_container => actual_x_data%store_forces%maxval_container(1)
857 : maxval_cache => actual_x_data%store_forces%maxval_cache(1)
858 :
859 : integral_containers => actual_x_data%store_forces%integral_containers(:, 1)
860 : integral_caches => actual_x_data%store_forces%integral_caches(:, 1)
861 : END IF
862 : END IF
863 :
864 : nblocks = load_balance_parameter%nblocks
865 :
866 : bins_left = .TRUE.
867 : do_it = .TRUE.
868 : bin = 0
869 : DO WHILE (bins_left)
870 : IF (.NOT. do_dynamic_load_balancing) THEN
871 : bin = bin + 1
872 : IF (bin > my_bin_size) THEN
873 : do_it = .FALSE.
874 : bins_left = .FALSE.
875 : ELSE
876 : do_it = .TRUE.
877 : bins_left = .TRUE.
878 : distribution_forces => actual_x_data%distribution_forces(bin)
879 : END IF
880 : ELSE
881 : !$OMP CRITICAL(hfxderivatives_critical)
882 : shm_task_counter = shm_task_counter + 1
883 : IF (shm_task_counter <= shm_total_bins) THEN
884 : my_thread_id = shm_task_list(shm_task_counter)%thread_id
885 : my_bin_id = shm_task_list(shm_task_counter)%bin_id
886 : IF (treat_forces_in_core) THEN
887 : maxval_cache => my_x_data(irep, my_thread_id)%store_forces%maxval_cache(my_bin_id)
888 : maxval_container => my_x_data(irep, my_thread_id)%store_forces%maxval_container(my_bin_id)
889 : integral_caches => my_x_data(irep, my_thread_id)%store_forces%integral_caches(:, my_bin_id)
890 : integral_containers => my_x_data(irep, my_thread_id)%store_forces%integral_containers(:, my_bin_id)
891 : END IF
892 : distribution_forces => my_x_data(irep, my_thread_id)%distribution_forces(my_bin_id)
893 : do_it = .TRUE.
894 : bins_left = .TRUE.
895 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
896 : distribution_forces%ram_counter = HUGE(distribution_forces%ram_counter)
897 : END IF
898 : counter = 0_Int_8
899 : ELSE
900 : do_it = .FALSE.
901 : bins_left = .FALSE.
902 : END IF
903 : !$OMP END CRITICAL(hfxderivatives_critical)
904 : END IF
905 :
906 : IF (.NOT. do_it) CYCLE
907 : !$OMP MASTER
908 : CALL timeset(routineN//"_bin", handle_bin)
909 : !$OMP END MASTER
910 : bintime_start = m_walltime()
911 : my_istart = distribution_forces%istart
912 : my_current_counter = 0
913 : IF (distribution_forces%number_of_atom_quartets == 0 .OR. &
914 : my_istart == -1_int_8) my_istart = nblocks**4
915 : atomic_blocks: DO atom_block = my_istart, nblocks**4 - 1, n_processes
916 : latom_block = INT(MODULO(atom_block, nblocks)) + 1
917 : tmp_block = atom_block/nblocks
918 : katom_block = INT(MODULO(tmp_block, nblocks)) + 1
919 : IF (latom_block < katom_block) CYCLE atomic_blocks
920 : tmp_block = tmp_block/nblocks
921 : jatom_block = INT(MODULO(tmp_block, nblocks)) + 1
922 : tmp_block = tmp_block/nblocks
923 : iatom_block = INT(MODULO(tmp_block, nblocks)) + 1
924 : IF (jatom_block < iatom_block) CYCLE atomic_blocks
925 : my_current_counter = my_current_counter + 1
926 : IF (my_current_counter > distribution_forces%number_of_atom_quartets) EXIT atomic_blocks
927 :
928 : iatom_start = actual_x_data%blocks(iatom_block)%istart
929 : iatom_end = actual_x_data%blocks(iatom_block)%iend
930 : jatom_start = actual_x_data%blocks(jatom_block)%istart
931 : jatom_end = actual_x_data%blocks(jatom_block)%iend
932 : katom_start = actual_x_data%blocks(katom_block)%istart
933 : katom_end = actual_x_data%blocks(katom_block)%iend
934 : latom_start = actual_x_data%blocks(latom_block)%istart
935 : latom_end = actual_x_data%blocks(latom_block)%iend
936 :
937 : pmax_blocks = MAX(shm_pmax_block(katom_block, iatom_block) + &
938 : shm_pmax_block(latom_block, jatom_block), &
939 : shm_pmax_block(latom_block, iatom_block) + &
940 : shm_pmax_block(katom_block, jatom_block))
941 :
942 : IF (2.0_dp*coeffs_kind_max0 + pmax_blocks < log10_eps_schwarz) CYCLE atomic_blocks
943 :
944 : CALL build_pair_list(natom, list_ij, set_list_ij, iatom_start, iatom_end, &
945 : jatom_start, jatom_end, &
946 : kind_of, basis_parameter, particle_set, &
947 : do_periodic, screen_coeffs_set, screen_coeffs_kind, coeffs_kind_max0, &
948 : log10_eps_schwarz, cell, pmax_blocks, shm_atomic_pair_list)
949 :
950 : CALL build_pair_list(natom, list_kl, set_list_kl, katom_start, katom_end, &
951 : latom_start, latom_end, &
952 : kind_of, basis_parameter, particle_set, &
953 : do_periodic, screen_coeffs_set, screen_coeffs_kind, coeffs_kind_max0, &
954 : log10_eps_schwarz, cell, pmax_blocks, shm_atomic_pair_list)
955 :
956 : DO i_list_ij = 1, list_ij%n_element
957 : iatom = list_ij%elements(i_list_ij)%pair(1)
958 : jatom = list_ij%elements(i_list_ij)%pair(2)
959 : i_set_list_ij_start = list_ij%elements(i_list_ij)%set_bounds(1)
960 : i_set_list_ij_stop = list_ij%elements(i_list_ij)%set_bounds(2)
961 : ikind = list_ij%elements(i_list_ij)%kind_pair(1)
962 : jkind = list_ij%elements(i_list_ij)%kind_pair(2)
963 : ra = list_ij%elements(i_list_ij)%r1
964 : rb = list_ij%elements(i_list_ij)%r2
965 : rab2 = list_ij%elements(i_list_ij)%dist2
966 :
967 : la_max => basis_parameter(ikind)%lmax
968 : la_min => basis_parameter(ikind)%lmin
969 : npgfa => basis_parameter(ikind)%npgf
970 : nseta = basis_parameter(ikind)%nset
971 : zeta => basis_parameter(ikind)%zet
972 : nsgfa => basis_parameter(ikind)%nsgf
973 : sphi_a_ext => basis_parameter(ikind)%sphi_ext(:, :, :, :)
974 : nsgfl_a => basis_parameter(ikind)%nsgfl
975 : sphi_a_u1 = UBOUND(sphi_a_ext, 1)
976 : sphi_a_u2 = UBOUND(sphi_a_ext, 2)
977 : sphi_a_u3 = UBOUND(sphi_a_ext, 3)
978 :
979 : lb_max => basis_parameter(jkind)%lmax
980 : lb_min => basis_parameter(jkind)%lmin
981 : npgfb => basis_parameter(jkind)%npgf
982 : nsetb = basis_parameter(jkind)%nset
983 : zetb => basis_parameter(jkind)%zet
984 : nsgfb => basis_parameter(jkind)%nsgf
985 : sphi_b_ext => basis_parameter(jkind)%sphi_ext(:, :, :, :)
986 : nsgfl_b => basis_parameter(jkind)%nsgfl
987 : sphi_b_u1 = UBOUND(sphi_b_ext, 1)
988 : sphi_b_u2 = UBOUND(sphi_b_ext, 2)
989 : sphi_b_u3 = UBOUND(sphi_b_ext, 3)
990 :
991 : i_atom = atom_of_kind(iatom)
992 : forces_map(1, 1) = ikind
993 : forces_map(1, 2) = i_atom
994 : j_atom = atom_of_kind(jatom)
995 : forces_map(2, 1) = jkind
996 : forces_map(2, 2) = j_atom
997 :
998 : DO i_list_kl = 1, list_kl%n_element
999 : katom = list_kl%elements(i_list_kl)%pair(1)
1000 : latom = list_kl%elements(i_list_kl)%pair(2)
1001 :
1002 : !All four centers equivalent => zero-contribution
1003 : !VIIRAL IF((iatom==jatom .AND. iatom==katom .AND. iatom==latom)) CYCLE
1004 : IF (.NOT. (katom + latom <= iatom + jatom)) CYCLE
1005 : IF (((iatom + jatom) == (katom + latom)) .AND. (katom < iatom)) CYCLE
1006 :
1007 : i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1)
1008 : i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2)
1009 : kkind = list_kl%elements(i_list_kl)%kind_pair(1)
1010 : lkind = list_kl%elements(i_list_kl)%kind_pair(2)
1011 : rc = list_kl%elements(i_list_kl)%r1
1012 : rd = list_kl%elements(i_list_kl)%r2
1013 : rcd2 = list_kl%elements(i_list_kl)%dist2
1014 :
1015 : IF (screen_pmat_forces) THEN
1016 : pmax_atom = MAX(shm_pmax_atom(katom, iatom) + &
1017 : shm_pmax_atom(latom, jatom), &
1018 : shm_pmax_atom(latom, iatom) + &
1019 : shm_pmax_atom(katom, jatom))
1020 : ELSE
1021 : pmax_atom = 0.0_dp
1022 : END IF
1023 :
1024 : IF ((screen_coeffs_kind(jkind, ikind)%x(1)*rab2 + &
1025 : screen_coeffs_kind(jkind, ikind)%x(2)) + &
1026 : (screen_coeffs_kind(lkind, kkind)%x(1)*rcd2 + &
1027 : screen_coeffs_kind(lkind, kkind)%x(2)) + pmax_atom < log10_eps_schwarz) CYCLE
1028 :
1029 : IF (.NOT. (shm_is_assoc_atomic_block(latom, iatom) >= 1 .AND. &
1030 : shm_is_assoc_atomic_block(katom, iatom) >= 1 .AND. &
1031 : shm_is_assoc_atomic_block(katom, jatom) >= 1 .AND. &
1032 : shm_is_assoc_atomic_block(latom, jatom) >= 1)) CYCLE
1033 : k_atom = atom_of_kind(katom)
1034 : forces_map(3, 1) = kkind
1035 : forces_map(3, 2) = k_atom
1036 :
1037 : l_atom = atom_of_kind(latom)
1038 : forces_map(4, 1) = lkind
1039 : forces_map(4, 2) = l_atom
1040 :
1041 : IF (my_nspins == 1) THEN
1042 : fac = 0.25_dp*hf_fraction
1043 : ELSE
1044 : fac = 0.5_dp*hf_fraction
1045 : END IF
1046 : !calculate symmetry_factor
1047 : symm_fac = 0.25_dp
1048 : IF (iatom == jatom) symm_fac = symm_fac*2.0_dp
1049 : IF (katom == latom) symm_fac = symm_fac*2.0_dp
1050 : IF (iatom == katom .AND. jatom == latom .AND. &
1051 : iatom /= jatom .AND. katom /= latom) symm_fac = symm_fac*2.0_dp
1052 : IF (iatom == katom .AND. iatom == jatom .AND. &
1053 : katom == latom) symm_fac = symm_fac*2.0_dp
1054 :
1055 : symm_fac = 1.0_dp/symm_fac
1056 : fac = fac*symm_fac
1057 :
1058 : lc_max => basis_parameter(kkind)%lmax
1059 : lc_min => basis_parameter(kkind)%lmin
1060 : npgfc => basis_parameter(kkind)%npgf
1061 : zetc => basis_parameter(kkind)%zet
1062 : nsgfc => basis_parameter(kkind)%nsgf
1063 : sphi_c_ext => basis_parameter(kkind)%sphi_ext(:, :, :, :)
1064 : nsgfl_c => basis_parameter(kkind)%nsgfl
1065 : sphi_c_u1 = UBOUND(sphi_c_ext, 1)
1066 : sphi_c_u2 = UBOUND(sphi_c_ext, 2)
1067 : sphi_c_u3 = UBOUND(sphi_c_ext, 3)
1068 :
1069 : ld_max => basis_parameter(lkind)%lmax
1070 : ld_min => basis_parameter(lkind)%lmin
1071 : npgfd => basis_parameter(lkind)%npgf
1072 : zetd => basis_parameter(lkind)%zet
1073 : nsgfd => basis_parameter(lkind)%nsgf
1074 : sphi_d_ext => basis_parameter(lkind)%sphi_ext(:, :, :, :)
1075 : nsgfl_d => basis_parameter(lkind)%nsgfl
1076 : sphi_d_u1 = UBOUND(sphi_d_ext, 1)
1077 : sphi_d_u2 = UBOUND(sphi_d_ext, 2)
1078 : sphi_d_u3 = UBOUND(sphi_d_ext, 3)
1079 :
1080 : atomic_offset_bd = shm_atomic_block_offset(jatom, latom)
1081 : atomic_offset_bc = shm_atomic_block_offset(jatom, katom)
1082 : atomic_offset_ad = shm_atomic_block_offset(iatom, latom)
1083 : atomic_offset_ac = shm_atomic_block_offset(iatom, katom)
1084 :
1085 : IF (jatom < latom) THEN
1086 : offset_bd_set => shm_set_offset(:, :, lkind, jkind)
1087 : ELSE
1088 : offset_bd_set => shm_set_offset(:, :, jkind, lkind)
1089 : END IF
1090 : IF (jatom < katom) THEN
1091 : offset_bc_set => shm_set_offset(:, :, kkind, jkind)
1092 : ELSE
1093 : offset_bc_set => shm_set_offset(:, :, jkind, kkind)
1094 : END IF
1095 : IF (iatom < latom) THEN
1096 : offset_ad_set => shm_set_offset(:, :, lkind, ikind)
1097 : ELSE
1098 : offset_ad_set => shm_set_offset(:, :, ikind, lkind)
1099 : END IF
1100 : IF (iatom < katom) THEN
1101 : offset_ac_set => shm_set_offset(:, :, kkind, ikind)
1102 : ELSE
1103 : offset_ac_set => shm_set_offset(:, :, ikind, kkind)
1104 : END IF
1105 :
1106 : IF (screen_pmat_forces) THEN
1107 : swap_id = 16
1108 : kind_kind_idx = INT(get_1D_idx(kkind, ikind, INT(nkind, int_8)))
1109 : IF (ikind >= kkind) THEN
1110 : ptr_p_1 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1111 : actual_x_data%map_atom_to_kind_atom(katom), &
1112 : actual_x_data%map_atom_to_kind_atom(iatom))
1113 : ELSE
1114 : ptr_p_1 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1115 : actual_x_data%map_atom_to_kind_atom(iatom), &
1116 : actual_x_data%map_atom_to_kind_atom(katom))
1117 : swap_id = swap_id + 1
1118 : END IF
1119 : kind_kind_idx = INT(get_1D_idx(lkind, jkind, INT(nkind, int_8)))
1120 : IF (jkind >= lkind) THEN
1121 : ptr_p_2 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1122 : actual_x_data%map_atom_to_kind_atom(latom), &
1123 : actual_x_data%map_atom_to_kind_atom(jatom))
1124 : ELSE
1125 : ptr_p_2 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1126 : actual_x_data%map_atom_to_kind_atom(jatom), &
1127 : actual_x_data%map_atom_to_kind_atom(latom))
1128 : swap_id = swap_id + 2
1129 : END IF
1130 : kind_kind_idx = INT(get_1D_idx(lkind, ikind, INT(nkind, int_8)))
1131 : IF (ikind >= lkind) THEN
1132 : ptr_p_3 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1133 : actual_x_data%map_atom_to_kind_atom(latom), &
1134 : actual_x_data%map_atom_to_kind_atom(iatom))
1135 : ELSE
1136 : ptr_p_3 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1137 : actual_x_data%map_atom_to_kind_atom(iatom), &
1138 : actual_x_data%map_atom_to_kind_atom(latom))
1139 : swap_id = swap_id + 4
1140 : END IF
1141 : kind_kind_idx = INT(get_1D_idx(kkind, jkind, INT(nkind, int_8)))
1142 : IF (jkind >= kkind) THEN
1143 : ptr_p_4 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1144 : actual_x_data%map_atom_to_kind_atom(katom), &
1145 : actual_x_data%map_atom_to_kind_atom(jatom))
1146 : ELSE
1147 : ptr_p_4 => shm_initial_p(kind_kind_idx)%p_kind(:, :, &
1148 : actual_x_data%map_atom_to_kind_atom(jatom), &
1149 : actual_x_data%map_atom_to_kind_atom(katom))
1150 : swap_id = swap_id + 8
1151 : END IF
1152 : END IF
1153 :
1154 : !! At this stage, check for memory used in compression
1155 :
1156 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
1157 : IF (treat_forces_in_core) THEN
1158 : ! ** We know the maximum amount of integrals that we can store per MPI-process
1159 : ! ** Now we need to sum the current memory usage among all openMP threads
1160 : ! ** We can just read what is currently stored on the corresponding x_data type
1161 : ! ** If this is thread i and it tries to read the data from thread j, that is
1162 : ! ** currently writing that data, we just dont care, because the possible error
1163 : ! ** is of the order of CACHE_SIZE
1164 : mem_compression_counter = 0
1165 : DO i = 1, n_threads
1166 : !$OMP ATOMIC READ
1167 : tmp_i4 = my_x_data(irep, i)%memory_parameter%actual_memory_usage
1168 : mem_compression_counter = mem_compression_counter + &
1169 : tmp_i4*memory_parameter%cache_size
1170 : END DO
1171 : IF (mem_compression_counter + memory_parameter%final_comp_counter_energy &
1172 : > memory_parameter%max_compression_counter) THEN
1173 : buffer_overflow = .TRUE.
1174 : IF (do_dynamic_load_balancing) THEN
1175 : distribution_forces%ram_counter = counter
1176 : ELSE
1177 : memory_parameter%ram_counter_forces = counter
1178 : END IF
1179 : ELSE
1180 : counter = counter + 1
1181 : buffer_overflow = .FALSE.
1182 : END IF
1183 : END IF
1184 : ELSE
1185 : IF (treat_forces_in_core) THEN
1186 : IF (do_dynamic_load_balancing) THEN
1187 : IF (distribution_forces%ram_counter == counter) THEN
1188 : buffer_overflow = .TRUE.
1189 : ELSE
1190 : counter = counter + 1
1191 : buffer_overflow = .FALSE.
1192 : END IF
1193 : ELSE
1194 : IF (memory_parameter%ram_counter_forces == counter) THEN
1195 : buffer_overflow = .TRUE.
1196 : ELSE
1197 : counter = counter + 1
1198 : buffer_overflow = .FALSE.
1199 : END IF
1200 : END IF
1201 : END IF
1202 : END IF
1203 :
1204 : DO i_set_list_ij = i_set_list_ij_start, i_set_list_ij_stop
1205 : iset = set_list_ij(i_set_list_ij)%pair(1)
1206 : jset = set_list_ij(i_set_list_ij)%pair(2)
1207 :
1208 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1209 : max_val1 = screen_coeffs_set(jset, iset, jkind, ikind)%x(1)*rab2 + &
1210 : screen_coeffs_set(jset, iset, jkind, ikind)%x(2)
1211 : !! Near field screening
1212 : IF (max_val1 + (screen_coeffs_kind(lkind, kkind)%x(1)*rcd2 + &
1213 : screen_coeffs_kind(lkind, kkind)%x(2)) + pmax_atom < log10_eps_schwarz) CYCLE
1214 :
1215 : DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop
1216 : kset = set_list_kl(i_set_list_kl)%pair(1)
1217 : lset = set_list_kl(i_set_list_kl)%pair(2)
1218 :
1219 : max_val2_set = (screen_coeffs_set(lset, kset, lkind, kkind)%x(1)*rcd2 + &
1220 : screen_coeffs_set(lset, kset, lkind, kkind)%x(2))
1221 : max_val2 = max_val1 + max_val2_set
1222 :
1223 : !! Near field screening
1224 : IF (max_val2 + pmax_atom < log10_eps_schwarz) CYCLE
1225 :
1226 : IF (screen_pmat_forces) THEN
1227 : CALL get_pmax_val(ptr_p_1, ptr_p_2, ptr_p_3, ptr_p_4, &
1228 : iset, jset, kset, lset, &
1229 : pmax_tmp, swap_id)
1230 :
1231 : log10_pmax = log_2 + pmax_tmp
1232 : ELSE
1233 : log10_pmax = 0.0_dp
1234 : END IF
1235 :
1236 : max_val2 = max_val2 + log10_pmax
1237 : IF (max_val2 < log10_eps_schwarz) CYCLE
1238 :
1239 : pmax_entry = EXP(log10_pmax*ln_10)
1240 : !! store current number of integrals, update total number and number of integrals in buffer
1241 : current_counter = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)*12
1242 : IF (buffer_overflow) THEN
1243 : neris_onthefly = neris_onthefly + current_counter
1244 : END IF
1245 :
1246 : !! Get integrals from buffer and update Kohn-Sham matrix
1247 : IF (.NOT. buffer_overflow .AND. .NOT. actual_x_data%memory_parameter%recalc_forces) THEN
1248 : nints = current_counter
1249 : CALL hfx_get_single_cache_element(estimate_to_store_int, 6, &
1250 : maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
1251 : use_disk_storage)
1252 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
1253 : ! IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE
1254 : IF (.NOT. use_virial) THEN
1255 : IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE
1256 : END IF
1257 : nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1
1258 : buffer_left = nints
1259 : buffer_start = 1
1260 : neris_incore = neris_incore + INT(nints, int_8)
1261 : DO WHILE (buffer_left > 0)
1262 : buffer_size = MIN(buffer_left, cache_size)
1263 : CALL hfx_get_mult_cache_elements(primitive_forces(buffer_start), &
1264 : buffer_size, nbits, &
1265 : integral_caches(nbits), &
1266 : integral_containers(nbits), &
1267 : eps_storage, pmax_entry, &
1268 : memory_parameter%actual_memory_usage, &
1269 : use_disk_storage)
1270 : buffer_left = buffer_left - buffer_size
1271 : buffer_start = buffer_start + buffer_size
1272 : END DO
1273 : END IF
1274 : !! Calculate integrals if we run out of buffer or the geometry did change
1275 : IF (actual_x_data%memory_parameter%recalc_forces .OR. buffer_overflow) THEN
1276 : sphi_a_ext_set => sphi_a_ext(:, :, :, iset)
1277 : sphi_b_ext_set => sphi_b_ext(:, :, :, jset)
1278 : sphi_c_ext_set => sphi_c_ext(:, :, :, kset)
1279 : sphi_d_ext_set => sphi_d_ext(:, :, :, lset)
1280 : max_contraction_val = max_contraction(iset, iatom)* &
1281 : max_contraction(jset, jatom)* &
1282 : max_contraction(kset, katom)* &
1283 : max_contraction(lset, latom)* &
1284 : pmax_entry
1285 : tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind)
1286 : tmp_R_2 => radii_pgf(:, :, lset, kset, lkind, kkind)
1287 : tmp_screen_pgf1 => screen_coeffs_pgf(:, :, jset, iset, jkind, ikind)
1288 : tmp_screen_pgf2 => screen_coeffs_pgf(:, :, lset, kset, lkind, kkind)
1289 :
1290 : CALL forces4(private_context, ra, rb, rc, rd, npgfa(iset), npgfb(jset), &
1291 : npgfc(kset), npgfd(lset), &
1292 : la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), &
1293 : lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), &
1294 : nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1295 : sphi_a_u1, sphi_a_u2, sphi_a_u3, &
1296 : sphi_b_u1, sphi_b_u2, sphi_b_u3, &
1297 : sphi_c_u1, sphi_c_u2, sphi_c_u3, &
1298 : sphi_d_u1, sphi_d_u2, sphi_d_u3, &
1299 : zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), &
1300 : zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), &
1301 : primitive_forces, &
1302 : potential_parameter, &
1303 : actual_x_data%neighbor_cells, &
1304 : screen_coeffs_set(jset, iset, jkind, ikind)%x, &
1305 : screen_coeffs_set(lset, kset, lkind, kkind)%x, eps_schwarz, &
1306 : max_contraction_val, cartesian_estimate, cell, neris_tmp, &
1307 : log10_pmax, log10_eps_schwarz, &
1308 : tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, &
1309 : pgf_list_ij, pgf_list_kl, pgf_product_list, &
1310 : nsgfl_a(:, iset), nsgfl_b(:, jset), &
1311 : nsgfl_c(:, kset), nsgfl_d(:, lset), &
1312 : sphi_a_ext_set, sphi_b_ext_set, sphi_c_ext_set, sphi_d_ext_set, &
1313 : ede_work, ede_work2, ede_work_forces, &
1314 : ede_buffer1, ede_buffer2, ede_primitives_tmp, &
1315 : nimages, do_periodic, use_virial, ede_work_virial, ede_work2_virial, &
1316 : ede_primitives_tmp_virial, primitive_forces_virial, &
1317 : cartesian_estimate_virial)
1318 :
1319 : nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset)*12
1320 : neris_total = neris_total + nints
1321 : nprim_ints = nprim_ints + neris_tmp
1322 : ! IF (cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate)
1323 : ! estimate_to_store_int = EXPONENT(cartesian_estimate)
1324 : ! estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
1325 : ! cartesian_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int+1)
1326 : ! IF (.NOT. buffer_overflow .AND. actual_x_data%memory_parameter%recalc_forces) THEN
1327 : ! IF (cartesian_estimate < eps_schwarz) THEN
1328 : ! CALL hfx_add_single_cache_element( &
1329 : ! estimate_to_store_int, 6, &
1330 : ! maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
1331 : ! use_disk_storage, max_val_memory)
1332 : ! END IF
1333 : ! END IF
1334 : !
1335 : ! IF (.NOT. use_virial) THEN
1336 : ! IF (cartesian_estimate < eps_schwarz) CYCLE
1337 : ! END IF
1338 :
1339 : !! Compress the array for storage
1340 : spherical_estimate = 0.0_dp
1341 : IF (cartesian_estimate /= 0.0_dp) THEN
1342 : DO i = 1, nints
1343 : spherical_estimate = MAX(spherical_estimate, ABS(primitive_forces(i)))
1344 : END DO
1345 : END IF
1346 :
1347 : spherical_estimate_virial = 0.0_dp
1348 : IF (use_virial .AND. cartesian_estimate_virial /= 0.0_dp) THEN
1349 : DO i = 1, 3*nints
1350 : spherical_estimate_virial = MAX(spherical_estimate_virial, ABS(primitive_forces_virial(i)))
1351 : END DO
1352 : END IF
1353 :
1354 : IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate)
1355 : estimate_to_store_int = EXPONENT(spherical_estimate)
1356 : estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
1357 : IF (.NOT. buffer_overflow .AND. actual_x_data%memory_parameter%recalc_forces) THEN
1358 : CALL hfx_add_single_cache_element(estimate_to_store_int, 6, &
1359 : maxval_cache, maxval_container, &
1360 : memory_parameter%actual_memory_usage, &
1361 : use_disk_storage, max_val_memory)
1362 : END IF
1363 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
1364 : IF (.NOT. use_virial) THEN
1365 : IF (spherical_estimate*pmax_entry < eps_schwarz) CYCLE
1366 : END IF
1367 : IF (.NOT. buffer_overflow) THEN
1368 : nbits = EXPONENT(ANINT(spherical_estimate*pmax_entry/eps_storage)) + 1
1369 : buffer_left = nints
1370 : buffer_start = 1
1371 : ! neris_incore = neris_incore+nints
1372 : neris_incore = neris_incore + INT(nints, int_8)
1373 :
1374 : DO WHILE (buffer_left > 0)
1375 : buffer_size = MIN(buffer_left, CACHE_SIZE)
1376 : CALL hfx_add_mult_cache_elements(primitive_forces(buffer_start), &
1377 : buffer_size, nbits, &
1378 : integral_caches(nbits), &
1379 : integral_containers(nbits), &
1380 : eps_storage, pmax_entry, &
1381 : memory_parameter%actual_memory_usage, &
1382 : use_disk_storage)
1383 : buffer_left = buffer_left - buffer_size
1384 : buffer_start = buffer_start + buffer_size
1385 : END DO
1386 : ELSE
1387 : !! In order to be consistent with in-core part, round all the eris wrt. eps_schwarz
1388 : !! but only if we treat forces in-core
1389 : IF (treat_forces_in_core .AND. cartesian_estimate /= 0.0_dp) THEN
1390 : DO i = 1, nints
1391 : primitive_forces(i) = primitive_forces(i)*pmax_entry
1392 : IF (ABS(primitive_forces(i)) > eps_storage) THEN
1393 : primitive_forces(i) = ANINT(primitive_forces(i)/eps_storage, dp)*eps_storage/pmax_entry
1394 : ELSE
1395 : primitive_forces(i) = 0.0_dp
1396 : END IF
1397 : END DO
1398 : END IF
1399 : END IF
1400 : IF (cartesian_estimate == 0.0_dp .AND. &
1401 : (.NOT. use_virial .OR. cartesian_estimate_virial == 0.0_dp)) CYCLE
1402 : END IF
1403 : IF (with_resp_density) THEN
1404 : CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1405 : full_density_alpha(:, 1), pbd_buf, pbc_buf, pad_buf, pac_buf, &
1406 : iatom, jatom, katom, latom, &
1407 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1408 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1409 : atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
1410 : CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1411 : full_density_resp, pbd_buf_resp, pbc_buf_resp, pad_buf_resp, pac_buf_resp, &
1412 : iatom, jatom, katom, latom, &
1413 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1414 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1415 : atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
1416 : ELSE
1417 : CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1418 : full_density_alpha(:, 1), pbd_buf, pbc_buf, pad_buf, pac_buf, &
1419 : iatom, jatom, katom, latom, &
1420 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1421 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1422 : atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
1423 : END IF
1424 : IF (my_nspins == 2) THEN
1425 : IF (with_resp_density) THEN
1426 : CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1427 : full_density_beta(:, 1), pbd_buf_beta, pbc_buf_beta, &
1428 : pad_buf_beta, pac_buf_beta, iatom, jatom, katom, latom, &
1429 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1430 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1431 : atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
1432 : CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1433 : full_density_resp_beta, pbd_buf_resp_beta, pbc_buf_resp_beta, &
1434 : pad_buf_resp_beta, pac_buf_resp_beta, iatom, jatom, katom, latom, &
1435 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1436 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1437 : atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
1438 : ELSE
1439 : CALL prefetch_density_matrix(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1440 : full_density_beta(:, 1), pbd_buf_beta, pbc_buf_beta, pad_buf_beta, &
1441 : pac_buf_beta, iatom, jatom, katom, latom, &
1442 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1443 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1444 : atomic_offset_ad, atomic_offset_ac, is_anti_symmetric)
1445 : END IF
1446 : END IF
1447 : IF (spherical_estimate*pmax_entry >= eps_schwarz) THEN
1448 : DO coord = 1, 12
1449 : T2 => primitive_forces((coord - 1)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) + 1: &
1450 : coord*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset))
1451 :
1452 : IF (with_resp_density) THEN
1453 : CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1454 : pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
1455 : T2, force, forces_map, coord, &
1456 : pbd_buf_resp, pbc_buf_resp, pad_buf_resp, pac_buf_resp, &
1457 : my_resp_only)
1458 : ELSE
1459 : CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1460 : pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
1461 : T2, force, forces_map, coord)
1462 : END IF
1463 : IF (my_nspins == 2) THEN
1464 : IF (with_resp_density) THEN
1465 : CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1466 : pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
1467 : fac, T2, force, forces_map, coord, &
1468 : pbd_buf_resp_beta, pbc_buf_resp_beta, &
1469 : pad_buf_resp_beta, pac_buf_resp_beta, &
1470 : my_resp_only)
1471 : ELSE
1472 : CALL update_forces(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1473 : pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, fac, &
1474 : T2, force, forces_map, coord)
1475 : END IF
1476 : END IF
1477 : END DO !coord
1478 : END IF
1479 : IF (use_virial .AND. spherical_estimate_virial*pmax_entry >= eps_schwarz) THEN
1480 : DO coord = 1, 12
1481 : DO i = 1, 3
1482 : T2 => primitive_forces_virial( &
1483 : ((i - 1)*12 + coord - 1)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) + 1: &
1484 : ((i - 1)*12 + coord)*nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset))
1485 : IF (with_resp_density) THEN
1486 : CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1487 : pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
1488 : T2, tmp_virial, coord, i, &
1489 : pbd_buf_resp, pbc_buf_resp, pad_buf_resp, &
1490 : pac_buf_resp, my_resp_only)
1491 : ELSE
1492 : CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1493 : pbd_buf, pbc_buf, pad_buf, pac_buf, fac, &
1494 : T2, tmp_virial, coord, i)
1495 : END IF
1496 : IF (my_nspins == 2) THEN
1497 : IF (with_resp_density) THEN
1498 : CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1499 : pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
1500 : fac, T2, tmp_virial, coord, i, &
1501 : pbd_buf_resp_beta, pbc_buf_resp_beta, &
1502 : pad_buf_resp_beta, pac_buf_resp_beta, &
1503 : my_resp_only)
1504 : ELSE
1505 : CALL update_virial(nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), &
1506 : pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta, &
1507 : fac, T2, tmp_virial, coord, i)
1508 : END IF
1509 : END IF
1510 : END DO
1511 : END DO !coord
1512 : END IF
1513 :
1514 : END DO !i_set_list_kl
1515 : END DO !i_set_list_ij
1516 : END DO !i_list_kl
1517 : END DO !i_list_ij
1518 : END DO atomic_blocks
1519 : bintime_stop = m_walltime()
1520 : distribution_forces%time_forces = bintime_stop - bintime_start
1521 : !$OMP MASTER
1522 : CALL timestop(handle_bin)
1523 : !$OMP END MASTER
1524 : END DO !bin
1525 :
1526 : !$OMP MASTER
1527 : logger => cp_get_default_logger()
1528 : do_print_load_balance_info = .FALSE.
1529 : do_print_load_balance_info = BTEST(cp_print_key_should_output(logger%iter_info, hfx_section, &
1530 : "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO"), cp_p_file)
1531 : !$OMP END MASTER
1532 : !$OMP BARRIER
1533 : IF (do_print_load_balance_info) THEN
1534 : iw = -1
1535 : !$OMP MASTER
1536 : iw = cp_print_key_unit_nr(logger, hfx_section, "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO", &
1537 : extension=".scfLog")
1538 : !$OMP END MASTER
1539 :
1540 : CALL collect_load_balance_info(para_env, actual_x_data, iw, n_threads, i_thread, &
1541 : hfx_do_eval_forces)
1542 :
1543 : !$OMP MASTER
1544 : CALL cp_print_key_finished_output(iw, logger, hfx_section, &
1545 : "LOAD_BALANCE%PRINT/LOAD_BALANCE_INFO")
1546 : !$OMP END MASTER
1547 : END IF
1548 :
1549 : IF (use_virial) THEN
1550 : DO i = 1, 3
1551 : DO j = 1, 3
1552 : DO k = 1, 3
1553 : !$OMP ATOMIC
1554 : virial%pv_fock_4c(i, j) = virial%pv_fock_4c(i, j) + tmp_virial(i, k)*cell%hmat(j, k)
1555 : END DO
1556 : END DO
1557 : END DO
1558 : END IF
1559 :
1560 : !$OMP BARRIER
1561 : !! Get some number about ERIS
1562 : !$OMP ATOMIC
1563 : shm_neris_total = shm_neris_total + neris_total
1564 : !$OMP ATOMIC
1565 : shm_neris_onthefly = shm_neris_onthefly + neris_onthefly
1566 : !$OMP ATOMIC
1567 : shm_nprim_ints = shm_nprim_ints + nprim_ints
1568 :
1569 : storage_counter_integrals = memory_parameter%actual_memory_usage* &
1570 : memory_parameter%cache_size
1571 : stor_count_max_val = max_val_memory*memory_parameter%cache_size
1572 : !$OMP ATOMIC
1573 : shm_storage_counter_integrals = shm_storage_counter_integrals + storage_counter_integrals
1574 : !$OMP ATOMIC
1575 : shm_neris_incore = shm_neris_incore + neris_incore
1576 : !$OMP ATOMIC
1577 : shm_stor_count_max_val = shm_stor_count_max_val + stor_count_max_val
1578 : !$OMP BARRIER
1579 :
1580 : ! IF(with_resp_density) THEN
1581 : ! ! restore original load_balance_parameter
1582 : ! actual_x_data%load_balance_parameter = load_balance_parameter_energy
1583 : ! DEALLOCATE(load_balance_parameter_energy)
1584 : ! END IF
1585 :
1586 : !$OMP MASTER
1587 : !! Print some memeory information if this is the first step
1588 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
1589 : tmp_i8(1:6) = [shm_storage_counter_integrals, shm_neris_onthefly, shm_neris_incore, &
1590 : shm_neris_total, shm_nprim_ints, shm_stor_count_max_val]
1591 : CALL para_env%sum(tmp_i8)
1592 : shm_storage_counter_integrals = tmp_i8(1)
1593 : shm_neris_onthefly = tmp_i8(2)
1594 : shm_neris_incore = tmp_i8(3)
1595 : shm_neris_total = tmp_i8(4)
1596 : shm_nprim_ints = tmp_i8(5)
1597 : shm_stor_count_max_val = tmp_i8(6)
1598 : mem_eris = (shm_storage_counter_integrals + 128*1024 - 1)/1024/128
1599 : compression_factor = REAL(shm_neris_incore, dp)/REAL(shm_storage_counter_integrals, dp)
1600 : mem_max_val = (shm_stor_count_max_val + 128*1024 - 1)/1024/128
1601 :
1602 : IF (shm_neris_incore == 0) THEN
1603 : mem_eris = 0
1604 : compression_factor = 0.0_dp
1605 : END IF
1606 :
1607 : logger => cp_get_default_logger()
1608 :
1609 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
1610 : extension=".scfLog")
1611 : IF (iw > 0) THEN
1612 : WRITE (UNIT=iw, FMT="(/,(T3,A,T65,I16))") &
1613 : "HFX_MEM_INFO| Number of cart. primitive DERIV's calculated: ", shm_nprim_ints
1614 :
1615 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
1616 : "HFX_MEM_INFO| Number of sph. DERIV's calculated: ", shm_neris_total
1617 :
1618 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
1619 : "HFX_MEM_INFO| Number of sph. DERIV's stored in-core: ", shm_neris_incore
1620 :
1621 : WRITE (UNIT=iw, FMT="((T3,A,T62,I19))") &
1622 : "HFX_MEM_INFO| Number of sph. DERIV's calculated on the fly: ", shm_neris_onthefly
1623 :
1624 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
1625 : "HFX_MEM_INFO| Total memory consumption DERIV's RAM [MiB]: ", mem_eris
1626 :
1627 : WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
1628 : "HFX_MEM_INFO| Whereof max-vals [MiB]: ", mem_max_val
1629 :
1630 : WRITE (UNIT=iw, FMT="((T3,A,T60,F21.2),/)") &
1631 : "HFX_MEM_INFO| Total compression factor DERIV's RAM: ", compression_factor
1632 : END IF
1633 :
1634 : CALL cp_print_key_finished_output(iw, logger, hfx_section, &
1635 : "HF_INFO")
1636 : END IF
1637 : !$OMP END MASTER
1638 :
1639 : !! flush caches if the geometry changed
1640 : IF (do_dynamic_load_balancing) THEN
1641 : my_bin_size = SIZE(actual_x_data%distribution_forces)
1642 : ELSE
1643 : my_bin_size = 1
1644 : END IF
1645 :
1646 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
1647 : IF (treat_forces_in_core) THEN
1648 : DO bin = 1, my_bin_size
1649 : maxval_cache => actual_x_data%store_forces%maxval_cache(bin)
1650 : maxval_container => actual_x_data%store_forces%maxval_container(bin)
1651 : integral_caches => actual_x_data%store_forces%integral_caches(:, bin)
1652 : integral_containers => actual_x_data%store_forces%integral_containers(:, bin)
1653 : CALL hfx_flush_last_cache(bits_max_val, maxval_cache, maxval_container, memory_parameter%actual_memory_usage, &
1654 : .FALSE.)
1655 : DO i = 1, 64
1656 : CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
1657 : memory_parameter%actual_memory_usage, .FALSE.)
1658 : END DO
1659 : END DO
1660 : END IF
1661 : END IF
1662 : !! reset all caches except we calculate all on the fly
1663 : IF (treat_forces_in_core) THEN
1664 : DO bin = 1, my_bin_size
1665 : maxval_cache => actual_x_data%store_forces%maxval_cache(bin)
1666 : maxval_container => actual_x_data%store_forces%maxval_container(bin)
1667 : integral_caches => actual_x_data%store_forces%integral_caches(:, bin)
1668 : integral_containers => actual_x_data%store_forces%integral_containers(:, bin)
1669 :
1670 : CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_parameter%actual_memory_usage, .FALSE.)
1671 : DO i = 1, 64
1672 : CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
1673 : memory_parameter%actual_memory_usage, &
1674 : .FALSE.)
1675 : END DO
1676 : END DO
1677 : END IF
1678 :
1679 : IF (actual_x_data%memory_parameter%recalc_forces) THEN
1680 : actual_x_data%memory_parameter%recalc_forces = .FALSE.
1681 : END IF
1682 :
1683 : !$OMP BARRIER
1684 : !$OMP MASTER
1685 : CALL timestop(handle_main)
1686 : !$OMP END MASTER
1687 : !$OMP BARRIER
1688 : DEALLOCATE (last_sgf_global)
1689 : !$OMP MASTER
1690 : DEALLOCATE (full_density_alpha)
1691 : IF (with_resp_density) THEN
1692 : DEALLOCATE (full_density_resp)
1693 : END IF
1694 : IF (my_nspins == 2) THEN
1695 : DEALLOCATE (full_density_beta)
1696 : IF (with_resp_density) THEN
1697 : DEALLOCATE (full_density_resp_beta)
1698 : END IF
1699 : END IF
1700 : IF (do_dynamic_load_balancing) THEN
1701 : DEALLOCATE (actual_x_data%task_list)
1702 : END IF
1703 : !$OMP END MASTER
1704 : DEALLOCATE (primitive_forces)
1705 : DEALLOCATE (atom_of_kind, kind_of)
1706 : DEALLOCATE (max_contraction)
1707 : DEALLOCATE (pbd_buf)
1708 : DEALLOCATE (pbc_buf)
1709 : DEALLOCATE (pad_buf)
1710 : DEALLOCATE (pac_buf)
1711 :
1712 : IF (with_resp_density) THEN
1713 : DEALLOCATE (pbd_buf_resp)
1714 : DEALLOCATE (pbc_buf_resp)
1715 : DEALLOCATE (pad_buf_resp)
1716 : DEALLOCATE (pac_buf_resp)
1717 : END IF
1718 :
1719 : DO i = 1, max_pgf**2
1720 : DEALLOCATE (pgf_list_ij(i)%image_list)
1721 : DEALLOCATE (pgf_list_kl(i)%image_list)
1722 : END DO
1723 :
1724 : DEALLOCATE (pgf_list_ij)
1725 : DEALLOCATE (pgf_list_kl)
1726 : DEALLOCATE (pgf_product_list)
1727 :
1728 : DEALLOCATE (set_list_ij, set_list_kl)
1729 :
1730 : DEALLOCATE (primitive_forces_virial)
1731 :
1732 : DEALLOCATE (ede_work, ede_work2, ede_work_forces, ede_buffer1, ede_buffer2, ede_primitives_tmp)
1733 :
1734 : DEALLOCATE (ede_work_virial, ede_work2_virial, ede_primitives_tmp_virial)
1735 :
1736 : DEALLOCATE (nimages)
1737 :
1738 : IF (my_nspins == 2) THEN
1739 : DEALLOCATE (pbd_buf_beta, pbc_buf_beta, pad_buf_beta, pac_buf_beta)
1740 : IF (with_resp_density) THEN
1741 : DEALLOCATE (pbd_buf_resp_beta)
1742 : DEALLOCATE (pbc_buf_resp_beta)
1743 : DEALLOCATE (pad_buf_resp_beta)
1744 : DEALLOCATE (pac_buf_resp_beta)
1745 : END IF
1746 : END IF
1747 :
1748 : !$OMP END PARALLEL
1749 :
1750 2008 : CALL timestop(handle)
1751 :
1752 4172888 : END SUBROUTINE derivatives_four_center
1753 :
1754 : ! **************************************************************************************************
1755 : !> \brief Given a 2d index pair, this function returns a 1d index pair for
1756 : !> a symmetric upper triangle NxN matrix
1757 : !> The compiler should inline this function, therefore it appears in
1758 : !> several modules
1759 : !> \param i 2d index
1760 : !> \param j 2d index
1761 : !> \param N matrix size
1762 : !> \return ...
1763 : !> \par History
1764 : !> 03.2009 created [Manuel Guidon]
1765 : !> \author Manuel Guidon
1766 : ! **************************************************************************************************
1767 63136 : PURE FUNCTION get_1D_idx(i, j, N)
1768 : INTEGER, INTENT(IN) :: i, j
1769 : INTEGER(int_8), INTENT(IN) :: N
1770 : INTEGER(int_8) :: get_1D_idx
1771 :
1772 : INTEGER(int_8) :: min_ij
1773 :
1774 63136 : min_ij = MIN(i, j)
1775 63136 : get_1D_idx = min_ij*N + MAX(i, j) - (min_ij - 1)*min_ij/2 - N
1776 :
1777 63136 : END FUNCTION get_1D_idx
1778 :
1779 : ! **************************************************************************************************
1780 : !> \brief This routine prefetches density matrix elements, i.e. reshuffles the
1781 : !> data in a way that they can be accessed later on in a cache friendly
1782 : !> way
1783 : !> \param ma_max Size of matrix blocks
1784 : !> \param mb_max Size of matrix blocks
1785 : !> \param mc_max Size of matrix blocks
1786 : !> \param md_max Size of matrix blocks
1787 : !> \param density ...
1788 : !> \param pbd buffer that will contain P(b,d)
1789 : !> \param pbc buffer that will contain P(b,c)
1790 : !> \param pad buffer that will contain P(a,d)
1791 : !> \param pac buffer that will contain P(a,c)
1792 : !> \param iatom ...
1793 : !> \param jatom ...
1794 : !> \param katom ...
1795 : !> \param latom ...
1796 : !> \param iset ...
1797 : !> \param jset ...
1798 : !> \param kset ...
1799 : !> \param lset ...
1800 : !> \param offset_bd_set ...
1801 : !> \param offset_bc_set ...
1802 : !> \param offset_ad_set ...
1803 : !> \param offset_ac_set ...
1804 : !> \param atomic_offset_bd ...
1805 : !> \param atomic_offset_bc ...
1806 : !> \param atomic_offset_ad ...
1807 : !> \param atomic_offset_ac ...
1808 : !> \param anti_symmetric ...
1809 : !> \par History
1810 : !> 03.2009 created [Manuel Guidon]
1811 : !> \author Manuel Guidon
1812 : ! **************************************************************************************************
1813 3193802 : SUBROUTINE prefetch_density_matrix(ma_max, mb_max, mc_max, md_max, &
1814 3193802 : density, pbd, pbc, pad, pac, &
1815 : iatom, jatom, katom, latom, &
1816 : iset, jset, kset, lset, offset_bd_set, offset_bc_set, &
1817 : offset_ad_set, offset_ac_set, atomic_offset_bd, atomic_offset_bc, &
1818 : atomic_offset_ad, atomic_offset_ac, anti_symmetric)
1819 :
1820 : INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
1821 : REAL(dp), DIMENSION(:), INTENT(IN) :: density
1822 : REAL(dp), DIMENSION(*), INTENT(INOUT) :: pbd, pbc, pad, pac
1823 : INTEGER, INTENT(IN) :: iatom, jatom, katom, latom, iset, jset, &
1824 : kset, lset
1825 : INTEGER, DIMENSION(:, :), POINTER :: offset_bd_set, offset_bc_set, &
1826 : offset_ad_set, offset_ac_set
1827 : INTEGER, INTENT(IN) :: atomic_offset_bd, atomic_offset_bc, &
1828 : atomic_offset_ad, atomic_offset_ac
1829 : LOGICAL :: anti_symmetric
1830 :
1831 : INTEGER :: i, j, ma, mb, mc, md, offset_ac, &
1832 : offset_ad, offset_bc, offset_bd
1833 : REAL(dp) :: fac
1834 :
1835 3193802 : fac = 1.0_dp
1836 3193802 : IF (anti_symmetric) fac = -1.0_dp
1837 3193802 : IF (jatom >= latom) THEN
1838 3187596 : i = 1
1839 3187596 : offset_bd = offset_bd_set(jset, lset) + atomic_offset_bd - 1
1840 3187596 : j = offset_bd
1841 9390202 : DO md = 1, md_max
1842 20309022 : DO mb = 1, mb_max
1843 10918820 : pbd(i) = density(j)*fac
1844 10918820 : i = i + 1
1845 17121426 : j = j + 1
1846 : END DO
1847 : END DO
1848 : ELSE
1849 6206 : i = 1
1850 6206 : offset_bd = offset_bd_set(lset, jset) + atomic_offset_bd - 1
1851 14840 : DO md = 1, md_max
1852 8634 : j = offset_bd + md - 1
1853 27686 : DO mb = 1, mb_max
1854 12846 : pbd(i) = density(j)
1855 12846 : i = i + 1
1856 21480 : j = j + md_max
1857 : END DO
1858 : END DO
1859 : END IF
1860 3193802 : IF (jatom >= katom) THEN
1861 3193802 : i = 1
1862 3193802 : offset_bc = offset_bc_set(jset, kset) + atomic_offset_bc - 1
1863 3193802 : j = offset_bc
1864 10478032 : DO mc = 1, mc_max
1865 22917310 : DO mb = 1, mb_max
1866 12439278 : pbc(i) = density(j)*fac
1867 12439278 : i = i + 1
1868 19723508 : j = j + 1
1869 : END DO
1870 : END DO
1871 : ELSE
1872 0 : i = 1
1873 0 : offset_bc = offset_bc_set(kset, jset) + atomic_offset_bc - 1
1874 0 : DO mc = 1, mc_max
1875 0 : j = offset_bc + mc - 1
1876 0 : DO mb = 1, mb_max
1877 0 : pbc(i) = density(j)
1878 0 : i = i + 1
1879 0 : j = j + mc_max
1880 : END DO
1881 : END DO
1882 : END IF
1883 3193802 : IF (iatom >= latom) THEN
1884 2457148 : i = 1
1885 2457148 : offset_ad = offset_ad_set(iset, lset) + atomic_offset_ad - 1
1886 2457148 : j = offset_ad
1887 7603951 : DO md = 1, md_max
1888 18448052 : DO ma = 1, ma_max
1889 10844101 : pad(i) = density(j)*fac
1890 10844101 : i = i + 1
1891 15990904 : j = j + 1
1892 : END DO
1893 : END DO
1894 : ELSE
1895 736654 : i = 1
1896 736654 : offset_ad = offset_ad_set(lset, iset) + atomic_offset_ad - 1
1897 1801091 : DO md = 1, md_max
1898 1064437 : j = offset_ad + md - 1
1899 4408196 : DO ma = 1, ma_max
1900 2607105 : pad(i) = density(j)
1901 2607105 : i = i + 1
1902 3671542 : j = j + md_max
1903 : END DO
1904 : END DO
1905 : END IF
1906 3193802 : IF (iatom >= katom) THEN
1907 3106373 : i = 1
1908 3106373 : offset_ac = offset_ac_set(iset, kset) + atomic_offset_ac - 1
1909 3106373 : j = offset_ac
1910 10263996 : DO mc = 1, mc_max
1911 25765115 : DO ma = 1, ma_max
1912 15501119 : pac(i) = density(j)*fac
1913 15501119 : i = i + 1
1914 22658742 : j = j + 1
1915 : END DO
1916 : END DO
1917 : ELSE
1918 87429 : i = 1
1919 87429 : offset_ac = offset_ac_set(kset, iset) + atomic_offset_ac - 1
1920 214036 : DO mc = 1, mc_max
1921 126607 : j = offset_ac + mc - 1
1922 596361 : DO ma = 1, ma_max
1923 382325 : pac(i) = density(j)
1924 382325 : i = i + 1
1925 508932 : j = j + mc_max
1926 : END DO
1927 : END DO
1928 : END IF
1929 3193802 : END SUBROUTINE prefetch_density_matrix
1930 :
1931 : ! **************************************************************************************************
1932 : !> \brief This routine updates the forces using buffered density matrices
1933 : !> \param ma_max Size of matrix blocks
1934 : !> \param mb_max Size of matrix blocks
1935 : !> \param mc_max Size of matrix blocks
1936 : !> \param md_max Size of matrix blocks
1937 : !> \param pbd buffer that will contain P(b,d)
1938 : !> \param pbc buffer that will contain P(b,c)
1939 : !> \param pad buffer that will contain P(a,d)
1940 : !> \param pac buffer that will contain P(a,c)
1941 : !> \param fac mulitplication factor (spin, symmetry)
1942 : !> \param prim primitive forces
1943 : !> \param force storage loacation for forces
1944 : !> \param forces_map index table
1945 : !> \param coord which of the 12 coords to be updated
1946 : !> \param pbd_resp ...
1947 : !> \param pbc_resp ...
1948 : !> \param pad_resp ...
1949 : !> \param pac_resp ...
1950 : !> \param resp_only ...
1951 : !> \par History
1952 : !> 03.2009 created [Manuel Guidon]
1953 : !> \author Manuel Guidon
1954 : ! **************************************************************************************************
1955 25764756 : SUBROUTINE update_forces(ma_max, mb_max, mc_max, md_max, &
1956 : pbd, pbc, pad, pac, fac, &
1957 25764756 : prim, force, forces_map, coord, &
1958 : pbd_resp, pbc_resp, pad_resp, pac_resp, resp_only)
1959 :
1960 : INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
1961 : REAL(dp), DIMENSION(*), INTENT(IN) :: pbd, pbc, pad, pac
1962 : REAL(dp), INTENT(IN) :: fac
1963 : REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
1964 : INTENT(IN) :: prim
1965 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1966 : INTEGER, INTENT(IN) :: forces_map(4, 2), coord
1967 : REAL(dp), DIMENSION(*), INTENT(IN), OPTIONAL :: pbd_resp, pbc_resp, pad_resp, pac_resp
1968 : LOGICAL, INTENT(IN), OPTIONAL :: resp_only
1969 :
1970 : INTEGER :: ma, mb, mc, md, p_index
1971 : LOGICAL :: full_force, with_resp_density
1972 : REAL(dp) :: temp1, temp1_resp, temp3, temp3_resp, &
1973 : temp4, teresp
1974 :
1975 25764756 : with_resp_density = .FALSE.
1976 : IF (PRESENT(pbd_resp) .AND. &
1977 : PRESENT(pbc_resp) .AND. &
1978 25764756 : PRESENT(pad_resp) .AND. &
1979 : PRESENT(pac_resp)) with_resp_density = .TRUE.
1980 :
1981 : IF (with_resp_density) THEN
1982 12534180 : full_force = .TRUE.
1983 12534180 : IF (PRESENT(resp_only)) full_force = .NOT. resp_only
1984 12534180 : p_index = 0
1985 12534180 : temp4 = 0.0_dp
1986 36648168 : DO md = 1, md_max
1987 95527968 : DO mc = 1, mc_max
1988 196567968 : DO mb = 1, mb_max
1989 113574180 : temp1 = pbc((mc - 1)*mb_max + mb)*fac
1990 113574180 : temp3 = pbd((md - 1)*mb_max + mb)*fac
1991 113574180 : temp1_resp = pbc_resp((mc - 1)*mb_max + mb)*fac
1992 113574180 : temp3_resp = pbd_resp((md - 1)*mb_max + mb)*fac
1993 491025648 : DO ma = 1, ma_max
1994 318571668 : p_index = p_index + 1
1995 : ! HF-SCF
1996 318571668 : IF (full_force) THEN
1997 : teresp = temp1*pad((md - 1)*ma_max + ma) + &
1998 102213732 : temp3*pac((mc - 1)*ma_max + ma)
1999 : ELSE
2000 : teresp = 0.0_dp
2001 : END IF
2002 : ! RESP+HF
2003 : teresp = teresp + &
2004 : pac((mc - 1)*ma_max + ma)*temp3_resp + &
2005 : pac_resp((mc - 1)*ma_max + ma)*temp3 + &
2006 : pad((md - 1)*ma_max + ma)*temp1_resp + &
2007 318571668 : pad_resp((md - 1)*ma_max + ma)*temp1
2008 432145848 : temp4 = temp4 + teresp*prim(p_index)
2009 : END DO !ma
2010 : END DO !mb
2011 : END DO !mc
2012 : END DO !md
2013 : ELSE
2014 : p_index = 0
2015 : temp4 = 0.0_dp
2016 39506952 : DO md = 1, md_max
2017 110089224 : DO mc = 1, mc_max
2018 237218544 : DO mb = 1, mb_max
2019 140359896 : temp1 = pbc((mc - 1)*mb_max + mb)*fac
2020 140359896 : temp3 = pbd((md - 1)*mb_max + mb)*fac
2021 679751700 : DO ma = 1, ma_max
2022 468809532 : p_index = p_index + 1
2023 : teresp = temp1*pad((md - 1)*ma_max + ma) + &
2024 468809532 : temp3*pac((mc - 1)*ma_max + ma)
2025 609169428 : temp4 = temp4 + teresp*prim(p_index)
2026 : END DO !ma
2027 : END DO !mb
2028 : END DO !mc
2029 : END DO !md
2030 : END IF
2031 :
2032 25764756 : !$OMP ATOMIC
2033 : force(forces_map((coord - 1)/3 + 1, 1))%fock_4c(MOD(coord - 1, 3) + 1, &
2034 : forces_map((coord - 1)/3 + 1, 2)) = &
2035 : force(forces_map((coord - 1)/3 + 1, 1))%fock_4c(MOD(coord - 1, 3) + 1, &
2036 : forces_map((coord - 1)/3 + 1, 2)) - &
2037 : temp4
2038 :
2039 25764756 : END SUBROUTINE update_forces
2040 :
2041 : ! **************************************************************************************************
2042 : !> \brief This routine updates the virial using buffered density matrices
2043 : !> \param ma_max Size of matrix blocks
2044 : !> \param mb_max Size of matrix blocks
2045 : !> \param mc_max Size of matrix blocks
2046 : !> \param md_max Size of matrix blocks
2047 : !> \param pbd buffer that will contain P(b,d)
2048 : !> \param pbc buffer that will contain P(b,c)
2049 : !> \param pad buffer that will contain P(a,d)
2050 : !> \param pac buffer that will contain P(a,c)
2051 : !> \param fac mulitplication factor (spin, symmetry)
2052 : !> \param prim primitive forces
2053 : !> \param tmp_virial ...
2054 : !> \param coord which of the 12 coords to be updated
2055 : !> \param l ...
2056 : !> \param pbd_resp ...
2057 : !> \param pbc_resp ...
2058 : !> \param pad_resp ...
2059 : !> \param pac_resp ...
2060 : !> \par History
2061 : !> 03.2009 created [Manuel Guidon]
2062 : !> \author Manuel Guidon
2063 : ! **************************************************************************************************
2064 6570468 : SUBROUTINE update_virial(ma_max, mb_max, mc_max, md_max, &
2065 : pbd, pbc, pad, pac, fac, &
2066 6570468 : prim, tmp_virial, coord, l, &
2067 : pbd_resp, pbc_resp, pad_resp, pac_resp, resp_only)
2068 :
2069 : INTEGER, INTENT(IN) :: ma_max, mb_max, mc_max, md_max
2070 : REAL(dp), DIMENSION(*), INTENT(IN) :: pbd, pbc, pad, pac
2071 : REAL(dp), INTENT(IN) :: fac
2072 : REAL(dp), DIMENSION(ma_max*mb_max*mc_max*md_max), &
2073 : INTENT(IN) :: prim
2074 : REAL(dp) :: tmp_virial(3, 3)
2075 : INTEGER, INTENT(IN) :: coord, l
2076 : REAL(dp), DIMENSION(*), INTENT(IN), OPTIONAL :: pbd_resp, pbc_resp, pad_resp, pac_resp
2077 : LOGICAL, INTENT(IN), OPTIONAL :: resp_only
2078 :
2079 : INTEGER :: i, j, ma, mb, mc, md, p_index
2080 : LOGICAL :: with_resp_density, full_force
2081 : REAL(dp) :: temp1, temp1_resp, temp3, temp3_resp, &
2082 : temp4, teresp
2083 :
2084 6570468 : with_resp_density = .FALSE.
2085 : IF (PRESENT(pbd_resp) .AND. &
2086 : PRESENT(pbc_resp) .AND. &
2087 6570468 : PRESENT(pad_resp) .AND. &
2088 : PRESENT(pac_resp)) with_resp_density = .TRUE.
2089 :
2090 : IF (with_resp_density) THEN
2091 5587056 : full_force = .TRUE.
2092 5587056 : IF (PRESENT(resp_only)) full_force = .NOT. resp_only
2093 5587056 : p_index = 0
2094 5587056 : temp4 = 0.0_dp
2095 16806024 : DO md = 1, md_max
2096 41978016 : DO mc = 1, mc_max
2097 91229400 : DO mb = 1, mb_max
2098 54838440 : temp1 = pbc((mc - 1)*mb_max + mb)*fac
2099 54838440 : temp3 = pbd((md - 1)*mb_max + mb)*fac
2100 54838440 : temp1_resp = pbc_resp((mc - 1)*mb_max + mb)*fac
2101 54838440 : temp3_resp = pbd_resp((md - 1)*mb_max + mb)*fac
2102 236614860 : DO ma = 1, ma_max
2103 156604428 : p_index = p_index + 1
2104 : ! HF-SCF
2105 156604428 : IF (full_force) THEN
2106 : teresp = temp1*pad((md - 1)*ma_max + ma) + &
2107 25168824 : temp3*pac((mc - 1)*ma_max + ma)
2108 : ELSE
2109 : teresp = 0.0_dp
2110 : END IF
2111 : ! RESP+HF
2112 : teresp = teresp + &
2113 : pac((mc - 1)*ma_max + ma)*temp3_resp + &
2114 : pac_resp((mc - 1)*ma_max + ma)*temp3 + &
2115 : pad((md - 1)*ma_max + ma)*temp1_resp + &
2116 156604428 : pad_resp((md - 1)*ma_max + ma)*temp1
2117 211442868 : temp4 = temp4 + teresp*prim(p_index)
2118 : END DO !ma
2119 : END DO !mb
2120 : END DO !mc
2121 : END DO !md
2122 : ELSE
2123 : p_index = 0
2124 : temp4 = 0.0_dp
2125 2505240 : DO md = 1, md_max
2126 5427684 : DO mc = 1, mc_max
2127 8390340 : DO mb = 1, mb_max
2128 3946068 : temp1 = pbc((mc - 1)*mb_max + mb)*fac
2129 3946068 : temp3 = pbd((md - 1)*mb_max + mb)*fac
2130 13873176 : DO ma = 1, ma_max
2131 7004664 : p_index = p_index + 1
2132 : teresp = temp1*pad((md - 1)*ma_max + ma) + &
2133 7004664 : temp3*pac((mc - 1)*ma_max + ma)
2134 10950732 : temp4 = temp4 + teresp*prim(p_index)
2135 : END DO !ma
2136 : END DO !mb
2137 : END DO !mc
2138 : END DO !md
2139 : END IF
2140 :
2141 6570468 : j = l
2142 6570468 : i = MOD(coord - 1, 3) + 1
2143 6570468 : tmp_virial(i, j) = tmp_virial(i, j) - temp4
2144 6570468 : END SUBROUTINE update_virial
2145 :
2146 : #:include 'hfx_get_pmax_val.fypp'
2147 :
2148 : END MODULE hfx_derivatives
|
