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