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