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