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