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 Types and set/get functions for HFX
10 : !> \par History
11 : !> 04.2008 created [Manuel Guidon]
12 : !> 05.2019 Moved erfc_cutoff to common/mathlib (A. Bussy)
13 : !> 10.2025 Added gcc from basis_parameter and hfx_library option
14 : !> \author Manuel Guidon
15 : ! **************************************************************************************************
16 : MODULE hfx_types
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind,&
19 : get_atomic_kind_set
20 : USE basis_set_types, ONLY: get_gto_basis_set,&
21 : gto_basis_set_p_type,&
22 : gto_basis_set_type
23 : USE bibliography, ONLY: bussy2023,&
24 : cite_reference,&
25 : guidon2008,&
26 : guidon2009
27 : USE cell_types, ONLY: cell_type,&
28 : get_cell,&
29 : plane_distance,&
30 : scaled_to_real
31 : USE cp_array_utils, ONLY: cp_1d_logical_p_type
32 : USE cp_control_types, ONLY: dft_control_type
33 : USE cp_dbcsr_api, ONLY: dbcsr_release,&
34 : dbcsr_type
35 : USE cp_files, ONLY: close_file,&
36 : file_exists,&
37 : open_file
38 : USE cp_log_handling, ONLY: cp_get_default_logger,&
39 : cp_logger_type
40 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
41 : cp_print_key_unit_nr
42 : USE cp_units, ONLY: cp_unit_from_cp2k
43 : USE dbt_api, ONLY: &
44 : dbt_create, dbt_default_distvec, dbt_destroy, dbt_distribution_destroy, &
45 : dbt_distribution_new, dbt_distribution_type, dbt_mp_dims_create, dbt_pgrid_create, &
46 : dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
47 : USE hfx_helpers, ONLY: count_cells_perd,&
48 : next_image_cell_perd
49 : USE input_constants, ONLY: &
50 : do_hfx_auto_shells, do_potential_coulomb, do_potential_gaussian, do_potential_id, &
51 : do_potential_long, do_potential_mix_cl, do_potential_mix_cl_trunc, do_potential_mix_lg, &
52 : do_potential_short, do_potential_truncated, hfx_ri_do_2c_diag, hfx_ri_do_2c_iter
53 : USE input_cp2k_hfx, ONLY: ri_mo,&
54 : ri_pmat
55 : USE input_section_types, ONLY: section_vals_get,&
56 : section_vals_get_subs_vals,&
57 : section_vals_type,&
58 : section_vals_val_get
59 : USE kinds, ONLY: default_path_length,&
60 : default_string_length,&
61 : dp,&
62 : int_8
63 : USE libint_2c_3c, ONLY: compare_potential_types,&
64 : libint_potential_type
65 : USE libint_wrapper, ONLY: &
66 : cp_libint_cleanup_eri, cp_libint_cleanup_eri1, cp_libint_init_eri, cp_libint_init_eri1, &
67 : cp_libint_set_contrdepth, cp_libint_static_cleanup, cp_libint_static_init, cp_libint_t, &
68 : prim_data_f_size
69 : USE machine, ONLY: m_chdir,&
70 : m_getcwd
71 : USE mathlib, ONLY: erfc_cutoff
72 : USE message_passing, ONLY: mp_cart_type,&
73 : mp_para_env_type
74 : USE orbital_pointers, ONLY: nco,&
75 : ncoset,&
76 : nso
77 : USE particle_methods, ONLY: get_particle_set
78 : USE particle_types, ONLY: particle_type
79 : USE physcon, ONLY: a_bohr
80 : USE qs_integral_utils, ONLY: basis_set_list_setup
81 : USE qs_kind_types, ONLY: get_qs_kind,&
82 : get_qs_kind_set,&
83 : qs_kind_type
84 : USE qs_tensors_types, ONLY: &
85 : create_2c_tensor, create_3c_tensor, create_tensor_batches, default_block_size, &
86 : distribution_3d_create, distribution_3d_destroy, distribution_3d_type, pgf_block_sizes, &
87 : split_block_sizes
88 : USE string_utilities, ONLY: compress
89 : USE t_c_g0, ONLY: free_C0
90 :
91 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
92 :
93 : #include "./base/base_uses.f90"
94 :
95 : IMPLICIT NONE
96 : PRIVATE
97 : PUBLIC :: hfx_type, hfx_create, hfx_release, &
98 : hfx_set_distr_energy, &
99 : hfx_set_distr_forces, &
100 : hfx_cell_type, hfx_distribution, &
101 : hfx_potential_type, hfx_screening_type, &
102 : hfx_memory_type, hfx_load_balance_type, hfx_general_type, &
103 : hfx_container_type, hfx_cache_type, &
104 : hfx_basis_type, parse_memory_section, &
105 : hfx_init_container, &
106 : hfx_basis_info_type, hfx_screen_coeff_type, &
107 : hfx_reset_memory_usage_counter, pair_list_type, pair_list_element_type, &
108 : pair_set_list_type, hfx_p_kind, hfx_2D_map, hfx_pgf_list, &
109 : hfx_pgf_product_list, hfx_block_range_type, &
110 : alloc_containers, dealloc_containers, hfx_task_list_type, init_t_c_g0_lmax, &
111 : hfx_create_neighbor_cells, hfx_create_basis_types, hfx_release_basis_types, &
112 : hfx_ri_type, hfx_compression_type, block_ind_type, hfx_ri_init, hfx_ri_release, &
113 : compare_hfx_sections
114 :
115 : #define CACHE_SIZE 1024
116 : #define BITS_MAX_VAL 6
117 :
118 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_types'
119 : INTEGER, PARAMETER, PUBLIC :: max_atom_block = 32
120 : INTEGER, PARAMETER, PUBLIC :: max_images = 27
121 : REAL(dp), PARAMETER, PUBLIC :: log_zero = -1000.0_dp
122 : REAL(dp), PARAMETER, PUBLIC :: powell_min_log = -20.0_dp
123 : REAL(KIND=dp), DIMENSION(0:10), &
124 : PARAMETER, PUBLIC :: mul_fact = [1.0_dp, &
125 : 1.1781_dp, &
126 : 1.3333_dp, &
127 : 1.4726_dp, &
128 : 1.6000_dp, &
129 : 1.7181_dp, &
130 : 1.8286_dp, &
131 : 1.9328_dp, &
132 : 2.0317_dp, &
133 : 2.1261_dp, &
134 : 2.2165_dp]
135 :
136 : INTEGER, SAVE :: init_t_c_g0_lmax = -1
137 :
138 : !***
139 :
140 : ! **************************************************************************************************
141 : TYPE hfx_potential_type
142 : INTEGER :: potential_type = do_potential_coulomb !! 1/r/ erfc(wr)/r ...
143 : REAL(dp) :: omega = 0.0_dp !! w
144 : REAL(dp) :: scale_coulomb = 0.0_dp !! scaling factor for mixed potential
145 : REAL(dp) :: scale_longrange = 0.0_dp !! scaling factor for mixed potential
146 : REAL(dp) :: scale_gaussian = 0.0_dp!! scaling factor for mixed potential
147 : REAL(dp) :: cutoff_radius = 0.0_dp!! cutoff radius if cutoff potential in use
148 : CHARACTER(default_path_length) :: filename = ""
149 : END TYPE hfx_potential_type
150 :
151 : ! **************************************************************************************************
152 : TYPE hfx_screening_type
153 : REAL(dp) :: eps_schwarz = 0.0_dp !! threshold
154 : REAL(dp) :: eps_schwarz_forces = 0.0_dp !! threshold
155 : LOGICAL :: do_p_screening_forces = .FALSE. !! screen on P^2 ?
156 : LOGICAL :: do_initial_p_screening = .FALSE. !! screen on initial guess?
157 : END TYPE hfx_screening_type
158 :
159 : ! **************************************************************************************************
160 : TYPE hfx_memory_type
161 : INTEGER :: max_memory = 0 !! user def max memory MiB
162 : INTEGER(int_8) :: max_compression_counter = 0_int_8 !! corresponding number of reals
163 : INTEGER(int_8) :: final_comp_counter_energy = 0_int_8
164 : LOGICAL :: do_all_on_the_fly = .FALSE. !! max mem == 0 ?
165 : REAL(dp) :: eps_storage_scaling = 0.0_dp
166 : INTEGER :: cache_size = 0
167 : INTEGER :: bits_max_val = 0
168 : INTEGER :: actual_memory_usage = 0
169 : INTEGER :: actual_memory_usage_disk = 0
170 : INTEGER(int_8) :: max_compression_counter_disk = 0_int_8
171 : LOGICAL :: do_disk_storage = .FALSE.
172 : CHARACTER(len=default_path_length) :: storage_location = ""
173 : INTEGER(int_8) :: ram_counter = 0_int_8
174 : INTEGER(int_8) :: ram_counter_forces = 0_int_8
175 : INTEGER(int_8) :: size_p_screen = 0_int_8
176 : LOGICAL :: treat_forces_in_core = .FALSE.
177 : LOGICAL :: recalc_forces = .FALSE.
178 : END TYPE hfx_memory_type
179 :
180 : ! **************************************************************************************************
181 : TYPE hfx_periodic_type
182 : INTEGER :: number_of_shells = -1 !! number of periodic image cells
183 : LOGICAL :: do_periodic = .FALSE. !! periodic ?
184 : INTEGER :: perd(3) = -1 !! x,xy,xyz,...
185 : INTEGER :: mode = -1
186 : REAL(dp) :: R_max_stress = 0.0_dp
187 : INTEGER :: number_of_shells_from_input = 0
188 : END TYPE hfx_periodic_type
189 :
190 : ! **************************************************************************************************
191 : TYPE hfx_load_balance_type
192 : INTEGER :: nbins = 0
193 : INTEGER :: block_size = 0
194 : INTEGER :: nblocks = 0
195 : LOGICAL :: rtp_redistribute = .FALSE.
196 : LOGICAL :: blocks_initialized = .FALSE.
197 : LOGICAL :: do_randomize = .FALSE.
198 : END TYPE hfx_load_balance_type
199 :
200 : ! **************************************************************************************************
201 : TYPE hfx_general_type
202 : REAL(dp) :: fraction = 0.0_dp !! for hybrids
203 : INTEGER :: hfx_library = 0
204 : LOGICAL :: treat_lsd_in_core = .FALSE.
205 : END TYPE hfx_general_type
206 :
207 : ! **************************************************************************************************
208 : TYPE hfx_cell_type
209 : REAL(dp) :: cell(3) = 0.0_dp
210 : REAL(dp) :: cell_r(3) = 0.0_dp
211 : END TYPE hfx_cell_type
212 :
213 : ! **************************************************************************************************
214 : TYPE hfx_distribution
215 : INTEGER(int_8) :: istart = 0_int_8
216 : INTEGER(int_8) :: number_of_atom_quartets = 0_int_8
217 : INTEGER(int_8) :: cost = 0_int_8
218 : REAL(KIND=dp) :: time_first_scf = 0.0_dp
219 : REAL(KIND=dp) :: time_other_scf = 0.0_dp
220 : REAL(KIND=dp) :: time_forces = 0.0_dp
221 : INTEGER(int_8) :: ram_counter = 0_int_8
222 : END TYPE hfx_distribution
223 :
224 : ! **************************************************************************************************
225 : TYPE pair_list_element_type
226 : INTEGER, DIMENSION(2) :: pair = 0
227 : INTEGER, DIMENSION(2) :: set_bounds = 0
228 : INTEGER, DIMENSION(2) :: kind_pair = 0
229 : REAL(KIND=dp) :: r1(3) = 0.0_dp, r2(3) = 0.0_dp
230 : REAL(KIND=dp) :: dist2 = 0.0_dp
231 : END TYPE pair_list_element_type
232 :
233 : ! **************************************************************************************************
234 : TYPE pair_set_list_type
235 : INTEGER, DIMENSION(2) :: pair = 0
236 : END TYPE pair_set_list_type
237 :
238 : ! **************************************************************************************************
239 : TYPE pair_list_type
240 : TYPE(pair_list_element_type), DIMENSION(max_atom_block**2) :: elements = pair_list_element_type()
241 : INTEGER :: n_element = 0
242 : END TYPE pair_list_type
243 :
244 : ! **************************************************************************************************
245 : TYPE hfx_cache_type
246 : INTEGER(int_8), DIMENSION(CACHE_SIZE) :: DATA = 0_int_8
247 : INTEGER :: element_counter = 0
248 : END TYPE hfx_cache_type
249 :
250 : ! **************************************************************************************************
251 : TYPE hfx_container_node
252 : TYPE(hfx_container_node), POINTER :: next => NULL(), prev => NULL()
253 : INTEGER(int_8), DIMENSION(CACHE_SIZE) :: DATA = 0_int_8
254 : END TYPE hfx_container_node
255 :
256 : ! **************************************************************************************************
257 : TYPE hfx_container_type
258 : TYPE(hfx_container_node), POINTER :: first => NULL(), current => NULL()
259 : INTEGER :: element_counter = 0
260 : INTEGER(int_8) :: file_counter = 0
261 : CHARACTER(LEN=5) :: desc = ""
262 : INTEGER :: unit = -1
263 : CHARACTER(default_path_length) :: filename = ""
264 : END TYPE hfx_container_type
265 :
266 : ! **************************************************************************************************
267 : TYPE hfx_basis_type
268 : INTEGER, DIMENSION(:), POINTER :: lmax => NULL()
269 : INTEGER, DIMENSION(:), POINTER :: lmin => NULL()
270 : INTEGER, DIMENSION(:), POINTER :: npgf => NULL()
271 : INTEGER :: nset = 0
272 : REAL(dp), DIMENSION(:, :), POINTER :: zet => NULL()
273 : INTEGER, DIMENSION(:), POINTER :: nsgf => NULL()
274 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf => NULL()
275 : REAL(dp), DIMENSION(:, :), POINTER :: sphi => NULL()
276 : INTEGER :: nsgf_total = 0
277 : INTEGER, DIMENSION(:, :), POINTER :: nl => NULL()
278 : INTEGER, DIMENSION(:, :), POINTER :: nsgfl => NULL()
279 : INTEGER, DIMENSION(:), POINTER :: nshell => NULL()
280 : REAL(dp), DIMENSION(:, :, :, :), POINTER &
281 : :: sphi_ext => NULL()
282 : REAL(dp), DIMENSION(:, :, :), POINTER :: gcc => NULL()
283 : REAL(dp), DIMENSION(:), POINTER :: set_radius => NULL()
284 : REAL(dp), DIMENSION(:, :), POINTER :: pgf_radius => NULL()
285 : REAL(dp) :: kind_radius = 0.0_dp
286 : END TYPE hfx_basis_type
287 :
288 : ! **************************************************************************************************
289 : TYPE hfx_basis_info_type
290 : INTEGER :: max_set = 0
291 : INTEGER :: max_sgf = 0
292 : INTEGER :: max_am = 0
293 : END TYPE hfx_basis_info_type
294 :
295 : ! **************************************************************************************************
296 : TYPE hfx_screen_coeff_type
297 : REAL(dp) :: x(2) = 0.0_dp
298 : END TYPE hfx_screen_coeff_type
299 :
300 : ! **************************************************************************************************
301 : TYPE hfx_p_kind
302 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: p_kind => NULL()
303 : END TYPE hfx_p_kind
304 :
305 : ! **************************************************************************************************
306 : TYPE hfx_2D_map
307 : INTEGER, DIMENSION(:), POINTER :: iatom_list => NULL()
308 : INTEGER, DIMENSION(:), POINTER :: jatom_list => NULL()
309 : END TYPE hfx_2D_map
310 :
311 : ! **************************************************************************************************
312 : TYPE hfx_pgf_image
313 : REAL(dp) :: ra(3) = 0.0_dp, rb(3) = 0.0_dp
314 : REAL(dp) :: rab2 = 0.0_dp
315 : REAL(dp) :: S1234 = 0.0_dp
316 : REAL(dp) :: P(3) = 0.0_dp
317 : REAL(dp) :: R = 0.0_dp
318 : REAL(dp) :: pgf_max = 0.0_dp
319 : REAL(dp), DIMENSION(3) :: bcell = 0.0_dp
320 : END TYPE hfx_pgf_image
321 :
322 : ! **************************************************************************************************
323 : TYPE hfx_pgf_list
324 : TYPE(hfx_pgf_image), DIMENSION(:), POINTER &
325 : :: image_list => NULL()
326 : INTEGER :: nimages = 0
327 : REAL(dp) :: zetapzetb = 0.0_dp
328 : REAL(dp) :: ZetaInv = 0.0_dp
329 : REAL(dp) :: zeta = 0.0_dp, zetb = 0.0_dp
330 : INTEGER :: ipgf = 0, jpgf = 0
331 : END TYPE hfx_pgf_list
332 :
333 : ! **************************************************************************************************
334 : TYPE hfx_pgf_product_list
335 : REAL(dp) :: ra(3) = 0.0_dp, rb(3) = 0.0_dp, rc(3) = 0.0_dp, rd(3) = 0.0_dp
336 : REAL(dp) :: ZetapEtaInv = 0.0_dp
337 : REAL(dp) :: Rho = 0.0_dp, RhoInv = 0.0_dp
338 : REAL(dp) :: P(3) = 0.0_dp, Q(3) = 0.0_dp, W(3) = 0.0_dp
339 : REAL(dp) :: AB(3) = 0.0_dp, CD(3) = 0.0_dp
340 : REAL(dp) :: Fm(prim_data_f_size) = 0.0_dp
341 : END TYPE hfx_pgf_product_list
342 :
343 : ! **************************************************************************************************
344 : TYPE hfx_block_range_type
345 : INTEGER :: istart = 0, iend = 0
346 : INTEGER(int_8) :: cost = 0_int_8
347 : END TYPE hfx_block_range_type
348 :
349 : ! **************************************************************************************************
350 : TYPE hfx_task_list_type
351 : INTEGER :: thread_id = 0
352 : INTEGER :: bin_id = 0
353 : INTEGER(int_8) :: cost = 0_int_8
354 : END TYPE hfx_task_list_type
355 :
356 : TYPE :: hfx_compression_type
357 : TYPE(hfx_container_type), DIMENSION(:), &
358 : POINTER :: maxval_container => NULL()
359 : TYPE(hfx_cache_type), DIMENSION(:), &
360 : POINTER :: maxval_cache => NULL()
361 : TYPE(hfx_container_type), DIMENSION(:, :), &
362 : POINTER :: integral_containers => NULL()
363 : TYPE(hfx_cache_type), DIMENSION(:, :), &
364 : POINTER :: integral_caches => NULL()
365 : TYPE(hfx_container_type), POINTER :: maxval_container_disk => NULL()
366 : TYPE(hfx_cache_type) :: maxval_cache_disk = hfx_cache_type()
367 : TYPE(hfx_cache_type) :: integral_caches_disk(64) = hfx_cache_type()
368 : TYPE(hfx_container_type), POINTER, &
369 : DIMENSION(:) :: integral_containers_disk => NULL()
370 : END TYPE hfx_compression_type
371 :
372 : TYPE :: block_ind_type
373 : INTEGER, DIMENSION(:, :), ALLOCATABLE :: ind
374 : END TYPE block_ind_type
375 :
376 : TYPE hfx_ri_type
377 : ! input parameters (see input_cp2k_hfx)
378 : REAL(KIND=dp) :: filter_eps = 0.0_dp, filter_eps_2c = 0.0_dp, filter_eps_storage = 0.0_dp, filter_eps_mo = 0.0_dp, &
379 : eps_lanczos = 0.0_dp, eps_pgf_orb = 0.0_dp, eps_eigval = 0.0_dp, kp_RI_range = 0.0_dp, &
380 : kp_image_range = 0.0_dp, kp_bump_rad = 0.0_dp
381 : INTEGER :: t2c_sqrt_order = 0, max_iter_lanczos = 0, flavor = 0, unit_nr_dbcsr = -1, unit_nr = -1, &
382 : min_bsize = 0, max_bsize_MO = 0, t2c_method = 0, nelectron_total = 0, input_flavor = 0, &
383 : ncell_RI = 0, nimg = 0, kp_stack_size = 0, nimg_nze = 0, kp_ngroups = 1
384 : LOGICAL :: check_2c_inv = .FALSE., calc_condnum = .FALSE.
385 :
386 : TYPE(libint_potential_type) :: ri_metric = libint_potential_type()
387 :
388 : ! input parameters from hfx
389 : TYPE(libint_potential_type) :: hfx_pot = libint_potential_type() ! interaction potential
390 : REAL(KIND=dp) :: eps_schwarz = 0.0_dp ! integral screening threshold
391 : REAL(KIND=dp) :: eps_schwarz_forces = 0.0_dp ! integral derivatives screening threshold
392 :
393 : LOGICAL :: same_op = .FALSE. ! whether RI operator is same as HF potential
394 :
395 : ! default process grid used for 3c tensors
396 : TYPE(dbt_pgrid_type), POINTER :: pgrid => NULL()
397 : TYPE(dbt_pgrid_type), POINTER :: pgrid_2d => NULL()
398 :
399 : ! distributions for (RI | AO AO) 3c integral tensor (non split)
400 : TYPE(distribution_3d_type) :: dist_3d = distribution_3d_type()
401 : TYPE(dbt_distribution_type) :: dist
402 :
403 : ! block sizes for RI and AO tensor dimensions (split)
404 : INTEGER, DIMENSION(:), ALLOCATABLE :: bsizes_RI, bsizes_AO, bsizes_RI_split, bsizes_AO_split, &
405 : bsizes_RI_fit, bsizes_AO_fit
406 :
407 : ! KP RI-HFX basis info
408 : INTEGER, DIMENSION(:), ALLOCATABLE :: img_to_RI_cell, present_images, idx_to_img, img_to_idx, &
409 : RI_cell_to_img
410 :
411 : ! KP RI-HFX cost information for a given atom pair i,j at a given cell b
412 : REAL(dp), DIMENSION(:, :, :), ALLOCATABLE :: kp_cost
413 :
414 : ! KP distribution of iatom (of i,j atom pairs) to subgroups
415 : TYPE(cp_1d_logical_p_type), DIMENSION(:), ALLOCATABLE :: iatom_to_subgroup
416 :
417 : ! KP 3c tensors replicated on the subgroups
418 : TYPE(dbt_type), DIMENSION(:), ALLOCATABLE :: kp_t_3c_int
419 :
420 : ! Note: changed static DIMENSION(1,1) of dbt_type to allocatables as workaround for gfortran 8.3.0,
421 : ! with static dimension gfortran gets stuck during compilation
422 :
423 : ! 2c tensors in (AO | AO) format
424 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: rho_ao_t, ks_t
425 :
426 : ! 2c tensors in (RI | RI) format for forces
427 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: t_2c_inv
428 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: t_2c_pot
429 :
430 : ! 2c tensor in matrix format for K-points RI-HFX
431 : TYPE(dbcsr_type), DIMENSION(:, :), ALLOCATABLE :: kp_mat_2c_pot
432 :
433 : ! 2c tensor in (RI | RI) format for contraction
434 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: t_2c_int
435 :
436 : ! 3c integral tensor in (AO RI | AO) format for contraction
437 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: t_3c_int_ctr_1
438 : TYPE(block_ind_type), DIMENSION(:, :), ALLOCATABLE :: blk_indices
439 : TYPE(dbt_pgrid_type), POINTER :: pgrid_1 => NULL()
440 :
441 : ! 3c integral tensor in ( AO | RI AO) (MO) or (AO RI | AO) (RHO) format for contraction
442 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: t_3c_int_ctr_2
443 : TYPE(dbt_pgrid_type), POINTER :: pgrid_2 => NULL()
444 :
445 : ! 3c integral tensor in ( RI | AO AO ) format for contraction
446 : TYPE(dbt_type), DIMENSION(:, :), ALLOCATABLE :: t_3c_int_ctr_3
447 :
448 : ! 3c integral tensor in (RI | MO AO ) format for contraction
449 : TYPE(dbt_type), DIMENSION(:, :, :), ALLOCATABLE :: t_3c_int_mo
450 : TYPE(dbt_type), DIMENSION(:, :, :), ALLOCATABLE :: t_3c_ctr_RI
451 : TYPE(dbt_type), DIMENSION(:, :, :), ALLOCATABLE :: t_3c_ctr_KS
452 : TYPE(dbt_type), DIMENSION(:, :, :), ALLOCATABLE :: t_3c_ctr_KS_copy
453 :
454 : ! optional: sections for output handling
455 : ! alternatively set unit_nr_dbcsr (for logging tensor operations) and unit_nr (for general
456 : ! output) directly
457 : TYPE(section_vals_type), POINTER :: ri_section => NULL(), hfx_section => NULL()
458 :
459 : ! types of primary and auxiliary basis
460 : CHARACTER(len=default_string_length) :: orb_basis_type = "", ri_basis_type = ""
461 :
462 : ! memory reduction factor
463 : INTEGER :: n_mem_input = 0, n_mem = 0, n_mem_RI = 0, n_mem_flavor_switch = 0
464 :
465 : ! offsets for memory batches
466 : INTEGER, DIMENSION(:), ALLOCATABLE :: starts_array_mem_block, ends_array_mem_block
467 : INTEGER, DIMENSION(:), ALLOCATABLE :: starts_array_mem, ends_array_mem
468 :
469 : INTEGER, DIMENSION(:), ALLOCATABLE :: starts_array_RI_mem_block, ends_array_RI_mem_block
470 : INTEGER, DIMENSION(:), ALLOCATABLE :: starts_array_RI_mem, ends_array_RI_mem
471 :
472 : INTEGER(int_8) :: dbcsr_nflop = 0_int_8
473 : REAL(dp) :: dbcsr_time = 0.0_dp
474 : INTEGER :: num_pe = 0
475 : TYPE(hfx_compression_type), DIMENSION(:, :), ALLOCATABLE :: store_3c
476 :
477 : END TYPE hfx_ri_type
478 :
479 : ! **************************************************************************************************
480 : !> \brief stores some data used in construction of Kohn-Sham matrix
481 : !> \param potential_parameter stores information on the potential (1/r, erfc(wr)/r
482 : !> \param screening_parameter stores screening infos such as epsilon
483 : !> \param memory_parameter stores infos on memory used for in-core calculations
484 : !> \param periodic_parameter stores information on how to apply pbc
485 : !> \param load_balance_parameter contains infos for Monte Carlo simulated annealing
486 : !> \param general_paramter at the moment stores the fraction of HF amount to be included
487 : !> \param maxval_container stores the maxvals in compressed form
488 : !> \param maxval_cache cache for maxvals in decompressed form
489 : !> \param integral_containers 64 containers for compressed integrals
490 : !> \param integral_caches 64 caches for decompressed integrals
491 : !> \param neighbor_cells manages handling of periodic cells
492 : !> \param distribution_energy stores information on parallelization of energy
493 : !> \param distribution_forces stores information on parallelization of forces
494 : !> \param initial_p stores the initial guess if requested
495 : !> \param is_assoc_atomic_block reflects KS sparsity
496 : !> \param number_of_p_entries Size of P matrix
497 : !> \param n_rep_hf Number of HFX replicas
498 : !> \param b_first_load_balance_x flag to indicate if it is enough just to update
499 : !> the distribution of the integrals
500 : !> \param full_ks_x full ks matrices
501 : !> \param lib libint type for eris
502 : !> \param basis_info contains information for basis sets
503 : !> \param screen_funct_coeffs_pgf pgf based near field screening coefficients
504 : !> \param pair_dist_radii_pgf pgf based radii coefficients of pair distributions
505 : !> \param screen_funct_coeffs_set set based near field screening coefficients
506 : !> \param screen_funct_coeffs_kind kind based near field screening coefficients
507 : !> \param screen_funct_is_initialized flag that indicates if the coefficients
508 : !> have already been fitted
509 : !> \par History
510 : !> 11.2006 created [Manuel Guidon]
511 : !> 02.2009 completely rewritten due to new screening
512 : !> \author Manuel Guidon
513 : ! **************************************************************************************************
514 : TYPE hfx_type
515 : TYPE(hfx_potential_type) :: potential_parameter = hfx_potential_type()
516 : TYPE(hfx_screening_type) :: screening_parameter = hfx_screening_type()
517 : TYPE(hfx_memory_type) :: memory_parameter = hfx_memory_type()
518 : TYPE(hfx_periodic_type) :: periodic_parameter = hfx_periodic_type()
519 : TYPE(hfx_load_balance_type) :: load_balance_parameter = hfx_load_balance_type()
520 : TYPE(hfx_general_type) :: general_parameter = hfx_general_type()
521 :
522 : TYPE(hfx_compression_type) :: store_ints = hfx_compression_type()
523 : TYPE(hfx_compression_type) :: store_forces = hfx_compression_type()
524 :
525 : TYPE(hfx_cell_type), DIMENSION(:), &
526 : POINTER :: neighbor_cells => NULL()
527 : TYPE(hfx_distribution), DIMENSION(:), &
528 : POINTER :: distribution_energy => NULL()
529 : TYPE(hfx_distribution), DIMENSION(:), &
530 : POINTER :: distribution_forces => NULL()
531 : INTEGER, DIMENSION(:, :), POINTER :: is_assoc_atomic_block => NULL()
532 : INTEGER :: number_of_p_entries = 0
533 : TYPE(hfx_basis_type), DIMENSION(:), &
534 : POINTER :: basis_parameter => NULL()
535 : INTEGER :: n_rep_hf = 0
536 : LOGICAL :: b_first_load_balance_energy = .FALSE., &
537 : b_first_load_balance_forces = .FALSE.
538 : REAL(dp), DIMENSION(:, :), POINTER :: full_ks_alpha => NULL()
539 : REAL(dp), DIMENSION(:, :), POINTER :: full_ks_beta => NULL()
540 : TYPE(cp_libint_t) :: lib
541 : TYPE(hfx_basis_info_type) :: basis_info = hfx_basis_info_type()
542 : TYPE(hfx_screen_coeff_type), &
543 : DIMENSION(:, :, :, :, :, :), POINTER :: screen_funct_coeffs_pgf => NULL(), &
544 : pair_dist_radii_pgf => NULL()
545 : TYPE(hfx_screen_coeff_type), &
546 : DIMENSION(:, :, :, :), POINTER :: screen_funct_coeffs_set => NULL()
547 : TYPE(hfx_screen_coeff_type), &
548 : DIMENSION(:, :), POINTER :: screen_funct_coeffs_kind => NULL()
549 : LOGICAL :: screen_funct_is_initialized = .FALSE.
550 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: initial_p => NULL()
551 : TYPE(hfx_p_kind), DIMENSION(:), POINTER :: initial_p_forces => NULL()
552 : INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom => NULL()
553 : TYPE(hfx_2D_map), DIMENSION(:), POINTER :: map_atoms_to_cpus => NULL()
554 : INTEGER, DIMENSION(:, :), POINTER :: atomic_block_offset => NULL()
555 : INTEGER, DIMENSION(:, :, :, :), POINTER :: set_offset => NULL()
556 : INTEGER, DIMENSION(:), POINTER :: block_offset => NULL()
557 : TYPE(hfx_block_range_type), DIMENSION(:), &
558 : POINTER :: blocks => NULL()
559 : TYPE(hfx_task_list_type), DIMENSION(:), &
560 : POINTER :: task_list => NULL()
561 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom => NULL(), pmax_atom_forces => NULL()
562 : TYPE(cp_libint_t) :: lib_deriv
563 : REAL(dp), DIMENSION(:, :), POINTER :: pmax_block => NULL()
564 : LOGICAL, DIMENSION(:, :), POINTER :: atomic_pair_list => NULL()
565 : LOGICAL, DIMENSION(:, :), POINTER :: atomic_pair_list_forces => NULL()
566 : LOGICAL :: do_hfx_ri = .FALSE.
567 : TYPE(hfx_ri_type), POINTER :: ri_data => NULL()
568 :
569 : ! ACE fields
570 : LOGICAL :: use_ace = .FALSE.
571 : INTEGER :: ace_rebuild_freq = 20
572 :
573 : ! ACE pprojectors will be declared in hfx_admm_utils.F
574 : LOGICAL :: ace_is_built = .FALSE.
575 : INTEGER :: ace_build_counter = 0
576 : END TYPE hfx_type
577 :
578 : CONTAINS
579 :
580 : ! **************************************************************************************************
581 : !> \brief - This routine allocates and initializes all types in hfx_data
582 : !> \param x_data contains all relevant data structures for hfx runs
583 : !> \param para_env ...
584 : !> \param hfx_section input section
585 : !> \param atomic_kind_set ...
586 : !> \param qs_kind_set ...
587 : !> \param particle_set ...
588 : !> \param dft_control ...
589 : !> \param cell ...
590 : !> \param orb_basis ...
591 : !> \param ri_basis ...
592 : !> \param nelectron_total ...
593 : !> \param nkp_grid ...
594 : !> \par History
595 : !> 09.2007 created [Manuel Guidon]
596 : !> 01.2024 pushed basis set decision outside of routine, keeps default as
597 : !> orb_basis = "ORB" and ri_basis = "AUX_FIT"
598 : !> No more ADMM references!
599 : !> \author Manuel Guidon
600 : !> \note
601 : !> - All POINTERS and ALLOCATABLES are allocated, even if their size is
602 : !> unknown at invocation time
603 : ! **************************************************************************************************
604 1478 : SUBROUTINE hfx_create(x_data, para_env, hfx_section, atomic_kind_set, qs_kind_set, &
605 : particle_set, dft_control, cell, orb_basis, ri_basis, &
606 : nelectron_total, nkp_grid)
607 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
608 : TYPE(mp_para_env_type) :: para_env
609 : TYPE(section_vals_type), POINTER :: hfx_section
610 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
611 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
612 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
613 : TYPE(dft_control_type), POINTER :: dft_control
614 : TYPE(cell_type), POINTER :: cell
615 : CHARACTER(LEN=*), OPTIONAL :: orb_basis, ri_basis
616 : INTEGER, OPTIONAL :: nelectron_total
617 : INTEGER, DIMENSION(3), OPTIONAL :: nkp_grid
618 :
619 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create'
620 :
621 : CHARACTER(LEN=512) :: error_msg
622 : CHARACTER(LEN=default_path_length) :: char_val
623 : CHARACTER(LEN=default_string_length) :: orb_basis_type, ri_basis_type
624 : INTEGER :: handle, i, i_thread, iatom, ikind, int_val, irep, jkind, max_set, n_rep_hf, &
625 : n_threads, natom, natom_a, natom_b, nkind, nseta, nsetb, pbc_shells, storage_id
626 1478 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2kind, kind_of
627 : LOGICAL :: do_ri, explicit, logic_val
628 : REAL(dp) :: real_val
629 : TYPE(hfx_type), POINTER :: actual_x_data
630 : TYPE(section_vals_type), POINTER :: hf_pbc_section, hf_sub_section, &
631 : hfx_ri_section
632 :
633 1478 : CALL timeset(routineN, handle)
634 :
635 1478 : CALL cite_reference(Guidon2008)
636 1478 : CALL cite_reference(Guidon2009)
637 :
638 1478 : natom = SIZE(particle_set)
639 :
640 : !! There might be 2 hf sections
641 1478 : CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
642 1478 : n_threads = 1
643 1478 : !$ n_threads = omp_get_max_threads()
644 :
645 1478 : CALL section_vals_val_get(hfx_section, "RI%_SECTION_PARAMETERS_", l_val=do_ri)
646 1478 : IF (do_ri) n_threads = 1 ! RI implementation does not use threads
647 :
648 1478 : IF (PRESENT(orb_basis)) THEN
649 1478 : orb_basis_type = orb_basis
650 : ELSE
651 0 : orb_basis_type = "ORB"
652 : END IF
653 1478 : IF (PRESENT(ri_basis)) THEN
654 0 : ri_basis_type = ri_basis
655 : ELSE
656 1478 : ri_basis_type = "RI_HFX"
657 : END IF
658 :
659 6257862 : ALLOCATE (x_data(n_rep_hf, n_threads))
660 2956 : DO i_thread = 1, n_threads
661 4444 : DO irep = 1, n_rep_hf
662 1488 : actual_x_data => x_data(irep, i_thread)
663 : !! Get data from input file
664 : !!
665 : !! GENERAL params
666 1488 : CALL section_vals_val_get(hfx_section, "FRACTION", r_val=real_val, i_rep_section=irep)
667 1488 : actual_x_data%general_parameter%fraction = real_val
668 1488 : actual_x_data%n_rep_hf = n_rep_hf
669 :
670 1488 : NULLIFY (actual_x_data%map_atoms_to_cpus)
671 :
672 1488 : CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=logic_val, i_rep_section=irep)
673 1488 : actual_x_data%general_parameter%treat_lsd_in_core = logic_val
674 :
675 1488 : CALL section_vals_val_get(hfx_section, "HFX_LIBRARY", i_val=int_val, i_rep_section=irep)
676 1488 : actual_x_data%general_parameter%hfx_library = int_val
677 :
678 1488 : hfx_ri_section => section_vals_get_subs_vals(hfx_section, "RI")
679 1488 : CALL section_vals_val_get(hfx_ri_section, "_SECTION_PARAMETERS_", l_val=actual_x_data%do_hfx_ri)
680 :
681 : !! MEMORY section
682 1488 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "MEMORY", i_rep_section=irep)
683 : CALL parse_memory_section(actual_x_data%memory_parameter, hf_sub_section, storage_id, i_thread, &
684 1488 : n_threads, para_env, irep, skip_disk=.FALSE., skip_in_core_forces=.FALSE.)
685 :
686 : !! PERIODIC section
687 1488 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "PERIODIC", i_rep_section=irep)
688 1488 : CALL section_vals_val_get(hf_sub_section, "NUMBER_OF_SHELLS", i_val=int_val)
689 1488 : actual_x_data%periodic_parameter%number_of_shells = int_val
690 1488 : actual_x_data%periodic_parameter%mode = int_val
691 1488 : CALL get_cell(cell=cell, periodic=actual_x_data%periodic_parameter%perd)
692 5952 : IF (SUM(actual_x_data%periodic_parameter%perd) == 0) THEN
693 1040 : actual_x_data%periodic_parameter%do_periodic = .FALSE.
694 : ELSE
695 448 : actual_x_data%periodic_parameter%do_periodic = .TRUE.
696 : END IF
697 :
698 : !! SCREENING section
699 1488 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "SCREENING", i_rep_section=irep)
700 1488 : CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ", r_val=real_val)
701 1488 : actual_x_data%screening_parameter%eps_schwarz = real_val
702 1488 : CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ_FORCES", r_val=real_val, explicit=explicit)
703 1488 : IF (explicit) THEN
704 204 : actual_x_data%screening_parameter%eps_schwarz_forces = real_val
705 : ELSE
706 : actual_x_data%screening_parameter%eps_schwarz_forces = &
707 1284 : 100._dp*actual_x_data%screening_parameter%eps_schwarz
708 : END IF
709 1488 : CALL section_vals_val_get(hf_sub_section, "SCREEN_P_FORCES", l_val=logic_val)
710 1488 : actual_x_data%screening_parameter%do_p_screening_forces = logic_val
711 1488 : CALL section_vals_val_get(hf_sub_section, "SCREEN_ON_INITIAL_P", l_val=logic_val)
712 1488 : actual_x_data%screening_parameter%do_initial_p_screening = logic_val
713 1488 : actual_x_data%screen_funct_is_initialized = .FALSE.
714 :
715 : !! INTERACTION_POTENTIAL section
716 1488 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL", i_rep_section=irep)
717 1488 : CALL section_vals_val_get(hf_sub_section, "POTENTIAL_TYPE", i_val=int_val)
718 1488 : actual_x_data%potential_parameter%potential_type = int_val
719 1488 : CALL section_vals_val_get(hf_sub_section, "OMEGA", r_val=real_val)
720 1488 : actual_x_data%potential_parameter%omega = real_val
721 1488 : CALL section_vals_val_get(hf_sub_section, "SCALE_COULOMB", r_val=real_val)
722 1488 : actual_x_data%potential_parameter%scale_coulomb = real_val
723 1488 : CALL section_vals_val_get(hf_sub_section, "SCALE_LONGRANGE", r_val=real_val)
724 1488 : actual_x_data%potential_parameter%scale_longrange = real_val
725 1488 : CALL section_vals_val_get(hf_sub_section, "SCALE_GAUSSIAN", r_val=real_val)
726 1488 : actual_x_data%potential_parameter%scale_gaussian = real_val
727 1488 : IF (actual_x_data%potential_parameter%potential_type == do_potential_truncated .OR. &
728 : actual_x_data%potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
729 376 : CALL section_vals_val_get(hf_sub_section, "CUTOFF_RADIUS", r_val=real_val)
730 376 : actual_x_data%potential_parameter%cutoff_radius = real_val
731 376 : CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val)
732 376 : CALL compress(char_val, .TRUE.)
733 : ! ** Check if file is there
734 376 : IF (.NOT. file_exists(char_val)) THEN
735 : WRITE (error_msg, '(A,A,A)') "Truncated hfx calculation requested. The file containing "// &
736 0 : "the data could not be found at ", TRIM(char_val), " Please check T_C_G_DATA "// &
737 0 : "in the INTERACTION_POTENTIAL section"
738 0 : CPABORT(error_msg)
739 : ELSE
740 376 : actual_x_data%potential_parameter%filename = char_val
741 : END IF
742 : END IF
743 1488 : IF (actual_x_data%potential_parameter%potential_type == do_potential_short) THEN
744 : CALL erfc_cutoff(actual_x_data%screening_parameter%eps_schwarz, &
745 : actual_x_data%potential_parameter%omega, &
746 48 : actual_x_data%potential_parameter%cutoff_radius)
747 48 : CALL section_vals_val_get(hf_sub_section, "CUTOFF_RADIUS", explicit=explicit)
748 48 : IF (explicit) THEN
749 0 : CALL section_vals_val_get(hf_sub_section, "CUTOFF_RADIUS", r_val=real_val)
750 : IF (real_val < actual_x_data%potential_parameter%cutoff_radius .AND. &
751 0 : i_thread == 1 .AND. irep == 1) THEN
752 : WRITE (error_msg, '(A,F6.3,A,ES8.1,A,F6.3,A,F6.3,A)') &
753 : "Periodic Hartree Fock calculation requested with the use "// &
754 0 : "of a shortrange potential erfc(omega*r)/r. Given omega = ", &
755 0 : actual_x_data%potential_parameter%omega, " and EPS_SCHWARZ = ", &
756 0 : actual_x_data%screening_parameter%eps_schwarz, ", the requested "// &
757 0 : "cutoff radius ", real_val*a_bohr*1e+10_dp, " A is smaller than "// &
758 0 : "what is necessary to satisfy erfc(omega*r)/r = EPS_SCHWARZ at r = ", &
759 0 : actual_x_data%potential_parameter%cutoff_radius*a_bohr*1e+10_dp, &
760 : " A. Increase input value (or omit keyword to use program default) "// &
761 0 : "to ensure accuracy."
762 0 : CPWARN(error_msg)
763 : END IF
764 0 : actual_x_data%potential_parameter%cutoff_radius = real_val
765 : END IF
766 : END IF
767 1488 : IF (actual_x_data%potential_parameter%potential_type == do_potential_id) THEN
768 28 : actual_x_data%potential_parameter%cutoff_radius = 0.0_dp
769 : END IF
770 :
771 : !! LOAD_BALANCE section
772 1488 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "LOAD_BALANCE", i_rep_section=irep)
773 1488 : CALL section_vals_val_get(hf_sub_section, "NBINS", i_val=int_val)
774 1488 : actual_x_data%load_balance_parameter%nbins = MAX(1, int_val)
775 1488 : actual_x_data%load_balance_parameter%blocks_initialized = .FALSE.
776 :
777 1488 : CALL section_vals_val_get(hf_sub_section, "RANDOMIZE", l_val=logic_val)
778 1488 : actual_x_data%load_balance_parameter%do_randomize = logic_val
779 :
780 1488 : actual_x_data%load_balance_parameter%rtp_redistribute = .FALSE.
781 1488 : IF (ASSOCIATED(dft_control%rtp_control)) THEN
782 36 : actual_x_data%load_balance_parameter%rtp_redistribute = dft_control%rtp_control%hfx_redistribute
783 : END IF
784 :
785 1488 : CALL section_vals_val_get(hf_sub_section, "BLOCK_SIZE", i_val=int_val)
786 : ! negative values ask for a computed default
787 1488 : IF (int_val <= 0) THEN
788 : ! this gives a reasonable number of blocks for binning, yet typically results in blocking.
789 : int_val = CEILING(0.1_dp*natom/ &
790 1488 : REAL(actual_x_data%load_balance_parameter%nbins*n_threads*para_env%num_pe, KIND=dp)**(0.25_dp))
791 : END IF
792 : ! at least 1 atom per block, and avoid overly large blocks
793 1488 : actual_x_data%load_balance_parameter%block_size = MIN(max_atom_block, MAX(1, int_val))
794 :
795 : CALL hfx_create_basis_types(actual_x_data%basis_parameter, actual_x_data%basis_info, qs_kind_set, &
796 1488 : orb_basis_type)
797 :
798 : !!**************************************************************************************************
799 : !! ** !! ** This code writes the contraction routines
800 : !! ** !! ** Very UGLY: BASIS_SET has to be 1 primitive and lmin=lmax=l. For g-functions
801 : !! ** !! **
802 : !! ** !! ** 1 4 4 1 1
803 : !! ** !! ** 1.0 1.0
804 : !! ** !! **
805 : !! ** k = max_am - 1
806 : !! ** write(filename,'(A,I0,A)') "sphi",k+1,"a"
807 : !! ** OPEN(UNIT=31415,FILE=filename)
808 : !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1)
809 : !! ** DO j=1,SIZE(sphi_a,2)
810 : !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN
811 : !! ** write(31415,'(A,I0,A,I0,A,I0,A,I0,A,I0,A)') "buffer1(i+imax*(",&
812 : !! ** j,&
813 : !! ** "-1)) = buffer1(i+imax*(",&
814 : !! ** j,&
815 : !! ** "-1)) + work(",&
816 : !! ** i-ncoset(k),&
817 : !! ** "+(i-1)*kmax) * sphi_a(",&
818 : !! ** i-ncoset(k),&
819 : !! ** ",",&
820 : !! ** j,&
821 : !! ** "+s_offset_a1)"
822 : !! ** END IF
823 : !! ** END DO
824 : !! ** END DO
825 : !! ** CLOSE(UNIT=31415)
826 : !! ** write(filename,'(A,I0,A)') "sphi",k+1,"b"
827 : !! ** OPEN(UNIT=31415,FILE=filename)
828 : !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1)
829 : !! ** DO j=1,SIZE(sphi_a,2)
830 : !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN
831 : !! ** write(31415,'(A,I0,A,I0,A,I0,A,I0,A,I0,A)') "buffer2(i+imax*(",&
832 : !! ** j,&
833 : !! ** "-1)) = buffer2(i+imax*(",&
834 : !! ** j,&
835 : !! ** "-1)) + buffer1(",&
836 : !! ** i-ncoset(k),&
837 : !! ** "+(i-1)*kmax) * sphi_b(",&
838 : !! ** i-ncoset(k),&
839 : !! ** ",",&
840 : !! ** j,&
841 : !! ** "+s_offset_b1)"
842 : !! **
843 : !! ** END IF
844 : !! ** END DO
845 : !! ** END DO
846 : !! ** CLOSE(UNIT=31415)
847 : !! ** write(filename,'(A,I0,A)') "sphi",k+1,"c"
848 : !! ** OPEN(UNIT=31415,FILE=filename)
849 : !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1)
850 : !! ** DO j=1,SIZE(sphi_a,2)
851 : !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN
852 : !! ** write(31415,'(A,I0,A,I0,A,I0,A,I0,A,I0,A)') "buffer1(i+imax*(",&
853 : !! ** j,&
854 : !! ** "-1)) = buffer1(i+imax*(",&
855 : !! ** j,&
856 : !! ** "-1)) + buffer2(",&
857 : !! ** i-ncoset(k),&
858 : !! ** "+(i-1)*kmax) * sphi_c(",&
859 : !! ** i-ncoset(k),&
860 : !! ** ",",&
861 : !! ** j,&
862 : !! ** "+s_offset_c1)"
863 : !! **
864 : !! ** END IF
865 : !! ** END DO
866 : !! ** END DO
867 : !! ** CLOSE(UNIT=31415)
868 : !! ** write(filename,'(A,I0,A)') "sphi",k+1,"d"
869 : !! ** OPEN(UNIT=31415,FILE=filename)
870 : !! ** DO i=ncoset(k)+1,SIZE(sphi_a,1)
871 : !! ** DO j=1,SIZE(sphi_a,2)
872 : !! ** IF( sphi_a(i,j) /= 0.0_dp) THEN
873 : !! **
874 : !! **
875 : !! ** write(31415,'(A,I0,A)') "primitives(s_offset_a1+i3, s_offset_b1+i2, s_offset_c1+i1, s_offset_d1+",&
876 : !! ** j,")= &"
877 : !! ** write(31415,'(A,I0,A)') "primitives(s_offset_a1+i3, s_offset_b1+i2, s_offset_c1+i1, s_offset_d1+",&
878 : !! ** j,")+ &"
879 : !! ** write(31415,'(A,I0,A,I0,A,I0,A)') "buffer1(",&
880 : !! ** i-ncoset(k),&
881 : !! ** "+(i-1)*kmax) * sphi_d(",&
882 : !! ** i-ncoset(k),&
883 : !! ** ",",&
884 : !! ** j,&
885 : !! ** "+s_offset_d1)"
886 : !! **
887 : !! **
888 : !! ** END IF
889 : !! ** END DO
890 : !! ** END DO
891 : !! ** CLOSE(UNIT=31415)
892 : !! ** stop
893 : !! *************************************************************************************************************************
894 :
895 1488 : IF (actual_x_data%periodic_parameter%do_periodic) THEN
896 448 : hf_pbc_section => section_vals_get_subs_vals(hfx_section, "PERIODIC", i_rep_section=irep)
897 448 : CALL section_vals_val_get(hf_pbc_section, "NUMBER_OF_SHELLS", i_val=pbc_shells)
898 448 : actual_x_data%periodic_parameter%number_of_shells_from_input = pbc_shells
899 3584 : ALLOCATE (actual_x_data%neighbor_cells(1))
900 896 : CALL hfx_create_neighbor_cells(actual_x_data, pbc_shells, cell, i_thread, nkp_grid=nkp_grid)
901 : ELSE
902 8320 : ALLOCATE (actual_x_data%neighbor_cells(1))
903 : ! ** Initialize this guy to enable non periodic stress regtests
904 1040 : actual_x_data%periodic_parameter%R_max_stress = 1.0_dp
905 : END IF
906 :
907 1488 : nkind = SIZE(qs_kind_set, 1)
908 1488 : max_set = actual_x_data%basis_info%max_set
909 :
910 : !! ** This guy is allocated on the master thread only
911 1488 : IF (i_thread == 1) THEN
912 5952 : ALLOCATE (actual_x_data%is_assoc_atomic_block(natom, natom))
913 4464 : ALLOCATE (actual_x_data%atomic_block_offset(natom, natom))
914 8928 : ALLOCATE (actual_x_data%set_offset(max_set, max_set, nkind, nkind))
915 4464 : ALLOCATE (actual_x_data%block_offset(para_env%num_pe + 1))
916 : END IF
917 :
918 2976 : ALLOCATE (actual_x_data%distribution_forces(1))
919 2976 : ALLOCATE (actual_x_data%distribution_energy(1))
920 :
921 1488 : actual_x_data%memory_parameter%size_p_screen = 0_int_8
922 1488 : IF (i_thread == 1) THEN
923 5952 : ALLOCATE (actual_x_data%atomic_pair_list(natom, natom))
924 4464 : ALLOCATE (actual_x_data%atomic_pair_list_forces(natom, natom))
925 : END IF
926 :
927 1488 : IF (actual_x_data%screening_parameter%do_initial_p_screening .OR. &
928 : actual_x_data%screening_parameter%do_p_screening_forces) THEN
929 : !! ** This guy is allocated on the master thread only
930 1458 : IF (i_thread == 1) THEN
931 5832 : ALLOCATE (actual_x_data%pmax_atom(natom, natom))
932 8768 : ALLOCATE (actual_x_data%initial_p(nkind*(nkind + 1)/2))
933 1458 : i = 1
934 4120 : DO ikind = 1, nkind
935 2662 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_a)
936 2662 : nseta = actual_x_data%basis_parameter(ikind)%nset
937 8514 : DO jkind = ikind, nkind
938 4394 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=natom_b)
939 4394 : nsetb = actual_x_data%basis_parameter(jkind)%nset
940 26364 : ALLOCATE (actual_x_data%initial_p(i)%p_kind(nseta, nsetb, natom_a, natom_b))
941 : actual_x_data%memory_parameter%size_p_screen = &
942 4394 : actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b
943 11450 : i = i + 1
944 : END DO
945 : END DO
946 :
947 4374 : ALLOCATE (actual_x_data%pmax_atom_forces(natom, natom))
948 7310 : ALLOCATE (actual_x_data%initial_p_forces(nkind*(nkind + 1)/2))
949 1458 : i = 1
950 4120 : DO ikind = 1, nkind
951 2662 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_a)
952 2662 : nseta = actual_x_data%basis_parameter(ikind)%nset
953 8514 : DO jkind = ikind, nkind
954 4394 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=natom_b)
955 4394 : nsetb = actual_x_data%basis_parameter(jkind)%nset
956 26364 : ALLOCATE (actual_x_data%initial_p_forces(i)%p_kind(nseta, nsetb, natom_a, natom_b))
957 : actual_x_data%memory_parameter%size_p_screen = &
958 4394 : actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b
959 11450 : i = i + 1
960 : END DO
961 : END DO
962 : END IF
963 4374 : ALLOCATE (actual_x_data%map_atom_to_kind_atom(natom))
964 1458 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
965 :
966 4374 : ALLOCATE (atom2kind(nkind))
967 1458 : atom2kind = 0
968 6138 : DO iatom = 1, natom
969 4680 : ikind = kind_of(iatom)
970 4680 : atom2kind(ikind) = atom2kind(ikind) + 1
971 6138 : actual_x_data%map_atom_to_kind_atom(iatom) = atom2kind(ikind)
972 : END DO
973 1458 : DEALLOCATE (kind_of, atom2kind)
974 : END IF
975 :
976 : ! ** Initialize libint type
977 1488 : CALL cp_libint_static_init()
978 1488 : CALL cp_libint_init_eri(actual_x_data%lib, actual_x_data%basis_info%max_am)
979 1488 : CALL cp_libint_init_eri1(actual_x_data%lib_deriv, actual_x_data%basis_info%max_am)
980 1488 : CALL cp_libint_set_contrdepth(actual_x_data%lib, 1)
981 1488 : CALL cp_libint_set_contrdepth(actual_x_data%lib_deriv, 1)
982 :
983 1488 : CALL alloc_containers(actual_x_data%store_ints, 1)
984 1488 : CALL alloc_containers(actual_x_data%store_forces, 1)
985 :
986 1488 : actual_x_data%store_ints%maxval_cache_disk%element_counter = 1
987 1488 : ALLOCATE (actual_x_data%store_ints%maxval_container_disk)
988 1525200 : ALLOCATE (actual_x_data%store_ints%maxval_container_disk%first)
989 1488 : actual_x_data%store_ints%maxval_container_disk%first%prev => NULL()
990 1488 : actual_x_data%store_ints%maxval_container_disk%first%next => NULL()
991 1488 : actual_x_data%store_ints%maxval_container_disk%current => actual_x_data%store_ints%maxval_container_disk%first
992 1525200 : actual_x_data%store_ints%maxval_container_disk%current%data = 0
993 1488 : actual_x_data%store_ints%maxval_container_disk%element_counter = 1
994 1488 : actual_x_data%store_ints%maxval_container_disk%file_counter = 1
995 1488 : actual_x_data%store_ints%maxval_container_disk%desc = 'Max_'
996 1488 : actual_x_data%store_ints%maxval_container_disk%unit = -1
997 : WRITE (actual_x_data%store_ints%maxval_container_disk%filename, '(A,I0,A,A,A)') &
998 1488 : TRIM(actual_x_data%memory_parameter%storage_location), &
999 2976 : storage_id, "_", actual_x_data%store_ints%maxval_container_disk%desc, "6"
1000 1488 : CALL compress(actual_x_data%store_ints%maxval_container_disk%filename, .TRUE.)
1001 96720 : ALLOCATE (actual_x_data%store_ints%integral_containers_disk(64))
1002 96720 : DO i = 1, 64
1003 95232 : actual_x_data%store_ints%integral_caches_disk(i)%element_counter = 1
1004 97612800 : actual_x_data%store_ints%integral_caches_disk(i)%data = 0
1005 97612800 : ALLOCATE (actual_x_data%store_ints%integral_containers_disk(i)%first)
1006 95232 : actual_x_data%store_ints%integral_containers_disk(i)%first%prev => NULL()
1007 95232 : actual_x_data%store_ints%integral_containers_disk(i)%first%next => NULL()
1008 : actual_x_data%store_ints%integral_containers_disk(i)%current => &
1009 95232 : actual_x_data%store_ints%integral_containers_disk(i)%first
1010 97612800 : actual_x_data%store_ints%integral_containers_disk(i)%current%data = 0
1011 95232 : actual_x_data%store_ints%integral_containers_disk(i)%element_counter = 1
1012 95232 : actual_x_data%store_ints%integral_containers_disk(i)%file_counter = 1
1013 95232 : actual_x_data%store_ints%integral_containers_disk(i)%desc = 'Int_'
1014 95232 : actual_x_data%store_ints%integral_containers_disk(i)%unit = -1
1015 : WRITE (actual_x_data%store_ints%integral_containers_disk(i)%filename, '(A,I0,A,A,I0)') &
1016 95232 : TRIM(actual_x_data%memory_parameter%storage_location), &
1017 190464 : storage_id, "_", actual_x_data%store_ints%integral_containers_disk(i)%desc, i
1018 96720 : CALL compress(actual_x_data%store_ints%integral_containers_disk(i)%filename, .TRUE.)
1019 : END DO
1020 :
1021 1488 : actual_x_data%b_first_load_balance_energy = .TRUE.
1022 1488 : actual_x_data%b_first_load_balance_forces = .TRUE.
1023 :
1024 1488 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "RI", i_rep_section=irep)
1025 1488 : IF (actual_x_data%do_hfx_ri) THEN
1026 114 : CPASSERT(PRESENT(nelectron_total))
1027 798 : ALLOCATE (actual_x_data%ri_data)
1028 : CALL hfx_ri_init_read_input_from_hfx(actual_x_data%ri_data, actual_x_data, hfx_section, &
1029 : hf_sub_section, qs_kind_set, &
1030 : particle_set, atomic_kind_set, dft_control, para_env, irep, &
1031 114 : nelectron_total, orb_basis_type, ri_basis_type)
1032 : END IF
1033 :
1034 : ! ACE section — read only on thread 1 to avoid redundant work
1035 13382 : IF (i_thread == 1) THEN
1036 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "ACE", &
1037 1488 : i_rep_section=irep)
1038 1488 : CALL section_vals_get(hf_sub_section, explicit=logic_val)
1039 1488 : IF (logic_val) THEN
1040 : CALL section_vals_val_get(hf_sub_section, "ACTIVE", &
1041 8 : l_val=actual_x_data%use_ace)
1042 : CALL section_vals_val_get(hf_sub_section, "REBUILD_FREQUENCY", &
1043 8 : i_val=actual_x_data%ace_rebuild_freq)
1044 : END IF
1045 : ! Sanity checks
1046 1488 : IF (actual_x_data%use_ace) THEN
1047 : ! ACE requires HFX to be meaningful
1048 8 : IF (actual_x_data%general_parameter%fraction <= 0.0_dp) THEN
1049 0 : CPABORT("ACE requires FRACTION > 0.")
1050 : END IF
1051 : ! If frequency is 1, it is full HFX
1052 8 : IF (actual_x_data%ace_rebuild_freq < 1) THEN
1053 0 : CPABORT("ACE: REBUILD_FREQUENCY must be >= 1")
1054 : END IF
1055 : END IF
1056 : END IF
1057 : END DO
1058 : END DO
1059 :
1060 2966 : DO irep = 1, n_rep_hf
1061 1488 : actual_x_data => x_data(irep, 1)
1062 2966 : CALL hfx_print_info(actual_x_data, hfx_section, irep)
1063 : END DO
1064 :
1065 1478 : CALL timestop(handle)
1066 :
1067 5912 : END SUBROUTINE hfx_create
1068 :
1069 : ! **************************************************************************************************
1070 : !> \brief Read RI input and initialize RI data for use within Hartree-Fock
1071 : !> \param ri_data ...
1072 : !> \param x_data ...
1073 : !> \param hfx_section ...
1074 : !> \param ri_section ...
1075 : !> \param qs_kind_set ...
1076 : !> \param particle_set ...
1077 : !> \param atomic_kind_set ...
1078 : !> \param dft_control ...
1079 : !> \param para_env ...
1080 : !> \param irep ...
1081 : !> \param nelectron_total ...
1082 : !> \param orb_basis_type ...
1083 : !> \param ri_basis_type ...
1084 : ! **************************************************************************************************
1085 114 : SUBROUTINE hfx_ri_init_read_input_from_hfx(ri_data, x_data, hfx_section, ri_section, qs_kind_set, &
1086 : particle_set, atomic_kind_set, dft_control, para_env, irep, &
1087 : nelectron_total, orb_basis_type, ri_basis_type)
1088 : TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1089 : TYPE(hfx_type), INTENT(INOUT) :: x_data
1090 : TYPE(section_vals_type), POINTER :: hfx_section, ri_section
1091 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1092 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1093 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1094 : TYPE(dft_control_type), POINTER :: dft_control
1095 : TYPE(mp_para_env_type) :: para_env
1096 : INTEGER, INTENT(IN) :: irep, nelectron_total
1097 : CHARACTER(LEN=*) :: orb_basis_type, ri_basis_type
1098 :
1099 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_init_read_input_from_hfx'
1100 :
1101 : CHARACTER(LEN=512) :: error_msg
1102 : CHARACTER(LEN=default_path_length) :: char_val, t_c_filename
1103 : INTEGER :: handle, unit_nr, unit_nr_dbcsr
1104 : TYPE(cp_logger_type), POINTER :: logger
1105 : TYPE(section_vals_type), POINTER :: hf_sub_section
1106 :
1107 114 : CALL timeset(routineN, handle)
1108 :
1109 114 : NULLIFY (hf_sub_section)
1110 :
1111 : ASSOCIATE (hfx_pot => ri_data%hfx_pot)
1112 114 : hfx_pot%potential_type = x_data%potential_parameter%potential_type
1113 114 : hfx_pot%omega = x_data%potential_parameter%omega
1114 114 : hfx_pot%cutoff_radius = x_data%potential_parameter%cutoff_radius
1115 114 : hfx_pot%scale_coulomb = x_data%potential_parameter%scale_coulomb
1116 114 : hfx_pot%scale_longrange = x_data%potential_parameter%scale_longrange
1117 : END ASSOCIATE
1118 114 : ri_data%ri_section => ri_section
1119 114 : ri_data%hfx_section => hfx_section
1120 114 : ri_data%eps_schwarz = x_data%screening_parameter%eps_schwarz
1121 114 : ri_data%eps_schwarz_forces = x_data%screening_parameter%eps_schwarz_forces
1122 :
1123 114 : logger => cp_get_default_logger()
1124 : unit_nr_dbcsr = cp_print_key_unit_nr(logger, ri_data%ri_section, "PRINT%RI_INFO", &
1125 114 : extension=".dbcsrLog")
1126 :
1127 : unit_nr = cp_print_key_unit_nr(logger, ri_data%hfx_section, "HF_INFO", &
1128 114 : extension=".scfLog")
1129 :
1130 114 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL", i_rep_section=irep)
1131 114 : CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val)
1132 114 : CALL compress(char_val, .TRUE.)
1133 :
1134 114 : IF (.NOT. file_exists(char_val)) THEN
1135 : WRITE (error_msg, '(A,A,A)') "File not found. Please check T_C_G_DATA "// &
1136 0 : "in the INTERACTION_POTENTIAL section"
1137 0 : CPABORT(error_msg)
1138 : ELSE
1139 114 : t_c_filename = char_val
1140 : END IF
1141 :
1142 : CALL hfx_ri_init_read_input(ri_data, ri_section, qs_kind_set, particle_set, atomic_kind_set, &
1143 : orb_basis_type, ri_basis_type, para_env, unit_nr, unit_nr_dbcsr, &
1144 114 : nelectron_total, t_c_filename=t_c_filename)
1145 :
1146 114 : IF (dft_control%smear .AND. ri_data%flavor == ri_mo) THEN
1147 0 : CPABORT("RI_FLAVOR MO is not consistent with smearing. Please use RI_FLAVOR RHO.")
1148 : END IF
1149 :
1150 114 : CALL timestop(handle)
1151 :
1152 114 : END SUBROUTINE hfx_ri_init_read_input_from_hfx
1153 :
1154 : ! **************************************************************************************************
1155 : !> \brief General routine for reading input of RI section and initializing RI data
1156 : !> \param ri_data ...
1157 : !> \param ri_section ...
1158 : !> \param qs_kind_set ...
1159 : !> \param particle_set ...
1160 : !> \param atomic_kind_set ...
1161 : !> \param orb_basis_type ...
1162 : !> \param ri_basis_type ...
1163 : !> \param para_env ...
1164 : !> \param unit_nr unit number of general output
1165 : !> \param unit_nr_dbcsr unit number for logging DBCSR tensor operations
1166 : !> \param nelectron_total ...
1167 : !> \param t_c_filename ...
1168 : ! **************************************************************************************************
1169 114 : SUBROUTINE hfx_ri_init_read_input(ri_data, ri_section, qs_kind_set, &
1170 : particle_set, atomic_kind_set, orb_basis_type, ri_basis_type, para_env, &
1171 : unit_nr, unit_nr_dbcsr, nelectron_total, t_c_filename)
1172 : TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1173 : TYPE(section_vals_type), POINTER :: ri_section
1174 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1175 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1176 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1177 : CHARACTER(LEN=*), INTENT(IN) :: orb_basis_type, ri_basis_type
1178 : TYPE(mp_para_env_type) :: para_env
1179 : INTEGER, INTENT(IN) :: unit_nr, unit_nr_dbcsr, nelectron_total
1180 : CHARACTER(len=*), INTENT(IN), OPTIONAL :: t_c_filename
1181 :
1182 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_init_read_input'
1183 :
1184 : INTEGER :: handle
1185 : LOGICAL :: explicit
1186 : REAL(dp) :: eps_storage_scaling
1187 :
1188 114 : CALL timeset(routineN, handle)
1189 :
1190 114 : CALL section_vals_val_get(ri_section, "EPS_FILTER", r_val=ri_data%filter_eps)
1191 114 : CALL section_vals_val_get(ri_section, "EPS_FILTER_2C", r_val=ri_data%filter_eps_2c)
1192 114 : CALL section_vals_val_get(ri_section, "EPS_STORAGE_SCALING", r_val=eps_storage_scaling)
1193 114 : ri_data%filter_eps_storage = ri_data%filter_eps*eps_storage_scaling
1194 114 : CALL section_vals_val_get(ri_section, "EPS_FILTER_MO", r_val=ri_data%filter_eps_mo)
1195 :
1196 : ASSOCIATE (ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
1197 114 : CALL section_vals_val_get(ri_section, "RI_METRIC", i_val=ri_metric%potential_type, explicit=explicit)
1198 114 : IF (.NOT. explicit .OR. ri_metric%potential_type == 0) THEN
1199 44 : ri_metric%potential_type = hfx_pot%potential_type
1200 : END IF
1201 :
1202 114 : CALL section_vals_val_get(ri_section, "OMEGA", r_val=ri_metric%omega, explicit=explicit)
1203 114 : IF (.NOT. explicit) THEN
1204 114 : ri_metric%omega = hfx_pot%omega
1205 : END IF
1206 :
1207 114 : CALL section_vals_val_get(ri_section, "CUTOFF_RADIUS", r_val=ri_metric%cutoff_radius, explicit=explicit)
1208 114 : IF (.NOT. explicit) THEN
1209 106 : ri_metric%cutoff_radius = hfx_pot%cutoff_radius
1210 : END IF
1211 :
1212 114 : CALL section_vals_val_get(ri_section, "SCALE_COULOMB", r_val=ri_metric%scale_coulomb, explicit=explicit)
1213 114 : IF (.NOT. explicit) THEN
1214 114 : ri_metric%scale_coulomb = hfx_pot%scale_coulomb
1215 : END IF
1216 :
1217 114 : CALL section_vals_val_get(ri_section, "SCALE_LONGRANGE", r_val=ri_metric%scale_longrange, explicit=explicit)
1218 114 : IF (.NOT. explicit) THEN
1219 114 : ri_metric%scale_longrange = hfx_pot%scale_longrange
1220 : END IF
1221 :
1222 114 : IF (ri_metric%potential_type == do_potential_short) THEN
1223 2 : CALL erfc_cutoff(ri_data%eps_schwarz, ri_metric%omega, ri_metric%cutoff_radius)
1224 : END IF
1225 114 : IF (ri_metric%potential_type == do_potential_id) ri_metric%cutoff_radius = 0.0_dp
1226 : END ASSOCIATE
1227 :
1228 114 : CALL section_vals_val_get(ri_section, "2C_MATRIX_FUNCTIONS", i_val=ri_data%t2c_method)
1229 114 : CALL section_vals_val_get(ri_section, "EPS_EIGVAL", r_val=ri_data%eps_eigval)
1230 114 : CALL section_vals_val_get(ri_section, "CHECK_2C_MATRIX", l_val=ri_data%check_2c_inv)
1231 114 : CALL section_vals_val_get(ri_section, "CALC_COND_NUM", l_val=ri_data%calc_condnum)
1232 114 : CALL section_vals_val_get(ri_section, "SQRT_ORDER", i_val=ri_data%t2c_sqrt_order)
1233 114 : CALL section_vals_val_get(ri_section, "EPS_LANCZOS", r_val=ri_data%eps_lanczos)
1234 114 : CALL section_vals_val_get(ri_section, "MAX_ITER_LANCZOS", i_val=ri_data%max_iter_lanczos)
1235 114 : CALL section_vals_val_get(ri_section, "RI_FLAVOR", i_val=ri_data%flavor)
1236 114 : CALL section_vals_val_get(ri_section, "EPS_PGF_ORB", r_val=ri_data%eps_pgf_orb)
1237 114 : CALL section_vals_val_get(ri_section, "MIN_BLOCK_SIZE", i_val=ri_data%min_bsize)
1238 114 : CALL section_vals_val_get(ri_section, "MAX_BLOCK_SIZE_MO", i_val=ri_data%max_bsize_MO)
1239 114 : CALL section_vals_val_get(ri_section, "MEMORY_CUT", i_val=ri_data%n_mem_input)
1240 114 : CALL section_vals_val_get(ri_section, "FLAVOR_SWITCH_MEMORY_CUT", i_val=ri_data%n_mem_flavor_switch)
1241 :
1242 114 : ri_data%orb_basis_type = orb_basis_type
1243 114 : ri_data%ri_basis_type = ri_basis_type
1244 114 : ri_data%nelectron_total = nelectron_total
1245 114 : ri_data%input_flavor = ri_data%flavor
1246 :
1247 114 : IF (PRESENT(t_c_filename)) THEN
1248 114 : ri_data%ri_metric%filename = t_c_filename
1249 114 : ri_data%hfx_pot%filename = t_c_filename
1250 : END IF
1251 :
1252 114 : ri_data%unit_nr_dbcsr = unit_nr_dbcsr
1253 114 : ri_data%unit_nr = unit_nr
1254 114 : ri_data%dbcsr_nflop = 0
1255 114 : ri_data%dbcsr_time = 0.0_dp
1256 :
1257 114 : CALL hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
1258 :
1259 114 : CALL timestop(handle)
1260 :
1261 798 : END SUBROUTINE hfx_ri_init_read_input
1262 :
1263 : ! **************************************************************************************************
1264 : !> \brief ...
1265 : !> \param ri_data ...
1266 : !> \param qs_kind_set ...
1267 : !> \param particle_set ...
1268 : !> \param atomic_kind_set ...
1269 : !> \param para_env ...
1270 : ! **************************************************************************************************
1271 136 : SUBROUTINE hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
1272 : TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1273 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1274 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1275 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1276 : TYPE(mp_para_env_type) :: para_env
1277 :
1278 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_init'
1279 :
1280 : INTEGER :: handle, i_mem, j_mem, MO_dim, natom, &
1281 : nkind, nproc
1282 136 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bsizes_AO_store, bsizes_RI_store, dist1, &
1283 136 : dist2, dist3, dist_AO_1, dist_AO_2, &
1284 : dist_RI
1285 : INTEGER, DIMENSION(2) :: pdims_2d
1286 : INTEGER, DIMENSION(3) :: pdims
1287 : LOGICAL :: same_op
1288 : TYPE(distribution_3d_type) :: dist_3d
1289 : TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1290 136 : DIMENSION(:) :: basis_set_AO, basis_set_RI
1291 136 : TYPE(mp_cart_type) :: mp_comm_3d
1292 :
1293 136 : CALL cite_reference(Bussy2023)
1294 :
1295 136 : CALL timeset(routineN, handle)
1296 :
1297 : ! initialize libint
1298 136 : CALL cp_libint_static_init()
1299 :
1300 136 : natom = SIZE(particle_set)
1301 136 : nkind = SIZE(qs_kind_set, 1)
1302 136 : nproc = para_env%num_pe
1303 :
1304 : ASSOCIATE (ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
1305 136 : IF (ri_metric%potential_type == do_potential_short) THEN
1306 2 : CALL erfc_cutoff(ri_data%eps_schwarz, ri_metric%omega, ri_metric%cutoff_radius)
1307 : END IF
1308 :
1309 136 : IF (hfx_pot%potential_type == do_potential_short) THEN
1310 : ! need a more accurate threshold for determining 2-center integral operator range
1311 : ! because stability of matrix inversion/sqrt is sensitive to this
1312 4 : CALL erfc_cutoff(ri_data%filter_eps_2c, hfx_pot%omega, hfx_pot%cutoff_radius)
1313 : END IF
1314 : ! determine whether RI metric is same operator as used in HFX
1315 136 : same_op = compare_potential_types(ri_metric, hfx_pot)
1316 : END ASSOCIATE
1317 :
1318 136 : ri_data%same_op = same_op
1319 :
1320 136 : pdims = 0
1321 136 : CALL mp_comm_3d%create(para_env, 3, pdims)
1322 :
1323 408 : ALLOCATE (ri_data%bsizes_RI(natom))
1324 272 : ALLOCATE (ri_data%bsizes_AO(natom))
1325 1016 : ALLOCATE (basis_set_RI(nkind), basis_set_AO(nkind))
1326 136 : CALL basis_set_list_setup(basis_set_RI, ri_data%ri_basis_type, qs_kind_set)
1327 136 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_data%bsizes_RI, basis=basis_set_RI)
1328 136 : CALL basis_set_list_setup(basis_set_AO, ri_data%orb_basis_type, qs_kind_set)
1329 136 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_data%bsizes_AO, basis=basis_set_AO)
1330 :
1331 272 : ALLOCATE (dist_RI(natom))
1332 272 : ALLOCATE (dist_AO_1(natom))
1333 272 : ALLOCATE (dist_AO_2(natom))
1334 136 : CALL dbt_default_distvec(natom, pdims(1), ri_data%bsizes_RI, dist_RI)
1335 136 : CALL dbt_default_distvec(natom, pdims(2), ri_data%bsizes_AO, dist_AO_1)
1336 136 : CALL dbt_default_distvec(natom, pdims(3), ri_data%bsizes_AO, dist_AO_2)
1337 : CALL distribution_3d_create(dist_3d, dist_RI, dist_ao_1, dist_ao_2, nkind, particle_set, &
1338 136 : mp_comm_3d, own_comm=.TRUE.)
1339 :
1340 408 : ALLOCATE (ri_data%pgrid)
1341 136 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid)
1342 :
1343 408 : ALLOCATE (ri_data%pgrid_2d)
1344 136 : pdims_2d = 0
1345 136 : CALL dbt_pgrid_create(para_env, pdims_2d, ri_data%pgrid_2d)
1346 :
1347 136 : ri_data%dist_3d = dist_3d
1348 :
1349 : CALL dbt_distribution_new(ri_data%dist, ri_data%pgrid, &
1350 136 : dist_RI, dist_AO_1, dist_AO_2)
1351 :
1352 136 : DEALLOCATE (dist_AO_1, dist_AO_2, dist_RI)
1353 :
1354 136 : ri_data%num_pe = para_env%num_pe
1355 :
1356 : ! initialize tensors expressed in basis representation
1357 136 : CALL pgf_block_sizes(atomic_kind_set, basis_set_AO, ri_data%min_bsize, ri_data%bsizes_AO_split)
1358 136 : CALL pgf_block_sizes(atomic_kind_set, basis_set_RI, ri_data%min_bsize, ri_data%bsizes_RI_split)
1359 :
1360 136 : CALL pgf_block_sizes(atomic_kind_set, basis_set_AO, 1, bsizes_AO_store)
1361 136 : CALL pgf_block_sizes(atomic_kind_set, basis_set_RI, 1, bsizes_RI_store)
1362 :
1363 666 : CALL split_block_sizes([SUM(ri_data%bsizes_AO)], ri_data%bsizes_AO_fit, default_block_size)
1364 666 : CALL split_block_sizes([SUM(ri_data%bsizes_RI)], ri_data%bsizes_RI_fit, default_block_size)
1365 :
1366 136 : IF (ri_data%flavor == ri_pmat) THEN
1367 :
1368 : !2 batching loops in RHO flavor SCF calculations => need to take the square root of MEMORY_CUT
1369 118 : ri_data%n_mem = ri_data%n_mem_input
1370 118 : ri_data%n_mem_RI = ri_data%n_mem_input
1371 :
1372 : CALL create_tensor_batches(ri_data%bsizes_AO_split, ri_data%n_mem, ri_data%starts_array_mem, &
1373 : ri_data%ends_array_mem, ri_data%starts_array_mem_block, &
1374 118 : ri_data%ends_array_mem_block)
1375 :
1376 : CALL create_tensor_batches(ri_data%bsizes_RI_split, ri_data%n_mem_RI, &
1377 : ri_data%starts_array_RI_mem, ri_data%ends_array_RI_mem, &
1378 118 : ri_data%starts_array_RI_mem_block, ri_data%ends_array_RI_mem_block)
1379 :
1380 354 : ALLOCATE (ri_data%pgrid_1)
1381 354 : ALLOCATE (ri_data%pgrid_2)
1382 118 : pdims = 0
1383 :
1384 : CALL dbt_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), &
1385 472 : SIZE(ri_data%bsizes_AO_split)])
1386 :
1387 118 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_1)
1388 :
1389 826 : pdims = pdims([2, 1, 3])
1390 118 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_2)
1391 :
1392 1062 : ALLOCATE (ri_data%t_3c_int_ctr_1(1, 1))
1393 : CALL create_3c_tensor(ri_data%t_3c_int_ctr_1(1, 1), dist1, dist2, dist3, &
1394 : ri_data%pgrid_1, ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, &
1395 118 : ri_data%bsizes_AO_split, [1, 2], [3], name="(AO RI | AO)")
1396 118 : DEALLOCATE (dist1, dist2, dist3)
1397 :
1398 1516 : ALLOCATE (ri_data%blk_indices(ri_data%n_mem, ri_data%n_mem_RI))
1399 250732 : ALLOCATE (ri_data%store_3c(ri_data%n_mem, ri_data%n_mem_RI))
1400 410 : DO i_mem = 1, ri_data%n_mem
1401 1162 : DO j_mem = 1, ri_data%n_mem_RI
1402 1044 : CALL alloc_containers(ri_data%store_3c(i_mem, j_mem), 1)
1403 : END DO
1404 : END DO
1405 :
1406 1062 : ALLOCATE (ri_data%t_3c_int_ctr_2(1, 1))
1407 : CALL create_3c_tensor(ri_data%t_3c_int_ctr_2(1, 1), dist1, dist2, dist3, &
1408 : ri_data%pgrid_1, ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, &
1409 118 : ri_data%bsizes_AO_split, [1, 2], [3], name="(AO RI | AO)")
1410 118 : DEALLOCATE (dist1, dist2, dist3)
1411 :
1412 1062 : ALLOCATE (ri_data%t_3c_int_ctr_3(1, 1))
1413 : CALL create_3c_tensor(ri_data%t_3c_int_ctr_3(1, 1), dist1, dist2, dist3, &
1414 : ri_data%pgrid_2, ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, &
1415 118 : ri_data%bsizes_AO_split, [1], [2, 3], name="(RI | AO AO)")
1416 118 : DEALLOCATE (dist1, dist2, dist3)
1417 :
1418 1062 : ALLOCATE (ri_data%t_2c_int(1, 1))
1419 : CALL create_2c_tensor(ri_data%t_2c_int(1, 1), dist1, dist2, ri_data%pgrid_2d, &
1420 : ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, &
1421 118 : name="(RI | RI)")
1422 118 : DEALLOCATE (dist1, dist2)
1423 :
1424 : !We store previous Pmat and KS mat, so that we can work with Delta P and gain sprasity as we go
1425 1180 : ALLOCATE (ri_data%rho_ao_t(2, 1))
1426 : CALL create_2c_tensor(ri_data%rho_ao_t(1, 1), dist1, dist2, ri_data%pgrid_2d, &
1427 : ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
1428 118 : name="(AO | AO)")
1429 118 : DEALLOCATE (dist1, dist2)
1430 118 : CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
1431 :
1432 1180 : ALLOCATE (ri_data%ks_t(2, 1))
1433 : CALL create_2c_tensor(ri_data%ks_t(1, 1), dist1, dist2, ri_data%pgrid_2d, &
1434 : ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
1435 118 : name="(AO | AO)")
1436 118 : DEALLOCATE (dist1, dist2)
1437 118 : CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
1438 :
1439 18 : ELSE IF (ri_data%flavor == ri_mo) THEN
1440 180 : ALLOCATE (ri_data%t_2c_int(2, 1))
1441 :
1442 : CALL create_2c_tensor(ri_data%t_2c_int(1, 1), dist1, dist2, ri_data%pgrid_2d, &
1443 : ri_data%bsizes_RI_fit, ri_data%bsizes_RI_fit, &
1444 18 : name="(RI | RI)")
1445 18 : CALL dbt_create(ri_data%t_2c_int(1, 1), ri_data%t_2c_int(2, 1))
1446 :
1447 18 : DEALLOCATE (dist1, dist2)
1448 :
1449 162 : ALLOCATE (ri_data%t_3c_int_ctr_1(1, 1))
1450 :
1451 54 : ALLOCATE (ri_data%pgrid_1)
1452 54 : ALLOCATE (ri_data%pgrid_2)
1453 : pdims = 0
1454 :
1455 18 : ri_data%n_mem = ri_data%n_mem_input**2
1456 18 : IF (ri_data%n_mem > ri_data%nelectron_total/2) ri_data%n_mem = MAX(ri_data%nelectron_total/2, 1)
1457 : ! Size of dimension corresponding to MOs is nelectron/2 and divided by the memory factor
1458 : ! we are using ceiling of that division to make sure that no MO dimension (after memory cut)
1459 : ! is larger than this (it is however not a problem for load balancing if actual MO dimension
1460 : ! is slightly smaller)
1461 18 : MO_dim = MAX((ri_data%nelectron_total/2 - 1)/ri_data%n_mem + 1, 1)
1462 18 : MO_dim = (MO_dim - 1)/ri_data%max_bsize_MO + 1
1463 :
1464 18 : pdims = 0
1465 72 : CALL dbt_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), MO_dim])
1466 :
1467 18 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_1)
1468 :
1469 126 : pdims = pdims([3, 2, 1])
1470 18 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_2)
1471 :
1472 : CALL create_3c_tensor(ri_data%t_3c_int_ctr_1(1, 1), dist1, dist2, dist3, &
1473 : ri_data%pgrid_1, ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, &
1474 18 : [1, 2], [3], name="(AO RI | AO)")
1475 18 : DEALLOCATE (dist1, dist2, dist3)
1476 :
1477 162 : ALLOCATE (ri_data%t_3c_int_ctr_2(1, 1))
1478 : CALL create_3c_tensor(ri_data%t_3c_int_ctr_2(1, 1), dist1, dist2, dist3, &
1479 : ri_data%pgrid_2, ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, &
1480 18 : [1], [2, 3], name="(AO | RI AO)")
1481 18 : DEALLOCATE (dist1, dist2, dist3)
1482 :
1483 : END IF
1484 :
1485 : !For forces
1486 1224 : ALLOCATE (ri_data%t_2c_inv(1, 1))
1487 : CALL create_2c_tensor(ri_data%t_2c_inv(1, 1), dist1, dist2, ri_data%pgrid_2d, &
1488 : ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, &
1489 136 : name="(RI | RI)")
1490 136 : DEALLOCATE (dist1, dist2)
1491 :
1492 1224 : ALLOCATE (ri_data%t_2c_pot(1, 1))
1493 : CALL create_2c_tensor(ri_data%t_2c_pot(1, 1), dist1, dist2, ri_data%pgrid_2d, &
1494 : ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, &
1495 136 : name="(RI | RI)")
1496 136 : DEALLOCATE (dist1, dist2)
1497 :
1498 136 : CALL timestop(handle)
1499 :
1500 272 : END SUBROUTINE hfx_ri_init
1501 :
1502 : ! **************************************************************************************************
1503 : !> \brief ...
1504 : !> \param ri_data ...
1505 : ! **************************************************************************************************
1506 114 : SUBROUTINE hfx_ri_write_stats(ri_data)
1507 : TYPE(hfx_ri_type), INTENT(IN) :: ri_data
1508 :
1509 : REAL(dp) :: my_flop_rate
1510 :
1511 : ASSOCIATE (unit_nr => ri_data%unit_nr, dbcsr_nflop => ri_data%dbcsr_nflop, &
1512 : dbcsr_time => ri_data%dbcsr_time, num_pe => ri_data%num_pe)
1513 114 : my_flop_rate = REAL(dbcsr_nflop, dp)/(1.0E09_dp*ri_data%dbcsr_time)
1514 114 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T2,A,T73,ES8.2)") &
1515 51 : "RI-HFX PERFORMANCE| DBT total number of flops:", REAL(dbcsr_nflop*num_pe, dp)
1516 114 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T2,A,T66,F15.2)") &
1517 51 : "RI-HFX PERFORMANCE| DBT total execution time:", dbcsr_time
1518 114 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T2,A,T66,F15.2)") &
1519 165 : "RI-HFX PERFORMANCE| DBT flop rate (Gflops / MPI rank):", my_flop_rate
1520 : END ASSOCIATE
1521 114 : END SUBROUTINE hfx_ri_write_stats
1522 :
1523 : ! **************************************************************************************************
1524 : !> \brief ...
1525 : !> \param ri_data ...
1526 : !> \param write_stats ...
1527 : ! **************************************************************************************************
1528 136 : SUBROUTINE hfx_ri_release(ri_data, write_stats)
1529 : TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1530 : LOGICAL, OPTIONAL :: write_stats
1531 :
1532 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_release'
1533 :
1534 : INTEGER :: handle, i, i_mem, ispin, j, j_mem, unused
1535 : LOGICAL :: my_write_stats
1536 :
1537 136 : CALL timeset(routineN, handle)
1538 :
1539 : ! cleanup libint
1540 136 : CALL cp_libint_static_cleanup()
1541 :
1542 136 : my_write_stats = .TRUE.
1543 136 : IF (PRESENT(write_stats)) my_write_stats = write_stats
1544 136 : IF (my_write_stats) CALL hfx_ri_write_stats(ri_data)
1545 :
1546 136 : IF (ASSOCIATED(ri_data%pgrid)) THEN
1547 136 : CALL dbt_pgrid_destroy(ri_data%pgrid)
1548 136 : DEALLOCATE (ri_data%pgrid)
1549 : END IF
1550 136 : IF (ASSOCIATED(ri_data%pgrid_1)) THEN
1551 136 : CALL dbt_pgrid_destroy(ri_data%pgrid_1)
1552 136 : DEALLOCATE (ri_data%pgrid_1)
1553 : END IF
1554 136 : IF (ASSOCIATED(ri_data%pgrid_2)) THEN
1555 136 : CALL dbt_pgrid_destroy(ri_data%pgrid_2)
1556 136 : DEALLOCATE (ri_data%pgrid_2)
1557 : END IF
1558 136 : IF (ASSOCIATED(ri_data%pgrid_2d)) THEN
1559 136 : CALL dbt_pgrid_destroy(ri_data%pgrid_2d)
1560 136 : DEALLOCATE (ri_data%pgrid_2d)
1561 : END IF
1562 :
1563 136 : CALL distribution_3d_destroy(ri_data%dist_3d)
1564 136 : CALL dbt_distribution_destroy(ri_data%dist)
1565 :
1566 136 : DEALLOCATE (ri_data%bsizes_RI)
1567 136 : DEALLOCATE (ri_data%bsizes_AO)
1568 136 : DEALLOCATE (ri_data%bsizes_AO_split)
1569 136 : DEALLOCATE (ri_data%bsizes_RI_split)
1570 136 : DEALLOCATE (ri_data%bsizes_AO_fit)
1571 136 : DEALLOCATE (ri_data%bsizes_RI_fit)
1572 :
1573 136 : IF (ri_data%flavor == ri_pmat) THEN
1574 410 : DO i_mem = 1, ri_data%n_mem
1575 1162 : DO j_mem = 1, ri_data%n_mem_RI
1576 1044 : CALL dealloc_containers(ri_data%store_3c(i_mem, j_mem), unused)
1577 : END DO
1578 : END DO
1579 :
1580 1688 : DO j = 1, SIZE(ri_data%t_3c_int_ctr_1, 2)
1581 3258 : DO i = 1, SIZE(ri_data%t_3c_int_ctr_1, 1)
1582 3140 : CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
1583 : END DO
1584 : END DO
1585 1688 : DEALLOCATE (ri_data%t_3c_int_ctr_1)
1586 :
1587 236 : DO j = 1, SIZE(ri_data%t_3c_int_ctr_2, 2)
1588 354 : DO i = 1, SIZE(ri_data%t_3c_int_ctr_2, 1)
1589 236 : CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
1590 : END DO
1591 : END DO
1592 236 : DEALLOCATE (ri_data%t_3c_int_ctr_2)
1593 :
1594 236 : DO j = 1, SIZE(ri_data%t_3c_int_ctr_3, 2)
1595 354 : DO i = 1, SIZE(ri_data%t_3c_int_ctr_3, 1)
1596 236 : CALL dbt_destroy(ri_data%t_3c_int_ctr_3(i, j))
1597 : END DO
1598 : END DO
1599 236 : DEALLOCATE (ri_data%t_3c_int_ctr_3)
1600 :
1601 296 : DO j = 1, SIZE(ri_data%t_2c_int, 2)
1602 474 : DO i = 1, SIZE(ri_data%t_2c_int, 1)
1603 356 : CALL dbt_destroy(ri_data%t_2c_int(i, j))
1604 : END DO
1605 : END DO
1606 296 : DEALLOCATE (ri_data%t_2c_int)
1607 :
1608 1688 : DO j = 1, SIZE(ri_data%rho_ao_t, 2)
1609 3520 : DO i = 1, SIZE(ri_data%rho_ao_t, 1)
1610 3402 : CALL dbt_destroy(ri_data%rho_ao_t(i, j))
1611 : END DO
1612 : END DO
1613 1950 : DEALLOCATE (ri_data%rho_ao_t)
1614 :
1615 1688 : DO j = 1, SIZE(ri_data%ks_t, 2)
1616 3520 : DO i = 1, SIZE(ri_data%ks_t, 1)
1617 3402 : CALL dbt_destroy(ri_data%ks_t(i, j))
1618 : END DO
1619 : END DO
1620 1950 : DEALLOCATE (ri_data%ks_t)
1621 :
1622 0 : DEALLOCATE (ri_data%starts_array_mem_block, ri_data%ends_array_mem_block, &
1623 118 : ri_data%starts_array_mem, ri_data%ends_array_mem)
1624 0 : DEALLOCATE (ri_data%starts_array_RI_mem_block, ri_data%ends_array_RI_mem_block, &
1625 118 : ri_data%starts_array_RI_mem, ri_data%ends_array_RI_mem)
1626 :
1627 870 : DEALLOCATE (ri_data%blk_indices)
1628 118 : DEALLOCATE (ri_data%store_3c)
1629 18 : ELSE IF (ri_data%flavor == ri_mo) THEN
1630 18 : CALL dbt_destroy(ri_data%t_3c_int_ctr_1(1, 1))
1631 18 : CALL dbt_destroy(ri_data%t_3c_int_ctr_2(1, 1))
1632 36 : DEALLOCATE (ri_data%t_3c_int_ctr_1)
1633 36 : DEALLOCATE (ri_data%t_3c_int_ctr_2)
1634 :
1635 40 : DO ispin = 1, SIZE(ri_data%t_3c_int_mo, 1)
1636 22 : CALL dbt_destroy(ri_data%t_3c_int_mo(ispin, 1, 1))
1637 22 : CALL dbt_destroy(ri_data%t_3c_ctr_RI(ispin, 1, 1))
1638 22 : CALL dbt_destroy(ri_data%t_3c_ctr_KS(ispin, 1, 1))
1639 40 : CALL dbt_destroy(ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
1640 : END DO
1641 54 : DO ispin = 1, 2
1642 54 : CALL dbt_destroy(ri_data%t_2c_int(ispin, 1))
1643 : END DO
1644 54 : DEALLOCATE (ri_data%t_2c_int)
1645 40 : DEALLOCATE (ri_data%t_3c_int_mo)
1646 40 : DEALLOCATE (ri_data%t_3c_ctr_RI)
1647 40 : DEALLOCATE (ri_data%t_3c_ctr_KS)
1648 40 : DEALLOCATE (ri_data%t_3c_ctr_KS_copy)
1649 : END IF
1650 :
1651 332 : DO j = 1, SIZE(ri_data%t_2c_inv, 2)
1652 528 : DO i = 1, SIZE(ri_data%t_2c_inv, 1)
1653 392 : CALL dbt_destroy(ri_data%t_2c_inv(i, j))
1654 : END DO
1655 : END DO
1656 332 : DEALLOCATE (ri_data%t_2c_inv)
1657 :
1658 332 : DO j = 1, SIZE(ri_data%t_2c_pot, 2)
1659 528 : DO i = 1, SIZE(ri_data%t_2c_pot, 1)
1660 392 : CALL dbt_destroy(ri_data%t_2c_pot(i, j))
1661 : END DO
1662 : END DO
1663 332 : DEALLOCATE (ri_data%t_2c_pot)
1664 :
1665 136 : IF (ALLOCATED(ri_data%kp_mat_2c_pot)) THEN
1666 1572 : DO j = 1, SIZE(ri_data%kp_mat_2c_pot, 2)
1667 3084 : DO i = 1, SIZE(ri_data%kp_mat_2c_pot, 1)
1668 3024 : CALL dbcsr_release(ri_data%kp_mat_2c_pot(i, j))
1669 : END DO
1670 : END DO
1671 60 : DEALLOCATE (ri_data%kp_mat_2c_pot)
1672 : END IF
1673 :
1674 136 : IF (ALLOCATED(ri_data%kp_t_3c_int)) THEN
1675 1572 : DO i = 1, SIZE(ri_data%kp_t_3c_int)
1676 1572 : CALL dbt_destroy(ri_data%kp_t_3c_int(i))
1677 : END DO
1678 1572 : DEALLOCATE (ri_data%kp_t_3c_int)
1679 : END IF
1680 :
1681 136 : IF (ALLOCATED(ri_data%rho_ao_t)) THEN
1682 0 : DO j = 1, SIZE(ri_data%rho_ao_t, 2)
1683 0 : DO i = 1, SIZE(ri_data%rho_ao_t, 1)
1684 0 : CALL dbt_destroy(ri_data%rho_ao_t(i, j))
1685 : END DO
1686 : END DO
1687 0 : DEALLOCATE (ri_data%rho_ao_t)
1688 : END IF
1689 :
1690 136 : IF (ALLOCATED(ri_data%ks_t)) THEN
1691 0 : DO j = 1, SIZE(ri_data%ks_t, 2)
1692 0 : DO i = 1, SIZE(ri_data%ks_t, 1)
1693 0 : CALL dbt_destroy(ri_data%ks_t(i, j))
1694 : END DO
1695 : END DO
1696 0 : DEALLOCATE (ri_data%ks_t)
1697 : END IF
1698 :
1699 136 : IF (ALLOCATED(ri_data%iatom_to_subgroup)) THEN
1700 180 : DO i = 1, SIZE(ri_data%iatom_to_subgroup)
1701 180 : DEALLOCATE (ri_data%iatom_to_subgroup(i)%array)
1702 : END DO
1703 60 : DEALLOCATE (ri_data%iatom_to_subgroup)
1704 : END IF
1705 :
1706 136 : CALL timestop(handle)
1707 136 : END SUBROUTINE hfx_ri_release
1708 :
1709 : ! **************************************************************************************************
1710 : !> \brief - This routine allocates and initializes the basis_info and basis_parameter types
1711 : !> \param basis_parameter ...
1712 : !> \param basis_info ...
1713 : !> \param qs_kind_set ...
1714 : !> \param basis_type ...
1715 : !> \par History
1716 : !> 07.2011 refactored
1717 : ! **************************************************************************************************
1718 2186 : SUBROUTINE hfx_create_basis_types(basis_parameter, basis_info, qs_kind_set, &
1719 : basis_type)
1720 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1721 : TYPE(hfx_basis_info_type) :: basis_info
1722 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1723 : CHARACTER(LEN=*) :: basis_type
1724 :
1725 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create_basis_types'
1726 :
1727 : INTEGER :: co_counter, handle, i, ikind, ipgf, iset, j, k, la, max_am_kind, max_coeff, &
1728 : max_nsgfl, max_pgf, max_pgf_kind, max_set, nkind, nl_count, nset, nseta, offset_a, &
1729 : offset_a1, s_offset_nl_a, sgfa, so_counter
1730 2186 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell
1731 2186 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl_a
1732 2186 : REAL(dp), DIMENSION(:, :), POINTER :: sphi_a
1733 : TYPE(gto_basis_set_type), POINTER :: orb_basis_a
1734 :
1735 2186 : CALL timeset(routineN, handle)
1736 :
1737 : ! BASIS parameter
1738 2186 : nkind = SIZE(qs_kind_set, 1)
1739 : !
1740 10556 : ALLOCATE (basis_parameter(nkind))
1741 6184 : max_set = 0
1742 6184 : DO ikind = 1, nkind
1743 3998 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_a, basis_type=basis_type)
1744 : CALL get_qs_kind_set(qs_kind_set, &
1745 : maxsgf=basis_info%max_sgf, &
1746 : maxnset=basis_info%max_set, &
1747 : maxlgto=basis_info%max_am, &
1748 3998 : basis_type=basis_type)
1749 3998 : IF (basis_info%max_set < max_set) CPABORT("UNEXPECTED MAX_SET")
1750 3998 : max_set = MAX(max_set, basis_info%max_set)
1751 : CALL get_gto_basis_set(gto_basis_set=orb_basis_a, &
1752 : lmax=basis_parameter(ikind)%lmax, &
1753 : lmin=basis_parameter(ikind)%lmin, &
1754 : npgf=basis_parameter(ikind)%npgf, &
1755 : nset=basis_parameter(ikind)%nset, &
1756 : zet=basis_parameter(ikind)%zet, &
1757 : nsgf_set=basis_parameter(ikind)%nsgf, &
1758 : first_sgf=basis_parameter(ikind)%first_sgf, &
1759 : sphi=basis_parameter(ikind)%sphi, &
1760 : gcc=basis_parameter(ikind)%gcc, &
1761 : nsgf=basis_parameter(ikind)%nsgf_total, &
1762 : l=basis_parameter(ikind)%nl, &
1763 : nshell=basis_parameter(ikind)%nshell, &
1764 : set_radius=basis_parameter(ikind)%set_radius, &
1765 : pgf_radius=basis_parameter(ikind)%pgf_radius, &
1766 6184 : kind_radius=basis_parameter(ikind)%kind_radius)
1767 : END DO
1768 6184 : DO ikind = 1, nkind
1769 15992 : ALLOCATE (basis_parameter(ikind)%nsgfl(0:basis_info%max_am, max_set))
1770 49604 : basis_parameter(ikind)%nsgfl = 0
1771 3998 : nset = basis_parameter(ikind)%nset
1772 3998 : nshell => basis_parameter(ikind)%nshell
1773 17548 : DO iset = 1, nset
1774 45580 : DO i = 0, basis_info%max_am
1775 30218 : nl_count = 0
1776 70570 : DO j = 1, nshell(iset)
1777 70570 : IF (basis_parameter(ikind)%nl(j, iset) == i) nl_count = nl_count + 1
1778 : END DO
1779 41582 : basis_parameter(ikind)%nsgfl(i, iset) = nl_count
1780 : END DO
1781 : END DO
1782 : END DO
1783 :
1784 : max_nsgfl = 0
1785 : max_pgf = 0
1786 6184 : DO ikind = 1, nkind
1787 3998 : max_coeff = 0
1788 3998 : max_am_kind = 0
1789 3998 : max_pgf_kind = 0
1790 3998 : npgfa => basis_parameter(ikind)%npgf
1791 3998 : nseta = basis_parameter(ikind)%nset
1792 3998 : nl_a => basis_parameter(ikind)%nsgfl
1793 3998 : la_max => basis_parameter(ikind)%lmax
1794 3998 : la_min => basis_parameter(ikind)%lmin
1795 15362 : DO iset = 1, nseta
1796 11364 : max_pgf_kind = MAX(max_pgf_kind, npgfa(iset))
1797 : max_pgf = MAX(max_pgf, npgfa(iset))
1798 29278 : DO la = la_min(iset), la_max(iset)
1799 13916 : max_nsgfl = MAX(max_nsgfl, nl_a(la, iset))
1800 13916 : max_coeff = MAX(max_coeff, nso(la)*nl_a(la, iset)*nco(la))
1801 25280 : max_am_kind = MAX(max_am_kind, la)
1802 : END DO
1803 : END DO
1804 23988 : ALLOCATE (basis_parameter(ikind)%sphi_ext(max_coeff, 0:max_am_kind, max_pgf_kind, nseta))
1805 2213272 : basis_parameter(ikind)%sphi_ext = 0.0_dp
1806 : END DO
1807 :
1808 6184 : DO ikind = 1, nkind
1809 3998 : sphi_a => basis_parameter(ikind)%sphi
1810 3998 : nseta = basis_parameter(ikind)%nset
1811 3998 : la_max => basis_parameter(ikind)%lmax
1812 3998 : la_min => basis_parameter(ikind)%lmin
1813 3998 : npgfa => basis_parameter(ikind)%npgf
1814 3998 : first_sgfa => basis_parameter(ikind)%first_sgf
1815 3998 : nl_a => basis_parameter(ikind)%nsgfl
1816 17548 : DO iset = 1, nseta
1817 11364 : sgfa = first_sgfa(1, iset)
1818 36866 : DO ipgf = 1, npgfa(iset)
1819 21504 : offset_a1 = (ipgf - 1)*ncoset(la_max(iset))
1820 21504 : s_offset_nl_a = 0
1821 60930 : DO la = la_min(iset), la_max(iset)
1822 28062 : offset_a = offset_a1 + ncoset(la - 1)
1823 : co_counter = 0
1824 28062 : co_counter = co_counter + 1
1825 28062 : so_counter = 0
1826 87312 : DO k = sgfa + s_offset_nl_a, sgfa + s_offset_nl_a + nso(la)*nl_a(la, iset) - 1
1827 249700 : DO i = offset_a + 1, offset_a + nco(la)
1828 162388 : so_counter = so_counter + 1
1829 221638 : basis_parameter(ikind)%sphi_ext(so_counter, la, ipgf, iset) = sphi_a(i, k)
1830 : END DO
1831 : END DO
1832 49566 : s_offset_nl_a = s_offset_nl_a + nso(la)*(nl_a(la, iset))
1833 : END DO
1834 : END DO
1835 : END DO
1836 : END DO
1837 :
1838 2186 : CALL timestop(handle)
1839 :
1840 2186 : END SUBROUTINE hfx_create_basis_types
1841 :
1842 : ! **************************************************************************************************
1843 : !> \brief ...
1844 : !> \param basis_parameter ...
1845 : ! **************************************************************************************************
1846 2186 : SUBROUTINE hfx_release_basis_types(basis_parameter)
1847 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
1848 :
1849 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_release_basis_types'
1850 :
1851 : INTEGER :: handle, i
1852 :
1853 2186 : CALL timeset(routineN, handle)
1854 :
1855 : !! BASIS parameter
1856 6184 : DO i = 1, SIZE(basis_parameter)
1857 3998 : DEALLOCATE (basis_parameter(i)%nsgfl)
1858 6184 : DEALLOCATE (basis_parameter(i)%sphi_ext)
1859 : END DO
1860 2186 : DEALLOCATE (basis_parameter)
1861 2186 : CALL timestop(handle)
1862 :
1863 2186 : END SUBROUTINE hfx_release_basis_types
1864 :
1865 : ! **************************************************************************************************
1866 : !> \brief - Parses the memory section
1867 : !> \param memory_parameter ...
1868 : !> \param hf_sub_section ...
1869 : !> \param storage_id ...
1870 : !> \param i_thread ...
1871 : !> \param n_threads ...
1872 : !> \param para_env ...
1873 : !> \param irep ...
1874 : !> \param skip_disk ...
1875 : !> \param skip_in_core_forces ...
1876 : ! **************************************************************************************************
1877 2488 : SUBROUTINE parse_memory_section(memory_parameter, hf_sub_section, storage_id, &
1878 : i_thread, n_threads, para_env, irep, skip_disk, skip_in_core_forces)
1879 : TYPE(hfx_memory_type) :: memory_parameter
1880 : TYPE(section_vals_type), POINTER :: hf_sub_section
1881 : INTEGER, INTENT(OUT), OPTIONAL :: storage_id
1882 : INTEGER, INTENT(IN), OPTIONAL :: i_thread, n_threads
1883 : TYPE(mp_para_env_type), OPTIONAL :: para_env
1884 : INTEGER, INTENT(IN), OPTIONAL :: irep
1885 : LOGICAL, INTENT(IN) :: skip_disk, skip_in_core_forces
1886 :
1887 : CHARACTER(LEN=512) :: error_msg
1888 : CHARACTER(LEN=default_path_length) :: char_val, filename, orig_wd
1889 : INTEGER :: int_val, stat
1890 : LOGICAL :: check, logic_val
1891 : REAL(dp) :: real_val
1892 :
1893 : check = (PRESENT(storage_id) .EQV. PRESENT(i_thread)) .AND. &
1894 : (PRESENT(storage_id) .EQV. PRESENT(n_threads)) .AND. &
1895 : (PRESENT(storage_id) .EQV. PRESENT(para_env)) .AND. &
1896 2488 : (PRESENT(storage_id) .EQV. PRESENT(irep))
1897 0 : CPASSERT(check)
1898 :
1899 : ! Memory Storage
1900 2488 : CALL section_vals_val_get(hf_sub_section, "MAX_MEMORY", i_val=int_val)
1901 2488 : memory_parameter%max_memory = int_val
1902 2488 : memory_parameter%max_compression_counter = int_val*1024_int_8*128_int_8
1903 2488 : CALL section_vals_val_get(hf_sub_section, "EPS_STORAGE", r_val=real_val)
1904 2488 : memory_parameter%eps_storage_scaling = real_val
1905 2488 : IF (int_val == 0) THEN
1906 20 : memory_parameter%do_all_on_the_fly = .TRUE.
1907 : ELSE
1908 2468 : memory_parameter%do_all_on_the_fly = .FALSE.
1909 : END IF
1910 2488 : memory_parameter%cache_size = CACHE_SIZE
1911 2488 : memory_parameter%bits_max_val = BITS_MAX_VAL
1912 2488 : memory_parameter%actual_memory_usage = 1
1913 2488 : IF (.NOT. skip_in_core_forces) THEN
1914 1488 : CALL section_vals_val_get(hf_sub_section, "TREAT_FORCES_IN_CORE", l_val=logic_val)
1915 1488 : memory_parameter%treat_forces_in_core = logic_val
1916 : END IF
1917 :
1918 : ! ** IF MAX_MEM == 0 overwrite this flag to false
1919 2488 : IF (memory_parameter%do_all_on_the_fly) memory_parameter%treat_forces_in_core = .FALSE.
1920 :
1921 : ! Disk Storage
1922 2488 : IF (.NOT. skip_disk) THEN
1923 1488 : memory_parameter%actual_memory_usage_disk = 1
1924 1488 : CALL section_vals_val_get(hf_sub_section, "MAX_DISK_SPACE", i_val=int_val)
1925 1488 : memory_parameter%max_compression_counter_disk = int_val*1024_int_8*128_int_8
1926 1488 : IF (int_val == 0) THEN
1927 1482 : memory_parameter%do_disk_storage = .FALSE.
1928 : ELSE
1929 6 : memory_parameter%do_disk_storage = .TRUE.
1930 : END IF
1931 1488 : CALL section_vals_val_get(hf_sub_section, "STORAGE_LOCATION", c_val=char_val)
1932 1488 : CALL compress(char_val, .TRUE.)
1933 : !! Add ending / if necessary
1934 :
1935 1488 : IF (SCAN(char_val, "/", .TRUE.) /= LEN_TRIM(char_val)) THEN
1936 1488 : WRITE (filename, '(A,A)') TRIM(char_val), "/"
1937 1488 : CALL compress(filename)
1938 : ELSE
1939 0 : filename = TRIM(char_val)
1940 : END IF
1941 1488 : CALL compress(filename, .TRUE.)
1942 :
1943 : !! quickly check if we can write on storage_location
1944 1488 : CALL m_getcwd(orig_wd)
1945 1488 : CALL m_chdir(TRIM(filename), stat)
1946 1488 : IF (stat /= 0) THEN
1947 0 : WRITE (error_msg, '(A,A,A)') "Request for disk storage failed due to unknown error while writing to ", &
1948 0 : TRIM(filename), ". Please check STORAGE_LOCATION"
1949 0 : CPABORT(error_msg)
1950 : END IF
1951 1488 : CALL m_chdir(orig_wd, stat)
1952 :
1953 1488 : memory_parameter%storage_location = filename
1954 1488 : CALL compress(memory_parameter%storage_location, .TRUE.)
1955 : ELSE
1956 1000 : memory_parameter%do_disk_storage = .FALSE.
1957 : END IF
1958 2488 : IF (PRESENT(storage_id)) THEN
1959 1488 : storage_id = (irep - 1)*para_env%num_pe*n_threads + para_env%mepos*n_threads + i_thread - 1
1960 : END IF
1961 2488 : END SUBROUTINE parse_memory_section
1962 :
1963 : ! **************************************************************************************************
1964 : !> \brief - This routine deallocates all data structures
1965 : !> \param x_data contains all relevant data structures for hfx runs
1966 : !> \par History
1967 : !> 09.2007 created [Manuel Guidon]
1968 : !> \author Manuel Guidon
1969 : ! **************************************************************************************************
1970 1478 : SUBROUTINE hfx_release(x_data)
1971 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1972 :
1973 : INTEGER :: i, i_thread, irep, n_rep_hf, n_threads
1974 : TYPE(cp_logger_type), POINTER :: logger
1975 : TYPE(hfx_type), POINTER :: actual_x_data
1976 :
1977 : !! There might be 2 hf sections
1978 :
1979 1478 : n_rep_hf = x_data(1, 1)%n_rep_hf
1980 1478 : n_threads = SIZE(x_data, 2)
1981 :
1982 1478 : IF (x_data(1, 1)%potential_parameter%potential_type == do_potential_truncated .OR. &
1983 : x_data(1, 1)%potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
1984 376 : init_t_c_g0_lmax = -1
1985 376 : CALL free_C0()
1986 : END IF
1987 2956 : DO i_thread = 1, n_threads
1988 4444 : DO irep = 1, n_rep_hf
1989 1488 : actual_x_data => x_data(irep, i_thread)
1990 1488 : DEALLOCATE (actual_x_data%neighbor_cells)
1991 1488 : DEALLOCATE (actual_x_data%distribution_energy)
1992 1488 : DEALLOCATE (actual_x_data%distribution_forces)
1993 :
1994 1488 : IF (actual_x_data%load_balance_parameter%blocks_initialized) THEN
1995 1366 : DEALLOCATE (actual_x_data%blocks)
1996 1366 : IF (i_thread == 1) THEN
1997 1366 : DEALLOCATE (actual_x_data%pmax_block)
1998 : END IF
1999 : END IF
2000 :
2001 1488 : IF (i_thread == 1) THEN
2002 1488 : DEALLOCATE (actual_x_data%atomic_pair_list)
2003 1488 : DEALLOCATE (actual_x_data%atomic_pair_list_forces)
2004 : END IF
2005 :
2006 1488 : IF (actual_x_data%screening_parameter%do_initial_p_screening .OR. &
2007 : actual_x_data%screening_parameter%do_p_screening_forces) THEN
2008 1458 : IF (i_thread == 1) THEN
2009 1458 : DEALLOCATE (actual_x_data%pmax_atom)
2010 5852 : DO i = 1, SIZE(actual_x_data%initial_p)
2011 5852 : DEALLOCATE (actual_x_data%initial_p(i)%p_kind)
2012 : END DO
2013 1458 : DEALLOCATE (actual_x_data%initial_p)
2014 :
2015 1458 : DEALLOCATE (actual_x_data%pmax_atom_forces)
2016 5852 : DO i = 1, SIZE(actual_x_data%initial_p_forces)
2017 5852 : DEALLOCATE (actual_x_data%initial_p_forces(i)%p_kind)
2018 : END DO
2019 1458 : DEALLOCATE (actual_x_data%initial_p_forces)
2020 : END IF
2021 1458 : DEALLOCATE (actual_x_data%map_atom_to_kind_atom)
2022 : END IF
2023 1488 : IF (i_thread == 1) THEN
2024 1488 : DEALLOCATE (actual_x_data%is_assoc_atomic_block)
2025 1488 : DEALLOCATE (actual_x_data%atomic_block_offset)
2026 1488 : DEALLOCATE (actual_x_data%set_offset)
2027 1488 : DEALLOCATE (actual_x_data%block_offset)
2028 : END IF
2029 :
2030 : !! BASIS parameter
2031 1488 : CALL hfx_release_basis_types(actual_x_data%basis_parameter)
2032 :
2033 : !MK Release libint and libderiv data structure
2034 1488 : CALL cp_libint_cleanup_eri(actual_x_data%lib)
2035 1488 : CALL cp_libint_cleanup_eri1(actual_x_data%lib_deriv)
2036 1488 : CALL cp_libint_static_cleanup()
2037 :
2038 : !! Deallocate containers
2039 1488 : CALL dealloc_containers(actual_x_data%store_ints, actual_x_data%memory_parameter%actual_memory_usage)
2040 1488 : CALL dealloc_containers(actual_x_data%store_forces, actual_x_data%memory_parameter%actual_memory_usage)
2041 :
2042 : !! Deallocate containers
2043 : CALL hfx_init_container(actual_x_data%store_ints%maxval_container_disk, &
2044 : actual_x_data%memory_parameter%actual_memory_usage_disk, &
2045 1488 : .FALSE.)
2046 1488 : IF (actual_x_data%memory_parameter%do_disk_storage) THEN
2047 6 : CALL close_file(unit_number=actual_x_data%store_ints%maxval_container_disk%unit, file_status="DELETE")
2048 : END IF
2049 1488 : DEALLOCATE (actual_x_data%store_ints%maxval_container_disk%first)
2050 1488 : DEALLOCATE (actual_x_data%store_ints%maxval_container_disk)
2051 :
2052 96720 : DO i = 1, 64
2053 : CALL hfx_init_container(actual_x_data%store_ints%integral_containers_disk(i), &
2054 : actual_x_data%memory_parameter%actual_memory_usage_disk, &
2055 95232 : .FALSE.)
2056 95232 : IF (actual_x_data%memory_parameter%do_disk_storage) THEN
2057 384 : CALL close_file(unit_number=actual_x_data%store_ints%integral_containers_disk(i)%unit, file_status="DELETE")
2058 : END IF
2059 96720 : DEALLOCATE (actual_x_data%store_ints%integral_containers_disk(i)%first)
2060 : END DO
2061 1488 : DEALLOCATE (actual_x_data%store_ints%integral_containers_disk)
2062 :
2063 : ! ** screening functions
2064 1488 : IF (actual_x_data%screen_funct_is_initialized) THEN
2065 1366 : DEALLOCATE (actual_x_data%screen_funct_coeffs_set)
2066 1366 : DEALLOCATE (actual_x_data%screen_funct_coeffs_kind)
2067 1366 : DEALLOCATE (actual_x_data%pair_dist_radii_pgf)
2068 1366 : DEALLOCATE (actual_x_data%screen_funct_coeffs_pgf)
2069 1366 : actual_x_data%screen_funct_is_initialized = .FALSE.
2070 : END IF
2071 :
2072 : ! ** maps
2073 1488 : IF (ASSOCIATED(actual_x_data%map_atoms_to_cpus)) THEN
2074 4096 : DO i = 1, SIZE(actual_x_data%map_atoms_to_cpus)
2075 2730 : DEALLOCATE (actual_x_data%map_atoms_to_cpus(i)%iatom_list)
2076 4096 : DEALLOCATE (actual_x_data%map_atoms_to_cpus(i)%jatom_list)
2077 : END DO
2078 1366 : DEALLOCATE (actual_x_data%map_atoms_to_cpus)
2079 : END IF
2080 :
2081 1488 : IF (actual_x_data%do_hfx_ri) THEN
2082 114 : CALL hfx_ri_release(actual_x_data%ri_data)
2083 114 : IF (ASSOCIATED(actual_x_data%ri_data%ri_section)) THEN
2084 114 : logger => cp_get_default_logger()
2085 : CALL cp_print_key_finished_output(actual_x_data%ri_data%unit_nr_dbcsr, logger, actual_x_data%ri_data%ri_section, &
2086 114 : "PRINT%RI_INFO")
2087 : END IF
2088 114 : IF (ASSOCIATED(actual_x_data%ri_data%hfx_section)) THEN
2089 114 : logger => cp_get_default_logger()
2090 : CALL cp_print_key_finished_output(actual_x_data%ri_data%unit_nr, logger, actual_x_data%ri_data%hfx_section, &
2091 114 : "HF_INFO")
2092 : END IF
2093 114 : DEALLOCATE (actual_x_data%ri_data)
2094 : END IF
2095 :
2096 : ! ACE cleanup — just reset scalars, ace_W is managed elsewhere
2097 1488 : actual_x_data%use_ace = .FALSE.
2098 1488 : actual_x_data%ace_is_built = .FALSE.
2099 2966 : actual_x_data%ace_build_counter = 0
2100 : END DO
2101 :
2102 : END DO
2103 :
2104 1478 : DEALLOCATE (x_data)
2105 1478 : END SUBROUTINE hfx_release
2106 :
2107 : ! **************************************************************************************************
2108 : !> \brief - This routine computes the neighbor cells that are taken into account
2109 : !> in periodic runs
2110 : !> \param x_data contains all relevant data structures for hfx runs
2111 : !> \param pbc_shells number of shells taken into account
2112 : !> \param cell cell
2113 : !> \param i_thread current thread ID
2114 : !> \param nkp_grid ...
2115 : !> \par History
2116 : !> 09.2007 created [Manuel Guidon]
2117 : !> \author Manuel Guidon
2118 : ! **************************************************************************************************
2119 10281 : SUBROUTINE hfx_create_neighbor_cells(x_data, pbc_shells, cell, i_thread, nkp_grid)
2120 : TYPE(hfx_type), POINTER :: x_data
2121 : INTEGER, INTENT(INOUT) :: pbc_shells
2122 : TYPE(cell_type), POINTER :: cell
2123 : INTEGER, INTENT(IN) :: i_thread
2124 : INTEGER, DIMENSION(3), OPTIONAL :: nkp_grid
2125 :
2126 : CHARACTER(LEN=512) :: error_msg
2127 : CHARACTER(LEN=64) :: char_nshells
2128 : INTEGER :: i, idx, ikind, ipgf, iset, ishell, j, jkind, jpgf, jset, jshell, k, kshell, l, &
2129 : m(3), max_shell, nkp(3), nseta, nsetb, perd(3), total_number_of_cells, ub, ub_max
2130 10281 : INTEGER, DIMENSION(:), POINTER :: la_max, lb_max, npgfa, npgfb
2131 : LOGICAL :: do_kpoints, image_cell_found, &
2132 : nothing_more_to_add
2133 : REAL(dp) :: cross_product(3), dist_min, distance(14), l_min, normal(3, 6), P(3, 14), &
2134 : plane_vector(3, 2), point_in_plane(3), r(3), R1, R_max, R_max_stress, s(3), x, y, z, Zeta1
2135 10281 : REAL(dp), DIMENSION(:, :), POINTER :: zeta, zetb
2136 10281 : TYPE(hfx_cell_type), ALLOCATABLE, DIMENSION(:) :: tmp_neighbor_cells
2137 :
2138 10281 : total_number_of_cells = 0
2139 :
2140 41124 : nkp = 1
2141 10281 : IF (PRESENT(nkp_grid)) nkp = nkp_grid
2142 40944 : do_kpoints = ANY(nkp > 1)
2143 :
2144 : ! ** Check some settings
2145 10281 : IF (i_thread == 1) THEN
2146 : IF (x_data%potential_parameter%potential_type /= do_potential_truncated .AND. &
2147 : x_data%potential_parameter%potential_type /= do_potential_short .AND. &
2148 448 : x_data%potential_parameter%potential_type /= do_potential_mix_cl_trunc .AND. &
2149 : x_data%potential_parameter%potential_type /= do_potential_id) THEN
2150 : CALL cp_warn(__LOCATION__, &
2151 : "Periodic Hartree Fock calculation requested without use "// &
2152 : "of a truncated or shortrange potential. This may lead to unphysical total energies. "// &
2153 104 : "Use a truncated potential to avoid possible problems.")
2154 344 : ELSE IF (x_data%potential_parameter%potential_type /= do_potential_id) THEN
2155 : !If k-points, use the Born-von Karman super cell as reference
2156 : l_min = MIN(REAL(nkp(1), dp)*plane_distance(1, 0, 0, cell), &
2157 : REAL(nkp(2), dp)*plane_distance(0, 1, 0, cell), &
2158 316 : REAL(nkp(3), dp)*plane_distance(0, 0, 1, cell))
2159 316 : l_min = 0.5_dp*l_min
2160 316 : IF (x_data%potential_parameter%cutoff_radius >= l_min) THEN
2161 38 : IF (.NOT. do_kpoints) THEN
2162 : WRITE (error_msg, "(A,F6.3,A,F6.3,A)") &
2163 : "Periodic Hartree Fock calculation requested with the use "// &
2164 : "of a truncated or shortrange potential. "// &
2165 38 : "The cutoff radius (", x_data%potential_parameter%cutoff_radius*a_bohr*1e+10_dp, &
2166 38 : " A) is larger than half the minimal cell dimension (", &
2167 38 : l_min*a_bohr*1e+10_dp, " A). This may lead to unphysical "// &
2168 : "total energies. Reduce the cutoff radius in order to avoid "// &
2169 76 : "possible problems."
2170 : ELSE
2171 : WRITE (error_msg, "(A,F6.3,A,F6.3,A)") &
2172 : "K-point Hartree-Fock calculation requested with the use of a "// &
2173 0 : "truncated or shortrange potential. The cutoff radius (", &
2174 0 : x_data%potential_parameter%cutoff_radius*a_bohr*1e+10_dp, &
2175 0 : " A) is larger than half the minimal Born-von Karman supercell dimension (", &
2176 0 : l_min*a_bohr*1e+10_dp, " A). This may lead "// &
2177 : "to unphysical total energies. Reduce the cutoff radius or increase "// &
2178 0 : "the number of K-points in order to avoid possible problems."
2179 : END IF
2180 38 : CALL cp_warn(__LOCATION__, error_msg)
2181 : END IF
2182 : END IF
2183 : END IF
2184 :
2185 17717 : SELECT CASE (x_data%potential_parameter%potential_type)
2186 : CASE (do_potential_truncated, do_potential_mix_cl_trunc, do_potential_short)
2187 7436 : R_max = 0.0_dp
2188 20546 : DO ikind = 1, SIZE(x_data%basis_parameter)
2189 13110 : la_max => x_data%basis_parameter(ikind)%lmax
2190 13110 : zeta => x_data%basis_parameter(ikind)%zet
2191 13110 : nseta = x_data%basis_parameter(ikind)%nset
2192 13110 : npgfa => x_data%basis_parameter(ikind)%npgf
2193 45124 : DO jkind = 1, SIZE(x_data%basis_parameter)
2194 24578 : lb_max => x_data%basis_parameter(jkind)%lmax
2195 24578 : zetb => x_data%basis_parameter(jkind)%zet
2196 24578 : nsetb = x_data%basis_parameter(jkind)%nset
2197 24578 : npgfb => x_data%basis_parameter(jkind)%npgf
2198 102590 : DO iset = 1, nseta
2199 280060 : DO jset = 1, nsetb
2200 582586 : DO ipgf = 1, npgfa(iset)
2201 1123098 : DO jpgf = 1, npgfb(jset)
2202 605414 : Zeta1 = zeta(ipgf, iset) + zetb(jpgf, jset)
2203 : R1 = 1.0_dp/SQRT(Zeta1)*mul_fact(la_max(iset) + lb_max(jset))* &
2204 605414 : SQRT(-LOG(x_data%screening_parameter%eps_schwarz))
2205 932518 : R_max = MAX(R1, R_max)
2206 : END DO
2207 : END DO
2208 : END DO
2209 : END DO
2210 : END DO
2211 : END DO
2212 :
2213 7436 : R_max = 2.0_dp*R_max + x_data%potential_parameter%cutoff_radius
2214 7436 : nothing_more_to_add = .FALSE.
2215 7436 : max_shell = 0
2216 7436 : total_number_of_cells = 0
2217 7436 : ub = 1
2218 7436 : DEALLOCATE (x_data%neighbor_cells)
2219 59488 : ALLOCATE (x_data%neighbor_cells(1))
2220 29744 : x_data%neighbor_cells(1)%cell = 0.0_dp
2221 29744 : x_data%neighbor_cells(1)%cell_r = 0.0_dp
2222 :
2223 : ! ** What follows is kind of a ray tracing algorithm
2224 : ! ** Given a image cell (ishell, jshell, kshell) we try to figure out the
2225 : ! ** shortest distance of this image cell to the basic unit cell (0,0,0), i.e. the point
2226 : ! ** (0.0, 0.0, 0.0)
2227 : ! ** This is achieved by checking the 8 Corners of the cell, and, in addition, the shortest distance
2228 : ! ** to all 6 faces. The faces are only taken into account if the penetration point of the normal
2229 : ! ** to the plane defined by a face lies within this face.
2230 : ! ** This is very fast, because no trigonometric functions are being used
2231 : ! ** The points are defined as follows
2232 : ! **
2233 : ! **
2234 : ! ** _________________________
2235 : ! ** /P4____________________P8/|
2236 : ! ** / / ___________________/ / |
2237 : ! ** / / /| | / / | z
2238 : ! ** / / / | | / / . | /|\ _ y
2239 : ! ** / / /| | | / / /| | | /|
2240 : ! ** / / / | | | / / / | | | /
2241 : ! ** / / / | | | / / /| | | | /
2242 : ! ** / /_/___| | |__________/ / / | | | |/
2243 : ! ** /P2______| | |_________P6/ / | | | ----------> x
2244 : ! ** | _______| | |_________| | | | | |
2245 : ! ** | | | | | |________________| | |
2246 : ! ** | | | |P3___________________P7 |
2247 : ! ** | | | / / _________________ / /
2248 : ! ** | | | / / / | | |/ / /
2249 : ! ** | | | / / / | | | / /
2250 : ! ** | | |/ / / | | |/ /
2251 : ! ** | | | / / | | ' /
2252 : ! ** | | |/_/_______________| | /
2253 : ! ** | |____________________| | /
2254 : ! ** |P1_____________________P5/
2255 : ! **
2256 : ! **
2257 :
2258 : DO WHILE (.NOT. nothing_more_to_add)
2259 : ! Calculate distances to the eight points P1 to P8
2260 30424 : image_cell_found = .FALSE.
2261 1307872 : ALLOCATE (tmp_neighbor_cells(1:ub))
2262 1034056 : DO i = 1, ub - 1
2263 1034056 : tmp_neighbor_cells(i) = x_data%neighbor_cells(i)
2264 : END DO
2265 30424 : ub_max = (2*max_shell + 1)**3
2266 30424 : DEALLOCATE (x_data%neighbor_cells)
2267 4734460 : ALLOCATE (x_data%neighbor_cells(1:ub_max))
2268 1034056 : DO i = 1, ub - 1
2269 1034056 : x_data%neighbor_cells(i) = tmp_neighbor_cells(i)
2270 : END DO
2271 3487436 : DO i = ub, ub_max
2272 13828048 : x_data%neighbor_cells(i)%cell = 0.0_dp
2273 13858472 : x_data%neighbor_cells(i)%cell_r = 0.0_dp
2274 : END DO
2275 :
2276 30424 : DEALLOCATE (tmp_neighbor_cells)
2277 :
2278 121696 : perd(1:3) = x_data%periodic_parameter%perd(1:3)
2279 :
2280 156092 : DO ishell = -max_shell*perd(1), max_shell*perd(1)
2281 850944 : DO jshell = -max_shell*perd(2), max_shell*perd(2)
2282 5164836 : DO kshell = -max_shell*perd(3), max_shell*perd(3)
2283 4344316 : IF (MAX(ABS(ishell), ABS(jshell), ABS(kshell)) /= max_shell) CYCLE
2284 : idx = 0
2285 8594784 : DO j = 0, 1
2286 5729856 : x = -1.0_dp/2.0_dp + j*1.0_dp
2287 20054496 : DO k = 0, 1
2288 11459712 : y = -1.0_dp/2.0_dp + k*1.0_dp
2289 40108992 : DO l = 0, 1
2290 22919424 : z = -1.0_dp/2.0_dp + l*1.0_dp
2291 22919424 : idx = idx + 1
2292 22919424 : P(1, idx) = x + ishell
2293 22919424 : P(2, idx) = y + jshell
2294 22919424 : P(3, idx) = z + kshell
2295 22919424 : CALL scaled_to_real(r, P(:, idx), cell)
2296 91677696 : distance(idx) = SQRT(SUM(r**2))
2297 103137408 : P(1:3, idx) = r
2298 : END DO
2299 : END DO
2300 : END DO
2301 : ! Now check distance to Faces and only take them into account if the base point lies within quadrilateral
2302 :
2303 : ! Face A (1342) 1 is the reference
2304 2864928 : idx = idx + 1
2305 11459712 : plane_vector(:, 1) = P(:, 3) - P(:, 1)
2306 11459712 : plane_vector(:, 2) = P(:, 2) - P(:, 1)
2307 2864928 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2308 2864928 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2309 2864928 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2310 20054496 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2311 11459712 : point_in_plane = -normal(:, 1)*(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2312 :
2313 2864928 : IF (point_is_in_quadrilateral(P(:, 1), P(:, 3), P(:, 4), P(:, 2), point_in_plane)) THEN
2314 54240 : distance(idx) = ABS(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2315 : ELSE
2316 2810688 : distance(idx) = HUGE(distance(idx))
2317 : END IF
2318 :
2319 : ! Face B (1562) 1 is the reference
2320 2864928 : idx = idx + 1
2321 11459712 : plane_vector(:, 1) = P(:, 2) - P(:, 1)
2322 11459712 : plane_vector(:, 2) = P(:, 5) - P(:, 1)
2323 2864928 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2324 2864928 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2325 2864928 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2326 20054496 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2327 11459712 : point_in_plane = -normal(:, 1)*(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2328 :
2329 2864928 : IF (point_is_in_quadrilateral(P(:, 1), P(:, 5), P(:, 6), P(:, 2), point_in_plane)) THEN
2330 54416 : distance(idx) = ABS(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2331 : ELSE
2332 2810512 : distance(idx) = HUGE(distance(idx))
2333 : END IF
2334 :
2335 : ! Face C (5786) 5 is the reference
2336 2864928 : idx = idx + 1
2337 11459712 : plane_vector(:, 1) = P(:, 7) - P(:, 5)
2338 11459712 : plane_vector(:, 2) = P(:, 6) - P(:, 5)
2339 2864928 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2340 2864928 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2341 2864928 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2342 20054496 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2343 11459712 : point_in_plane = -normal(:, 1)*(normal(1, 1)*P(1, 5) + normal(2, 1)*P(2, 5) + normal(3, 1)*P(3, 5))
2344 :
2345 2864928 : IF (point_is_in_quadrilateral(P(:, 5), P(:, 7), P(:, 8), P(:, 6), point_in_plane)) THEN
2346 54240 : distance(idx) = ABS(normal(1, 1)*P(1, 5) + normal(2, 1)*P(2, 5) + normal(3, 1)*P(3, 5))
2347 : ELSE
2348 2810688 : distance(idx) = HUGE(distance(idx))
2349 : END IF
2350 :
2351 : ! Face D (3784) 3 is the reference
2352 2864928 : idx = idx + 1
2353 11459712 : plane_vector(:, 1) = P(:, 7) - P(:, 3)
2354 11459712 : plane_vector(:, 2) = P(:, 4) - P(:, 3)
2355 2864928 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2356 2864928 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2357 2864928 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2358 20054496 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2359 11459712 : point_in_plane = -normal(:, 1)*(normal(1, 1)*P(1, 3) + normal(2, 1)*P(2, 3) + normal(3, 1)*P(3, 3))
2360 :
2361 2864928 : IF (point_is_in_quadrilateral(P(:, 3), P(:, 7), P(:, 8), P(:, 4), point_in_plane)) THEN
2362 54416 : distance(idx) = ABS(normal(1, 1)*P(1, 3) + normal(2, 1)*P(2, 3) + normal(3, 1)*P(3, 3))
2363 : ELSE
2364 2810512 : distance(idx) = HUGE(distance(idx))
2365 : END IF
2366 :
2367 : ! Face E (2684) 2 is the reference
2368 2864928 : idx = idx + 1
2369 11459712 : plane_vector(:, 1) = P(:, 6) - P(:, 2)
2370 11459712 : plane_vector(:, 2) = P(:, 4) - P(:, 2)
2371 2864928 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2372 2864928 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2373 2864928 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2374 20054496 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2375 11459712 : point_in_plane = -normal(:, 1)*(normal(1, 1)*P(1, 2) + normal(2, 1)*P(2, 2) + normal(3, 1)*P(3, 2))
2376 :
2377 2864928 : IF (point_is_in_quadrilateral(P(:, 2), P(:, 6), P(:, 8), P(:, 4), point_in_plane)) THEN
2378 54220 : distance(idx) = ABS(normal(1, 1)*P(1, 2) + normal(2, 1)*P(2, 2) + normal(3, 1)*P(3, 2))
2379 : ELSE
2380 2810708 : distance(idx) = HUGE(distance(idx))
2381 : END IF
2382 :
2383 : ! Face F (1573) 1 is the reference
2384 2864928 : idx = idx + 1
2385 11459712 : plane_vector(:, 1) = P(:, 5) - P(:, 1)
2386 11459712 : plane_vector(:, 2) = P(:, 3) - P(:, 1)
2387 2864928 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2388 2864928 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2389 2864928 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2390 20054496 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2391 11459712 : point_in_plane = -normal(:, 1)*(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2392 :
2393 2864928 : IF (point_is_in_quadrilateral(P(:, 1), P(:, 5), P(:, 7), P(:, 3), point_in_plane)) THEN
2394 54220 : distance(idx) = ABS(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2395 : ELSE
2396 2810708 : distance(idx) = HUGE(distance(idx))
2397 : END IF
2398 :
2399 42973920 : dist_min = MINVAL(distance)
2400 2864928 : IF (max_shell == 0) THEN
2401 7436 : image_cell_found = .TRUE.
2402 : END IF
2403 3559780 : IF (dist_min < R_max) THEN
2404 665848 : total_number_of_cells = total_number_of_cells + 1
2405 2663392 : x_data%neighbor_cells(ub)%cell = REAL([ishell, jshell, kshell], dp)
2406 665848 : ub = ub + 1
2407 665848 : image_cell_found = .TRUE.
2408 : END IF
2409 :
2410 : END DO
2411 : END DO
2412 : END DO
2413 30424 : IF (image_cell_found) THEN
2414 22988 : max_shell = max_shell + 1
2415 : ELSE
2416 : nothing_more_to_add = .TRUE.
2417 : END IF
2418 : END DO
2419 : ! now remove what is not needed
2420 732772 : ALLOCATE (tmp_neighbor_cells(total_number_of_cells))
2421 673284 : DO i = 1, ub - 1
2422 673284 : tmp_neighbor_cells(i) = x_data%neighbor_cells(i)
2423 : END DO
2424 7436 : DEALLOCATE (x_data%neighbor_cells)
2425 : ! If we only need the supercell, total_number_of_cells is still 0, repair
2426 7436 : IF (total_number_of_cells == 0) THEN
2427 0 : total_number_of_cells = 1
2428 0 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2429 0 : DO i = 1, total_number_of_cells
2430 0 : x_data%neighbor_cells(i)%cell = 0.0_dp
2431 0 : x_data%neighbor_cells(i)%cell_r = 0.0_dp
2432 : END DO
2433 : ELSE
2434 725336 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2435 673284 : DO i = 1, total_number_of_cells
2436 673284 : x_data%neighbor_cells(i) = tmp_neighbor_cells(i)
2437 : END DO
2438 : END IF
2439 7436 : DEALLOCATE (tmp_neighbor_cells)
2440 :
2441 7436 : IF (x_data%periodic_parameter%number_of_shells == do_hfx_auto_shells) THEN
2442 : ! Do nothing
2443 : ELSE
2444 60 : total_number_of_cells = 0
2445 206 : DO i = 0, x_data%periodic_parameter%number_of_shells
2446 206 : total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2447 : END DO
2448 60 : IF (total_number_of_cells < SIZE(x_data%neighbor_cells)) THEN
2449 60 : IF (i_thread == 1) THEN
2450 4 : WRITE (char_nshells, '(I3)') SIZE(x_data%neighbor_cells)
2451 : WRITE (error_msg, '(A,A,A)') "Periodic Hartree Fock calculation requested with use "// &
2452 : "of a truncated potential. The number of shells to be considered "// &
2453 : "might be too small. CP2K conservatively estimates to need "//TRIM(char_nshells)//" periodic images "// &
2454 4 : "Please carefully check if you get converged results."
2455 4 : CPWARN(error_msg)
2456 : END IF
2457 : END IF
2458 60 : total_number_of_cells = 0
2459 206 : DO i = 0, x_data%periodic_parameter%number_of_shells
2460 206 : total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2461 : END DO
2462 60 : DEALLOCATE (x_data%neighbor_cells)
2463 :
2464 1272 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2465 60 : m = 0
2466 60 : i = 1
2467 3168 : DO WHILE (SUM(m**2) <= x_data%periodic_parameter%number_of_shells)
2468 2928 : x_data%neighbor_cells(i)%cell = REAL(m, dp)
2469 732 : CALL next_image_cell_perd(m, x_data%periodic_parameter%perd)
2470 732 : i = i + 1
2471 : END DO
2472 : END IF
2473 : CASE DEFAULT
2474 2845 : total_number_of_cells = 0
2475 2845 : IF (pbc_shells == -1) pbc_shells = 0
2476 5690 : DO i = 0, pbc_shells
2477 5690 : total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2478 : END DO
2479 2845 : DEALLOCATE (x_data%neighbor_cells)
2480 :
2481 28450 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2482 :
2483 2845 : m = 0
2484 2845 : i = 1
2485 33041 : DO WHILE (SUM(m**2) <= pbc_shells)
2486 11380 : x_data%neighbor_cells(i)%cell = REAL(m, dp)
2487 2845 : CALL next_image_cell_perd(m, x_data%periodic_parameter%perd)
2488 5690 : i = i + 1
2489 : END DO
2490 : END SELECT
2491 :
2492 : ! ** Transform into real coord
2493 674846 : DO i = 1, SIZE(x_data%neighbor_cells)
2494 : r = 0.0_dp
2495 2658260 : x_data%neighbor_cells(i)%cell_r(:) = 0.0_dp
2496 2658260 : s = x_data%neighbor_cells(i)%cell(:)
2497 674846 : CALL scaled_to_real(x_data%neighbor_cells(i)%cell_r, s, cell)
2498 : END DO
2499 10281 : x_data%periodic_parameter%number_of_shells = pbc_shells
2500 :
2501 10281 : R_max_stress = 0.0_dp
2502 674846 : DO i = 1, SIZE(x_data%neighbor_cells)
2503 2668541 : R_max_stress = MAX(R_max_stress, MAXVAL(ABS(x_data%neighbor_cells(i)%cell_r(:))))
2504 : END DO
2505 133653 : R_max_stress = R_max_stress + ABS(MAXVAL(cell%hmat(:, :)))
2506 10281 : x_data%periodic_parameter%R_max_stress = R_max_stress
2507 :
2508 10281 : END SUBROUTINE hfx_create_neighbor_cells
2509 :
2510 : ! performs a fuzzy check of being in a quadrilateral
2511 : ! **************************************************************************************************
2512 : !> \brief ...
2513 : !> \param A ...
2514 : !> \param B ...
2515 : !> \param C ...
2516 : !> \param D ...
2517 : !> \param P ...
2518 : !> \return ...
2519 : ! **************************************************************************************************
2520 17189568 : FUNCTION point_is_in_quadrilateral(A, B, C, D, P)
2521 : REAL(dp) :: A(3), B(3), C(3), D(3), P(3)
2522 : LOGICAL :: point_is_in_quadrilateral
2523 :
2524 : REAL(dp), PARAMETER :: fuzzy = 1000.0_dp*EPSILON(1.0_dp)
2525 :
2526 : REAL(dp) :: dot00, dot01, dot02, dot11, dot12, &
2527 : invDenom, u, v, v0(3), v1(3), v2(3)
2528 :
2529 17189568 : point_is_in_quadrilateral = .FALSE.
2530 :
2531 : ! ** Check for both triangles ABC and ACD
2532 : ! **
2533 : ! ** D -------------- C
2534 : ! ** / /
2535 : ! ** / /
2536 : ! ** A----------------B
2537 : ! **
2538 : ! **
2539 : ! **
2540 :
2541 : ! ** ABC
2542 :
2543 68758272 : v0 = D - A
2544 68758272 : v1 = C - A
2545 68758272 : v2 = P - A
2546 :
2547 : ! ** Compute dot products
2548 68758272 : dot00 = DOT_PRODUCT(v0, v0)
2549 68758272 : dot01 = DOT_PRODUCT(v0, v1)
2550 68758272 : dot02 = DOT_PRODUCT(v0, v2)
2551 68758272 : dot11 = DOT_PRODUCT(v1, v1)
2552 68758272 : dot12 = DOT_PRODUCT(v1, v2)
2553 :
2554 : ! ** Compute barycentric coordinates
2555 17189568 : invDenom = 1/(dot00*dot11 - dot01*dot01)
2556 17189568 : u = (dot11*dot02 - dot01*dot12)*invDenom
2557 17189568 : v = (dot00*dot12 - dot01*dot02)*invDenom
2558 : ! ** Check if point is in triangle
2559 17189568 : IF ((u >= 0 - fuzzy) .AND. (v >= 0 - fuzzy) .AND. (u + v <= 1 + fuzzy)) THEN
2560 17189568 : point_is_in_quadrilateral = .TRUE.
2561 : RETURN
2562 : END IF
2563 67479288 : v0 = C - A
2564 67479288 : v1 = B - A
2565 67479288 : v2 = P - A
2566 :
2567 : ! ** Compute dot products
2568 67479288 : dot00 = DOT_PRODUCT(v0, v0)
2569 67479288 : dot01 = DOT_PRODUCT(v0, v1)
2570 67479288 : dot02 = DOT_PRODUCT(v0, v2)
2571 67479288 : dot11 = DOT_PRODUCT(v1, v1)
2572 67479288 : dot12 = DOT_PRODUCT(v1, v2)
2573 :
2574 : ! ** Compute barycentric coordinates
2575 16869822 : invDenom = 1/(dot00*dot11 - dot01*dot01)
2576 16869822 : u = (dot11*dot02 - dot01*dot12)*invDenom
2577 16869822 : v = (dot00*dot12 - dot01*dot02)*invDenom
2578 :
2579 : ! ** Check if point is in triangle
2580 16869822 : IF ((u >= 0 - fuzzy) .AND. (v >= 0 - fuzzy) .AND. (u + v <= 1 + fuzzy)) THEN
2581 6006 : point_is_in_quadrilateral = .TRUE.
2582 6006 : RETURN
2583 : END IF
2584 :
2585 : END FUNCTION point_is_in_quadrilateral
2586 :
2587 : ! **************************************************************************************************
2588 : !> \brief - This routine deletes all list entries in a container in order to
2589 : !> deallocate the memory.
2590 : !> \param container container that contains the compressed elements
2591 : !> \param memory_usage ...
2592 : !> \param do_disk_storage ...
2593 : !> \par History
2594 : !> 10.2007 created [Manuel Guidon]
2595 : !> \author Manuel Guidon
2596 : ! **************************************************************************************************
2597 3871780 : SUBROUTINE hfx_init_container(container, memory_usage, do_disk_storage)
2598 : TYPE(hfx_container_type) :: container
2599 : INTEGER :: memory_usage
2600 : LOGICAL :: do_disk_storage
2601 :
2602 : TYPE(hfx_container_node), POINTER :: current, next
2603 :
2604 : !! DEALLOCATE memory
2605 :
2606 3871780 : current => container%first
2607 7924277 : DO WHILE (ASSOCIATED(current))
2608 4052497 : next => current%next
2609 4052497 : DEALLOCATE (current)
2610 4052497 : current => next
2611 : END DO
2612 :
2613 : !! Allocate first list entry, init members
2614 3972446280 : ALLOCATE (container%first)
2615 : container%first%prev => NULL()
2616 : container%first%next => NULL()
2617 3871780 : container%current => container%first
2618 3968574500 : container%current%data = 0
2619 3871780 : container%element_counter = 1
2620 3871780 : memory_usage = 1
2621 :
2622 3871780 : IF (do_disk_storage) THEN
2623 : !! close the file, if this is no the first time
2624 390 : IF (container%unit /= -1) THEN
2625 0 : CALL close_file(unit_number=container%unit)
2626 : END IF
2627 : CALL open_file(file_name=TRIM(container%filename), file_status="UNKNOWN", file_form="UNFORMATTED", file_action="WRITE", &
2628 390 : unit_number=container%unit)
2629 : END IF
2630 :
2631 3871780 : END SUBROUTINE hfx_init_container
2632 :
2633 : ! **************************************************************************************************
2634 : !> \brief - This routine stores the data obtained from the load balance routine
2635 : !> for the energy
2636 : !> \param ptr_to_distr contains data to store
2637 : !> \param x_data contains all relevant data structures for hfx runs
2638 : !> \par History
2639 : !> 09.2007 created [Manuel Guidon]
2640 : !> \author Manuel Guidon
2641 : ! **************************************************************************************************
2642 2448 : SUBROUTINE hfx_set_distr_energy(ptr_to_distr, x_data)
2643 : TYPE(hfx_distribution), DIMENSION(:), POINTER :: ptr_to_distr
2644 : TYPE(hfx_type), POINTER :: x_data
2645 :
2646 2448 : DEALLOCATE (x_data%distribution_energy)
2647 :
2648 163890 : ALLOCATE (x_data%distribution_energy(SIZE(ptr_to_distr)))
2649 317988 : x_data%distribution_energy = ptr_to_distr
2650 :
2651 2448 : END SUBROUTINE hfx_set_distr_energy
2652 :
2653 : ! **************************************************************************************************
2654 : !> \brief - This routine stores the data obtained from the load balance routine
2655 : !> for the forces
2656 : !> \param ptr_to_distr contains data to store
2657 : !> \param x_data contains all relevant data structures for hfx runs
2658 : !> \par History
2659 : !> 09.2007 created [Manuel Guidon]
2660 : !> \author Manuel Guidon
2661 : ! **************************************************************************************************
2662 1498 : SUBROUTINE hfx_set_distr_forces(ptr_to_distr, x_data)
2663 : TYPE(hfx_distribution), DIMENSION(:), POINTER :: ptr_to_distr
2664 : TYPE(hfx_type), POINTER :: x_data
2665 :
2666 1498 : DEALLOCATE (x_data%distribution_forces)
2667 :
2668 100360 : ALLOCATE (x_data%distribution_forces(SIZE(ptr_to_distr)))
2669 194740 : x_data%distribution_forces = ptr_to_distr
2670 :
2671 1498 : END SUBROUTINE hfx_set_distr_forces
2672 :
2673 : ! **************************************************************************************************
2674 : !> \brief - resets the maximum memory usage for a HFX calculation subtracting
2675 : !> all relevant buffers from the input MAX_MEM value and add 10% of
2676 : !> safety margin
2677 : !> \param memory_parameter Memory information
2678 : !> \param subtr_size_mb size of buffers in MiB
2679 : !> \par History
2680 : !> 02.2009 created [Manuel Guidon]
2681 : !> \author Manuel Guidon
2682 : ! **************************************************************************************************
2683 41035 : SUBROUTINE hfx_reset_memory_usage_counter(memory_parameter, subtr_size_mb)
2684 :
2685 : TYPE(hfx_memory_type) :: memory_parameter
2686 : INTEGER(int_8), INTENT(IN) :: subtr_size_mb
2687 :
2688 : INTEGER(int_8) :: max_memory
2689 :
2690 41035 : max_memory = memory_parameter%max_memory
2691 41035 : max_memory = max_memory - subtr_size_mb
2692 41035 : IF (max_memory <= 0) THEN
2693 38 : memory_parameter%do_all_on_the_fly = .TRUE.
2694 38 : memory_parameter%max_compression_counter = 0
2695 : ELSE
2696 40997 : memory_parameter%do_all_on_the_fly = .FALSE.
2697 40997 : memory_parameter%max_compression_counter = max_memory*1024_int_8*128_int_8
2698 : END IF
2699 41035 : END SUBROUTINE hfx_reset_memory_usage_counter
2700 :
2701 : ! **************************************************************************************************
2702 : !> \brief - This routine prints some information on HFX
2703 : !> \param x_data contains all relevant data structures for hfx runs
2704 : !> \param hfx_section HFX input section
2705 : !> \par History
2706 : !> 03.2008 created [Manuel Guidon]
2707 : !> \author Manuel Guidon
2708 : ! **************************************************************************************************
2709 1374 : SUBROUTINE hfx_print_std_info(x_data, hfx_section)
2710 : TYPE(hfx_type), POINTER :: x_data
2711 : TYPE(section_vals_type), POINTER :: hfx_section
2712 :
2713 : INTEGER :: iw
2714 : TYPE(cp_logger_type), POINTER :: logger
2715 :
2716 1374 : NULLIFY (logger)
2717 1374 : logger => cp_get_default_logger()
2718 :
2719 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2720 1374 : extension=".scfLog")
2721 :
2722 1374 : IF (iw > 0) THEN
2723 : WRITE (UNIT=iw, FMT="((T3,A,T73,ES8.1))") &
2724 338 : "HFX_INFO| EPS_SCHWARZ: ", x_data%screening_parameter%eps_schwarz
2725 : WRITE (UNIT=iw, FMT="((T3,A,T73,ES8.1))") &
2726 338 : "HFX_INFO| EPS_SCHWARZ_FORCES ", x_data%screening_parameter%eps_schwarz_forces
2727 : WRITE (UNIT=iw, FMT="((T3,A,T73,ES8.1))") &
2728 338 : "HFX_INFO| EPS_STORAGE_SCALING: ", x_data%memory_parameter%eps_storage_scaling
2729 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2730 338 : "HFX_INFO| NBINS: ", x_data%load_balance_parameter%nbins
2731 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2732 338 : "HFX_INFO| BLOCK_SIZE: ", x_data%load_balance_parameter%block_size
2733 338 : IF (x_data%periodic_parameter%do_periodic) THEN
2734 98 : IF (x_data%periodic_parameter%mode == -1) THEN
2735 : WRITE (UNIT=iw, FMT="((T3,A,T77,A))") &
2736 96 : "HFX_INFO| NUMBER_OF_SHELLS: ", "AUTO"
2737 : ELSE
2738 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2739 2 : "HFX_INFO| NUMBER_OF_SHELLS: ", x_data%periodic_parameter%mode
2740 : END IF
2741 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2742 98 : "HFX_INFO| Number of periodic shells considered: ", x_data%periodic_parameter%number_of_shells
2743 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20),/)") &
2744 98 : "HFX_INFO| Number of periodic cells considered: ", SIZE(x_data%neighbor_cells)
2745 : ELSE
2746 : WRITE (UNIT=iw, FMT="((T3,A,T77,A))") &
2747 240 : "HFX_INFO| Number of periodic shells considered: ", "NONE"
2748 : WRITE (UNIT=iw, FMT="((T3,A,T77,A),/)") &
2749 240 : "HFX_INFO| Number of periodic cells considered: ", "NONE"
2750 : END IF
2751 : END IF
2752 1374 : END SUBROUTINE hfx_print_std_info
2753 :
2754 : ! **************************************************************************************************
2755 : !> \brief ...
2756 : !> \param ri_data ...
2757 : !> \param hfx_section ...
2758 : ! **************************************************************************************************
2759 114 : SUBROUTINE hfx_print_ri_info(ri_data, hfx_section)
2760 : TYPE(hfx_ri_type), POINTER :: ri_data
2761 : TYPE(section_vals_type), POINTER :: hfx_section
2762 :
2763 : INTEGER :: iw
2764 : REAL(dp) :: rc_ang
2765 : TYPE(cp_logger_type), POINTER :: logger
2766 : TYPE(section_vals_type), POINTER :: ri_section
2767 :
2768 114 : NULLIFY (logger, ri_section)
2769 114 : logger => cp_get_default_logger()
2770 :
2771 114 : ri_section => ri_data%ri_section
2772 :
2773 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2774 114 : extension=".scfLog")
2775 :
2776 114 : IF (iw > 0) THEN
2777 :
2778 : ASSOCIATE (ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
2779 62 : SELECT CASE (ri_metric%potential_type)
2780 : CASE (do_potential_coulomb)
2781 : WRITE (UNIT=iw, FMT="(/T3,A,T74,A)") &
2782 11 : "HFX_RI_INFO| RI metric: ", "COULOMB"
2783 : CASE (do_potential_short)
2784 : WRITE (UNIT=iw, FMT="(T3,A,T71,A)") &
2785 1 : "HFX_RI_INFO| RI metric: ", "SHORTRANGE"
2786 : WRITE (iw, '(T3,A,T61,F20.10)') &
2787 1 : "HFX_RI_INFO| Omega: ", ri_metric%omega
2788 1 : rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
2789 : WRITE (iw, '(T3,A,T61,F20.10)') &
2790 1 : "HFX_RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
2791 : CASE (do_potential_long)
2792 : WRITE (UNIT=iw, FMT="(T3,A,T72,A)") &
2793 0 : "HFX_RI_INFO| RI metric: ", "LONGRANGE"
2794 : WRITE (iw, '(T3,A,T61,F20.10)') &
2795 0 : "HFX_RI_INFO| Omega: ", ri_metric%omega
2796 : CASE (do_potential_id)
2797 : WRITE (UNIT=iw, FMT="(T3,A,T74,A)") &
2798 33 : "HFX_RI_INFO| RI metric: ", "OVERLAP"
2799 : CASE (do_potential_truncated)
2800 : WRITE (UNIT=iw, FMT="(T3,A,T64,A)") &
2801 5 : "HFX_RI_INFO| RI metric: ", "TRUNCATED COULOMB"
2802 5 : rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
2803 : WRITE (iw, '(T3,A,T61,F20.10)') &
2804 56 : "HFX_RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
2805 : END SELECT
2806 :
2807 : END ASSOCIATE
2808 54 : SELECT CASE (ri_data%flavor)
2809 : CASE (ri_mo)
2810 : WRITE (UNIT=iw, FMT="(T3, A, T79, A)") &
2811 3 : "HFX_RI_INFO| RI flavor: ", "MO"
2812 : CASE (ri_pmat)
2813 : WRITE (UNIT=iw, FMT="(T3, A, T78, A)") &
2814 51 : "HFX_RI_INFO| RI flavor: ", "RHO"
2815 : END SELECT
2816 51 : SELECT CASE (ri_data%t2c_method)
2817 : CASE (hfx_ri_do_2c_iter)
2818 : WRITE (UNIT=iw, FMT="(T3, A, T69, A)") &
2819 0 : "HFX_RI_INFO| Matrix SQRT/INV", "DBCSR / iter"
2820 : CASE (hfx_ri_do_2c_diag)
2821 : WRITE (UNIT=iw, FMT="(T3, A, T65, A)") &
2822 51 : "HFX_RI_INFO| Matrix SQRT/INV", "Dense / diag"
2823 : END SELECT
2824 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2825 51 : "HFX_RI_INFO| EPS_FILTER", ri_data%filter_eps
2826 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2827 51 : "HFX_RI_INFO| EPS_FILTER 2-center", ri_data%filter_eps_2c
2828 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2829 51 : "HFX_RI_INFO| EPS_FILTER storage", ri_data%filter_eps_storage
2830 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2831 51 : "HFX_RI_INFO| EPS_FILTER MO", ri_data%filter_eps_mo
2832 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2833 51 : "HFX_RI_INFO| EPS_PGF_ORB", ri_data%eps_pgf_orb
2834 : WRITE (UNIT=iw, FMT="((T3, A, T73, ES8.1))") &
2835 51 : "HFX_RI_INFO| EPS_SCHWARZ: ", ri_data%eps_schwarz
2836 : WRITE (UNIT=iw, FMT="((T3, A, T73, ES8.1))") &
2837 51 : "HFX_RI_INFO| EPS_SCHWARZ_FORCES: ", ri_data%eps_schwarz_forces
2838 : WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
2839 51 : "HFX_RI_INFO| Minimum block size", ri_data%min_bsize
2840 : WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
2841 51 : "HFX_RI_INFO| MO block size", ri_data%max_bsize_MO
2842 : WRITE (UNIT=iw, FMT="(T3, A, T79, I2)") &
2843 51 : "HFX_RI_INFO| Memory reduction factor", ri_data%n_mem_input
2844 : END IF
2845 :
2846 114 : END SUBROUTINE hfx_print_ri_info
2847 :
2848 : ! **************************************************************************************************
2849 : !> \brief ...
2850 : !> \param x_data ...
2851 : !> \param hfx_section ...
2852 : !> \param i_rep ...
2853 : ! **************************************************************************************************
2854 1488 : SUBROUTINE hfx_print_info(x_data, hfx_section, i_rep)
2855 : TYPE(hfx_type), POINTER :: x_data
2856 : TYPE(section_vals_type), POINTER :: hfx_section
2857 : INTEGER, INTENT(IN) :: i_rep
2858 :
2859 : INTEGER :: iw
2860 : REAL(dp) :: rc_ang
2861 : TYPE(cp_logger_type), POINTER :: logger
2862 :
2863 1488 : NULLIFY (logger)
2864 1488 : logger => cp_get_default_logger()
2865 :
2866 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2867 1488 : extension=".scfLog")
2868 :
2869 1488 : IF (iw > 0) THEN
2870 : WRITE (UNIT=iw, FMT="(/,(T3,A,T61,I20))") &
2871 389 : "HFX_INFO| Replica ID: ", i_rep
2872 :
2873 : WRITE (iw, '(T3,A,T61,F20.10)') &
2874 389 : "HFX_INFO| FRACTION: ", x_data%general_parameter%fraction
2875 636 : SELECT CASE (x_data%potential_parameter%potential_type)
2876 : CASE (do_potential_coulomb)
2877 : WRITE (UNIT=iw, FMT="((T3,A,T74,A))") &
2878 247 : "HFX_INFO| Interaction Potential: ", "COULOMB"
2879 : CASE (do_potential_short)
2880 : WRITE (UNIT=iw, FMT="((T3,A,T71,A))") &
2881 12 : "HFX_INFO| Interaction Potential: ", "SHORTRANGE"
2882 : WRITE (iw, '(T3,A,T61,F20.10)') &
2883 12 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2884 12 : rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2885 : WRITE (iw, '(T3,A,T61,F20.10)') &
2886 12 : "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2887 : CASE (do_potential_long)
2888 : WRITE (UNIT=iw, FMT="((T3,A,T72,A))") &
2889 4 : "HFX_INFO| Interaction Potential: ", "LONGRANGE"
2890 : WRITE (iw, '(T3,A,T61,F20.10)') &
2891 4 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2892 : CASE (do_potential_mix_cl)
2893 : WRITE (UNIT=iw, FMT="((T3,A,T75,A))") &
2894 7 : "HFX_INFO| Interaction Potential: ", "MIX_CL"
2895 : WRITE (iw, '(T3,A,T61,F20.10)') &
2896 7 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2897 : WRITE (iw, '(T3,A,T61,F20.10)') &
2898 7 : "HFX_INFO| SCALE_COULOMB: ", x_data%potential_parameter%scale_coulomb
2899 : WRITE (iw, '(T3,A,T61,F20.10)') &
2900 7 : "HFX_INFO| SCALE_LONGRANGE: ", x_data%potential_parameter%scale_longrange
2901 : CASE (do_potential_gaussian)
2902 : WRITE (UNIT=iw, FMT="((T3,A,T73,A))") &
2903 0 : "HFX_INFO| Interaction Potential: ", "GAUSSIAN"
2904 : WRITE (iw, '(T3,A,T61,F20.10)') &
2905 0 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2906 : CASE (do_potential_mix_lg)
2907 : WRITE (UNIT=iw, FMT="((T3,A,T75,A))") &
2908 2 : "HFX_INFO| Interaction Potential: ", "MIX_LG"
2909 : WRITE (iw, '(T3,A,T61,F20.10)') &
2910 2 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2911 : WRITE (iw, '(T3,A,T61,F20.10)') &
2912 2 : "HFX_INFO| SCALE_LONGRANGE: ", x_data%potential_parameter%scale_longrange
2913 : WRITE (iw, '(T3,A,T61,F20.10)') &
2914 2 : "HFX_INFO| SCALE_GAUSSIAN: ", x_data%potential_parameter%scale_gaussian
2915 : CASE (do_potential_id)
2916 : WRITE (UNIT=iw, FMT="((T3,A,T73,A))") &
2917 14 : "HFX_INFO| Interaction Potential: ", "IDENTITY"
2918 : CASE (do_potential_truncated)
2919 : WRITE (UNIT=iw, FMT="((T3,A,T72,A))") &
2920 94 : "HFX_INFO| Interaction Potential: ", "TRUNCATED"
2921 94 : rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2922 : WRITE (iw, '(T3,A,T61,F20.10)') &
2923 94 : "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2924 : CASE (do_potential_mix_cl_trunc)
2925 : WRITE (UNIT=iw, FMT="((T3,A,T65,A))") &
2926 9 : "HFX_INFO| Interaction Potential: ", "TRUNCATED MIX_CL"
2927 9 : rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2928 : WRITE (iw, '(T3,A,T61,F20.10)') &
2929 398 : "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2930 : END SELECT
2931 :
2932 : END IF
2933 1488 : IF (x_data%do_hfx_ri) THEN
2934 114 : CALL hfx_print_ri_info(x_data%ri_data, hfx_section)
2935 : ELSE
2936 1374 : CALL hfx_print_std_info(x_data, hfx_section)
2937 : END IF
2938 :
2939 : ! ACE section
2940 1488 : IF (x_data%use_ace .AND. iw > 0) THEN
2941 : WRITE (UNIT=iw, FMT="(/,T3,A)") &
2942 4 : "HFX_INFO| ACE (Adaptively Compressed Exchange): ACTIVE"
2943 : WRITE (UNIT=iw, FMT="(T3,A,T61,I20)") &
2944 4 : "HFX_INFO| ACE rebuild frequency: ", x_data%ace_rebuild_freq
2945 : END IF
2946 :
2947 : CALL cp_print_key_finished_output(iw, logger, hfx_section, &
2948 1488 : "HF_INFO")
2949 1488 : END SUBROUTINE hfx_print_info
2950 :
2951 : ! **************************************************************************************************
2952 : !> \brief ...
2953 : !> \param DATA ...
2954 : !> \param memory_usage ...
2955 : ! **************************************************************************************************
2956 32960 : SUBROUTINE dealloc_containers(DATA, memory_usage)
2957 : TYPE(hfx_compression_type) :: data
2958 : INTEGER :: memory_usage
2959 :
2960 : INTEGER :: bin, i
2961 :
2962 65920 : DO bin = 1, SIZE(data%maxval_container)
2963 : CALL hfx_init_container(data%maxval_container(bin), memory_usage, &
2964 32960 : .FALSE.)
2965 65920 : DEALLOCATE (data%maxval_container(bin)%first)
2966 : END DO
2967 32960 : DEALLOCATE (data%maxval_container)
2968 32960 : DEALLOCATE (data%maxval_cache)
2969 :
2970 65920 : DO bin = 1, SIZE(data%integral_containers, 2)
2971 2175360 : DO i = 1, 64
2972 : CALL hfx_init_container(data%integral_containers(i, bin), memory_usage, &
2973 2109440 : .FALSE.)
2974 2142400 : DEALLOCATE (data%integral_containers(i, bin)%first)
2975 : END DO
2976 : END DO
2977 32960 : DEALLOCATE (data%integral_containers)
2978 :
2979 32960 : DEALLOCATE (data%integral_caches)
2980 :
2981 32960 : END SUBROUTINE dealloc_containers
2982 :
2983 : ! **************************************************************************************************
2984 : !> \brief ...
2985 : !> \param DATA ...
2986 : !> \param bin_size ...
2987 : ! **************************************************************************************************
2988 32960 : SUBROUTINE alloc_containers(DATA, bin_size)
2989 : TYPE(hfx_compression_type) :: data
2990 : INTEGER, INTENT(IN) :: bin_size
2991 :
2992 : INTEGER :: bin, i
2993 :
2994 33882880 : ALLOCATE (data%maxval_cache(bin_size))
2995 65920 : DO bin = 1, bin_size
2996 65920 : data%maxval_cache(bin)%element_counter = 1
2997 : END DO
2998 131840 : ALLOCATE (data%maxval_container(bin_size))
2999 65920 : DO bin = 1, bin_size
3000 33816960 : ALLOCATE (data%maxval_container(bin)%first)
3001 : data%maxval_container(bin)%first%prev => NULL()
3002 : data%maxval_container(bin)%first%next => NULL()
3003 32960 : data%maxval_container(bin)%current => data%maxval_container(bin)%first
3004 33784000 : data%maxval_container(bin)%current%data = 0
3005 65920 : data%maxval_container(bin)%element_counter = 1
3006 : END DO
3007 :
3008 2241280 : ALLOCATE (data%integral_containers(64, bin_size))
3009 35992320 : ALLOCATE (data%integral_caches(64, bin_size))
3010 :
3011 65920 : DO bin = 1, bin_size
3012 2175360 : DO i = 1, 64
3013 2109440 : data%integral_caches(i, bin)%element_counter = 1
3014 2162176000 : data%integral_caches(i, bin)%data = 0
3015 2164285440 : ALLOCATE (data%integral_containers(i, bin)%first)
3016 : data%integral_containers(i, bin)%first%prev => NULL()
3017 : data%integral_containers(i, bin)%first%next => NULL()
3018 2109440 : data%integral_containers(i, bin)%current => data%integral_containers(i, bin)%first
3019 2162176000 : data%integral_containers(i, bin)%current%data = 0
3020 2142400 : data%integral_containers(i, bin)%element_counter = 1
3021 : END DO
3022 : END DO
3023 :
3024 32960 : END SUBROUTINE alloc_containers
3025 :
3026 : ! **************************************************************************************************
3027 : !> \brief Compares the non-technical parts of two HFX input section and check whether they are the same
3028 : !> Ignore things that would not change results (MEMORY, LOAD_BALANCE)
3029 : !> \param hfx_section1 ...
3030 : !> \param hfx_section2 ...
3031 : !> \param is_identical ...
3032 : !> \param same_except_frac ...
3033 : !> \return ...
3034 : ! **************************************************************************************************
3035 582 : SUBROUTINE compare_hfx_sections(hfx_section1, hfx_section2, is_identical, same_except_frac)
3036 :
3037 : TYPE(section_vals_type), POINTER :: hfx_section1, hfx_section2
3038 : LOGICAL, INTENT(OUT) :: is_identical
3039 : LOGICAL, INTENT(OUT), OPTIONAL :: same_except_frac
3040 :
3041 : CHARACTER(LEN=default_path_length) :: cval1, cval2
3042 : INTEGER :: irep, ival1, ival2, n_rep_hf1, n_rep_hf2
3043 : LOGICAL :: lval1, lval2
3044 : REAL(dp) :: rval1, rval2
3045 : TYPE(section_vals_type), POINTER :: hfx_sub_section1, hfx_sub_section2
3046 :
3047 194 : is_identical = .TRUE.
3048 194 : IF (PRESENT(same_except_frac)) same_except_frac = .FALSE.
3049 :
3050 194 : CALL section_vals_get(hfx_section1, n_repetition=n_rep_hf1)
3051 194 : CALL section_vals_get(hfx_section2, n_repetition=n_rep_hf2)
3052 194 : is_identical = n_rep_hf1 == n_rep_hf2
3053 200 : IF (.NOT. is_identical) RETURN
3054 :
3055 134 : DO irep = 1, n_rep_hf1
3056 70 : CALL section_vals_val_get(hfx_section1, "PW_HFX", l_val=lval1, i_rep_section=irep)
3057 70 : CALL section_vals_val_get(hfx_section2, "PW_HFX", l_val=lval2, i_rep_section=irep)
3058 70 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3059 :
3060 70 : CALL section_vals_val_get(hfx_section1, "PW_HFX_BLOCKSIZE", i_val=ival1, i_rep_section=irep)
3061 70 : CALL section_vals_val_get(hfx_section2, "PW_HFX_BLOCKSIZE", i_val=ival2, i_rep_section=irep)
3062 70 : IF (ival1 /= ival2) is_identical = .FALSE.
3063 :
3064 70 : CALL section_vals_val_get(hfx_section1, "TREAT_LSD_IN_CORE", l_val=lval1, i_rep_section=irep)
3065 70 : CALL section_vals_val_get(hfx_section2, "TREAT_LSD_IN_CORE", l_val=lval2, i_rep_section=irep)
3066 70 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3067 :
3068 70 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "INTERACTION_POTENTIAL", i_rep_section=irep)
3069 70 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "INTERACTION_POTENTIAL", i_rep_section=irep)
3070 :
3071 70 : CALL section_vals_val_get(hfx_sub_section1, "OMEGA", r_val=rval1, i_rep_section=irep)
3072 70 : CALL section_vals_val_get(hfx_sub_section2, "OMEGA", r_val=rval2, i_rep_section=irep)
3073 70 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3074 :
3075 70 : CALL section_vals_val_get(hfx_sub_section1, "POTENTIAL_TYPE", i_val=ival1, i_rep_section=irep)
3076 70 : CALL section_vals_val_get(hfx_sub_section2, "POTENTIAL_TYPE", i_val=ival2, i_rep_section=irep)
3077 70 : IF (ival1 /= ival2) is_identical = .FALSE.
3078 70 : IF (.NOT. is_identical) RETURN
3079 :
3080 64 : IF (ival1 == do_potential_truncated .OR. ival1 == do_potential_mix_cl_trunc) THEN
3081 6 : CALL section_vals_val_get(hfx_sub_section1, "CUTOFF_RADIUS", r_val=rval1, i_rep_section=irep)
3082 6 : CALL section_vals_val_get(hfx_sub_section2, "CUTOFF_RADIUS", r_val=rval2, i_rep_section=irep)
3083 6 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3084 :
3085 6 : CALL section_vals_val_get(hfx_sub_section1, "T_C_G_DATA", c_val=cval1, i_rep_section=irep)
3086 6 : CALL section_vals_val_get(hfx_sub_section2, "T_C_G_DATA", c_val=cval2, i_rep_section=irep)
3087 6 : IF (cval1 /= cval2) is_identical = .FALSE.
3088 : END IF
3089 :
3090 64 : CALL section_vals_val_get(hfx_sub_section1, "SCALE_COULOMB", r_val=rval1, i_rep_section=irep)
3091 64 : CALL section_vals_val_get(hfx_sub_section2, "SCALE_COULOMB", r_val=rval2, i_rep_section=irep)
3092 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3093 :
3094 64 : CALL section_vals_val_get(hfx_sub_section1, "SCALE_GAUSSIAN", r_val=rval1, i_rep_section=irep)
3095 64 : CALL section_vals_val_get(hfx_sub_section2, "SCALE_GAUSSIAN", r_val=rval2, i_rep_section=irep)
3096 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3097 :
3098 64 : CALL section_vals_val_get(hfx_sub_section1, "SCALE_LONGRANGE", r_val=rval1, i_rep_section=irep)
3099 64 : CALL section_vals_val_get(hfx_sub_section2, "SCALE_LONGRANGE", r_val=rval2, i_rep_section=irep)
3100 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3101 :
3102 64 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "PERIODIC", i_rep_section=irep)
3103 64 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "PERIODIC", i_rep_section=irep)
3104 :
3105 64 : CALL section_vals_val_get(hfx_sub_section1, "NUMBER_OF_SHELLS", i_val=ival1, i_rep_section=irep)
3106 64 : CALL section_vals_val_get(hfx_sub_section2, "NUMBER_OF_SHELLS", i_val=ival2, i_rep_section=irep)
3107 64 : IF (ival1 /= ival2) is_identical = .FALSE.
3108 :
3109 64 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "RI", i_rep_section=irep)
3110 64 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "RI", i_rep_section=irep)
3111 :
3112 64 : CALL section_vals_val_get(hfx_sub_section1, "_SECTION_PARAMETERS_", l_val=lval1, i_rep_section=irep)
3113 64 : CALL section_vals_val_get(hfx_sub_section2, "_SECTION_PARAMETERS_", l_val=lval2, i_rep_section=irep)
3114 64 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3115 :
3116 64 : CALL section_vals_val_get(hfx_sub_section1, "CUTOFF_RADIUS", r_val=rval1, i_rep_section=irep)
3117 64 : CALL section_vals_val_get(hfx_sub_section2, "CUTOFF_RADIUS", r_val=rval2, i_rep_section=irep)
3118 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3119 :
3120 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_EIGVAL", r_val=rval1, i_rep_section=irep)
3121 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_EIGVAL", r_val=rval2, i_rep_section=irep)
3122 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3123 :
3124 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER", r_val=rval1, i_rep_section=irep)
3125 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER", r_val=rval2, i_rep_section=irep)
3126 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3127 :
3128 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER_2C", r_val=rval1, i_rep_section=irep)
3129 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER_2C", r_val=rval2, i_rep_section=irep)
3130 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3131 :
3132 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER_MO", r_val=rval1, i_rep_section=irep)
3133 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER_MO", r_val=rval2, i_rep_section=irep)
3134 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3135 :
3136 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_PGF_ORB", r_val=rval1, i_rep_section=irep)
3137 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_PGF_ORB", r_val=rval2, i_rep_section=irep)
3138 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3139 :
3140 64 : CALL section_vals_val_get(hfx_sub_section1, "MAX_BLOCK_SIZE_MO", i_val=ival1, i_rep_section=irep)
3141 64 : CALL section_vals_val_get(hfx_sub_section2, "MAX_BLOCK_SIZE_MO", i_val=ival2, i_rep_section=irep)
3142 64 : IF (ival1 /= ival2) is_identical = .FALSE.
3143 :
3144 64 : CALL section_vals_val_get(hfx_sub_section1, "MIN_BLOCK_SIZE", i_val=ival1, i_rep_section=irep)
3145 64 : CALL section_vals_val_get(hfx_sub_section2, "MIN_BLOCK_SIZE", i_val=ival2, i_rep_section=irep)
3146 64 : IF (ival1 /= ival2) is_identical = .FALSE.
3147 :
3148 64 : CALL section_vals_val_get(hfx_sub_section1, "OMEGA", r_val=rval1, i_rep_section=irep)
3149 64 : CALL section_vals_val_get(hfx_sub_section2, "OMEGA", r_val=rval2, i_rep_section=irep)
3150 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3151 :
3152 64 : CALL section_vals_val_get(hfx_sub_section1, "RI_FLAVOR", i_val=ival1, i_rep_section=irep)
3153 64 : CALL section_vals_val_get(hfx_sub_section2, "RI_FLAVOR", i_val=ival2, i_rep_section=irep)
3154 64 : IF (ival1 /= ival2) is_identical = .FALSE.
3155 :
3156 64 : CALL section_vals_val_get(hfx_sub_section1, "RI_METRIC", i_val=ival1, i_rep_section=irep)
3157 64 : CALL section_vals_val_get(hfx_sub_section2, "RI_METRIC", i_val=ival2, i_rep_section=irep)
3158 64 : IF (ival1 /= ival2) is_identical = .FALSE.
3159 :
3160 64 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "SCREENING", i_rep_section=irep)
3161 64 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "SCREENING", i_rep_section=irep)
3162 :
3163 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_SCHWARZ", r_val=rval1, i_rep_section=irep)
3164 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_SCHWARZ", r_val=rval2, i_rep_section=irep)
3165 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3166 :
3167 64 : CALL section_vals_val_get(hfx_sub_section1, "EPS_SCHWARZ_FORCES", r_val=rval1, i_rep_section=irep)
3168 64 : CALL section_vals_val_get(hfx_sub_section2, "EPS_SCHWARZ_FORCES", r_val=rval2, i_rep_section=irep)
3169 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3170 :
3171 64 : CALL section_vals_val_get(hfx_sub_section1, "P_SCREEN_CORRECTION_FACTOR", r_val=rval1, i_rep_section=irep)
3172 64 : CALL section_vals_val_get(hfx_sub_section2, "P_SCREEN_CORRECTION_FACTOR", r_val=rval2, i_rep_section=irep)
3173 64 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3174 :
3175 64 : CALL section_vals_val_get(hfx_sub_section1, "SCREEN_ON_INITIAL_P", l_val=lval1, i_rep_section=irep)
3176 64 : CALL section_vals_val_get(hfx_sub_section2, "SCREEN_ON_INITIAL_P", l_val=lval2, i_rep_section=irep)
3177 64 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3178 :
3179 64 : CALL section_vals_val_get(hfx_sub_section1, "SCREEN_P_FORCES", l_val=lval1, i_rep_section=irep)
3180 64 : CALL section_vals_val_get(hfx_sub_section2, "SCREEN_P_FORCES", l_val=lval2, i_rep_section=irep)
3181 1758 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3182 :
3183 : END DO
3184 :
3185 : !Test of the fraction
3186 64 : IF (is_identical) THEN
3187 120 : DO irep = 1, n_rep_hf1
3188 60 : CALL section_vals_val_get(hfx_section1, "FRACTION", r_val=rval1, i_rep_section=irep)
3189 60 : CALL section_vals_val_get(hfx_section2, "FRACTION", r_val=rval2, i_rep_section=irep)
3190 120 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3191 : END DO
3192 :
3193 60 : IF (PRESENT(same_except_frac)) THEN
3194 36 : IF (.NOT. is_identical) same_except_frac = .TRUE.
3195 : END IF
3196 : END IF
3197 :
3198 : END SUBROUTINE compare_hfx_sections
3199 :
3200 0 : END MODULE hfx_types
3201 :
|
