Line data Source code
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 1310 : 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 1310 : 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 1310 : CALL timeset(routineN, handle)
621 :
622 1310 : CALL cite_reference(Guidon2008)
623 1310 : CALL cite_reference(Guidon2009)
624 :
625 1310 : natom = SIZE(particle_set)
626 :
627 : !! There might be 2 hf sections
628 1310 : CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
629 1310 : n_threads = 1
630 1310 : !$ n_threads = omp_get_max_threads()
631 :
632 1310 : CALL section_vals_val_get(hfx_section, "RI%_SECTION_PARAMETERS_", l_val=do_ri)
633 1310 : IF (do_ri) n_threads = 1 ! RI implementation does not use threads
634 :
635 1310 : IF (PRESENT(orb_basis)) THEN
636 1310 : orb_basis_type = orb_basis
637 : ELSE
638 0 : orb_basis_type = "ORB"
639 : END IF
640 1310 : IF (PRESENT(ri_basis)) THEN
641 0 : ri_basis_type = ri_basis
642 : ELSE
643 1310 : ri_basis_type = "RI_HFX"
644 : END IF
645 :
646 5546550 : ALLOCATE (x_data(n_rep_hf, n_threads))
647 2620 : DO i_thread = 1, n_threads
648 3940 : DO irep = 1, n_rep_hf
649 1320 : actual_x_data => x_data(irep, i_thread)
650 : !! Get data from input file
651 : !!
652 : !! GENERAL params
653 1320 : CALL section_vals_val_get(hfx_section, "FRACTION", r_val=real_val, i_rep_section=irep)
654 1320 : actual_x_data%general_parameter%fraction = real_val
655 1320 : actual_x_data%n_rep_hf = n_rep_hf
656 :
657 1320 : NULLIFY (actual_x_data%map_atoms_to_cpus)
658 :
659 1320 : CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=logic_val, i_rep_section=irep)
660 1320 : actual_x_data%general_parameter%treat_lsd_in_core = logic_val
661 :
662 1320 : hfx_ri_section => section_vals_get_subs_vals(hfx_section, "RI")
663 1320 : CALL section_vals_val_get(hfx_ri_section, "_SECTION_PARAMETERS_", l_val=actual_x_data%do_hfx_ri)
664 :
665 : !! MEMORY section
666 1320 : 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 1320 : n_threads, para_env, irep, skip_disk=.FALSE., skip_in_core_forces=.FALSE.)
669 :
670 : !! PERIODIC section
671 1320 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "PERIODIC", i_rep_section=irep)
672 1320 : CALL section_vals_val_get(hf_sub_section, "NUMBER_OF_SHELLS", i_val=int_val)
673 1320 : actual_x_data%periodic_parameter%number_of_shells = int_val
674 1320 : actual_x_data%periodic_parameter%mode = int_val
675 1320 : CALL get_cell(cell=cell, periodic=actual_x_data%periodic_parameter%perd)
676 5280 : IF (SUM(actual_x_data%periodic_parameter%perd) == 0) THEN
677 924 : actual_x_data%periodic_parameter%do_periodic = .FALSE.
678 : ELSE
679 396 : actual_x_data%periodic_parameter%do_periodic = .TRUE.
680 : END IF
681 :
682 : !! SCREENING section
683 1320 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "SCREENING", i_rep_section=irep)
684 1320 : CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ", r_val=real_val)
685 1320 : actual_x_data%screening_parameter%eps_schwarz = real_val
686 1320 : CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ_FORCES", r_val=real_val, explicit=explicit)
687 1320 : IF (explicit) THEN
688 186 : actual_x_data%screening_parameter%eps_schwarz_forces = real_val
689 : ELSE
690 : actual_x_data%screening_parameter%eps_schwarz_forces = &
691 1134 : 100._dp*actual_x_data%screening_parameter%eps_schwarz
692 : END IF
693 1320 : CALL section_vals_val_get(hf_sub_section, "SCREEN_P_FORCES", l_val=logic_val)
694 1320 : actual_x_data%screening_parameter%do_p_screening_forces = logic_val
695 1320 : CALL section_vals_val_get(hf_sub_section, "SCREEN_ON_INITIAL_P", l_val=logic_val)
696 1320 : actual_x_data%screening_parameter%do_initial_p_screening = logic_val
697 1320 : actual_x_data%screen_funct_is_initialized = .FALSE.
698 :
699 : !! INTERACTION_POTENTIAL section
700 1320 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL", i_rep_section=irep)
701 1320 : CALL section_vals_val_get(hf_sub_section, "POTENTIAL_TYPE", i_val=int_val)
702 1320 : actual_x_data%potential_parameter%potential_type = int_val
703 1320 : CALL section_vals_val_get(hf_sub_section, "OMEGA", r_val=real_val)
704 1320 : actual_x_data%potential_parameter%omega = real_val
705 1320 : CALL section_vals_val_get(hf_sub_section, "SCALE_COULOMB", r_val=real_val)
706 1320 : actual_x_data%potential_parameter%scale_coulomb = real_val
707 1320 : CALL section_vals_val_get(hf_sub_section, "SCALE_LONGRANGE", r_val=real_val)
708 1320 : actual_x_data%potential_parameter%scale_longrange = real_val
709 1320 : CALL section_vals_val_get(hf_sub_section, "SCALE_GAUSSIAN", r_val=real_val)
710 1320 : actual_x_data%potential_parameter%scale_gaussian = real_val
711 1320 : 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 342 : CALL section_vals_val_get(hf_sub_section, "CUTOFF_RADIUS", r_val=real_val)
714 342 : actual_x_data%potential_parameter%cutoff_radius = real_val
715 342 : CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val)
716 342 : CALL compress(char_val, .TRUE.)
717 : ! ** Check if file is there
718 342 : IF (.NOT. file_exists(char_val)) THEN
719 : WRITE (error_msg, '(A,A,A)') "Truncated hfx calculation requested. The file containing "// &
720 0 : "the data could not be found at ", TRIM(char_val), " Please check T_C_G_DATA "// &
721 0 : "in the INTERACTION_POTENTIAL section"
722 0 : CPABORT(error_msg)
723 : ELSE
724 342 : actual_x_data%potential_parameter%filename = char_val
725 : END IF
726 : END IF
727 1320 : 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 46 : actual_x_data%potential_parameter%cutoff_radius)
731 : END IF
732 1320 : IF (actual_x_data%potential_parameter%potential_type == do_potential_id) THEN
733 22 : actual_x_data%potential_parameter%cutoff_radius = 0.0_dp
734 : END IF
735 :
736 : !! LOAD_BALANCE section
737 1320 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "LOAD_BALANCE", i_rep_section=irep)
738 1320 : CALL section_vals_val_get(hf_sub_section, "NBINS", i_val=int_val)
739 1320 : actual_x_data%load_balance_parameter%nbins = MAX(1, int_val)
740 1320 : actual_x_data%load_balance_parameter%blocks_initialized = .FALSE.
741 :
742 1320 : CALL section_vals_val_get(hf_sub_section, "RANDOMIZE", l_val=logic_val)
743 1320 : actual_x_data%load_balance_parameter%do_randomize = logic_val
744 :
745 1320 : actual_x_data%load_balance_parameter%rtp_redistribute = .FALSE.
746 1320 : IF (ASSOCIATED(dft_control%rtp_control)) &
747 34 : actual_x_data%load_balance_parameter%rtp_redistribute = dft_control%rtp_control%hfx_redistribute
748 :
749 1320 : CALL section_vals_val_get(hf_sub_section, "BLOCK_SIZE", i_val=int_val)
750 : ! negative values ask for a computed default
751 1320 : 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 1320 : 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 1320 : 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 1320 : 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 1320 : IF (actual_x_data%periodic_parameter%do_periodic) THEN
860 396 : hf_pbc_section => section_vals_get_subs_vals(hfx_section, "PERIODIC", i_rep_section=irep)
861 396 : CALL section_vals_val_get(hf_pbc_section, "NUMBER_OF_SHELLS", i_val=pbc_shells)
862 396 : actual_x_data%periodic_parameter%number_of_shells_from_input = pbc_shells
863 3168 : ALLOCATE (actual_x_data%neighbor_cells(1))
864 792 : CALL hfx_create_neighbor_cells(actual_x_data, pbc_shells, cell, i_thread, nkp_grid=nkp_grid)
865 : ELSE
866 7392 : ALLOCATE (actual_x_data%neighbor_cells(1))
867 : ! ** Initialize this guy to enable non periodic stress regtests
868 924 : actual_x_data%periodic_parameter%R_max_stress = 1.0_dp
869 : END IF
870 :
871 1320 : nkind = SIZE(qs_kind_set, 1)
872 1320 : max_set = actual_x_data%basis_info%max_set
873 :
874 : !! ** This guy is allocated on the master thread only
875 1320 : IF (i_thread == 1) THEN
876 5280 : ALLOCATE (actual_x_data%is_assoc_atomic_block(natom, natom))
877 5280 : ALLOCATE (actual_x_data%atomic_block_offset(natom, natom))
878 7920 : ALLOCATE (actual_x_data%set_offset(max_set, max_set, nkind, nkind))
879 3960 : ALLOCATE (actual_x_data%block_offset(para_env%num_pe + 1))
880 : END IF
881 :
882 2640 : ALLOCATE (actual_x_data%distribution_forces(1))
883 2640 : ALLOCATE (actual_x_data%distribution_energy(1))
884 :
885 1320 : actual_x_data%memory_parameter%size_p_screen = 0_int_8
886 1320 : IF (i_thread == 1) THEN
887 5280 : ALLOCATE (actual_x_data%atomic_pair_list(natom, natom))
888 5280 : ALLOCATE (actual_x_data%atomic_pair_list_forces(natom, natom))
889 : END IF
890 :
891 1320 : 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 1298 : IF (i_thread == 1) THEN
895 5192 : ALLOCATE (actual_x_data%pmax_atom(natom, natom))
896 7848 : ALLOCATE (actual_x_data%initial_p(nkind*(nkind + 1)/2))
897 1298 : i = 1
898 3662 : DO ikind = 1, nkind
899 2364 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_a)
900 2364 : nseta = actual_x_data%basis_parameter(ikind)%nset
901 7616 : DO jkind = ikind, nkind
902 3954 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=natom_b)
903 3954 : nsetb = actual_x_data%basis_parameter(jkind)%nset
904 23724 : 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 3954 : actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b
907 10272 : i = i + 1
908 : END DO
909 : END DO
910 :
911 5192 : ALLOCATE (actual_x_data%pmax_atom_forces(natom, natom))
912 7848 : ALLOCATE (actual_x_data%initial_p_forces(nkind*(nkind + 1)/2))
913 1298 : i = 1
914 3662 : DO ikind = 1, nkind
915 2364 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_a)
916 2364 : nseta = actual_x_data%basis_parameter(ikind)%nset
917 7616 : DO jkind = ikind, nkind
918 3954 : CALL get_atomic_kind(atomic_kind_set(jkind), natom=natom_b)
919 3954 : nsetb = actual_x_data%basis_parameter(jkind)%nset
920 23724 : 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 3954 : actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b
923 10272 : i = i + 1
924 : END DO
925 : END DO
926 : END IF
927 3894 : ALLOCATE (actual_x_data%map_atom_to_kind_atom(natom))
928 1298 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
929 :
930 3894 : ALLOCATE (atom2kind(nkind))
931 3662 : atom2kind = 0
932 5376 : DO iatom = 1, natom
933 4078 : ikind = kind_of(iatom)
934 4078 : atom2kind(ikind) = atom2kind(ikind) + 1
935 5376 : actual_x_data%map_atom_to_kind_atom(iatom) = atom2kind(ikind)
936 : END DO
937 1298 : DEALLOCATE (kind_of, atom2kind)
938 : END IF
939 :
940 : ! ** Initialize libint type
941 1320 : CALL cp_libint_static_init()
942 1320 : CALL cp_libint_init_eri(actual_x_data%lib, actual_x_data%basis_info%max_am)
943 1320 : CALL cp_libint_init_eri1(actual_x_data%lib_deriv, actual_x_data%basis_info%max_am)
944 1320 : CALL cp_libint_set_contrdepth(actual_x_data%lib, 1)
945 1320 : CALL cp_libint_set_contrdepth(actual_x_data%lib_deriv, 1)
946 :
947 1320 : CALL alloc_containers(actual_x_data%store_ints, 1)
948 1320 : CALL alloc_containers(actual_x_data%store_forces, 1)
949 :
950 1320 : actual_x_data%store_ints%maxval_cache_disk%element_counter = 1
951 1320 : ALLOCATE (actual_x_data%store_ints%maxval_container_disk)
952 1353000 : ALLOCATE (actual_x_data%store_ints%maxval_container_disk%first)
953 1320 : actual_x_data%store_ints%maxval_container_disk%first%prev => NULL()
954 1320 : actual_x_data%store_ints%maxval_container_disk%first%next => NULL()
955 1320 : actual_x_data%store_ints%maxval_container_disk%current => actual_x_data%store_ints%maxval_container_disk%first
956 1353000 : actual_x_data%store_ints%maxval_container_disk%current%data = 0
957 1320 : actual_x_data%store_ints%maxval_container_disk%element_counter = 1
958 1320 : actual_x_data%store_ints%maxval_container_disk%file_counter = 1
959 1320 : actual_x_data%store_ints%maxval_container_disk%desc = 'Max_'
960 1320 : 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 1320 : TRIM(actual_x_data%memory_parameter%storage_location), &
963 2640 : storage_id, "_", actual_x_data%store_ints%maxval_container_disk%desc, "6"
964 1320 : CALL compress(actual_x_data%store_ints%maxval_container_disk%filename, .TRUE.)
965 85800 : ALLOCATE (actual_x_data%store_ints%integral_containers_disk(64))
966 85800 : DO i = 1, 64
967 84480 : actual_x_data%store_ints%integral_caches_disk(i)%element_counter = 1
968 86592000 : actual_x_data%store_ints%integral_caches_disk(i)%data = 0
969 86592000 : ALLOCATE (actual_x_data%store_ints%integral_containers_disk(i)%first)
970 84480 : actual_x_data%store_ints%integral_containers_disk(i)%first%prev => NULL()
971 84480 : actual_x_data%store_ints%integral_containers_disk(i)%first%next => NULL()
972 : actual_x_data%store_ints%integral_containers_disk(i)%current => &
973 84480 : actual_x_data%store_ints%integral_containers_disk(i)%first
974 86592000 : actual_x_data%store_ints%integral_containers_disk(i)%current%data = 0
975 84480 : actual_x_data%store_ints%integral_containers_disk(i)%element_counter = 1
976 84480 : actual_x_data%store_ints%integral_containers_disk(i)%file_counter = 1
977 84480 : actual_x_data%store_ints%integral_containers_disk(i)%desc = 'Int_'
978 84480 : 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 84480 : TRIM(actual_x_data%memory_parameter%storage_location), &
981 168960 : storage_id, "_", actual_x_data%store_ints%integral_containers_disk(i)%desc, i
982 85800 : CALL compress(actual_x_data%store_ints%integral_containers_disk(i)%filename, .TRUE.)
983 : END DO
984 :
985 1320 : actual_x_data%b_first_load_balance_energy = .TRUE.
986 1320 : actual_x_data%b_first_load_balance_forces = .TRUE.
987 :
988 1320 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "RI", i_rep_section=irep)
989 11870 : IF (actual_x_data%do_hfx_ri) THEN
990 100 : CPASSERT(PRESENT(nelectron_total))
991 700 : 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 100 : nelectron_total, orb_basis_type, ri_basis_type)
996 : END IF
997 : END DO
998 : END DO
999 :
1000 2630 : DO irep = 1, n_rep_hf
1001 1320 : actual_x_data => x_data(irep, 1)
1002 2630 : CALL hfx_print_info(actual_x_data, hfx_section, irep)
1003 : END DO
1004 :
1005 1310 : CALL timestop(handle)
1006 :
1007 5240 : 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 100 : 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 100 : CALL timeset(routineN, handle)
1048 :
1049 100 : NULLIFY (hf_sub_section)
1050 :
1051 : ASSOCIATE (hfx_pot => ri_data%hfx_pot)
1052 100 : hfx_pot%potential_type = x_data%potential_parameter%potential_type
1053 100 : hfx_pot%omega = x_data%potential_parameter%omega
1054 100 : hfx_pot%cutoff_radius = x_data%potential_parameter%cutoff_radius
1055 : END ASSOCIATE
1056 100 : ri_data%ri_section => ri_section
1057 100 : ri_data%hfx_section => hfx_section
1058 100 : ri_data%eps_schwarz = x_data%screening_parameter%eps_schwarz
1059 100 : ri_data%eps_schwarz_forces = x_data%screening_parameter%eps_schwarz_forces
1060 :
1061 100 : logger => cp_get_default_logger()
1062 : unit_nr_dbcsr = cp_print_key_unit_nr(logger, ri_data%ri_section, "PRINT%RI_INFO", &
1063 100 : extension=".dbcsrLog")
1064 :
1065 : unit_nr = cp_print_key_unit_nr(logger, ri_data%hfx_section, "HF_INFO", &
1066 100 : extension=".scfLog")
1067 :
1068 100 : hf_sub_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL", i_rep_section=irep)
1069 100 : CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val)
1070 100 : CALL compress(char_val, .TRUE.)
1071 :
1072 100 : IF (.NOT. file_exists(char_val)) THEN
1073 : WRITE (error_msg, '(A,A,A)') "File not found. Please check T_C_G_DATA "// &
1074 0 : "in the INTERACTION_POTENTIAL section"
1075 0 : CPABORT(error_msg)
1076 : ELSE
1077 100 : 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 100 : nelectron_total, t_c_filename=t_c_filename)
1083 :
1084 100 : IF (dft_control%smear .AND. ri_data%flavor == ri_mo) THEN
1085 0 : CPABORT("RI_FLAVOR MO is not consistent with smearing. Please use RI_FLAVOR RHO.")
1086 : END IF
1087 :
1088 100 : CALL timestop(handle)
1089 :
1090 100 : 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 100 : 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 100 : CALL timeset(routineN, handle)
1127 :
1128 100 : CALL section_vals_val_get(ri_section, "EPS_FILTER", r_val=ri_data%filter_eps)
1129 100 : CALL section_vals_val_get(ri_section, "EPS_FILTER_2C", r_val=ri_data%filter_eps_2c)
1130 100 : CALL section_vals_val_get(ri_section, "EPS_STORAGE_SCALING", r_val=eps_storage_scaling)
1131 100 : ri_data%filter_eps_storage = ri_data%filter_eps*eps_storage_scaling
1132 100 : 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 100 : CALL section_vals_val_get(ri_section, "RI_METRIC", i_val=ri_metric%potential_type, explicit=explicit)
1136 100 : IF (.NOT. explicit .OR. ri_metric%potential_type == 0) THEN
1137 42 : ri_metric%potential_type = hfx_pot%potential_type
1138 : END IF
1139 :
1140 100 : CALL section_vals_val_get(ri_section, "OMEGA", r_val=ri_metric%omega, explicit=explicit)
1141 100 : IF (.NOT. explicit) THEN
1142 100 : ri_metric%omega = hfx_pot%omega
1143 : END IF
1144 :
1145 100 : CALL section_vals_val_get(ri_section, "CUTOFF_RADIUS", r_val=ri_metric%cutoff_radius, explicit=explicit)
1146 100 : IF (.NOT. explicit) THEN
1147 94 : ri_metric%cutoff_radius = hfx_pot%cutoff_radius
1148 : END IF
1149 :
1150 100 : IF (ri_metric%potential_type == do_potential_short) &
1151 2 : CALL erfc_cutoff(ri_data%eps_schwarz, ri_metric%omega, ri_metric%cutoff_radius)
1152 100 : IF (ri_metric%potential_type == do_potential_id) ri_metric%cutoff_radius = 0.0_dp
1153 : END ASSOCIATE
1154 :
1155 100 : CALL section_vals_val_get(ri_section, "2C_MATRIX_FUNCTIONS", i_val=ri_data%t2c_method)
1156 100 : CALL section_vals_val_get(ri_section, "EPS_EIGVAL", r_val=ri_data%eps_eigval)
1157 100 : CALL section_vals_val_get(ri_section, "CHECK_2C_MATRIX", l_val=ri_data%check_2c_inv)
1158 100 : CALL section_vals_val_get(ri_section, "CALC_COND_NUM", l_val=ri_data%calc_condnum)
1159 100 : CALL section_vals_val_get(ri_section, "SQRT_ORDER", i_val=ri_data%t2c_sqrt_order)
1160 100 : CALL section_vals_val_get(ri_section, "EPS_LANCZOS", r_val=ri_data%eps_lanczos)
1161 100 : CALL section_vals_val_get(ri_section, "MAX_ITER_LANCZOS", i_val=ri_data%max_iter_lanczos)
1162 100 : CALL section_vals_val_get(ri_section, "RI_FLAVOR", i_val=ri_data%flavor)
1163 100 : CALL section_vals_val_get(ri_section, "EPS_PGF_ORB", r_val=ri_data%eps_pgf_orb)
1164 100 : CALL section_vals_val_get(ri_section, "MIN_BLOCK_SIZE", i_val=ri_data%min_bsize)
1165 100 : CALL section_vals_val_get(ri_section, "MAX_BLOCK_SIZE_MO", i_val=ri_data%max_bsize_MO)
1166 100 : CALL section_vals_val_get(ri_section, "MEMORY_CUT", i_val=ri_data%n_mem_input)
1167 100 : CALL section_vals_val_get(ri_section, "FLAVOR_SWITCH_MEMORY_CUT", i_val=ri_data%n_mem_flavor_switch)
1168 :
1169 100 : ri_data%orb_basis_type = orb_basis_type
1170 100 : ri_data%ri_basis_type = ri_basis_type
1171 100 : ri_data%nelectron_total = nelectron_total
1172 100 : ri_data%input_flavor = ri_data%flavor
1173 :
1174 100 : IF (PRESENT(t_c_filename)) THEN
1175 100 : ri_data%ri_metric%filename = t_c_filename
1176 100 : ri_data%hfx_pot%filename = t_c_filename
1177 : END IF
1178 :
1179 100 : ri_data%unit_nr_dbcsr = unit_nr_dbcsr
1180 100 : ri_data%unit_nr = unit_nr
1181 100 : ri_data%dbcsr_nflop = 0
1182 100 : ri_data%dbcsr_time = 0.0_dp
1183 :
1184 100 : CALL hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
1185 :
1186 100 : CALL timestop(handle)
1187 :
1188 500 : 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 122 : 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 122 : INTEGER, ALLOCATABLE, DIMENSION(:) :: bsizes_AO_store, bsizes_RI_store, dist1, &
1210 122 : 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 122 : DIMENSION(:) :: basis_set_AO, basis_set_RI
1218 122 : TYPE(mp_cart_type) :: mp_comm_3d
1219 :
1220 122 : CALL cite_reference(Bussy2023)
1221 :
1222 122 : CALL timeset(routineN, handle)
1223 :
1224 : ! initialize libint
1225 122 : CALL cp_libint_static_init()
1226 :
1227 122 : natom = SIZE(particle_set)
1228 122 : nkind = SIZE(qs_kind_set, 1)
1229 122 : nproc = para_env%num_pe
1230 :
1231 : ASSOCIATE (ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
1232 122 : IF (ri_metric%potential_type == do_potential_short) THEN
1233 2 : CALL erfc_cutoff(ri_data%eps_schwarz, ri_metric%omega, ri_metric%cutoff_radius)
1234 : END IF
1235 :
1236 122 : 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 4 : 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 122 : same_op = ri_metric%potential_type == hfx_pot%potential_type
1243 :
1244 122 : IF (same_op .AND. hfx_pot%potential_type == do_potential_truncated) THEN
1245 14 : same_op = ABS(ri_metric%cutoff_radius - hfx_pot%cutoff_radius) < 1.0E-16_dp
1246 : END IF
1247 :
1248 244 : IF (same_op .AND. hfx_pot%potential_type == do_potential_short) THEN
1249 2 : same_op = ABS(ri_metric%omega - hfx_pot%omega) < 1.0E-16_dp
1250 : END IF
1251 : END ASSOCIATE
1252 :
1253 122 : ri_data%same_op = same_op
1254 :
1255 122 : pdims = 0
1256 122 : CALL mp_comm_3d%create(para_env, 3, pdims)
1257 :
1258 366 : ALLOCATE (ri_data%bsizes_RI(natom))
1259 366 : ALLOCATE (ri_data%bsizes_AO(natom))
1260 1038 : ALLOCATE (basis_set_RI(nkind), basis_set_AO(nkind))
1261 122 : CALL basis_set_list_setup(basis_set_RI, ri_data%ri_basis_type, qs_kind_set)
1262 122 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_data%bsizes_RI, basis=basis_set_RI)
1263 122 : CALL basis_set_list_setup(basis_set_AO, ri_data%orb_basis_type, qs_kind_set)
1264 122 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_data%bsizes_AO, basis=basis_set_AO)
1265 :
1266 366 : ALLOCATE (dist_RI(natom))
1267 366 : ALLOCATE (dist_AO_1(natom))
1268 366 : ALLOCATE (dist_AO_2(natom))
1269 122 : CALL dbt_default_distvec(natom, pdims(1), ri_data%bsizes_RI, dist_RI)
1270 122 : CALL dbt_default_distvec(natom, pdims(2), ri_data%bsizes_AO, dist_AO_1)
1271 122 : 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 122 : mp_comm_3d, own_comm=.TRUE.)
1274 :
1275 488 : ALLOCATE (ri_data%pgrid)
1276 122 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid)
1277 :
1278 488 : ALLOCATE (ri_data%pgrid_2d)
1279 122 : pdims_2d = 0
1280 122 : CALL dbt_pgrid_create(para_env, pdims_2d, ri_data%pgrid_2d)
1281 :
1282 122 : ri_data%dist_3d = dist_3d
1283 :
1284 : CALL dbt_distribution_new(ri_data%dist, ri_data%pgrid, &
1285 122 : dist_RI, dist_AO_1, dist_AO_2)
1286 :
1287 122 : DEALLOCATE (dist_AO_1, dist_AO_2, dist_RI)
1288 :
1289 122 : ri_data%num_pe = para_env%num_pe
1290 :
1291 : ! initialize tensors expressed in basis representation
1292 122 : CALL pgf_block_sizes(atomic_kind_set, basis_set_AO, ri_data%min_bsize, ri_data%bsizes_AO_split)
1293 122 : CALL pgf_block_sizes(atomic_kind_set, basis_set_RI, ri_data%min_bsize, ri_data%bsizes_RI_split)
1294 :
1295 122 : CALL pgf_block_sizes(atomic_kind_set, basis_set_AO, 1, bsizes_AO_store)
1296 122 : CALL pgf_block_sizes(atomic_kind_set, basis_set_RI, 1, bsizes_RI_store)
1297 :
1298 590 : CALL split_block_sizes([SUM(ri_data%bsizes_AO)], ri_data%bsizes_AO_fit, default_block_size)
1299 590 : CALL split_block_sizes([SUM(ri_data%bsizes_RI)], ri_data%bsizes_RI_fit, default_block_size)
1300 :
1301 122 : 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 104 : ri_data%n_mem = ri_data%n_mem_input
1305 104 : 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 104 : 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 104 : ri_data%starts_array_RI_mem_block, ri_data%ends_array_RI_mem_block)
1314 :
1315 416 : ALLOCATE (ri_data%pgrid_1)
1316 416 : ALLOCATE (ri_data%pgrid_2)
1317 104 : pdims = 0
1318 :
1319 : CALL dbt_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), &
1320 416 : SIZE(ri_data%bsizes_AO_split)])
1321 :
1322 104 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_1)
1323 :
1324 832 : pdims = pdims([2, 1, 3])
1325 104 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_2)
1326 :
1327 1144 : 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 104 : ri_data%bsizes_AO_split, [1, 2], [3], name="(AO RI | AO)")
1331 104 : DEALLOCATE (dist1, dist2, dist3)
1332 :
1333 1364 : ALLOCATE (ri_data%blk_indices(ri_data%n_mem, ri_data%n_mem_RI))
1334 221116 : ALLOCATE (ri_data%store_3c(ri_data%n_mem, ri_data%n_mem_RI))
1335 366 : DO i_mem = 1, ri_data%n_mem
1336 1052 : DO j_mem = 1, ri_data%n_mem_RI
1337 948 : CALL alloc_containers(ri_data%store_3c(i_mem, j_mem), 1)
1338 : END DO
1339 : END DO
1340 :
1341 1144 : 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 104 : ri_data%bsizes_AO_split, [1, 2], [3], name="(AO RI | AO)")
1345 104 : DEALLOCATE (dist1, dist2, dist3)
1346 :
1347 1144 : 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 104 : ri_data%bsizes_AO_split, [1], [2, 3], name="(RI | AO AO)")
1351 104 : DEALLOCATE (dist1, dist2, dist3)
1352 :
1353 1144 : 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 104 : name="(RI | RI)")
1357 104 : 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 1248 : 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 104 : name="(AO | AO)")
1364 104 : DEALLOCATE (dist1, dist2)
1365 104 : CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
1366 :
1367 1248 : 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 104 : name="(AO | AO)")
1371 104 : DEALLOCATE (dist1, dist2)
1372 624 : CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
1373 :
1374 18 : ELSEIF (ri_data%flavor == ri_mo) THEN
1375 216 : 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 18 : name="(RI | RI)")
1380 18 : CALL dbt_create(ri_data%t_2c_int(1, 1), ri_data%t_2c_int(2, 1))
1381 :
1382 18 : DEALLOCATE (dist1, dist2)
1383 :
1384 198 : ALLOCATE (ri_data%t_3c_int_ctr_1(1, 1))
1385 :
1386 72 : ALLOCATE (ri_data%pgrid_1)
1387 72 : ALLOCATE (ri_data%pgrid_2)
1388 : pdims = 0
1389 :
1390 18 : ri_data%n_mem = ri_data%n_mem_input**2
1391 18 : 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 18 : MO_dim = MAX((ri_data%nelectron_total/2 - 1)/ri_data%n_mem + 1, 1)
1397 18 : MO_dim = (MO_dim - 1)/ri_data%max_bsize_MO + 1
1398 :
1399 18 : pdims = 0
1400 72 : CALL dbt_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), MO_dim])
1401 :
1402 18 : CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_1)
1403 :
1404 144 : pdims = pdims([3, 2, 1])
1405 18 : 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 18 : [1, 2], [3], name="(AO RI | AO)")
1410 18 : DEALLOCATE (dist1, dist2, dist3)
1411 :
1412 198 : 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 18 : [1], [2, 3], name="(AO | RI AO)")
1416 36 : DEALLOCATE (dist1, dist2, dist3)
1417 :
1418 : END IF
1419 :
1420 : !For forces
1421 1342 : 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 122 : name="(RI | RI)")
1425 122 : DEALLOCATE (dist1, dist2)
1426 :
1427 1342 : 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 122 : name="(RI | RI)")
1431 122 : DEALLOCATE (dist1, dist2)
1432 :
1433 122 : CALL timestop(handle)
1434 :
1435 732 : END SUBROUTINE
1436 :
1437 : ! **************************************************************************************************
1438 : !> \brief ...
1439 : !> \param ri_data ...
1440 : ! **************************************************************************************************
1441 100 : 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 100 : my_flop_rate = REAL(dbcsr_nflop, dp)/(1.0E09_dp*ri_data%dbcsr_time)
1449 100 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(/T2,A,T73,ES8.2)") &
1450 44 : "RI-HFX PERFORMANCE| DBT total number of flops:", REAL(dbcsr_nflop*num_pe, dp)
1451 100 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T2,A,T66,F15.2)") &
1452 44 : "RI-HFX PERFORMANCE| DBT total execution time:", dbcsr_time
1453 100 : IF (unit_nr > 0) WRITE (UNIT=unit_nr, FMT="(T2,A,T66,F15.2)") &
1454 144 : "RI-HFX PERFORMANCE| DBT flop rate (Gflops / MPI rank):", my_flop_rate
1455 : END ASSOCIATE
1456 100 : END SUBROUTINE
1457 :
1458 : ! **************************************************************************************************
1459 : !> \brief ...
1460 : !> \param ri_data ...
1461 : !> \param write_stats ...
1462 : ! **************************************************************************************************
1463 122 : 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 122 : CALL timeset(routineN, handle)
1473 :
1474 : ! cleanup libint
1475 122 : CALL cp_libint_static_cleanup()
1476 :
1477 122 : my_write_stats = .TRUE.
1478 122 : IF (PRESENT(write_stats)) my_write_stats = write_stats
1479 122 : IF (my_write_stats) CALL hfx_ri_write_stats(ri_data)
1480 :
1481 122 : IF (ASSOCIATED(ri_data%pgrid)) THEN
1482 122 : CALL dbt_pgrid_destroy(ri_data%pgrid)
1483 122 : DEALLOCATE (ri_data%pgrid)
1484 : END IF
1485 122 : IF (ASSOCIATED(ri_data%pgrid_1)) THEN
1486 122 : CALL dbt_pgrid_destroy(ri_data%pgrid_1)
1487 122 : DEALLOCATE (ri_data%pgrid_1)
1488 : END IF
1489 122 : IF (ASSOCIATED(ri_data%pgrid_2)) THEN
1490 122 : CALL dbt_pgrid_destroy(ri_data%pgrid_2)
1491 122 : DEALLOCATE (ri_data%pgrid_2)
1492 : END IF
1493 122 : IF (ASSOCIATED(ri_data%pgrid_2d)) THEN
1494 122 : CALL dbt_pgrid_destroy(ri_data%pgrid_2d)
1495 122 : DEALLOCATE (ri_data%pgrid_2d)
1496 : END IF
1497 :
1498 122 : CALL distribution_3d_destroy(ri_data%dist_3d)
1499 122 : CALL dbt_distribution_destroy(ri_data%dist)
1500 :
1501 122 : DEALLOCATE (ri_data%bsizes_RI)
1502 122 : DEALLOCATE (ri_data%bsizes_AO)
1503 122 : DEALLOCATE (ri_data%bsizes_AO_split)
1504 122 : DEALLOCATE (ri_data%bsizes_RI_split)
1505 122 : DEALLOCATE (ri_data%bsizes_AO_fit)
1506 122 : DEALLOCATE (ri_data%bsizes_RI_fit)
1507 :
1508 122 : IF (ri_data%flavor == ri_pmat) THEN
1509 366 : DO i_mem = 1, ri_data%n_mem
1510 1052 : DO j_mem = 1, ri_data%n_mem_RI
1511 948 : CALL dealloc_containers(ri_data%store_3c(i_mem, j_mem), unused)
1512 : END DO
1513 : END DO
1514 :
1515 1282 : DO j = 1, SIZE(ri_data%t_3c_int_ctr_1, 2)
1516 2460 : DO i = 1, SIZE(ri_data%t_3c_int_ctr_1, 1)
1517 2356 : CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
1518 : END DO
1519 : END DO
1520 1282 : DEALLOCATE (ri_data%t_3c_int_ctr_1)
1521 :
1522 208 : DO j = 1, SIZE(ri_data%t_3c_int_ctr_2, 2)
1523 312 : DO i = 1, SIZE(ri_data%t_3c_int_ctr_2, 1)
1524 208 : CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
1525 : END DO
1526 : END DO
1527 208 : DEALLOCATE (ri_data%t_3c_int_ctr_2)
1528 :
1529 208 : DO j = 1, SIZE(ri_data%t_3c_int_ctr_3, 2)
1530 312 : DO i = 1, SIZE(ri_data%t_3c_int_ctr_3, 1)
1531 208 : CALL dbt_destroy(ri_data%t_3c_int_ctr_3(i, j))
1532 : END DO
1533 : END DO
1534 208 : DEALLOCATE (ri_data%t_3c_int_ctr_3)
1535 :
1536 258 : DO j = 1, SIZE(ri_data%t_2c_int, 2)
1537 412 : DO i = 1, SIZE(ri_data%t_2c_int, 1)
1538 308 : CALL dbt_destroy(ri_data%t_2c_int(i, j))
1539 : END DO
1540 : END DO
1541 258 : DEALLOCATE (ri_data%t_2c_int)
1542 :
1543 1282 : DO j = 1, SIZE(ri_data%rho_ao_t, 2)
1544 2718 : DO i = 1, SIZE(ri_data%rho_ao_t, 1)
1545 2614 : CALL dbt_destroy(ri_data%rho_ao_t(i, j))
1546 : END DO
1547 : END DO
1548 1540 : DEALLOCATE (ri_data%rho_ao_t)
1549 :
1550 1282 : DO j = 1, SIZE(ri_data%ks_t, 2)
1551 2718 : DO i = 1, SIZE(ri_data%ks_t, 1)
1552 2614 : CALL dbt_destroy(ri_data%ks_t(i, j))
1553 : END DO
1554 : END DO
1555 1540 : DEALLOCATE (ri_data%ks_t)
1556 :
1557 0 : DEALLOCATE (ri_data%starts_array_mem_block, ri_data%ends_array_mem_block, &
1558 104 : ri_data%starts_array_mem, ri_data%ends_array_mem)
1559 0 : DEALLOCATE (ri_data%starts_array_RI_mem_block, ri_data%ends_array_RI_mem_block, &
1560 104 : ri_data%starts_array_RI_mem, ri_data%ends_array_RI_mem)
1561 :
1562 790 : DEALLOCATE (ri_data%blk_indices)
1563 104 : DEALLOCATE (ri_data%store_3c)
1564 18 : ELSEIF (ri_data%flavor == ri_mo) THEN
1565 18 : CALL dbt_destroy(ri_data%t_3c_int_ctr_1(1, 1))
1566 18 : CALL dbt_destroy(ri_data%t_3c_int_ctr_2(1, 1))
1567 36 : DEALLOCATE (ri_data%t_3c_int_ctr_1)
1568 36 : DEALLOCATE (ri_data%t_3c_int_ctr_2)
1569 :
1570 40 : DO ispin = 1, SIZE(ri_data%t_3c_int_mo, 1)
1571 22 : CALL dbt_destroy(ri_data%t_3c_int_mo(ispin, 1, 1))
1572 22 : CALL dbt_destroy(ri_data%t_3c_ctr_RI(ispin, 1, 1))
1573 22 : CALL dbt_destroy(ri_data%t_3c_ctr_KS(ispin, 1, 1))
1574 40 : CALL dbt_destroy(ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
1575 : END DO
1576 54 : DO ispin = 1, 2
1577 54 : CALL dbt_destroy(ri_data%t_2c_int(ispin, 1))
1578 : END DO
1579 54 : DEALLOCATE (ri_data%t_2c_int)
1580 40 : DEALLOCATE (ri_data%t_3c_int_mo)
1581 40 : DEALLOCATE (ri_data%t_3c_ctr_RI)
1582 40 : DEALLOCATE (ri_data%t_3c_ctr_KS)
1583 40 : DEALLOCATE (ri_data%t_3c_ctr_KS_copy)
1584 : END IF
1585 :
1586 294 : DO j = 1, SIZE(ri_data%t_2c_inv, 2)
1587 466 : DO i = 1, SIZE(ri_data%t_2c_inv, 1)
1588 344 : CALL dbt_destroy(ri_data%t_2c_inv(i, j))
1589 : END DO
1590 : END DO
1591 294 : DEALLOCATE (ri_data%t_2c_inv)
1592 :
1593 294 : DO j = 1, SIZE(ri_data%t_2c_pot, 2)
1594 466 : DO i = 1, SIZE(ri_data%t_2c_pot, 1)
1595 344 : CALL dbt_destroy(ri_data%t_2c_pot(i, j))
1596 : END DO
1597 : END DO
1598 294 : DEALLOCATE (ri_data%t_2c_pot)
1599 :
1600 122 : IF (ALLOCATED(ri_data%kp_mat_2c_pot)) THEN
1601 1174 : DO j = 1, SIZE(ri_data%kp_mat_2c_pot, 2)
1602 2298 : DO i = 1, SIZE(ri_data%kp_mat_2c_pot, 1)
1603 2248 : CALL dbcsr_release(ri_data%kp_mat_2c_pot(i, j))
1604 : END DO
1605 : END DO
1606 50 : DEALLOCATE (ri_data%kp_mat_2c_pot)
1607 : END IF
1608 :
1609 122 : IF (ALLOCATED(ri_data%kp_t_3c_int)) THEN
1610 1174 : DO i = 1, SIZE(ri_data%kp_t_3c_int)
1611 1174 : CALL dbt_destroy(ri_data%kp_t_3c_int(i))
1612 : END DO
1613 1174 : DEALLOCATE (ri_data%kp_t_3c_int)
1614 : END IF
1615 :
1616 122 : IF (ALLOCATED(ri_data%rho_ao_t)) THEN
1617 0 : DO j = 1, SIZE(ri_data%rho_ao_t, 2)
1618 0 : DO i = 1, SIZE(ri_data%rho_ao_t, 1)
1619 0 : CALL dbt_destroy(ri_data%rho_ao_t(i, j))
1620 : END DO
1621 : END DO
1622 0 : DEALLOCATE (ri_data%rho_ao_t)
1623 : END IF
1624 :
1625 122 : IF (ALLOCATED(ri_data%ks_t)) THEN
1626 0 : DO j = 1, SIZE(ri_data%ks_t, 2)
1627 0 : DO i = 1, SIZE(ri_data%ks_t, 1)
1628 0 : CALL dbt_destroy(ri_data%ks_t(i, j))
1629 : END DO
1630 : END DO
1631 0 : DEALLOCATE (ri_data%ks_t)
1632 : END IF
1633 :
1634 122 : IF (ALLOCATED(ri_data%iatom_to_subgroup)) THEN
1635 150 : DO i = 1, SIZE(ri_data%iatom_to_subgroup)
1636 150 : DEALLOCATE (ri_data%iatom_to_subgroup(i)%array)
1637 : END DO
1638 50 : DEALLOCATE (ri_data%iatom_to_subgroup)
1639 : END IF
1640 :
1641 122 : CALL timestop(handle)
1642 122 : 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 1950 : 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 1950 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell
1666 1950 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl_a
1667 1950 : REAL(dp), DIMENSION(:, :), POINTER :: sphi_a
1668 : TYPE(gto_basis_set_type), POINTER :: orb_basis_a
1669 :
1670 1950 : CALL timeset(routineN, handle)
1671 :
1672 : ! BASIS parameter
1673 1950 : nkind = SIZE(qs_kind_set, 1)
1674 : !
1675 9398 : ALLOCATE (basis_parameter(nkind))
1676 1950 : max_set = 0
1677 5498 : DO ikind = 1, nkind
1678 3548 : 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 3548 : basis_type=basis_type)
1684 3548 : IF (basis_info%max_set < max_set) CPABORT("UNEXPECTED MAX_SET")
1685 3548 : 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 5498 : kind_radius=basis_parameter(ikind)%kind_radius)
1701 : END DO
1702 5498 : DO ikind = 1, nkind
1703 14192 : ALLOCATE (basis_parameter(ikind)%nsgfl(0:basis_info%max_am, max_set))
1704 44196 : basis_parameter(ikind)%nsgfl = 0
1705 3548 : nset = basis_parameter(ikind)%nset
1706 3548 : nshell => basis_parameter(ikind)%nshell
1707 15622 : DO iset = 1, nset
1708 40602 : DO i = 0, basis_info%max_am
1709 26930 : nl_count = 0
1710 62844 : DO j = 1, nshell(iset)
1711 62844 : IF (basis_parameter(ikind)%nl(j, iset) == i) nl_count = nl_count + 1
1712 : END DO
1713 37054 : 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 5498 : DO ikind = 1, nkind
1721 3548 : max_coeff = 0
1722 3548 : max_am_kind = 0
1723 3548 : max_pgf_kind = 0
1724 3548 : npgfa => basis_parameter(ikind)%npgf
1725 3548 : nseta = basis_parameter(ikind)%nset
1726 3548 : nl_a => basis_parameter(ikind)%nsgfl
1727 3548 : la_max => basis_parameter(ikind)%lmax
1728 3548 : la_min => basis_parameter(ikind)%lmin
1729 13672 : DO iset = 1, nseta
1730 10124 : max_pgf_kind = MAX(max_pgf_kind, npgfa(iset))
1731 : max_pgf = MAX(max_pgf, npgfa(iset))
1732 26086 : DO la = la_min(iset), la_max(iset)
1733 12414 : max_nsgfl = MAX(max_nsgfl, nl_a(la, iset))
1734 12414 : max_coeff = MAX(max_coeff, nso(la)*nl_a(la, iset)*nco(la))
1735 22538 : max_am_kind = MAX(max_am_kind, la)
1736 : END DO
1737 : END DO
1738 21288 : ALLOCATE (basis_parameter(ikind)%sphi_ext(max_coeff, 0:max_am_kind, max_pgf_kind, nseta))
1739 2032464 : basis_parameter(ikind)%sphi_ext = 0.0_dp
1740 : END DO
1741 :
1742 5498 : DO ikind = 1, nkind
1743 3548 : sphi_a => basis_parameter(ikind)%sphi
1744 3548 : nseta = basis_parameter(ikind)%nset
1745 3548 : la_max => basis_parameter(ikind)%lmax
1746 3548 : la_min => basis_parameter(ikind)%lmin
1747 3548 : npgfa => basis_parameter(ikind)%npgf
1748 3548 : first_sgfa => basis_parameter(ikind)%first_sgf
1749 3548 : nl_a => basis_parameter(ikind)%nsgfl
1750 15622 : DO iset = 1, nseta
1751 10124 : sgfa = first_sgfa(1, iset)
1752 32876 : DO ipgf = 1, npgfa(iset)
1753 19204 : offset_a1 = (ipgf - 1)*ncoset(la_max(iset))
1754 19204 : s_offset_nl_a = 0
1755 54400 : DO la = la_min(iset), la_max(iset)
1756 25072 : offset_a = offset_a1 + ncoset(la - 1)
1757 : co_counter = 0
1758 25072 : co_counter = co_counter + 1
1759 25072 : so_counter = 0
1760 77476 : DO k = sgfa + s_offset_nl_a, sgfa + s_offset_nl_a + nso(la)*nl_a(la, iset) - 1
1761 221010 : DO i = offset_a + 1, offset_a + nco(la)
1762 143534 : so_counter = so_counter + 1
1763 195938 : basis_parameter(ikind)%sphi_ext(so_counter, la, ipgf, iset) = sphi_a(i, k)
1764 : END DO
1765 : END DO
1766 44276 : 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 1950 : CALL timestop(handle)
1773 :
1774 1950 : END SUBROUTINE hfx_create_basis_types
1775 :
1776 : ! **************************************************************************************************
1777 : !> \brief ...
1778 : !> \param basis_parameter ...
1779 : ! **************************************************************************************************
1780 1950 : 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 1950 : CALL timeset(routineN, handle)
1788 :
1789 : !! BASIS parameter
1790 5498 : DO i = 1, SIZE(basis_parameter)
1791 3548 : DEALLOCATE (basis_parameter(i)%nsgfl)
1792 5498 : DEALLOCATE (basis_parameter(i)%sphi_ext)
1793 : END DO
1794 1950 : DEALLOCATE (basis_parameter)
1795 1950 : CALL timestop(handle)
1796 :
1797 1950 : 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 2318 : 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 2318 : (PRESENT(storage_id) .EQV. PRESENT(irep))
1831 0 : CPASSERT(check)
1832 :
1833 : ! Memory Storage
1834 2318 : CALL section_vals_val_get(hf_sub_section, "MAX_MEMORY", i_val=int_val)
1835 2318 : memory_parameter%max_memory = int_val
1836 2318 : memory_parameter%max_compression_counter = int_val*1024_int_8*128_int_8
1837 2318 : CALL section_vals_val_get(hf_sub_section, "EPS_STORAGE", r_val=real_val)
1838 2318 : memory_parameter%eps_storage_scaling = real_val
1839 2318 : IF (int_val == 0) THEN
1840 20 : memory_parameter%do_all_on_the_fly = .TRUE.
1841 : ELSE
1842 2298 : memory_parameter%do_all_on_the_fly = .FALSE.
1843 : END IF
1844 2318 : memory_parameter%cache_size = CACHE_SIZE
1845 2318 : memory_parameter%bits_max_val = BITS_MAX_VAL
1846 2318 : memory_parameter%actual_memory_usage = 1
1847 2318 : IF (.NOT. skip_in_core_forces) THEN
1848 1320 : CALL section_vals_val_get(hf_sub_section, "TREAT_FORCES_IN_CORE", l_val=logic_val)
1849 1320 : memory_parameter%treat_forces_in_core = logic_val
1850 : END IF
1851 :
1852 : ! ** IF MAX_MEM == 0 overwrite this flag to false
1853 2318 : IF (memory_parameter%do_all_on_the_fly) memory_parameter%treat_forces_in_core = .FALSE.
1854 :
1855 : ! Disk Storage
1856 2318 : IF (.NOT. skip_disk) THEN
1857 1320 : memory_parameter%actual_memory_usage_disk = 1
1858 1320 : CALL section_vals_val_get(hf_sub_section, "MAX_DISK_SPACE", i_val=int_val)
1859 1320 : memory_parameter%max_compression_counter_disk = int_val*1024_int_8*128_int_8
1860 1320 : IF (int_val == 0) THEN
1861 1314 : memory_parameter%do_disk_storage = .FALSE.
1862 : ELSE
1863 6 : memory_parameter%do_disk_storage = .TRUE.
1864 : END IF
1865 1320 : CALL section_vals_val_get(hf_sub_section, "STORAGE_LOCATION", c_val=char_val)
1866 1320 : CALL compress(char_val, .TRUE.)
1867 : !! Add ending / if necessary
1868 :
1869 1320 : IF (SCAN(char_val, "/", .TRUE.) /= LEN_TRIM(char_val)) THEN
1870 1320 : WRITE (filename, '(A,A)') TRIM(char_val), "/"
1871 1320 : CALL compress(filename)
1872 : ELSE
1873 0 : filename = TRIM(char_val)
1874 : END IF
1875 1320 : CALL compress(filename, .TRUE.)
1876 :
1877 : !! quickly check if we can write on storage_location
1878 1320 : CALL m_getcwd(orig_wd)
1879 1320 : CALL m_chdir(TRIM(filename), stat)
1880 1320 : IF (stat /= 0) THEN
1881 0 : WRITE (error_msg, '(A,A,A)') "Request for disk storage failed due to unknown error while writing to ", &
1882 0 : TRIM(filename), ". Please check STORAGE_LOCATION"
1883 0 : CPABORT(error_msg)
1884 : END IF
1885 1320 : CALL m_chdir(orig_wd, stat)
1886 :
1887 1320 : memory_parameter%storage_location = filename
1888 1320 : CALL compress(memory_parameter%storage_location, .TRUE.)
1889 : ELSE
1890 998 : memory_parameter%do_disk_storage = .FALSE.
1891 : END IF
1892 2318 : IF (PRESENT(storage_id)) THEN
1893 1320 : storage_id = (irep - 1)*para_env%num_pe*n_threads + para_env%mepos*n_threads + i_thread - 1
1894 : END IF
1895 2318 : 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 1310 : 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 1310 : n_rep_hf = x_data(1, 1)%n_rep_hf
1914 1310 : n_threads = SIZE(x_data, 2)
1915 :
1916 1310 : 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 342 : init_t_c_g0_lmax = -1
1919 342 : CALL free_C0()
1920 : END IF
1921 2620 : DO i_thread = 1, n_threads
1922 3940 : DO irep = 1, n_rep_hf
1923 1320 : actual_x_data => x_data(irep, i_thread)
1924 1320 : DEALLOCATE (actual_x_data%neighbor_cells)
1925 1320 : DEALLOCATE (actual_x_data%distribution_energy)
1926 1320 : DEALLOCATE (actual_x_data%distribution_forces)
1927 :
1928 1320 : IF (actual_x_data%load_balance_parameter%blocks_initialized) THEN
1929 1216 : DEALLOCATE (actual_x_data%blocks)
1930 1216 : IF (i_thread == 1) THEN
1931 1216 : DEALLOCATE (actual_x_data%pmax_block)
1932 : END IF
1933 : END IF
1934 :
1935 1320 : IF (i_thread == 1) THEN
1936 1320 : DEALLOCATE (actual_x_data%atomic_pair_list)
1937 1320 : DEALLOCATE (actual_x_data%atomic_pair_list_forces)
1938 : END IF
1939 :
1940 1320 : IF (actual_x_data%screening_parameter%do_initial_p_screening .OR. &
1941 : actual_x_data%screening_parameter%do_p_screening_forces) THEN
1942 1298 : IF (i_thread == 1) THEN
1943 1298 : DEALLOCATE (actual_x_data%pmax_atom)
1944 5252 : DO i = 1, SIZE(actual_x_data%initial_p)
1945 5252 : DEALLOCATE (actual_x_data%initial_p(i)%p_kind)
1946 : END DO
1947 1298 : DEALLOCATE (actual_x_data%initial_p)
1948 :
1949 1298 : DEALLOCATE (actual_x_data%pmax_atom_forces)
1950 5252 : DO i = 1, SIZE(actual_x_data%initial_p_forces)
1951 5252 : DEALLOCATE (actual_x_data%initial_p_forces(i)%p_kind)
1952 : END DO
1953 1298 : DEALLOCATE (actual_x_data%initial_p_forces)
1954 : END IF
1955 1298 : DEALLOCATE (actual_x_data%map_atom_to_kind_atom)
1956 : END IF
1957 1320 : IF (i_thread == 1) THEN
1958 1320 : DEALLOCATE (actual_x_data%is_assoc_atomic_block)
1959 1320 : DEALLOCATE (actual_x_data%atomic_block_offset)
1960 1320 : DEALLOCATE (actual_x_data%set_offset)
1961 1320 : DEALLOCATE (actual_x_data%block_offset)
1962 : END IF
1963 :
1964 : !! BASIS parameter
1965 1320 : CALL hfx_release_basis_types(actual_x_data%basis_parameter)
1966 :
1967 : !MK Release libint and libderiv data structure
1968 1320 : CALL cp_libint_cleanup_eri(actual_x_data%lib)
1969 1320 : CALL cp_libint_cleanup_eri1(actual_x_data%lib_deriv)
1970 1320 : CALL cp_libint_static_cleanup()
1971 :
1972 : !! Deallocate containers
1973 1320 : CALL dealloc_containers(actual_x_data%store_ints, actual_x_data%memory_parameter%actual_memory_usage)
1974 1320 : 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 1320 : .FALSE.)
1980 1320 : IF (actual_x_data%memory_parameter%do_disk_storage) THEN
1981 6 : CALL close_file(unit_number=actual_x_data%store_ints%maxval_container_disk%unit, file_status="DELETE")
1982 : END IF
1983 1320 : DEALLOCATE (actual_x_data%store_ints%maxval_container_disk%first)
1984 1320 : DEALLOCATE (actual_x_data%store_ints%maxval_container_disk)
1985 :
1986 85800 : 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 84480 : .FALSE.)
1990 84480 : IF (actual_x_data%memory_parameter%do_disk_storage) THEN
1991 384 : CALL close_file(unit_number=actual_x_data%store_ints%integral_containers_disk(i)%unit, file_status="DELETE")
1992 : END IF
1993 85800 : DEALLOCATE (actual_x_data%store_ints%integral_containers_disk(i)%first)
1994 : END DO
1995 1320 : DEALLOCATE (actual_x_data%store_ints%integral_containers_disk)
1996 :
1997 : ! ** screening functions
1998 1320 : IF (actual_x_data%screen_funct_is_initialized) THEN
1999 1216 : DEALLOCATE (actual_x_data%screen_funct_coeffs_set)
2000 1216 : DEALLOCATE (actual_x_data%screen_funct_coeffs_kind)
2001 1216 : DEALLOCATE (actual_x_data%pair_dist_radii_pgf)
2002 1216 : DEALLOCATE (actual_x_data%screen_funct_coeffs_pgf)
2003 1216 : actual_x_data%screen_funct_is_initialized = .FALSE.
2004 : END IF
2005 :
2006 : ! ** maps
2007 1320 : IF (ASSOCIATED(actual_x_data%map_atoms_to_cpus)) THEN
2008 3634 : DO i = 1, SIZE(actual_x_data%map_atoms_to_cpus)
2009 2418 : DEALLOCATE (actual_x_data%map_atoms_to_cpus(i)%iatom_list)
2010 3634 : DEALLOCATE (actual_x_data%map_atoms_to_cpus(i)%jatom_list)
2011 : END DO
2012 1216 : DEALLOCATE (actual_x_data%map_atoms_to_cpus)
2013 : END IF
2014 :
2015 2630 : IF (actual_x_data%do_hfx_ri) THEN
2016 100 : CALL hfx_ri_release(actual_x_data%ri_data)
2017 100 : IF (ASSOCIATED(actual_x_data%ri_data%ri_section)) THEN
2018 100 : 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 100 : "PRINT%RI_INFO")
2021 : END IF
2022 100 : IF (ASSOCIATED(actual_x_data%ri_data%hfx_section)) THEN
2023 100 : 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 100 : "HF_INFO")
2026 : END IF
2027 100 : DEALLOCATE (actual_x_data%ri_data)
2028 : END IF
2029 : END DO
2030 :
2031 : END DO
2032 :
2033 1310 : DEALLOCATE (x_data)
2034 1310 : 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 9166 : 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 9166 : 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 9166 : REAL(dp), DIMENSION(:, :), POINTER :: zeta, zetb
2065 9166 : TYPE(hfx_cell_type), ALLOCATABLE, DIMENSION(:) :: tmp_neighbor_cells
2066 :
2067 9166 : total_number_of_cells = 0
2068 :
2069 36664 : nkp = 1
2070 9166 : IF (PRESENT(nkp_grid)) nkp = nkp_grid
2071 36514 : do_kpoints = ANY(nkp > 1)
2072 :
2073 : ! ** Check some settings
2074 9166 : 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 396 : 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 92 : "Use a truncated potential to avoid possible problems.")
2083 304 : 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 282 : REAL(nkp(3), dp)*plane_distance(0, 0, 1, cell))
2088 282 : l_min = 0.5_dp*l_min
2089 282 : IF (x_data%potential_parameter%cutoff_radius >= l_min) THEN
2090 28 : 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 28 : "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 0 : "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 15970 : SELECT CASE (x_data%potential_parameter%potential_type)
2110 : CASE (do_potential_truncated, do_potential_mix_cl_trunc, do_potential_short)
2111 6804 : R_max = 0.0_dp
2112 18788 : DO ikind = 1, SIZE(x_data%basis_parameter)
2113 11984 : la_max => x_data%basis_parameter(ikind)%lmax
2114 11984 : zeta => x_data%basis_parameter(ikind)%zet
2115 11984 : nseta = x_data%basis_parameter(ikind)%nset
2116 11984 : npgfa => x_data%basis_parameter(ikind)%npgf
2117 41252 : DO jkind = 1, SIZE(x_data%basis_parameter)
2118 22464 : lb_max => x_data%basis_parameter(jkind)%lmax
2119 22464 : zetb => x_data%basis_parameter(jkind)%zet
2120 22464 : nsetb = x_data%basis_parameter(jkind)%nset
2121 22464 : npgfb => x_data%basis_parameter(jkind)%npgf
2122 93276 : DO iset = 1, nseta
2123 253808 : DO jset = 1, nsetb
2124 529396 : DO ipgf = 1, npgfa(iset)
2125 1027996 : DO jpgf = 1, npgfb(jset)
2126 557428 : Zeta1 = zeta(ipgf, iset) + zetb(jpgf, jset)
2127 : R1 = 1.0_dp/SQRT(Zeta1)*mul_fact(la_max(iset) + lb_max(jset))* &
2128 557428 : SQRT(-LOG(x_data%screening_parameter%eps_schwarz))
2129 855480 : 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 6804 : R_max = 2.0_dp*R_max + x_data%potential_parameter%cutoff_radius
2138 6804 : nothing_more_to_add = .FALSE.
2139 6804 : max_shell = 0
2140 6804 : total_number_of_cells = 0
2141 6804 : ub = 1
2142 6804 : DEALLOCATE (x_data%neighbor_cells)
2143 54432 : ALLOCATE (x_data%neighbor_cells(1))
2144 27216 : x_data%neighbor_cells(1)%cell = 0.0_dp
2145 27216 : 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 34412 : DO WHILE (.NOT. nothing_more_to_add)
2183 : ! Calculate distances to the eight points P1 to P8
2184 27608 : image_cell_found = .FALSE.
2185 1114464 : ALLOCATE (tmp_neighbor_cells(1:ub))
2186 865992 : DO i = 1, ub - 1
2187 865992 : tmp_neighbor_cells(i) = x_data%neighbor_cells(i)
2188 : END DO
2189 27608 : ub_max = (2*max_shell + 1)**3
2190 27608 : DEALLOCATE (x_data%neighbor_cells)
2191 4118168 : ALLOCATE (x_data%neighbor_cells(1:ub_max))
2192 865992 : DO i = 1, ub - 1
2193 865992 : x_data%neighbor_cells(i) = tmp_neighbor_cells(i)
2194 : END DO
2195 3058920 : DO i = ub, ub_max
2196 12125248 : x_data%neighbor_cells(i)%cell = 0.0_dp
2197 12152856 : x_data%neighbor_cells(i)%cell_r = 0.0_dp
2198 : END DO
2199 :
2200 27608 : DEALLOCATE (tmp_neighbor_cells)
2201 :
2202 110432 : perd(1:3) = x_data%periodic_parameter%perd(1:3)
2203 :
2204 140296 : DO ishell = -max_shell*perd(1), max_shell*perd(1)
2205 753788 : DO jshell = -max_shell*perd(2), max_shell*perd(2)
2206 4479548 : DO kshell = -max_shell*perd(3), max_shell*perd(3)
2207 3753368 : IF (MAX(ABS(ishell), ABS(jshell), ABS(kshell)) /= max_shell) CYCLE
2208 : idx = 0
2209 7531896 : DO j = 0, 1
2210 5021264 : x = -1.0_dp/2.0_dp + j*1.0_dp
2211 17574424 : DO k = 0, 1
2212 10042528 : y = -1.0_dp/2.0_dp + k*1.0_dp
2213 35148848 : DO l = 0, 1
2214 20085056 : z = -1.0_dp/2.0_dp + l*1.0_dp
2215 20085056 : idx = idx + 1
2216 20085056 : P(1, idx) = x + ishell
2217 20085056 : P(2, idx) = y + jshell
2218 20085056 : P(3, idx) = z + kshell
2219 20085056 : CALL scaled_to_real(r, P(:, idx), cell)
2220 80340224 : distance(idx) = SQRT(SUM(r**2))
2221 90382752 : 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 2510632 : idx = idx + 1
2229 10042528 : plane_vector(:, 1) = P(:, 3) - P(:, 1)
2230 10042528 : plane_vector(:, 2) = P(:, 2) - P(:, 1)
2231 2510632 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2232 2510632 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2233 2510632 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2234 17574424 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2235 10042528 : 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 2510632 : IF (point_is_in_quadrilateral(P(:, 1), P(:, 3), P(:, 4), P(:, 2), point_in_plane)) THEN
2238 49240 : distance(idx) = ABS(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2239 : ELSE
2240 2461392 : distance(idx) = HUGE(distance(idx))
2241 : END IF
2242 :
2243 : ! Face B (1562) 1 is the reference
2244 2510632 : idx = idx + 1
2245 10042528 : plane_vector(:, 1) = P(:, 2) - P(:, 1)
2246 10042528 : plane_vector(:, 2) = P(:, 5) - P(:, 1)
2247 2510632 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2248 2510632 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2249 2510632 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2250 17574424 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2251 10042528 : 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 2510632 : IF (point_is_in_quadrilateral(P(:, 1), P(:, 5), P(:, 6), P(:, 2), point_in_plane)) THEN
2254 49416 : distance(idx) = ABS(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2255 : ELSE
2256 2461216 : distance(idx) = HUGE(distance(idx))
2257 : END IF
2258 :
2259 : ! Face C (5786) 5 is the reference
2260 2510632 : idx = idx + 1
2261 10042528 : plane_vector(:, 1) = P(:, 7) - P(:, 5)
2262 10042528 : plane_vector(:, 2) = P(:, 6) - P(:, 5)
2263 2510632 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2264 2510632 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2265 2510632 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2266 17574424 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2267 10042528 : 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 2510632 : IF (point_is_in_quadrilateral(P(:, 5), P(:, 7), P(:, 8), P(:, 6), point_in_plane)) THEN
2270 49240 : distance(idx) = ABS(normal(1, 1)*P(1, 5) + normal(2, 1)*P(2, 5) + normal(3, 1)*P(3, 5))
2271 : ELSE
2272 2461392 : distance(idx) = HUGE(distance(idx))
2273 : END IF
2274 :
2275 : ! Face D (3784) 3 is the reference
2276 2510632 : idx = idx + 1
2277 10042528 : plane_vector(:, 1) = P(:, 7) - P(:, 3)
2278 10042528 : plane_vector(:, 2) = P(:, 4) - P(:, 3)
2279 2510632 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2280 2510632 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2281 2510632 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2282 17574424 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2283 10042528 : 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 2510632 : IF (point_is_in_quadrilateral(P(:, 3), P(:, 7), P(:, 8), P(:, 4), point_in_plane)) THEN
2286 49416 : distance(idx) = ABS(normal(1, 1)*P(1, 3) + normal(2, 1)*P(2, 3) + normal(3, 1)*P(3, 3))
2287 : ELSE
2288 2461216 : distance(idx) = HUGE(distance(idx))
2289 : END IF
2290 :
2291 : ! Face E (2684) 2 is the reference
2292 2510632 : idx = idx + 1
2293 10042528 : plane_vector(:, 1) = P(:, 6) - P(:, 2)
2294 10042528 : plane_vector(:, 2) = P(:, 4) - P(:, 2)
2295 2510632 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2296 2510632 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2297 2510632 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2298 17574424 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2299 10042528 : 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 2510632 : IF (point_is_in_quadrilateral(P(:, 2), P(:, 6), P(:, 8), P(:, 4), point_in_plane)) THEN
2302 49220 : distance(idx) = ABS(normal(1, 1)*P(1, 2) + normal(2, 1)*P(2, 2) + normal(3, 1)*P(3, 2))
2303 : ELSE
2304 2461412 : distance(idx) = HUGE(distance(idx))
2305 : END IF
2306 :
2307 : ! Face F (1573) 1 is the reference
2308 2510632 : idx = idx + 1
2309 10042528 : plane_vector(:, 1) = P(:, 5) - P(:, 1)
2310 10042528 : plane_vector(:, 2) = P(:, 3) - P(:, 1)
2311 2510632 : cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2312 2510632 : cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2313 2510632 : cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2314 17574424 : normal(:, 1) = cross_product/SQRT(SUM(cross_product**2))
2315 10042528 : 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 2510632 : IF (point_is_in_quadrilateral(P(:, 1), P(:, 5), P(:, 7), P(:, 3), point_in_plane)) THEN
2318 49220 : distance(idx) = ABS(normal(1, 1)*P(1, 1) + normal(2, 1)*P(2, 1) + normal(3, 1)*P(3, 1))
2319 : ELSE
2320 2461412 : distance(idx) = HUGE(distance(idx))
2321 : END IF
2322 :
2323 40170112 : dist_min = MINVAL(distance)
2324 2510632 : IF (max_shell == 0) THEN
2325 6804 : image_cell_found = .TRUE.
2326 : END IF
2327 3124124 : IF (dist_min < R_max) THEN
2328 573644 : total_number_of_cells = total_number_of_cells + 1
2329 2294576 : x_data%neighbor_cells(ub)%cell = REAL((/ishell, jshell, kshell/), dp)
2330 573644 : ub = ub + 1
2331 573644 : image_cell_found = .TRUE.
2332 : END IF
2333 :
2334 : END DO
2335 : END DO
2336 : END DO
2337 34412 : IF (image_cell_found) THEN
2338 20804 : 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 634880 : ALLOCATE (tmp_neighbor_cells(total_number_of_cells))
2345 580448 : DO i = 1, ub - 1
2346 580448 : tmp_neighbor_cells(i) = x_data%neighbor_cells(i)
2347 : END DO
2348 6804 : DEALLOCATE (x_data%neighbor_cells)
2349 : ! If we only need the supercell, total_number_of_cells is still 0, repair
2350 6804 : IF (total_number_of_cells == 0) THEN
2351 0 : total_number_of_cells = 1
2352 0 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2353 0 : DO i = 1, total_number_of_cells
2354 0 : x_data%neighbor_cells(i)%cell = 0.0_dp
2355 0 : x_data%neighbor_cells(i)%cell_r = 0.0_dp
2356 : END DO
2357 : ELSE
2358 634880 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2359 580448 : DO i = 1, total_number_of_cells
2360 580448 : x_data%neighbor_cells(i) = tmp_neighbor_cells(i)
2361 : END DO
2362 : END IF
2363 6804 : DEALLOCATE (tmp_neighbor_cells)
2364 :
2365 6804 : IF (x_data%periodic_parameter%number_of_shells == do_hfx_auto_shells) THEN
2366 : ! Do nothing
2367 : ELSE
2368 60 : total_number_of_cells = 0
2369 206 : DO i = 0, x_data%periodic_parameter%number_of_shells
2370 206 : total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2371 : END DO
2372 60 : IF (total_number_of_cells < SIZE(x_data%neighbor_cells)) THEN
2373 60 : IF (i_thread == 1) THEN
2374 :
2375 4 : 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 4 : "Please carefully check if you get converged results."
2380 4 : CPWARN(error_msg)
2381 : END IF
2382 : END IF
2383 60 : total_number_of_cells = 0
2384 206 : DO i = 0, x_data%periodic_parameter%number_of_shells
2385 206 : total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2386 : END DO
2387 60 : DEALLOCATE (x_data%neighbor_cells)
2388 :
2389 1272 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2390 60 : m = 0
2391 60 : i = 1
2392 3168 : DO WHILE (SUM(m**2) <= x_data%periodic_parameter%number_of_shells)
2393 2928 : x_data%neighbor_cells(i)%cell = REAL(m, dp)
2394 732 : CALL next_image_cell_perd(m, x_data%periodic_parameter%perd)
2395 732 : i = i + 1
2396 : END DO
2397 : END IF
2398 : CASE DEFAULT
2399 2362 : total_number_of_cells = 0
2400 2362 : IF (pbc_shells == -1) pbc_shells = 0
2401 4724 : DO i = 0, pbc_shells
2402 4724 : total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2403 : END DO
2404 2362 : DEALLOCATE (x_data%neighbor_cells)
2405 :
2406 23620 : ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2407 :
2408 2362 : m = 0
2409 2362 : i = 1
2410 28062 : DO WHILE (SUM(m**2) <= pbc_shells)
2411 9448 : x_data%neighbor_cells(i)%cell = REAL(m, dp)
2412 2362 : CALL next_image_cell_perd(m, x_data%periodic_parameter%perd)
2413 4724 : i = i + 1
2414 : END DO
2415 : END SELECT
2416 :
2417 : ! ** Transform into real coord
2418 581044 : DO i = 1, SIZE(x_data%neighbor_cells)
2419 : r = 0.0_dp
2420 2287512 : x_data%neighbor_cells(i)%cell_r(:) = 0.0_dp
2421 2287512 : s = x_data%neighbor_cells(i)%cell(:)
2422 581044 : CALL scaled_to_real(x_data%neighbor_cells(i)%cell_r, s, cell)
2423 : END DO
2424 9166 : x_data%periodic_parameter%number_of_shells = pbc_shells
2425 :
2426 9166 : R_max_stress = 0.0_dp
2427 581044 : DO i = 1, SIZE(x_data%neighbor_cells)
2428 2868556 : R_max_stress = MAX(R_max_stress, MAXVAL(ABS(x_data%neighbor_cells(i)%cell_r(:))))
2429 : END DO
2430 119158 : R_max_stress = R_max_stress + ABS(MAXVAL(cell%hmat(:, :)))
2431 9166 : x_data%periodic_parameter%R_max_stress = R_max_stress
2432 :
2433 9166 : 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 15063792 : 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 15063792 : 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 60255168 : v0 = D - A
2469 60255168 : v1 = C - A
2470 60255168 : v2 = P - A
2471 :
2472 : ! ** Compute dot products
2473 60255168 : dot00 = DOT_PRODUCT(v0, v0)
2474 60255168 : dot01 = DOT_PRODUCT(v0, v1)
2475 60255168 : dot02 = DOT_PRODUCT(v0, v2)
2476 60255168 : dot11 = DOT_PRODUCT(v1, v1)
2477 60255168 : dot12 = DOT_PRODUCT(v1, v2)
2478 :
2479 : ! ** Compute barycentric coordinates
2480 15063792 : invDenom = 1/(dot00*dot11 - dot01*dot01)
2481 15063792 : u = (dot11*dot02 - dot01*dot12)*invDenom
2482 15063792 : v = (dot00*dot12 - dot01*dot02)*invDenom
2483 : ! ** Check if point is in triangle
2484 15063792 : IF ((u >= 0 - fuzzy) .AND. (v >= 0 - fuzzy) .AND. (u + v <= 1 + fuzzy)) THEN
2485 15063792 : point_is_in_quadrilateral = .TRUE.
2486 : RETURN
2487 : END IF
2488 59092704 : v0 = C - A
2489 59092704 : v1 = B - A
2490 59092704 : v2 = P - A
2491 :
2492 : ! ** Compute dot products
2493 59092704 : dot00 = DOT_PRODUCT(v0, v0)
2494 59092704 : dot01 = DOT_PRODUCT(v0, v1)
2495 59092704 : dot02 = DOT_PRODUCT(v0, v2)
2496 59092704 : dot11 = DOT_PRODUCT(v1, v1)
2497 59092704 : dot12 = DOT_PRODUCT(v1, v2)
2498 :
2499 : ! ** Compute barycentric coordinates
2500 14773176 : invDenom = 1/(dot00*dot11 - dot01*dot01)
2501 14773176 : u = (dot11*dot02 - dot01*dot12)*invDenom
2502 14773176 : v = (dot00*dot12 - dot01*dot02)*invDenom
2503 :
2504 : ! ** Check if point is in triangle
2505 14773176 : IF ((u >= 0 - fuzzy) .AND. (v >= 0 - fuzzy) .AND. (u + v <= 1 + fuzzy)) THEN
2506 5136 : point_is_in_quadrilateral = .TRUE.
2507 5136 : 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 3545350 : 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 3545350 : current => container%first
2532 7223083 : DO WHILE (ASSOCIATED(current))
2533 3677733 : next => current%next
2534 3677733 : DEALLOCATE (current)
2535 3677733 : current => next
2536 : END DO
2537 :
2538 : !! Allocate first list entry, init members
2539 3637529100 : ALLOCATE (container%first)
2540 : container%first%prev => NULL()
2541 : container%first%next => NULL()
2542 3545350 : container%current => container%first
2543 3633983750 : container%current%data = 0
2544 3545350 : container%element_counter = 1
2545 3545350 : memory_usage = 1
2546 :
2547 3545350 : IF (do_disk_storage) THEN
2548 : !! close the file, if this is no the first time
2549 390 : IF (container%unit /= -1) THEN
2550 0 : 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 390 : unit_number=container%unit)
2554 : END IF
2555 :
2556 3545350 : 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 2126 : 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 2126 : DEALLOCATE (x_data%distribution_energy)
2572 :
2573 141560 : ALLOCATE (x_data%distribution_energy(SIZE(ptr_to_distr)))
2574 274616 : x_data%distribution_energy = ptr_to_distr
2575 :
2576 2126 : 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 1330 : 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 1330 : DEALLOCATE (x_data%distribution_forces)
2592 :
2593 89104 : ALLOCATE (x_data%distribution_forces(SIZE(ptr_to_distr)))
2594 172900 : x_data%distribution_forces = ptr_to_distr
2595 :
2596 1330 : 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 34546 : 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 34546 : max_memory = memory_parameter%max_memory
2616 34546 : max_memory = max_memory - subtr_size_mb
2617 34546 : IF (max_memory <= 0) THEN
2618 30 : memory_parameter%do_all_on_the_fly = .TRUE.
2619 30 : memory_parameter%max_compression_counter = 0
2620 : ELSE
2621 34516 : memory_parameter%do_all_on_the_fly = .FALSE.
2622 34516 : memory_parameter%max_compression_counter = max_memory*1024_int_8*128_int_8
2623 : END IF
2624 34546 : 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 1220 : 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 1220 : NULLIFY (logger)
2642 1220 : logger => cp_get_default_logger()
2643 :
2644 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2645 1220 : extension=".scfLog")
2646 :
2647 1220 : IF (iw > 0) THEN
2648 : WRITE (UNIT=iw, FMT="((T3,A,T73,ES8.1))") &
2649 296 : "HFX_INFO| EPS_SCHWARZ: ", x_data%screening_parameter%eps_schwarz
2650 : WRITE (UNIT=iw, FMT="((T3,A,T73,ES8.1))") &
2651 296 : "HFX_INFO| EPS_SCHWARZ_FORCES ", x_data%screening_parameter%eps_schwarz_forces
2652 : WRITE (UNIT=iw, FMT="((T3,A,T73,ES8.1))") &
2653 296 : "HFX_INFO| EPS_STORAGE_SCALING: ", x_data%memory_parameter%eps_storage_scaling
2654 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2655 296 : "HFX_INFO| NBINS: ", x_data%load_balance_parameter%nbins
2656 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2657 296 : "HFX_INFO| BLOCK_SIZE: ", x_data%load_balance_parameter%block_size
2658 296 : IF (x_data%periodic_parameter%do_periodic) THEN
2659 90 : IF (x_data%periodic_parameter%mode == -1) THEN
2660 : WRITE (UNIT=iw, FMT="((T3,A,T77,A))") &
2661 88 : "HFX_INFO| NUMBER_OF_SHELLS: ", "AUTO"
2662 : ELSE
2663 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2664 2 : "HFX_INFO| NUMBER_OF_SHELLS: ", x_data%periodic_parameter%mode
2665 : END IF
2666 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20))") &
2667 90 : "HFX_INFO| Number of periodic shells considered: ", x_data%periodic_parameter%number_of_shells
2668 : WRITE (UNIT=iw, FMT="((T3,A,T61,I20),/)") &
2669 90 : "HFX_INFO| Number of periodic cells considered: ", SIZE(x_data%neighbor_cells)
2670 : ELSE
2671 : WRITE (UNIT=iw, FMT="((T3,A,T77,A))") &
2672 206 : "HFX_INFO| Number of periodic shells considered: ", "NONE"
2673 : WRITE (UNIT=iw, FMT="((T3,A,T77,A),/)") &
2674 206 : "HFX_INFO| Number of periodic cells considered: ", "NONE"
2675 : END IF
2676 : END IF
2677 1220 : END SUBROUTINE hfx_print_std_info
2678 :
2679 : ! **************************************************************************************************
2680 : !> \brief ...
2681 : !> \param ri_data ...
2682 : !> \param hfx_section ...
2683 : ! **************************************************************************************************
2684 100 : 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 100 : NULLIFY (logger, ri_section)
2694 100 : logger => cp_get_default_logger()
2695 :
2696 100 : ri_section => ri_data%ri_section
2697 :
2698 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2699 100 : extension=".scfLog")
2700 :
2701 100 : IF (iw > 0) THEN
2702 :
2703 : ASSOCIATE (ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
2704 55 : SELECT CASE (ri_metric%potential_type)
2705 : CASE (do_potential_coulomb)
2706 : WRITE (UNIT=iw, FMT="(/T3,A,T74,A)") &
2707 11 : "HFX_RI_INFO| RI metric: ", "COULOMB"
2708 : CASE (do_potential_short)
2709 : WRITE (UNIT=iw, FMT="(T3,A,T71,A)") &
2710 1 : "HFX_RI_INFO| RI metric: ", "SHORTRANGE"
2711 : WRITE (iw, '(T3,A,T61,F20.10)') &
2712 1 : "HFX_RI_INFO| Omega: ", ri_metric%omega
2713 1 : rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
2714 : WRITE (iw, '(T3,A,T61,F20.10)') &
2715 1 : "HFX_RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
2716 : CASE (do_potential_long)
2717 : WRITE (UNIT=iw, FMT="(T3,A,T72,A)") &
2718 0 : "HFX_RI_INFO| RI metric: ", "LONGRANGE"
2719 : WRITE (iw, '(T3,A,T61,F20.10)') &
2720 0 : "HFX_RI_INFO| Omega: ", ri_metric%omega
2721 : CASE (do_potential_id)
2722 : WRITE (UNIT=iw, FMT="(T3,A,T73,A)") &
2723 27 : "HFX_RI_INFO| RI metric: ", "OVERLAP"
2724 : CASE (do_potential_truncated)
2725 : WRITE (UNIT=iw, FMT="(T3,A,T72,A)") &
2726 5 : "HFX_RI_INFO| RI metric: ", "TRUNCATED COULOMB"
2727 5 : rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
2728 : WRITE (iw, '(T3,A,T61,F20.10)') &
2729 49 : "HFX_RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
2730 : END SELECT
2731 :
2732 : END ASSOCIATE
2733 47 : SELECT CASE (ri_data%flavor)
2734 : CASE (ri_mo)
2735 : WRITE (UNIT=iw, FMT="(T3, A, T79, A)") &
2736 3 : "HFX_RI_INFO| RI flavor: ", "MO"
2737 : CASE (ri_pmat)
2738 : WRITE (UNIT=iw, FMT="(T3, A, T78, A)") &
2739 44 : "HFX_RI_INFO| RI flavor: ", "RHO"
2740 : END SELECT
2741 44 : SELECT CASE (ri_data%t2c_method)
2742 : CASE (hfx_ri_do_2c_iter)
2743 : WRITE (UNIT=iw, FMT="(T3, A, T69, A)") &
2744 0 : "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 44 : "HFX_RI_INFO| Matrix SQRT/INV", "Dense / diag"
2748 : END SELECT
2749 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2750 44 : "HFX_RI_INFO| EPS_FILTER", ri_data%filter_eps
2751 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2752 44 : "HFX_RI_INFO| EPS_FILTER 2-center", ri_data%filter_eps_2c
2753 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2754 44 : "HFX_RI_INFO| EPS_FILTER storage", ri_data%filter_eps_storage
2755 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2756 44 : "HFX_RI_INFO| EPS_FILTER MO", ri_data%filter_eps_mo
2757 : WRITE (UNIT=iw, FMT="(T3, A, T73, ES8.1)") &
2758 44 : "HFX_RI_INFO| EPS_PGF_ORB", ri_data%eps_pgf_orb
2759 : WRITE (UNIT=iw, FMT="((T3, A, T73, ES8.1))") &
2760 44 : "HFX_RI_INFO| EPS_SCHWARZ: ", ri_data%eps_schwarz
2761 : WRITE (UNIT=iw, FMT="((T3, A, T73, ES8.1))") &
2762 44 : "HFX_RI_INFO| EPS_SCHWARZ_FORCES: ", ri_data%eps_schwarz_forces
2763 : WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
2764 44 : "HFX_RI_INFO| Minimum block size", ri_data%min_bsize
2765 : WRITE (UNIT=iw, FMT="(T3, A, T78, I3)") &
2766 44 : "HFX_RI_INFO| MO block size", ri_data%max_bsize_MO
2767 : WRITE (UNIT=iw, FMT="(T3, A, T79, I2)") &
2768 44 : "HFX_RI_INFO| Memory reduction factor", ri_data%n_mem_input
2769 : END IF
2770 :
2771 100 : END SUBROUTINE
2772 :
2773 : ! **************************************************************************************************
2774 : !> \brief ...
2775 : !> \param x_data ...
2776 : !> \param hfx_section ...
2777 : !> \param i_rep ...
2778 : ! **************************************************************************************************
2779 1320 : 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 1320 : NULLIFY (logger)
2789 1320 : logger => cp_get_default_logger()
2790 :
2791 : iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2792 1320 : extension=".scfLog")
2793 :
2794 1320 : IF (iw > 0) THEN
2795 : WRITE (UNIT=iw, FMT="(/,(T3,A,T61,I20))") &
2796 340 : "HFX_INFO| Replica ID: ", i_rep
2797 :
2798 : WRITE (iw, '(T3,A,T61,F20.10)') &
2799 340 : "HFX_INFO| FRACTION: ", x_data%general_parameter%fraction
2800 547 : SELECT CASE (x_data%potential_parameter%potential_type)
2801 : CASE (do_potential_coulomb)
2802 : WRITE (UNIT=iw, FMT="((T3,A,T74,A))") &
2803 207 : "HFX_INFO| Interaction Potential: ", "COULOMB"
2804 : CASE (do_potential_short)
2805 : WRITE (UNIT=iw, FMT="((T3,A,T71,A))") &
2806 12 : "HFX_INFO| Interaction Potential: ", "SHORTRANGE"
2807 : WRITE (iw, '(T3,A,T61,F20.10)') &
2808 12 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2809 12 : rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2810 : WRITE (iw, '(T3,A,T61,F20.10)') &
2811 12 : "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2812 : CASE (do_potential_long)
2813 : WRITE (UNIT=iw, FMT="((T3,A,T72,A))") &
2814 4 : "HFX_INFO| Interaction Potential: ", "LONGRANGE"
2815 : WRITE (iw, '(T3,A,T61,F20.10)') &
2816 4 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2817 : CASE (do_potential_mix_cl)
2818 : WRITE (UNIT=iw, FMT="((T3,A,T75,A))") &
2819 7 : "HFX_INFO| Interaction Potential: ", "MIX_CL"
2820 : WRITE (iw, '(T3,A,T61,F20.10)') &
2821 7 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2822 : WRITE (iw, '(T3,A,T61,F20.10)') &
2823 7 : "HFX_INFO| SCALE_COULOMB: ", x_data%potential_parameter%scale_coulomb
2824 : WRITE (iw, '(T3,A,T61,F20.10)') &
2825 7 : "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 0 : "HFX_INFO| Interaction Potential: ", "GAUSSIAN"
2829 : WRITE (iw, '(T3,A,T61,F20.10)') &
2830 0 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2831 : CASE (do_potential_mix_lg)
2832 : WRITE (UNIT=iw, FMT="((T3,A,T75,A))") &
2833 2 : "HFX_INFO| Interaction Potential: ", "MIX_LG"
2834 : WRITE (iw, '(T3,A,T61,F20.10)') &
2835 2 : "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2836 : WRITE (iw, '(T3,A,T61,F20.10)') &
2837 2 : "HFX_INFO| SCALE_LONGRANGE: ", x_data%potential_parameter%scale_longrange
2838 : WRITE (iw, '(T3,A,T61,F20.10)') &
2839 2 : "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 11 : "HFX_INFO| Interaction Potential: ", "IDENTITY"
2843 : CASE (do_potential_truncated)
2844 : WRITE (UNIT=iw, FMT="((T3,A,T72,A))") &
2845 92 : "HFX_INFO| Interaction Potential: ", "TRUNCATED"
2846 92 : rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2847 : WRITE (iw, '(T3,A,T61,F20.10)') &
2848 92 : "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2849 : CASE (do_potential_mix_cl_trunc)
2850 : WRITE (UNIT=iw, FMT="((T3,A,T65,A))") &
2851 5 : "HFX_INFO| Interaction Potential: ", "TRUNCATED MIX_CL"
2852 5 : rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2853 : WRITE (iw, '(T3,A,T61,F20.10)') &
2854 345 : "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2855 : END SELECT
2856 :
2857 : END IF
2858 1320 : IF (x_data%do_hfx_ri) THEN
2859 100 : CALL hfx_print_ri_info(x_data%ri_data, hfx_section)
2860 : ELSE
2861 1220 : 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 1320 : "HF_INFO")
2866 1320 : END SUBROUTINE
2867 :
2868 : ! **************************************************************************************************
2869 : !> \brief ...
2870 : !> \param DATA ...
2871 : !> \param memory_usage ...
2872 : ! **************************************************************************************************
2873 30296 : SUBROUTINE dealloc_containers(DATA, memory_usage)
2874 : TYPE(hfx_compression_type) :: data
2875 : INTEGER :: memory_usage
2876 :
2877 : INTEGER :: bin, i
2878 :
2879 60592 : DO bin = 1, SIZE(data%maxval_container)
2880 : CALL hfx_init_container(data%maxval_container(bin), memory_usage, &
2881 30296 : .FALSE.)
2882 60592 : DEALLOCATE (data%maxval_container(bin)%first)
2883 : END DO
2884 30296 : DEALLOCATE (data%maxval_container)
2885 30296 : DEALLOCATE (data%maxval_cache)
2886 :
2887 60592 : DO bin = 1, SIZE(data%integral_containers, 2)
2888 1999536 : DO i = 1, 64
2889 : CALL hfx_init_container(data%integral_containers(i, bin), memory_usage, &
2890 1938944 : .FALSE.)
2891 1969240 : DEALLOCATE (data%integral_containers(i, bin)%first)
2892 : END DO
2893 : END DO
2894 30296 : DEALLOCATE (data%integral_containers)
2895 :
2896 30296 : DEALLOCATE (data%integral_caches)
2897 :
2898 30296 : END SUBROUTINE dealloc_containers
2899 :
2900 : ! **************************************************************************************************
2901 : !> \brief ...
2902 : !> \param DATA ...
2903 : !> \param bin_size ...
2904 : ! **************************************************************************************************
2905 30296 : 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 31144288 : ALLOCATE (data%maxval_cache(bin_size))
2912 60592 : DO bin = 1, bin_size
2913 60592 : data%maxval_cache(bin)%element_counter = 1
2914 : END DO
2915 121184 : ALLOCATE (data%maxval_container(bin_size))
2916 60592 : DO bin = 1, bin_size
2917 31083696 : ALLOCATE (data%maxval_container(bin)%first)
2918 : data%maxval_container(bin)%first%prev => NULL()
2919 : data%maxval_container(bin)%first%next => NULL()
2920 30296 : data%maxval_container(bin)%current => data%maxval_container(bin)%first
2921 31053400 : data%maxval_container(bin)%current%data = 0
2922 60592 : data%maxval_container(bin)%element_counter = 1
2923 : END DO
2924 :
2925 2060128 : ALLOCATE (data%integral_containers(64, bin_size))
2926 33083232 : ALLOCATE (data%integral_caches(64, bin_size))
2927 :
2928 60592 : DO bin = 1, bin_size
2929 1999536 : DO i = 1, 64
2930 1938944 : data%integral_caches(i, bin)%element_counter = 1
2931 1987417600 : data%integral_caches(i, bin)%data = 0
2932 1989356544 : 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 1938944 : data%integral_containers(i, bin)%current => data%integral_containers(i, bin)%first
2936 1987417600 : data%integral_containers(i, bin)%current%data = 0
2937 1969240 : data%integral_containers(i, bin)%element_counter = 1
2938 : END DO
2939 : END DO
2940 :
2941 30296 : 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 450 : 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 150 : is_identical = .TRUE.
2965 150 : IF (PRESENT(same_except_frac)) same_except_frac = .FALSE.
2966 :
2967 150 : CALL section_vals_get(hfx_section1, n_repetition=n_rep_hf1)
2968 150 : CALL section_vals_get(hfx_section2, n_repetition=n_rep_hf2)
2969 150 : is_identical = n_rep_hf1 == n_rep_hf2
2970 156 : IF (.NOT. is_identical) RETURN
2971 :
2972 98 : DO irep = 1, n_rep_hf1
2973 52 : CALL section_vals_val_get(hfx_section1, "PW_HFX", l_val=lval1, i_rep_section=irep)
2974 52 : CALL section_vals_val_get(hfx_section2, "PW_HFX", l_val=lval2, i_rep_section=irep)
2975 52 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
2976 :
2977 52 : CALL section_vals_val_get(hfx_section1, "PW_HFX_BLOCKSIZE", i_val=ival1, i_rep_section=irep)
2978 52 : CALL section_vals_val_get(hfx_section2, "PW_HFX_BLOCKSIZE", i_val=ival2, i_rep_section=irep)
2979 52 : IF (ival1 .NE. ival2) is_identical = .FALSE.
2980 :
2981 52 : CALL section_vals_val_get(hfx_section1, "TREAT_LSD_IN_CORE", l_val=lval1, i_rep_section=irep)
2982 52 : CALL section_vals_val_get(hfx_section2, "TREAT_LSD_IN_CORE", l_val=lval2, i_rep_section=irep)
2983 52 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
2984 :
2985 52 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "INTERACTION_POTENTIAL", i_rep_section=irep)
2986 52 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "INTERACTION_POTENTIAL", i_rep_section=irep)
2987 :
2988 52 : CALL section_vals_val_get(hfx_sub_section1, "OMEGA", r_val=rval1, i_rep_section=irep)
2989 52 : CALL section_vals_val_get(hfx_sub_section2, "OMEGA", r_val=rval2, i_rep_section=irep)
2990 52 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
2991 :
2992 52 : CALL section_vals_val_get(hfx_sub_section1, "POTENTIAL_TYPE", i_val=ival1, i_rep_section=irep)
2993 52 : CALL section_vals_val_get(hfx_sub_section2, "POTENTIAL_TYPE", i_val=ival2, i_rep_section=irep)
2994 52 : IF (ival1 .NE. ival2) is_identical = .FALSE.
2995 52 : IF (.NOT. is_identical) RETURN
2996 :
2997 46 : IF (ival1 == do_potential_truncated .OR. ival1 == do_potential_mix_cl_trunc) THEN
2998 6 : CALL section_vals_val_get(hfx_sub_section1, "CUTOFF_RADIUS", r_val=rval1, i_rep_section=irep)
2999 6 : CALL section_vals_val_get(hfx_sub_section2, "CUTOFF_RADIUS", r_val=rval2, i_rep_section=irep)
3000 6 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3001 :
3002 6 : CALL section_vals_val_get(hfx_sub_section1, "T_C_G_DATA", c_val=cval1, i_rep_section=irep)
3003 6 : CALL section_vals_val_get(hfx_sub_section2, "T_C_G_DATA", c_val=cval2, i_rep_section=irep)
3004 6 : IF (cval1 .NE. cval2) is_identical = .FALSE.
3005 : END IF
3006 :
3007 46 : CALL section_vals_val_get(hfx_sub_section1, "SCALE_COULOMB", r_val=rval1, i_rep_section=irep)
3008 46 : CALL section_vals_val_get(hfx_sub_section2, "SCALE_COULOMB", r_val=rval2, i_rep_section=irep)
3009 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3010 :
3011 46 : CALL section_vals_val_get(hfx_sub_section1, "SCALE_GAUSSIAN", r_val=rval1, i_rep_section=irep)
3012 46 : CALL section_vals_val_get(hfx_sub_section2, "SCALE_GAUSSIAN", r_val=rval2, i_rep_section=irep)
3013 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3014 :
3015 46 : CALL section_vals_val_get(hfx_sub_section1, "SCALE_LONGRANGE", r_val=rval1, i_rep_section=irep)
3016 46 : CALL section_vals_val_get(hfx_sub_section2, "SCALE_LONGRANGE", r_val=rval2, i_rep_section=irep)
3017 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3018 :
3019 46 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "PERIODIC", i_rep_section=irep)
3020 46 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "PERIODIC", i_rep_section=irep)
3021 :
3022 46 : CALL section_vals_val_get(hfx_sub_section1, "NUMBER_OF_SHELLS", i_val=ival1, i_rep_section=irep)
3023 46 : CALL section_vals_val_get(hfx_sub_section2, "NUMBER_OF_SHELLS", i_val=ival2, i_rep_section=irep)
3024 46 : IF (ival1 .NE. ival2) is_identical = .FALSE.
3025 :
3026 46 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "RI", i_rep_section=irep)
3027 46 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "RI", i_rep_section=irep)
3028 :
3029 46 : CALL section_vals_val_get(hfx_sub_section1, "_SECTION_PARAMETERS_", l_val=lval1, i_rep_section=irep)
3030 46 : CALL section_vals_val_get(hfx_sub_section2, "_SECTION_PARAMETERS_", l_val=lval2, i_rep_section=irep)
3031 46 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3032 :
3033 46 : CALL section_vals_val_get(hfx_sub_section1, "CUTOFF_RADIUS", r_val=rval1, i_rep_section=irep)
3034 46 : CALL section_vals_val_get(hfx_sub_section2, "CUTOFF_RADIUS", r_val=rval2, i_rep_section=irep)
3035 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3036 :
3037 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_EIGVAL", r_val=rval1, i_rep_section=irep)
3038 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_EIGVAL", r_val=rval2, i_rep_section=irep)
3039 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3040 :
3041 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER", r_val=rval1, i_rep_section=irep)
3042 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER", r_val=rval2, i_rep_section=irep)
3043 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3044 :
3045 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER_2C", r_val=rval1, i_rep_section=irep)
3046 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER_2C", r_val=rval2, i_rep_section=irep)
3047 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3048 :
3049 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER_MO", r_val=rval1, i_rep_section=irep)
3050 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER_MO", r_val=rval2, i_rep_section=irep)
3051 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3052 :
3053 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_PGF_ORB", r_val=rval1, i_rep_section=irep)
3054 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_PGF_ORB", r_val=rval2, i_rep_section=irep)
3055 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3056 :
3057 46 : CALL section_vals_val_get(hfx_sub_section1, "MAX_BLOCK_SIZE_MO", i_val=ival1, i_rep_section=irep)
3058 46 : CALL section_vals_val_get(hfx_sub_section2, "MAX_BLOCK_SIZE_MO", i_val=ival2, i_rep_section=irep)
3059 46 : IF (ival1 .NE. ival2) is_identical = .FALSE.
3060 :
3061 46 : CALL section_vals_val_get(hfx_sub_section1, "MIN_BLOCK_SIZE", i_val=ival1, i_rep_section=irep)
3062 46 : CALL section_vals_val_get(hfx_sub_section2, "MIN_BLOCK_SIZE", i_val=ival2, i_rep_section=irep)
3063 46 : IF (ival1 .NE. ival2) is_identical = .FALSE.
3064 :
3065 46 : CALL section_vals_val_get(hfx_sub_section1, "OMEGA", r_val=rval1, i_rep_section=irep)
3066 46 : CALL section_vals_val_get(hfx_sub_section2, "OMEGA", r_val=rval2, i_rep_section=irep)
3067 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3068 :
3069 46 : CALL section_vals_val_get(hfx_sub_section1, "RI_FLAVOR", i_val=ival1, i_rep_section=irep)
3070 46 : CALL section_vals_val_get(hfx_sub_section2, "RI_FLAVOR", i_val=ival2, i_rep_section=irep)
3071 46 : IF (ival1 .NE. ival2) is_identical = .FALSE.
3072 :
3073 46 : CALL section_vals_val_get(hfx_sub_section1, "RI_METRIC", i_val=ival1, i_rep_section=irep)
3074 46 : CALL section_vals_val_get(hfx_sub_section2, "RI_METRIC", i_val=ival2, i_rep_section=irep)
3075 46 : IF (ival1 .NE. ival2) is_identical = .FALSE.
3076 :
3077 46 : hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "SCREENING", i_rep_section=irep)
3078 46 : hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "SCREENING", i_rep_section=irep)
3079 :
3080 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_SCHWARZ", r_val=rval1, i_rep_section=irep)
3081 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_SCHWARZ", r_val=rval2, i_rep_section=irep)
3082 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3083 :
3084 46 : CALL section_vals_val_get(hfx_sub_section1, "EPS_SCHWARZ_FORCES", r_val=rval1, i_rep_section=irep)
3085 46 : CALL section_vals_val_get(hfx_sub_section2, "EPS_SCHWARZ_FORCES", r_val=rval2, i_rep_section=irep)
3086 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3087 :
3088 46 : CALL section_vals_val_get(hfx_sub_section1, "P_SCREEN_CORRECTION_FACTOR", r_val=rval1, i_rep_section=irep)
3089 46 : CALL section_vals_val_get(hfx_sub_section2, "P_SCREEN_CORRECTION_FACTOR", r_val=rval2, i_rep_section=irep)
3090 46 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3091 :
3092 46 : CALL section_vals_val_get(hfx_sub_section1, "SCREEN_ON_INITIAL_P", l_val=lval1, i_rep_section=irep)
3093 46 : CALL section_vals_val_get(hfx_sub_section2, "SCREEN_ON_INITIAL_P", l_val=lval2, i_rep_section=irep)
3094 46 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3095 :
3096 46 : CALL section_vals_val_get(hfx_sub_section1, "SCREEN_P_FORCES", l_val=lval1, i_rep_section=irep)
3097 46 : CALL section_vals_val_get(hfx_sub_section2, "SCREEN_P_FORCES", l_val=lval2, i_rep_section=irep)
3098 1272 : IF (lval1 .NEQV. lval2) is_identical = .FALSE.
3099 :
3100 : END DO
3101 :
3102 : !Test of the fraction
3103 46 : IF (is_identical) THEN
3104 84 : DO irep = 1, n_rep_hf1
3105 42 : CALL section_vals_val_get(hfx_section1, "FRACTION", r_val=rval1, i_rep_section=irep)
3106 42 : CALL section_vals_val_get(hfx_section2, "FRACTION", r_val=rval2, i_rep_section=irep)
3107 84 : IF (ABS(rval1 - rval2) > EPSILON(1.0_dp)) is_identical = .FALSE.
3108 : END DO
3109 :
3110 42 : IF (PRESENT(same_except_frac)) THEN
3111 30 : IF (.NOT. is_identical) same_except_frac = .TRUE.
3112 : END IF
3113 : END IF
3114 :
3115 : END SUBROUTINE compare_hfx_sections
3116 :
3117 0 : END MODULE hfx_types
3118 :
|
