(git:97501a3)
Loading...
Searching...
No Matches
hfx_types.F
Go to the documentation of this file.
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! **************************************************************************************************
15MODULE 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
566CONTAINS
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 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 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 CALL timeset(routinen, handle)
622
623 CALL cite_reference(guidon2008)
624 CALL cite_reference(guidon2009)
625
626 natom = SIZE(particle_set)
627
628 !! There might be 2 hf sections
629 CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
630 n_threads = 1
631!$ n_threads = omp_get_max_threads()
632
633 CALL section_vals_val_get(hfx_section, "RI%_SECTION_PARAMETERS_", l_val=do_ri)
634 IF (do_ri) n_threads = 1 ! RI implementation does not use threads
635
636 IF (PRESENT(orb_basis)) THEN
637 orb_basis_type = orb_basis
638 ELSE
639 orb_basis_type = "ORB"
640 END IF
641 IF (PRESENT(ri_basis)) THEN
642 ri_basis_type = ri_basis
643 ELSE
644 ri_basis_type = "RI_HFX"
645 END IF
646
647 ALLOCATE (x_data(n_rep_hf, n_threads))
648 DO i_thread = 1, n_threads
649 DO irep = 1, n_rep_hf
650 actual_x_data => x_data(irep, i_thread)
651 !! Get data from input file
652 !!
653 !! GENERAL params
654 CALL section_vals_val_get(hfx_section, "FRACTION", r_val=real_val, i_rep_section=irep)
655 actual_x_data%general_parameter%fraction = real_val
656 actual_x_data%n_rep_hf = n_rep_hf
657
658 NULLIFY (actual_x_data%map_atoms_to_cpus)
659
660 CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=logic_val, i_rep_section=irep)
661 actual_x_data%general_parameter%treat_lsd_in_core = logic_val
662
663 hfx_ri_section => section_vals_get_subs_vals(hfx_section, "RI")
664 CALL section_vals_val_get(hfx_ri_section, "_SECTION_PARAMETERS_", l_val=actual_x_data%do_hfx_ri)
665
666 !! MEMORY section
667 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 n_threads, para_env, irep, skip_disk=.false., skip_in_core_forces=.false.)
670
671 !! PERIODIC section
672 hf_sub_section => section_vals_get_subs_vals(hfx_section, "PERIODIC", i_rep_section=irep)
673 CALL section_vals_val_get(hf_sub_section, "NUMBER_OF_SHELLS", i_val=int_val)
674 actual_x_data%periodic_parameter%number_of_shells = int_val
675 actual_x_data%periodic_parameter%mode = int_val
676 CALL get_cell(cell=cell, periodic=actual_x_data%periodic_parameter%perd)
677 IF (sum(actual_x_data%periodic_parameter%perd) == 0) THEN
678 actual_x_data%periodic_parameter%do_periodic = .false.
679 ELSE
680 actual_x_data%periodic_parameter%do_periodic = .true.
681 END IF
682
683 !! SCREENING section
684 hf_sub_section => section_vals_get_subs_vals(hfx_section, "SCREENING", i_rep_section=irep)
685 CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ", r_val=real_val)
686 actual_x_data%screening_parameter%eps_schwarz = real_val
687 CALL section_vals_val_get(hf_sub_section, "EPS_SCHWARZ_FORCES", r_val=real_val, explicit=explicit)
688 IF (explicit) THEN
689 actual_x_data%screening_parameter%eps_schwarz_forces = real_val
690 ELSE
691 actual_x_data%screening_parameter%eps_schwarz_forces = &
692 100._dp*actual_x_data%screening_parameter%eps_schwarz
693 END IF
694 CALL section_vals_val_get(hf_sub_section, "SCREEN_P_FORCES", l_val=logic_val)
695 actual_x_data%screening_parameter%do_p_screening_forces = logic_val
696 CALL section_vals_val_get(hf_sub_section, "SCREEN_ON_INITIAL_P", l_val=logic_val)
697 actual_x_data%screening_parameter%do_initial_p_screening = logic_val
698 actual_x_data%screen_funct_is_initialized = .false.
699
700 !! INTERACTION_POTENTIAL section
701 hf_sub_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL", i_rep_section=irep)
702 CALL section_vals_val_get(hf_sub_section, "POTENTIAL_TYPE", i_val=int_val)
703 actual_x_data%potential_parameter%potential_type = int_val
704 CALL section_vals_val_get(hf_sub_section, "OMEGA", r_val=real_val)
705 actual_x_data%potential_parameter%omega = real_val
706 CALL section_vals_val_get(hf_sub_section, "SCALE_COULOMB", r_val=real_val)
707 actual_x_data%potential_parameter%scale_coulomb = real_val
708 CALL section_vals_val_get(hf_sub_section, "SCALE_LONGRANGE", r_val=real_val)
709 actual_x_data%potential_parameter%scale_longrange = real_val
710 CALL section_vals_val_get(hf_sub_section, "SCALE_GAUSSIAN", r_val=real_val)
711 actual_x_data%potential_parameter%scale_gaussian = real_val
712 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 CALL section_vals_val_get(hf_sub_section, "CUTOFF_RADIUS", r_val=real_val)
715 actual_x_data%potential_parameter%cutoff_radius = real_val
716 CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val)
717 CALL compress(char_val, .true.)
718 ! ** Check if file is there
719 IF (.NOT. file_exists(char_val)) THEN
720 WRITE (error_msg, '(A,A,A)') "Truncated hfx calculation requested. The file containing "// &
721 "the data could not be found at ", trim(char_val), " Please check T_C_G_DATA "// &
722 "in the INTERACTION_POTENTIAL section"
723 cpabort(error_msg)
724 ELSE
725 actual_x_data%potential_parameter%filename = char_val
726 END IF
727 END IF
728 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 actual_x_data%potential_parameter%cutoff_radius)
732 END IF
733 IF (actual_x_data%potential_parameter%potential_type == do_potential_id) THEN
734 actual_x_data%potential_parameter%cutoff_radius = 0.0_dp
735 END IF
736
737 !! LOAD_BALANCE section
738 hf_sub_section => section_vals_get_subs_vals(hfx_section, "LOAD_BALANCE", i_rep_section=irep)
739 CALL section_vals_val_get(hf_sub_section, "NBINS", i_val=int_val)
740 actual_x_data%load_balance_parameter%nbins = max(1, int_val)
741 actual_x_data%load_balance_parameter%blocks_initialized = .false.
742
743 CALL section_vals_val_get(hf_sub_section, "RANDOMIZE", l_val=logic_val)
744 actual_x_data%load_balance_parameter%do_randomize = logic_val
745
746 actual_x_data%load_balance_parameter%rtp_redistribute = .false.
747 IF (ASSOCIATED(dft_control%rtp_control)) &
748 actual_x_data%load_balance_parameter%rtp_redistribute = dft_control%rtp_control%hfx_redistribute
749
750 CALL section_vals_val_get(hf_sub_section, "BLOCK_SIZE", i_val=int_val)
751 ! negative values ask for a computed default
752 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 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 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 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 IF (actual_x_data%periodic_parameter%do_periodic) THEN
861 hf_pbc_section => section_vals_get_subs_vals(hfx_section, "PERIODIC", i_rep_section=irep)
862 CALL section_vals_val_get(hf_pbc_section, "NUMBER_OF_SHELLS", i_val=pbc_shells)
863 actual_x_data%periodic_parameter%number_of_shells_from_input = pbc_shells
864 ALLOCATE (actual_x_data%neighbor_cells(1))
865 CALL hfx_create_neighbor_cells(actual_x_data, pbc_shells, cell, i_thread, nkp_grid=nkp_grid)
866 ELSE
867 ALLOCATE (actual_x_data%neighbor_cells(1))
868 ! ** Initialize this guy to enable non periodic stress regtests
869 actual_x_data%periodic_parameter%R_max_stress = 1.0_dp
870 END IF
871
872 nkind = SIZE(qs_kind_set, 1)
873 max_set = actual_x_data%basis_info%max_set
874
875 !! ** This guy is allocated on the master thread only
876 IF (i_thread == 1) THEN
877 ALLOCATE (actual_x_data%is_assoc_atomic_block(natom, natom))
878 ALLOCATE (actual_x_data%atomic_block_offset(natom, natom))
879 ALLOCATE (actual_x_data%set_offset(max_set, max_set, nkind, nkind))
880 ALLOCATE (actual_x_data%block_offset(para_env%num_pe + 1))
881 END IF
882
883 ALLOCATE (actual_x_data%distribution_forces(1))
884 ALLOCATE (actual_x_data%distribution_energy(1))
885
886 actual_x_data%memory_parameter%size_p_screen = 0_int_8
887 IF (i_thread == 1) THEN
888 ALLOCATE (actual_x_data%atomic_pair_list(natom, natom))
889 ALLOCATE (actual_x_data%atomic_pair_list_forces(natom, natom))
890 END IF
891
892 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 IF (i_thread == 1) THEN
896 ALLOCATE (actual_x_data%pmax_atom(natom, natom))
897 ALLOCATE (actual_x_data%initial_p(nkind*(nkind + 1)/2))
898 i = 1
899 DO ikind = 1, nkind
900 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_a)
901 nseta = actual_x_data%basis_parameter(ikind)%nset
902 DO jkind = ikind, nkind
903 CALL get_atomic_kind(atomic_kind_set(jkind), natom=natom_b)
904 nsetb = actual_x_data%basis_parameter(jkind)%nset
905 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 actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b
908 i = i + 1
909 END DO
910 END DO
911
912 ALLOCATE (actual_x_data%pmax_atom_forces(natom, natom))
913 ALLOCATE (actual_x_data%initial_p_forces(nkind*(nkind + 1)/2))
914 i = 1
915 DO ikind = 1, nkind
916 CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_a)
917 nseta = actual_x_data%basis_parameter(ikind)%nset
918 DO jkind = ikind, nkind
919 CALL get_atomic_kind(atomic_kind_set(jkind), natom=natom_b)
920 nsetb = actual_x_data%basis_parameter(jkind)%nset
921 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 actual_x_data%memory_parameter%size_p_screen + nseta*nsetb*natom_a*natom_b
924 i = i + 1
925 END DO
926 END DO
927 END IF
928 ALLOCATE (actual_x_data%map_atom_to_kind_atom(natom))
929 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
930
931 ALLOCATE (atom2kind(nkind))
932 atom2kind = 0
933 DO iatom = 1, natom
934 ikind = kind_of(iatom)
935 atom2kind(ikind) = atom2kind(ikind) + 1
936 actual_x_data%map_atom_to_kind_atom(iatom) = atom2kind(ikind)
937 END DO
938 DEALLOCATE (kind_of, atom2kind)
939 END IF
940
941 ! ** Initialize libint type
942 CALL cp_libint_static_init()
943 CALL cp_libint_init_eri(actual_x_data%lib, actual_x_data%basis_info%max_am)
944 CALL cp_libint_init_eri1(actual_x_data%lib_deriv, actual_x_data%basis_info%max_am)
945 CALL cp_libint_set_contrdepth(actual_x_data%lib, 1)
946 CALL cp_libint_set_contrdepth(actual_x_data%lib_deriv, 1)
947
948 CALL alloc_containers(actual_x_data%store_ints, 1)
949 CALL alloc_containers(actual_x_data%store_forces, 1)
950
951 actual_x_data%store_ints%maxval_cache_disk%element_counter = 1
952 ALLOCATE (actual_x_data%store_ints%maxval_container_disk)
953 ALLOCATE (actual_x_data%store_ints%maxval_container_disk%first)
954 actual_x_data%store_ints%maxval_container_disk%first%prev => null()
955 actual_x_data%store_ints%maxval_container_disk%first%next => null()
956 actual_x_data%store_ints%maxval_container_disk%current => actual_x_data%store_ints%maxval_container_disk%first
957 actual_x_data%store_ints%maxval_container_disk%current%data = 0
958 actual_x_data%store_ints%maxval_container_disk%element_counter = 1
959 actual_x_data%store_ints%maxval_container_disk%file_counter = 1
960 actual_x_data%store_ints%maxval_container_disk%desc = 'Max_'
961 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 trim(actual_x_data%memory_parameter%storage_location), &
964 storage_id, "_", actual_x_data%store_ints%maxval_container_disk%desc, "6"
965 CALL compress(actual_x_data%store_ints%maxval_container_disk%filename, .true.)
966 ALLOCATE (actual_x_data%store_ints%integral_containers_disk(64))
967 DO i = 1, 64
968 actual_x_data%store_ints%integral_caches_disk(i)%element_counter = 1
969 actual_x_data%store_ints%integral_caches_disk(i)%data = 0
970 ALLOCATE (actual_x_data%store_ints%integral_containers_disk(i)%first)
971 actual_x_data%store_ints%integral_containers_disk(i)%first%prev => null()
972 actual_x_data%store_ints%integral_containers_disk(i)%first%next => null()
973 actual_x_data%store_ints%integral_containers_disk(i)%current => &
974 actual_x_data%store_ints%integral_containers_disk(i)%first
975 actual_x_data%store_ints%integral_containers_disk(i)%current%data = 0
976 actual_x_data%store_ints%integral_containers_disk(i)%element_counter = 1
977 actual_x_data%store_ints%integral_containers_disk(i)%file_counter = 1
978 actual_x_data%store_ints%integral_containers_disk(i)%desc = 'Int_'
979 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 trim(actual_x_data%memory_parameter%storage_location), &
982 storage_id, "_", actual_x_data%store_ints%integral_containers_disk(i)%desc, i
983 CALL compress(actual_x_data%store_ints%integral_containers_disk(i)%filename, .true.)
984 END DO
985
986 actual_x_data%b_first_load_balance_energy = .true.
987 actual_x_data%b_first_load_balance_forces = .true.
988
989 hf_sub_section => section_vals_get_subs_vals(hfx_section, "RI", i_rep_section=irep)
990 IF (actual_x_data%do_hfx_ri) THEN
991 cpassert(PRESENT(nelectron_total))
992 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 nelectron_total, orb_basis_type, ri_basis_type)
997 END IF
998 END DO
999 END DO
1000
1001 DO irep = 1, n_rep_hf
1002 actual_x_data => x_data(irep, 1)
1003 CALL hfx_print_info(actual_x_data, hfx_section, irep)
1004 END DO
1005
1006 CALL timestop(handle)
1007
1008 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 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 CALL timeset(routinen, handle)
1049
1050 NULLIFY (hf_sub_section)
1051
1052 associate(hfx_pot => ri_data%hfx_pot)
1053 hfx_pot%potential_type = x_data%potential_parameter%potential_type
1054 hfx_pot%omega = x_data%potential_parameter%omega
1055 hfx_pot%cutoff_radius = x_data%potential_parameter%cutoff_radius
1056 hfx_pot%scale_coulomb = x_data%potential_parameter%scale_coulomb
1057 hfx_pot%scale_longrange = x_data%potential_parameter%scale_longrange
1058 END associate
1059 ri_data%ri_section => ri_section
1060 ri_data%hfx_section => hfx_section
1061 ri_data%eps_schwarz = x_data%screening_parameter%eps_schwarz
1062 ri_data%eps_schwarz_forces = x_data%screening_parameter%eps_schwarz_forces
1063
1064 logger => cp_get_default_logger()
1065 unit_nr_dbcsr = cp_print_key_unit_nr(logger, ri_data%ri_section, "PRINT%RI_INFO", &
1066 extension=".dbcsrLog")
1067
1068 unit_nr = cp_print_key_unit_nr(logger, ri_data%hfx_section, "HF_INFO", &
1069 extension=".scfLog")
1070
1071 hf_sub_section => section_vals_get_subs_vals(hfx_section, "INTERACTION_POTENTIAL", i_rep_section=irep)
1072 CALL section_vals_val_get(hf_sub_section, "T_C_G_DATA", c_val=char_val)
1073 CALL compress(char_val, .true.)
1074
1075 IF (.NOT. file_exists(char_val)) THEN
1076 WRITE (error_msg, '(A,A,A)') "File not found. Please check T_C_G_DATA "// &
1077 "in the INTERACTION_POTENTIAL section"
1078 cpabort(error_msg)
1079 ELSE
1080 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 nelectron_total, t_c_filename=t_c_filename)
1086
1087 IF (dft_control%smear .AND. ri_data%flavor == ri_mo) THEN
1088 cpabort("RI_FLAVOR MO is not consistent with smearing. Please use RI_FLAVOR RHO.")
1089 END IF
1090
1091 CALL timestop(handle)
1092
1093 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 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 CALL timeset(routinen, handle)
1130
1131 CALL section_vals_val_get(ri_section, "EPS_FILTER", r_val=ri_data%filter_eps)
1132 CALL section_vals_val_get(ri_section, "EPS_FILTER_2C", r_val=ri_data%filter_eps_2c)
1133 CALL section_vals_val_get(ri_section, "EPS_STORAGE_SCALING", r_val=eps_storage_scaling)
1134 ri_data%filter_eps_storage = ri_data%filter_eps*eps_storage_scaling
1135 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 CALL section_vals_val_get(ri_section, "RI_METRIC", i_val=ri_metric%potential_type, explicit=explicit)
1139 IF (.NOT. explicit .OR. ri_metric%potential_type == 0) THEN
1140 ri_metric%potential_type = hfx_pot%potential_type
1141 END IF
1142
1143 CALL section_vals_val_get(ri_section, "OMEGA", r_val=ri_metric%omega, explicit=explicit)
1144 IF (.NOT. explicit) THEN
1145 ri_metric%omega = hfx_pot%omega
1146 END IF
1147
1148 CALL section_vals_val_get(ri_section, "CUTOFF_RADIUS", r_val=ri_metric%cutoff_radius, explicit=explicit)
1149 IF (.NOT. explicit) THEN
1150 ri_metric%cutoff_radius = hfx_pot%cutoff_radius
1151 END IF
1152
1153 CALL section_vals_val_get(ri_section, "SCALE_COULOMB", r_val=ri_metric%scale_coulomb, explicit=explicit)
1154 IF (.NOT. explicit) THEN
1155 ri_metric%scale_coulomb = hfx_pot%scale_coulomb
1156 END IF
1157
1158 CALL section_vals_val_get(ri_section, "SCALE_LONGRANGE", r_val=ri_metric%scale_longrange, explicit=explicit)
1159 IF (.NOT. explicit) THEN
1160 ri_metric%scale_longrange = hfx_pot%scale_longrange
1161 END IF
1162
1163 IF (ri_metric%potential_type == do_potential_short) &
1164 CALL erfc_cutoff(ri_data%eps_schwarz, ri_metric%omega, ri_metric%cutoff_radius)
1165 IF (ri_metric%potential_type == do_potential_id) ri_metric%cutoff_radius = 0.0_dp
1166 END associate
1167
1168 CALL section_vals_val_get(ri_section, "2C_MATRIX_FUNCTIONS", i_val=ri_data%t2c_method)
1169 CALL section_vals_val_get(ri_section, "EPS_EIGVAL", r_val=ri_data%eps_eigval)
1170 CALL section_vals_val_get(ri_section, "CHECK_2C_MATRIX", l_val=ri_data%check_2c_inv)
1171 CALL section_vals_val_get(ri_section, "CALC_COND_NUM", l_val=ri_data%calc_condnum)
1172 CALL section_vals_val_get(ri_section, "SQRT_ORDER", i_val=ri_data%t2c_sqrt_order)
1173 CALL section_vals_val_get(ri_section, "EPS_LANCZOS", r_val=ri_data%eps_lanczos)
1174 CALL section_vals_val_get(ri_section, "MAX_ITER_LANCZOS", i_val=ri_data%max_iter_lanczos)
1175 CALL section_vals_val_get(ri_section, "RI_FLAVOR", i_val=ri_data%flavor)
1176 CALL section_vals_val_get(ri_section, "EPS_PGF_ORB", r_val=ri_data%eps_pgf_orb)
1177 CALL section_vals_val_get(ri_section, "MIN_BLOCK_SIZE", i_val=ri_data%min_bsize)
1178 CALL section_vals_val_get(ri_section, "MAX_BLOCK_SIZE_MO", i_val=ri_data%max_bsize_MO)
1179 CALL section_vals_val_get(ri_section, "MEMORY_CUT", i_val=ri_data%n_mem_input)
1180 CALL section_vals_val_get(ri_section, "FLAVOR_SWITCH_MEMORY_CUT", i_val=ri_data%n_mem_flavor_switch)
1181
1182 ri_data%orb_basis_type = orb_basis_type
1183 ri_data%ri_basis_type = ri_basis_type
1184 ri_data%nelectron_total = nelectron_total
1185 ri_data%input_flavor = ri_data%flavor
1186
1187 IF (PRESENT(t_c_filename)) THEN
1188 ri_data%ri_metric%filename = t_c_filename
1189 ri_data%hfx_pot%filename = t_c_filename
1190 END IF
1191
1192 ri_data%unit_nr_dbcsr = unit_nr_dbcsr
1193 ri_data%unit_nr = unit_nr
1194 ri_data%dbcsr_nflop = 0
1195 ri_data%dbcsr_time = 0.0_dp
1196
1197 CALL hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
1198
1199 CALL timestop(handle)
1200
1201 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: bsizes_ao_store, bsizes_ri_store, dist1, &
1223 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 DIMENSION(:) :: basis_set_ao, basis_set_ri
1231 TYPE(mp_cart_type) :: mp_comm_3d
1232
1233 CALL cite_reference(bussy2023)
1234
1235 CALL timeset(routinen, handle)
1236
1237 ! initialize libint
1238 CALL cp_libint_static_init()
1239
1240 natom = SIZE(particle_set)
1241 nkind = SIZE(qs_kind_set, 1)
1242 nproc = para_env%num_pe
1243
1244 associate(ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
1245 IF (ri_metric%potential_type == do_potential_short) THEN
1246 CALL erfc_cutoff(ri_data%eps_schwarz, ri_metric%omega, ri_metric%cutoff_radius)
1247 END IF
1248
1249 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 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 same_op = compare_potential_types(ri_metric, hfx_pot)
1256 END associate
1257
1258 ri_data%same_op = same_op
1259
1260 pdims = 0
1261 CALL mp_comm_3d%create(para_env, 3, pdims)
1262
1263 ALLOCATE (ri_data%bsizes_RI(natom))
1264 ALLOCATE (ri_data%bsizes_AO(natom))
1265 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
1266 CALL basis_set_list_setup(basis_set_ri, ri_data%ri_basis_type, qs_kind_set)
1267 CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_data%bsizes_RI, basis=basis_set_ri)
1268 CALL basis_set_list_setup(basis_set_ao, ri_data%orb_basis_type, qs_kind_set)
1269 CALL get_particle_set(particle_set, qs_kind_set, nsgf=ri_data%bsizes_AO, basis=basis_set_ao)
1270
1271 ALLOCATE (dist_ri(natom))
1272 ALLOCATE (dist_ao_1(natom))
1273 ALLOCATE (dist_ao_2(natom))
1274 CALL dbt_default_distvec(natom, pdims(1), ri_data%bsizes_RI, dist_ri)
1275 CALL dbt_default_distvec(natom, pdims(2), ri_data%bsizes_AO, dist_ao_1)
1276 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 mp_comm_3d, own_comm=.true.)
1279
1280 ALLOCATE (ri_data%pgrid)
1281 CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid)
1282
1283 ALLOCATE (ri_data%pgrid_2d)
1284 pdims_2d = 0
1285 CALL dbt_pgrid_create(para_env, pdims_2d, ri_data%pgrid_2d)
1286
1287 ri_data%dist_3d = dist_3d
1288
1289 CALL dbt_distribution_new(ri_data%dist, ri_data%pgrid, &
1290 dist_ri, dist_ao_1, dist_ao_2)
1291
1292 DEALLOCATE (dist_ao_1, dist_ao_2, dist_ri)
1293
1294 ri_data%num_pe = para_env%num_pe
1295
1296 ! initialize tensors expressed in basis representation
1297 CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, ri_data%min_bsize, ri_data%bsizes_AO_split)
1298 CALL pgf_block_sizes(atomic_kind_set, basis_set_ri, ri_data%min_bsize, ri_data%bsizes_RI_split)
1299
1300 CALL pgf_block_sizes(atomic_kind_set, basis_set_ao, 1, bsizes_ao_store)
1301 CALL pgf_block_sizes(atomic_kind_set, basis_set_ri, 1, bsizes_ri_store)
1302
1303 CALL split_block_sizes([sum(ri_data%bsizes_AO)], ri_data%bsizes_AO_fit, default_block_size)
1304 CALL split_block_sizes([sum(ri_data%bsizes_RI)], ri_data%bsizes_RI_fit, default_block_size)
1305
1306 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 ri_data%n_mem = ri_data%n_mem_input
1310 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 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 ri_data%starts_array_RI_mem_block, ri_data%ends_array_RI_mem_block)
1319
1320 ALLOCATE (ri_data%pgrid_1)
1321 ALLOCATE (ri_data%pgrid_2)
1322 pdims = 0
1323
1324 CALL dbt_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), &
1325 SIZE(ri_data%bsizes_AO_split)])
1326
1327 CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_1)
1328
1329 pdims = pdims([2, 1, 3])
1330 CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_2)
1331
1332 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 ri_data%bsizes_AO_split, [1, 2], [3], name="(AO RI | AO)")
1336 DEALLOCATE (dist1, dist2, dist3)
1337
1338 ALLOCATE (ri_data%blk_indices(ri_data%n_mem, ri_data%n_mem_RI))
1339 ALLOCATE (ri_data%store_3c(ri_data%n_mem, ri_data%n_mem_RI))
1340 DO i_mem = 1, ri_data%n_mem
1341 DO j_mem = 1, ri_data%n_mem_RI
1342 CALL alloc_containers(ri_data%store_3c(i_mem, j_mem), 1)
1343 END DO
1344 END DO
1345
1346 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 ri_data%bsizes_AO_split, [1, 2], [3], name="(AO RI | AO)")
1350 DEALLOCATE (dist1, dist2, dist3)
1351
1352 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 ri_data%bsizes_AO_split, [1], [2, 3], name="(RI | AO AO)")
1356 DEALLOCATE (dist1, dist2, dist3)
1357
1358 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 name="(RI | RI)")
1362 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 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 name="(AO | AO)")
1369 DEALLOCATE (dist1, dist2)
1370 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
1371
1372 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 name="(AO | AO)")
1376 DEALLOCATE (dist1, dist2)
1377 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
1378
1379 ELSEIF (ri_data%flavor == ri_mo) THEN
1380 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 name="(RI | RI)")
1385 CALL dbt_create(ri_data%t_2c_int(1, 1), ri_data%t_2c_int(2, 1))
1386
1387 DEALLOCATE (dist1, dist2)
1388
1389 ALLOCATE (ri_data%t_3c_int_ctr_1(1, 1))
1390
1391 ALLOCATE (ri_data%pgrid_1)
1392 ALLOCATE (ri_data%pgrid_2)
1393 pdims = 0
1394
1395 ri_data%n_mem = ri_data%n_mem_input**2
1396 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 mo_dim = max((ri_data%nelectron_total/2 - 1)/ri_data%n_mem + 1, 1)
1402 mo_dim = (mo_dim - 1)/ri_data%max_bsize_MO + 1
1403
1404 pdims = 0
1405 CALL dbt_mp_dims_create(nproc, pdims, [SIZE(ri_data%bsizes_AO_split), SIZE(ri_data%bsizes_RI_split), mo_dim])
1406
1407 CALL dbt_pgrid_create(para_env, pdims, ri_data%pgrid_1)
1408
1409 pdims = pdims([3, 2, 1])
1410 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 [1, 2], [3], name="(AO RI | AO)")
1415 DEALLOCATE (dist1, dist2, dist3)
1416
1417 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 [1], [2, 3], name="(AO | RI AO)")
1421 DEALLOCATE (dist1, dist2, dist3)
1422
1423 END IF
1424
1425 !For forces
1426 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 name="(RI | RI)")
1430 DEALLOCATE (dist1, dist2)
1431
1432 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 name="(RI | RI)")
1436 DEALLOCATE (dist1, dist2)
1437
1438 CALL timestop(handle)
1439
1440 END SUBROUTINE
1441
1442! **************************************************************************************************
1443!> \brief ...
1444!> \param ri_data ...
1445! **************************************************************************************************
1446 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 my_flop_rate = real(dbcsr_nflop, dp)/(1.0e09_dp*ri_data%dbcsr_time)
1454 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(/T2,A,T73,ES8.2)") &
1455 "RI-HFX PERFORMANCE| DBT total number of flops:", real(dbcsr_nflop*num_pe, dp)
1456 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T2,A,T66,F15.2)") &
1457 "RI-HFX PERFORMANCE| DBT total execution time:", dbcsr_time
1458 IF (unit_nr > 0) WRITE (unit=unit_nr, fmt="(T2,A,T66,F15.2)") &
1459 "RI-HFX PERFORMANCE| DBT flop rate (Gflops / MPI rank):", my_flop_rate
1460 END associate
1461 END SUBROUTINE
1462
1463! **************************************************************************************************
1464!> \brief ...
1465!> \param ri_data ...
1466!> \param write_stats ...
1467! **************************************************************************************************
1468 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 CALL timeset(routinen, handle)
1478
1479 ! cleanup libint
1480 CALL cp_libint_static_cleanup()
1481
1482 my_write_stats = .true.
1483 IF (PRESENT(write_stats)) my_write_stats = write_stats
1484 IF (my_write_stats) CALL hfx_ri_write_stats(ri_data)
1485
1486 IF (ASSOCIATED(ri_data%pgrid)) THEN
1487 CALL dbt_pgrid_destroy(ri_data%pgrid)
1488 DEALLOCATE (ri_data%pgrid)
1489 END IF
1490 IF (ASSOCIATED(ri_data%pgrid_1)) THEN
1491 CALL dbt_pgrid_destroy(ri_data%pgrid_1)
1492 DEALLOCATE (ri_data%pgrid_1)
1493 END IF
1494 IF (ASSOCIATED(ri_data%pgrid_2)) THEN
1495 CALL dbt_pgrid_destroy(ri_data%pgrid_2)
1496 DEALLOCATE (ri_data%pgrid_2)
1497 END IF
1498 IF (ASSOCIATED(ri_data%pgrid_2d)) THEN
1499 CALL dbt_pgrid_destroy(ri_data%pgrid_2d)
1500 DEALLOCATE (ri_data%pgrid_2d)
1501 END IF
1502
1503 CALL distribution_3d_destroy(ri_data%dist_3d)
1504 CALL dbt_distribution_destroy(ri_data%dist)
1505
1506 DEALLOCATE (ri_data%bsizes_RI)
1507 DEALLOCATE (ri_data%bsizes_AO)
1508 DEALLOCATE (ri_data%bsizes_AO_split)
1509 DEALLOCATE (ri_data%bsizes_RI_split)
1510 DEALLOCATE (ri_data%bsizes_AO_fit)
1511 DEALLOCATE (ri_data%bsizes_RI_fit)
1512
1513 IF (ri_data%flavor == ri_pmat) THEN
1514 DO i_mem = 1, ri_data%n_mem
1515 DO j_mem = 1, ri_data%n_mem_RI
1516 CALL dealloc_containers(ri_data%store_3c(i_mem, j_mem), unused)
1517 END DO
1518 END DO
1519
1520 DO j = 1, SIZE(ri_data%t_3c_int_ctr_1, 2)
1521 DO i = 1, SIZE(ri_data%t_3c_int_ctr_1, 1)
1522 CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
1523 END DO
1524 END DO
1525 DEALLOCATE (ri_data%t_3c_int_ctr_1)
1526
1527 DO j = 1, SIZE(ri_data%t_3c_int_ctr_2, 2)
1528 DO i = 1, SIZE(ri_data%t_3c_int_ctr_2, 1)
1529 CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
1530 END DO
1531 END DO
1532 DEALLOCATE (ri_data%t_3c_int_ctr_2)
1533
1534 DO j = 1, SIZE(ri_data%t_3c_int_ctr_3, 2)
1535 DO i = 1, SIZE(ri_data%t_3c_int_ctr_3, 1)
1536 CALL dbt_destroy(ri_data%t_3c_int_ctr_3(i, j))
1537 END DO
1538 END DO
1539 DEALLOCATE (ri_data%t_3c_int_ctr_3)
1540
1541 DO j = 1, SIZE(ri_data%t_2c_int, 2)
1542 DO i = 1, SIZE(ri_data%t_2c_int, 1)
1543 CALL dbt_destroy(ri_data%t_2c_int(i, j))
1544 END DO
1545 END DO
1546 DEALLOCATE (ri_data%t_2c_int)
1547
1548 DO j = 1, SIZE(ri_data%rho_ao_t, 2)
1549 DO i = 1, SIZE(ri_data%rho_ao_t, 1)
1550 CALL dbt_destroy(ri_data%rho_ao_t(i, j))
1551 END DO
1552 END DO
1553 DEALLOCATE (ri_data%rho_ao_t)
1554
1555 DO j = 1, SIZE(ri_data%ks_t, 2)
1556 DO i = 1, SIZE(ri_data%ks_t, 1)
1557 CALL dbt_destroy(ri_data%ks_t(i, j))
1558 END DO
1559 END DO
1560 DEALLOCATE (ri_data%ks_t)
1561
1562 DEALLOCATE (ri_data%starts_array_mem_block, ri_data%ends_array_mem_block, &
1563 ri_data%starts_array_mem, ri_data%ends_array_mem)
1564 DEALLOCATE (ri_data%starts_array_RI_mem_block, ri_data%ends_array_RI_mem_block, &
1565 ri_data%starts_array_RI_mem, ri_data%ends_array_RI_mem)
1566
1567 DEALLOCATE (ri_data%blk_indices)
1568 DEALLOCATE (ri_data%store_3c)
1569 ELSEIF (ri_data%flavor == ri_mo) THEN
1570 CALL dbt_destroy(ri_data%t_3c_int_ctr_1(1, 1))
1571 CALL dbt_destroy(ri_data%t_3c_int_ctr_2(1, 1))
1572 DEALLOCATE (ri_data%t_3c_int_ctr_1)
1573 DEALLOCATE (ri_data%t_3c_int_ctr_2)
1574
1575 DO ispin = 1, SIZE(ri_data%t_3c_int_mo, 1)
1576 CALL dbt_destroy(ri_data%t_3c_int_mo(ispin, 1, 1))
1577 CALL dbt_destroy(ri_data%t_3c_ctr_RI(ispin, 1, 1))
1578 CALL dbt_destroy(ri_data%t_3c_ctr_KS(ispin, 1, 1))
1579 CALL dbt_destroy(ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
1580 END DO
1581 DO ispin = 1, 2
1582 CALL dbt_destroy(ri_data%t_2c_int(ispin, 1))
1583 END DO
1584 DEALLOCATE (ri_data%t_2c_int)
1585 DEALLOCATE (ri_data%t_3c_int_mo)
1586 DEALLOCATE (ri_data%t_3c_ctr_RI)
1587 DEALLOCATE (ri_data%t_3c_ctr_KS)
1588 DEALLOCATE (ri_data%t_3c_ctr_KS_copy)
1589 END IF
1590
1591 DO j = 1, SIZE(ri_data%t_2c_inv, 2)
1592 DO i = 1, SIZE(ri_data%t_2c_inv, 1)
1593 CALL dbt_destroy(ri_data%t_2c_inv(i, j))
1594 END DO
1595 END DO
1596 DEALLOCATE (ri_data%t_2c_inv)
1597
1598 DO j = 1, SIZE(ri_data%t_2c_pot, 2)
1599 DO i = 1, SIZE(ri_data%t_2c_pot, 1)
1600 CALL dbt_destroy(ri_data%t_2c_pot(i, j))
1601 END DO
1602 END DO
1603 DEALLOCATE (ri_data%t_2c_pot)
1604
1605 IF (ALLOCATED(ri_data%kp_mat_2c_pot)) THEN
1606 DO j = 1, SIZE(ri_data%kp_mat_2c_pot, 2)
1607 DO i = 1, SIZE(ri_data%kp_mat_2c_pot, 1)
1608 CALL dbcsr_release(ri_data%kp_mat_2c_pot(i, j))
1609 END DO
1610 END DO
1611 DEALLOCATE (ri_data%kp_mat_2c_pot)
1612 END IF
1613
1614 IF (ALLOCATED(ri_data%kp_t_3c_int)) THEN
1615 DO i = 1, SIZE(ri_data%kp_t_3c_int)
1616 CALL dbt_destroy(ri_data%kp_t_3c_int(i))
1617 END DO
1618 DEALLOCATE (ri_data%kp_t_3c_int)
1619 END IF
1620
1621 IF (ALLOCATED(ri_data%rho_ao_t)) THEN
1622 DO j = 1, SIZE(ri_data%rho_ao_t, 2)
1623 DO i = 1, SIZE(ri_data%rho_ao_t, 1)
1624 CALL dbt_destroy(ri_data%rho_ao_t(i, j))
1625 END DO
1626 END DO
1627 DEALLOCATE (ri_data%rho_ao_t)
1628 END IF
1629
1630 IF (ALLOCATED(ri_data%ks_t)) THEN
1631 DO j = 1, SIZE(ri_data%ks_t, 2)
1632 DO i = 1, SIZE(ri_data%ks_t, 1)
1633 CALL dbt_destroy(ri_data%ks_t(i, j))
1634 END DO
1635 END DO
1636 DEALLOCATE (ri_data%ks_t)
1637 END IF
1638
1639 IF (ALLOCATED(ri_data%iatom_to_subgroup)) THEN
1640 DO i = 1, SIZE(ri_data%iatom_to_subgroup)
1641 DEALLOCATE (ri_data%iatom_to_subgroup(i)%array)
1642 END DO
1643 DEALLOCATE (ri_data%iatom_to_subgroup)
1644 END IF
1645
1646 CALL timestop(handle)
1647 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 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 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell
1671 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl_a
1672 REAL(dp), DIMENSION(:, :), POINTER :: sphi_a
1673 TYPE(gto_basis_set_type), POINTER :: orb_basis_a
1674
1675 CALL timeset(routinen, handle)
1676
1677 ! BASIS parameter
1678 nkind = SIZE(qs_kind_set, 1)
1679 !
1680 ALLOCATE (basis_parameter(nkind))
1681 max_set = 0
1682 DO ikind = 1, nkind
1683 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 basis_type=basis_type)
1689 IF (basis_info%max_set < max_set) cpabort("UNEXPECTED MAX_SET")
1690 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 kind_radius=basis_parameter(ikind)%kind_radius)
1706 END DO
1707 DO ikind = 1, nkind
1708 ALLOCATE (basis_parameter(ikind)%nsgfl(0:basis_info%max_am, max_set))
1709 basis_parameter(ikind)%nsgfl = 0
1710 nset = basis_parameter(ikind)%nset
1711 nshell => basis_parameter(ikind)%nshell
1712 DO iset = 1, nset
1713 DO i = 0, basis_info%max_am
1714 nl_count = 0
1715 DO j = 1, nshell(iset)
1716 IF (basis_parameter(ikind)%nl(j, iset) == i) nl_count = nl_count + 1
1717 END DO
1718 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 DO ikind = 1, nkind
1726 max_coeff = 0
1727 max_am_kind = 0
1728 max_pgf_kind = 0
1729 npgfa => basis_parameter(ikind)%npgf
1730 nseta = basis_parameter(ikind)%nset
1731 nl_a => basis_parameter(ikind)%nsgfl
1732 la_max => basis_parameter(ikind)%lmax
1733 la_min => basis_parameter(ikind)%lmin
1734 DO iset = 1, nseta
1735 max_pgf_kind = max(max_pgf_kind, npgfa(iset))
1736 max_pgf = max(max_pgf, npgfa(iset))
1737 DO la = la_min(iset), la_max(iset)
1738 max_nsgfl = max(max_nsgfl, nl_a(la, iset))
1739 max_coeff = max(max_coeff, nso(la)*nl_a(la, iset)*nco(la))
1740 max_am_kind = max(max_am_kind, la)
1741 END DO
1742 END DO
1743 ALLOCATE (basis_parameter(ikind)%sphi_ext(max_coeff, 0:max_am_kind, max_pgf_kind, nseta))
1744 basis_parameter(ikind)%sphi_ext = 0.0_dp
1745 END DO
1746
1747 DO ikind = 1, nkind
1748 sphi_a => basis_parameter(ikind)%sphi
1749 nseta = basis_parameter(ikind)%nset
1750 la_max => basis_parameter(ikind)%lmax
1751 la_min => basis_parameter(ikind)%lmin
1752 npgfa => basis_parameter(ikind)%npgf
1753 first_sgfa => basis_parameter(ikind)%first_sgf
1754 nl_a => basis_parameter(ikind)%nsgfl
1755 DO iset = 1, nseta
1756 sgfa = first_sgfa(1, iset)
1757 DO ipgf = 1, npgfa(iset)
1758 offset_a1 = (ipgf - 1)*ncoset(la_max(iset))
1759 s_offset_nl_a = 0
1760 DO la = la_min(iset), la_max(iset)
1761 offset_a = offset_a1 + ncoset(la - 1)
1762 co_counter = 0
1763 co_counter = co_counter + 1
1764 so_counter = 0
1765 DO k = sgfa + s_offset_nl_a, sgfa + s_offset_nl_a + nso(la)*nl_a(la, iset) - 1
1766 DO i = offset_a + 1, offset_a + nco(la)
1767 so_counter = so_counter + 1
1768 basis_parameter(ikind)%sphi_ext(so_counter, la, ipgf, iset) = sphi_a(i, k)
1769 END DO
1770 END DO
1771 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 CALL timestop(handle)
1778
1779 END SUBROUTINE hfx_create_basis_types
1780
1781! **************************************************************************************************
1782!> \brief ...
1783!> \param basis_parameter ...
1784! **************************************************************************************************
1785 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 CALL timeset(routinen, handle)
1793
1794 !! BASIS parameter
1795 DO i = 1, SIZE(basis_parameter)
1796 DEALLOCATE (basis_parameter(i)%nsgfl)
1797 DEALLOCATE (basis_parameter(i)%sphi_ext)
1798 END DO
1799 DEALLOCATE (basis_parameter)
1800 CALL timestop(handle)
1801
1802 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 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 (PRESENT(storage_id) .EQV. PRESENT(irep))
1836 cpassert(check)
1837
1838 ! Memory Storage
1839 CALL section_vals_val_get(hf_sub_section, "MAX_MEMORY", i_val=int_val)
1840 memory_parameter%max_memory = int_val
1841 memory_parameter%max_compression_counter = int_val*1024_int_8*128_int_8
1842 CALL section_vals_val_get(hf_sub_section, "EPS_STORAGE", r_val=real_val)
1843 memory_parameter%eps_storage_scaling = real_val
1844 IF (int_val == 0) THEN
1845 memory_parameter%do_all_on_the_fly = .true.
1846 ELSE
1847 memory_parameter%do_all_on_the_fly = .false.
1848 END IF
1849 memory_parameter%cache_size = cache_size
1850 memory_parameter%bits_max_val = bits_max_val
1851 memory_parameter%actual_memory_usage = 1
1852 IF (.NOT. skip_in_core_forces) THEN
1853 CALL section_vals_val_get(hf_sub_section, "TREAT_FORCES_IN_CORE", l_val=logic_val)
1854 memory_parameter%treat_forces_in_core = logic_val
1855 END IF
1856
1857 ! ** IF MAX_MEM == 0 overwrite this flag to false
1858 IF (memory_parameter%do_all_on_the_fly) memory_parameter%treat_forces_in_core = .false.
1859
1860 ! Disk Storage
1861 IF (.NOT. skip_disk) THEN
1862 memory_parameter%actual_memory_usage_disk = 1
1863 CALL section_vals_val_get(hf_sub_section, "MAX_DISK_SPACE", i_val=int_val)
1864 memory_parameter%max_compression_counter_disk = int_val*1024_int_8*128_int_8
1865 IF (int_val == 0) THEN
1866 memory_parameter%do_disk_storage = .false.
1867 ELSE
1868 memory_parameter%do_disk_storage = .true.
1869 END IF
1870 CALL section_vals_val_get(hf_sub_section, "STORAGE_LOCATION", c_val=char_val)
1871 CALL compress(char_val, .true.)
1872 !! Add ending / if necessary
1873
1874 IF (scan(char_val, "/", .true.) /= len_trim(char_val)) THEN
1875 WRITE (filename, '(A,A)') trim(char_val), "/"
1876 CALL compress(filename)
1877 ELSE
1878 filename = trim(char_val)
1879 END IF
1880 CALL compress(filename, .true.)
1881
1882 !! quickly check if we can write on storage_location
1883 CALL m_getcwd(orig_wd)
1884 CALL m_chdir(trim(filename), stat)
1885 IF (stat /= 0) THEN
1886 WRITE (error_msg, '(A,A,A)') "Request for disk storage failed due to unknown error while writing to ", &
1887 trim(filename), ". Please check STORAGE_LOCATION"
1888 cpabort(error_msg)
1889 END IF
1890 CALL m_chdir(orig_wd, stat)
1891
1892 memory_parameter%storage_location = filename
1893 CALL compress(memory_parameter%storage_location, .true.)
1894 ELSE
1895 memory_parameter%do_disk_storage = .false.
1896 END IF
1897 IF (PRESENT(storage_id)) THEN
1898 storage_id = (irep - 1)*para_env%num_pe*n_threads + para_env%mepos*n_threads + i_thread - 1
1899 END IF
1900 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 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 n_rep_hf = x_data(1, 1)%n_rep_hf
1919 n_threads = SIZE(x_data, 2)
1920
1921 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 init_t_c_g0_lmax = -1
1924 CALL free_c0()
1925 END IF
1926 DO i_thread = 1, n_threads
1927 DO irep = 1, n_rep_hf
1928 actual_x_data => x_data(irep, i_thread)
1929 DEALLOCATE (actual_x_data%neighbor_cells)
1930 DEALLOCATE (actual_x_data%distribution_energy)
1931 DEALLOCATE (actual_x_data%distribution_forces)
1932
1933 IF (actual_x_data%load_balance_parameter%blocks_initialized) THEN
1934 DEALLOCATE (actual_x_data%blocks)
1935 IF (i_thread == 1) THEN
1936 DEALLOCATE (actual_x_data%pmax_block)
1937 END IF
1938 END IF
1939
1940 IF (i_thread == 1) THEN
1941 DEALLOCATE (actual_x_data%atomic_pair_list)
1942 DEALLOCATE (actual_x_data%atomic_pair_list_forces)
1943 END IF
1944
1945 IF (actual_x_data%screening_parameter%do_initial_p_screening .OR. &
1946 actual_x_data%screening_parameter%do_p_screening_forces) THEN
1947 IF (i_thread == 1) THEN
1948 DEALLOCATE (actual_x_data%pmax_atom)
1949 DO i = 1, SIZE(actual_x_data%initial_p)
1950 DEALLOCATE (actual_x_data%initial_p(i)%p_kind)
1951 END DO
1952 DEALLOCATE (actual_x_data%initial_p)
1953
1954 DEALLOCATE (actual_x_data%pmax_atom_forces)
1955 DO i = 1, SIZE(actual_x_data%initial_p_forces)
1956 DEALLOCATE (actual_x_data%initial_p_forces(i)%p_kind)
1957 END DO
1958 DEALLOCATE (actual_x_data%initial_p_forces)
1959 END IF
1960 DEALLOCATE (actual_x_data%map_atom_to_kind_atom)
1961 END IF
1962 IF (i_thread == 1) THEN
1963 DEALLOCATE (actual_x_data%is_assoc_atomic_block)
1964 DEALLOCATE (actual_x_data%atomic_block_offset)
1965 DEALLOCATE (actual_x_data%set_offset)
1966 DEALLOCATE (actual_x_data%block_offset)
1967 END IF
1968
1969 !! BASIS parameter
1970 CALL hfx_release_basis_types(actual_x_data%basis_parameter)
1971
1972 !MK Release libint and libderiv data structure
1973 CALL cp_libint_cleanup_eri(actual_x_data%lib)
1974 CALL cp_libint_cleanup_eri1(actual_x_data%lib_deriv)
1975 CALL cp_libint_static_cleanup()
1976
1977 !! Deallocate containers
1978 CALL dealloc_containers(actual_x_data%store_ints, actual_x_data%memory_parameter%actual_memory_usage)
1979 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 .false.)
1985 IF (actual_x_data%memory_parameter%do_disk_storage) THEN
1986 CALL close_file(unit_number=actual_x_data%store_ints%maxval_container_disk%unit, file_status="DELETE")
1987 END IF
1988 DEALLOCATE (actual_x_data%store_ints%maxval_container_disk%first)
1989 DEALLOCATE (actual_x_data%store_ints%maxval_container_disk)
1990
1991 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 .false.)
1995 IF (actual_x_data%memory_parameter%do_disk_storage) THEN
1996 CALL close_file(unit_number=actual_x_data%store_ints%integral_containers_disk(i)%unit, file_status="DELETE")
1997 END IF
1998 DEALLOCATE (actual_x_data%store_ints%integral_containers_disk(i)%first)
1999 END DO
2000 DEALLOCATE (actual_x_data%store_ints%integral_containers_disk)
2001
2002 ! ** screening functions
2003 IF (actual_x_data%screen_funct_is_initialized) THEN
2004 DEALLOCATE (actual_x_data%screen_funct_coeffs_set)
2005 DEALLOCATE (actual_x_data%screen_funct_coeffs_kind)
2006 DEALLOCATE (actual_x_data%pair_dist_radii_pgf)
2007 DEALLOCATE (actual_x_data%screen_funct_coeffs_pgf)
2008 actual_x_data%screen_funct_is_initialized = .false.
2009 END IF
2010
2011 ! ** maps
2012 IF (ASSOCIATED(actual_x_data%map_atoms_to_cpus)) THEN
2013 DO i = 1, SIZE(actual_x_data%map_atoms_to_cpus)
2014 DEALLOCATE (actual_x_data%map_atoms_to_cpus(i)%iatom_list)
2015 DEALLOCATE (actual_x_data%map_atoms_to_cpus(i)%jatom_list)
2016 END DO
2017 DEALLOCATE (actual_x_data%map_atoms_to_cpus)
2018 END IF
2019
2020 IF (actual_x_data%do_hfx_ri) THEN
2021 CALL hfx_ri_release(actual_x_data%ri_data)
2022 IF (ASSOCIATED(actual_x_data%ri_data%ri_section)) THEN
2023 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 "PRINT%RI_INFO")
2026 END IF
2027 IF (ASSOCIATED(actual_x_data%ri_data%hfx_section)) THEN
2028 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 "HF_INFO")
2031 END IF
2032 DEALLOCATE (actual_x_data%ri_data)
2033 END IF
2034 END DO
2035
2036 END DO
2037
2038 DEALLOCATE (x_data)
2039 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 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 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 REAL(dp), DIMENSION(:, :), POINTER :: zeta, zetb
2070 TYPE(hfx_cell_type), ALLOCATABLE, DIMENSION(:) :: tmp_neighbor_cells
2071
2072 total_number_of_cells = 0
2073
2074 nkp = 1
2075 IF (PRESENT(nkp_grid)) nkp = nkp_grid
2076 do_kpoints = any(nkp > 1)
2077
2078 ! ** Check some settings
2079 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 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 "Use a truncated potential to avoid possible problems.")
2088 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 REAL(nkp(3), dp)*plane_distance(0, 0, 1, cell))
2093 l_min = 0.5_dp*l_min
2094 IF (x_data%potential_parameter%cutoff_radius >= l_min) THEN
2095 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 "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 "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 SELECT CASE (x_data%potential_parameter%potential_type)
2115 CASE (do_potential_truncated, do_potential_mix_cl_trunc, do_potential_short)
2116 r_max = 0.0_dp
2117 DO ikind = 1, SIZE(x_data%basis_parameter)
2118 la_max => x_data%basis_parameter(ikind)%lmax
2119 zeta => x_data%basis_parameter(ikind)%zet
2120 nseta = x_data%basis_parameter(ikind)%nset
2121 npgfa => x_data%basis_parameter(ikind)%npgf
2122 DO jkind = 1, SIZE(x_data%basis_parameter)
2123 lb_max => x_data%basis_parameter(jkind)%lmax
2124 zetb => x_data%basis_parameter(jkind)%zet
2125 nsetb = x_data%basis_parameter(jkind)%nset
2126 npgfb => x_data%basis_parameter(jkind)%npgf
2127 DO iset = 1, nseta
2128 DO jset = 1, nsetb
2129 DO ipgf = 1, npgfa(iset)
2130 DO jpgf = 1, npgfb(jset)
2131 zeta1 = zeta(ipgf, iset) + zetb(jpgf, jset)
2132 r1 = 1.0_dp/sqrt(zeta1)*mul_fact(la_max(iset) + lb_max(jset))* &
2133 sqrt(-log(x_data%screening_parameter%eps_schwarz))
2134 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 r_max = 2.0_dp*r_max + x_data%potential_parameter%cutoff_radius
2143 nothing_more_to_add = .false.
2144 max_shell = 0
2145 total_number_of_cells = 0
2146 ub = 1
2147 DEALLOCATE (x_data%neighbor_cells)
2148 ALLOCATE (x_data%neighbor_cells(1))
2149 x_data%neighbor_cells(1)%cell = 0.0_dp
2150 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 DO WHILE (.NOT. nothing_more_to_add)
2188 ! Calculate distances to the eight points P1 to P8
2189 image_cell_found = .false.
2190 ALLOCATE (tmp_neighbor_cells(1:ub))
2191 DO i = 1, ub - 1
2192 tmp_neighbor_cells(i) = x_data%neighbor_cells(i)
2193 END DO
2194 ub_max = (2*max_shell + 1)**3
2195 DEALLOCATE (x_data%neighbor_cells)
2196 ALLOCATE (x_data%neighbor_cells(1:ub_max))
2197 DO i = 1, ub - 1
2198 x_data%neighbor_cells(i) = tmp_neighbor_cells(i)
2199 END DO
2200 DO i = ub, ub_max
2201 x_data%neighbor_cells(i)%cell = 0.0_dp
2202 x_data%neighbor_cells(i)%cell_r = 0.0_dp
2203 END DO
2204
2205 DEALLOCATE (tmp_neighbor_cells)
2206
2207 perd(1:3) = x_data%periodic_parameter%perd(1:3)
2208
2209 DO ishell = -max_shell*perd(1), max_shell*perd(1)
2210 DO jshell = -max_shell*perd(2), max_shell*perd(2)
2211 DO kshell = -max_shell*perd(3), max_shell*perd(3)
2212 IF (max(abs(ishell), abs(jshell), abs(kshell)) /= max_shell) cycle
2213 idx = 0
2214 DO j = 0, 1
2215 x = -1.0_dp/2.0_dp + j*1.0_dp
2216 DO k = 0, 1
2217 y = -1.0_dp/2.0_dp + k*1.0_dp
2218 DO l = 0, 1
2219 z = -1.0_dp/2.0_dp + l*1.0_dp
2220 idx = idx + 1
2221 p(1, idx) = x + ishell
2222 p(2, idx) = y + jshell
2223 p(3, idx) = z + kshell
2224 CALL scaled_to_real(r, p(:, idx), cell)
2225 distance(idx) = sqrt(sum(r**2))
2226 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 idx = idx + 1
2234 plane_vector(:, 1) = p(:, 3) - p(:, 1)
2235 plane_vector(:, 2) = p(:, 2) - p(:, 1)
2236 cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2237 cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2238 cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2239 normal(:, 1) = cross_product/sqrt(sum(cross_product**2))
2240 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 IF (point_is_in_quadrilateral(p(:, 1), p(:, 3), p(:, 4), p(:, 2), point_in_plane)) THEN
2243 distance(idx) = abs(normal(1, 1)*p(1, 1) + normal(2, 1)*p(2, 1) + normal(3, 1)*p(3, 1))
2244 ELSE
2245 distance(idx) = huge(distance(idx))
2246 END IF
2247
2248 ! Face B (1562) 1 is the reference
2249 idx = idx + 1
2250 plane_vector(:, 1) = p(:, 2) - p(:, 1)
2251 plane_vector(:, 2) = p(:, 5) - p(:, 1)
2252 cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2253 cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2254 cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2255 normal(:, 1) = cross_product/sqrt(sum(cross_product**2))
2256 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 IF (point_is_in_quadrilateral(p(:, 1), p(:, 5), p(:, 6), p(:, 2), point_in_plane)) THEN
2259 distance(idx) = abs(normal(1, 1)*p(1, 1) + normal(2, 1)*p(2, 1) + normal(3, 1)*p(3, 1))
2260 ELSE
2261 distance(idx) = huge(distance(idx))
2262 END IF
2263
2264 ! Face C (5786) 5 is the reference
2265 idx = idx + 1
2266 plane_vector(:, 1) = p(:, 7) - p(:, 5)
2267 plane_vector(:, 2) = p(:, 6) - p(:, 5)
2268 cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2269 cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2270 cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2271 normal(:, 1) = cross_product/sqrt(sum(cross_product**2))
2272 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 IF (point_is_in_quadrilateral(p(:, 5), p(:, 7), p(:, 8), p(:, 6), point_in_plane)) THEN
2275 distance(idx) = abs(normal(1, 1)*p(1, 5) + normal(2, 1)*p(2, 5) + normal(3, 1)*p(3, 5))
2276 ELSE
2277 distance(idx) = huge(distance(idx))
2278 END IF
2279
2280 ! Face D (3784) 3 is the reference
2281 idx = idx + 1
2282 plane_vector(:, 1) = p(:, 7) - p(:, 3)
2283 plane_vector(:, 2) = p(:, 4) - p(:, 3)
2284 cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2285 cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2286 cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2287 normal(:, 1) = cross_product/sqrt(sum(cross_product**2))
2288 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 IF (point_is_in_quadrilateral(p(:, 3), p(:, 7), p(:, 8), p(:, 4), point_in_plane)) THEN
2291 distance(idx) = abs(normal(1, 1)*p(1, 3) + normal(2, 1)*p(2, 3) + normal(3, 1)*p(3, 3))
2292 ELSE
2293 distance(idx) = huge(distance(idx))
2294 END IF
2295
2296 ! Face E (2684) 2 is the reference
2297 idx = idx + 1
2298 plane_vector(:, 1) = p(:, 6) - p(:, 2)
2299 plane_vector(:, 2) = p(:, 4) - p(:, 2)
2300 cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2301 cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2302 cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2303 normal(:, 1) = cross_product/sqrt(sum(cross_product**2))
2304 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 IF (point_is_in_quadrilateral(p(:, 2), p(:, 6), p(:, 8), p(:, 4), point_in_plane)) THEN
2307 distance(idx) = abs(normal(1, 1)*p(1, 2) + normal(2, 1)*p(2, 2) + normal(3, 1)*p(3, 2))
2308 ELSE
2309 distance(idx) = huge(distance(idx))
2310 END IF
2311
2312 ! Face F (1573) 1 is the reference
2313 idx = idx + 1
2314 plane_vector(:, 1) = p(:, 5) - p(:, 1)
2315 plane_vector(:, 2) = p(:, 3) - p(:, 1)
2316 cross_product(1) = plane_vector(2, 1)*plane_vector(3, 2) - plane_vector(3, 1)*plane_vector(2, 2)
2317 cross_product(2) = plane_vector(3, 1)*plane_vector(1, 2) - plane_vector(1, 1)*plane_vector(3, 2)
2318 cross_product(3) = plane_vector(1, 1)*plane_vector(2, 2) - plane_vector(2, 1)*plane_vector(1, 2)
2319 normal(:, 1) = cross_product/sqrt(sum(cross_product**2))
2320 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 IF (point_is_in_quadrilateral(p(:, 1), p(:, 5), p(:, 7), p(:, 3), point_in_plane)) THEN
2323 distance(idx) = abs(normal(1, 1)*p(1, 1) + normal(2, 1)*p(2, 1) + normal(3, 1)*p(3, 1))
2324 ELSE
2325 distance(idx) = huge(distance(idx))
2326 END IF
2327
2328 dist_min = minval(distance)
2329 IF (max_shell == 0) THEN
2330 image_cell_found = .true.
2331 END IF
2332 IF (dist_min < r_max) THEN
2333 total_number_of_cells = total_number_of_cells + 1
2334 x_data%neighbor_cells(ub)%cell = real((/ishell, jshell, kshell/), dp)
2335 ub = ub + 1
2336 image_cell_found = .true.
2337 END IF
2338
2339 END DO
2340 END DO
2341 END DO
2342 IF (image_cell_found) THEN
2343 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 ALLOCATE (tmp_neighbor_cells(total_number_of_cells))
2350 DO i = 1, ub - 1
2351 tmp_neighbor_cells(i) = x_data%neighbor_cells(i)
2352 END DO
2353 DEALLOCATE (x_data%neighbor_cells)
2354 ! If we only need the supercell, total_number_of_cells is still 0, repair
2355 IF (total_number_of_cells == 0) THEN
2356 total_number_of_cells = 1
2357 ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2358 DO i = 1, total_number_of_cells
2359 x_data%neighbor_cells(i)%cell = 0.0_dp
2360 x_data%neighbor_cells(i)%cell_r = 0.0_dp
2361 END DO
2362 ELSE
2363 ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2364 DO i = 1, total_number_of_cells
2365 x_data%neighbor_cells(i) = tmp_neighbor_cells(i)
2366 END DO
2367 END IF
2368 DEALLOCATE (tmp_neighbor_cells)
2369
2370 IF (x_data%periodic_parameter%number_of_shells == do_hfx_auto_shells) THEN
2371 ! Do nothing
2372 ELSE
2373 total_number_of_cells = 0
2374 DO i = 0, x_data%periodic_parameter%number_of_shells
2375 total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2376 END DO
2377 IF (total_number_of_cells < SIZE(x_data%neighbor_cells)) THEN
2378 IF (i_thread == 1) THEN
2379 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 "Please carefully check if you get converged results."
2384 cpwarn(error_msg)
2385 END IF
2386 END IF
2387 total_number_of_cells = 0
2388 DO i = 0, x_data%periodic_parameter%number_of_shells
2389 total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2390 END DO
2391 DEALLOCATE (x_data%neighbor_cells)
2392
2393 ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2394 m = 0
2395 i = 1
2396 DO WHILE (sum(m**2) <= x_data%periodic_parameter%number_of_shells)
2397 x_data%neighbor_cells(i)%cell = real(m, dp)
2398 CALL next_image_cell_perd(m, x_data%periodic_parameter%perd)
2399 i = i + 1
2400 END DO
2401 END IF
2402 CASE DEFAULT
2403 total_number_of_cells = 0
2404 IF (pbc_shells == -1) pbc_shells = 0
2405 DO i = 0, pbc_shells
2406 total_number_of_cells = total_number_of_cells + count_cells_perd(i, x_data%periodic_parameter%perd)
2407 END DO
2408 DEALLOCATE (x_data%neighbor_cells)
2409
2410 ALLOCATE (x_data%neighbor_cells(total_number_of_cells))
2411
2412 m = 0
2413 i = 1
2414 DO WHILE (sum(m**2) <= pbc_shells)
2415 x_data%neighbor_cells(i)%cell = real(m, dp)
2416 CALL next_image_cell_perd(m, x_data%periodic_parameter%perd)
2417 i = i + 1
2418 END DO
2419 END SELECT
2420
2421 ! ** Transform into real coord
2422 DO i = 1, SIZE(x_data%neighbor_cells)
2423 r = 0.0_dp
2424 x_data%neighbor_cells(i)%cell_r(:) = 0.0_dp
2425 s = x_data%neighbor_cells(i)%cell(:)
2426 CALL scaled_to_real(x_data%neighbor_cells(i)%cell_r, s, cell)
2427 END DO
2428 x_data%periodic_parameter%number_of_shells = pbc_shells
2429
2430 r_max_stress = 0.0_dp
2431 DO i = 1, SIZE(x_data%neighbor_cells)
2432 r_max_stress = max(r_max_stress, maxval(abs(x_data%neighbor_cells(i)%cell_r(:))))
2433 END DO
2434 r_max_stress = r_max_stress + abs(maxval(cell%hmat(:, :)))
2435 x_data%periodic_parameter%R_max_stress = r_max_stress
2436
2437 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 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 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 v0 = d - a
2473 v1 = c - a
2474 v2 = p - a
2475
2476 ! ** Compute dot products
2477 dot00 = dot_product(v0, v0)
2478 dot01 = dot_product(v0, v1)
2479 dot02 = dot_product(v0, v2)
2480 dot11 = dot_product(v1, v1)
2481 dot12 = dot_product(v1, v2)
2482
2483 ! ** Compute barycentric coordinates
2484 invdenom = 1/(dot00*dot11 - dot01*dot01)
2485 u = (dot11*dot02 - dot01*dot12)*invdenom
2486 v = (dot00*dot12 - dot01*dot02)*invdenom
2487 ! ** Check if point is in triangle
2488 IF ((u >= 0 - fuzzy) .AND. (v >= 0 - fuzzy) .AND. (u + v <= 1 + fuzzy)) THEN
2489 point_is_in_quadrilateral = .true.
2490 RETURN
2491 END IF
2492 v0 = c - a
2493 v1 = b - a
2494 v2 = p - a
2495
2496 ! ** Compute dot products
2497 dot00 = dot_product(v0, v0)
2498 dot01 = dot_product(v0, v1)
2499 dot02 = dot_product(v0, v2)
2500 dot11 = dot_product(v1, v1)
2501 dot12 = dot_product(v1, v2)
2502
2503 ! ** Compute barycentric coordinates
2504 invdenom = 1/(dot00*dot11 - dot01*dot01)
2505 u = (dot11*dot02 - dot01*dot12)*invdenom
2506 v = (dot00*dot12 - dot01*dot02)*invdenom
2507
2508 ! ** Check if point is in triangle
2509 IF ((u >= 0 - fuzzy) .AND. (v >= 0 - fuzzy) .AND. (u + v <= 1 + fuzzy)) THEN
2510 point_is_in_quadrilateral = .true.
2511 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 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 current => container%first
2536 DO WHILE (ASSOCIATED(current))
2537 next => current%next
2538 DEALLOCATE (current)
2539 current => next
2540 END DO
2541
2542 !! Allocate first list entry, init members
2543 ALLOCATE (container%first)
2544 container%first%prev => null()
2545 container%first%next => null()
2546 container%current => container%first
2547 container%current%data = 0
2548 container%element_counter = 1
2549 memory_usage = 1
2550
2551 IF (do_disk_storage) THEN
2552 !! close the file, if this is no the first time
2553 IF (container%unit /= -1) THEN
2554 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 unit_number=container%unit)
2558 END IF
2559
2560 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 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 DEALLOCATE (x_data%distribution_energy)
2576
2577 ALLOCATE (x_data%distribution_energy(SIZE(ptr_to_distr)))
2578 x_data%distribution_energy = ptr_to_distr
2579
2580 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 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 DEALLOCATE (x_data%distribution_forces)
2596
2597 ALLOCATE (x_data%distribution_forces(SIZE(ptr_to_distr)))
2598 x_data%distribution_forces = ptr_to_distr
2599
2600 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 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 max_memory = memory_parameter%max_memory
2620 max_memory = max_memory - subtr_size_mb
2621 IF (max_memory <= 0) THEN
2622 memory_parameter%do_all_on_the_fly = .true.
2623 memory_parameter%max_compression_counter = 0
2624 ELSE
2625 memory_parameter%do_all_on_the_fly = .false.
2626 memory_parameter%max_compression_counter = max_memory*1024_int_8*128_int_8
2627 END IF
2628 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 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 NULLIFY (logger)
2646 logger => cp_get_default_logger()
2647
2648 iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2649 extension=".scfLog")
2650
2651 IF (iw > 0) THEN
2652 WRITE (unit=iw, fmt="((T3,A,T73,ES8.1))") &
2653 "HFX_INFO| EPS_SCHWARZ: ", x_data%screening_parameter%eps_schwarz
2654 WRITE (unit=iw, fmt="((T3,A,T73,ES8.1))") &
2655 "HFX_INFO| EPS_SCHWARZ_FORCES ", x_data%screening_parameter%eps_schwarz_forces
2656 WRITE (unit=iw, fmt="((T3,A,T73,ES8.1))") &
2657 "HFX_INFO| EPS_STORAGE_SCALING: ", x_data%memory_parameter%eps_storage_scaling
2658 WRITE (unit=iw, fmt="((T3,A,T61,I20))") &
2659 "HFX_INFO| NBINS: ", x_data%load_balance_parameter%nbins
2660 WRITE (unit=iw, fmt="((T3,A,T61,I20))") &
2661 "HFX_INFO| BLOCK_SIZE: ", x_data%load_balance_parameter%block_size
2662 IF (x_data%periodic_parameter%do_periodic) THEN
2663 IF (x_data%periodic_parameter%mode == -1) THEN
2664 WRITE (unit=iw, fmt="((T3,A,T77,A))") &
2665 "HFX_INFO| NUMBER_OF_SHELLS: ", "AUTO"
2666 ELSE
2667 WRITE (unit=iw, fmt="((T3,A,T61,I20))") &
2668 "HFX_INFO| NUMBER_OF_SHELLS: ", x_data%periodic_parameter%mode
2669 END IF
2670 WRITE (unit=iw, fmt="((T3,A,T61,I20))") &
2671 "HFX_INFO| Number of periodic shells considered: ", x_data%periodic_parameter%number_of_shells
2672 WRITE (unit=iw, fmt="((T3,A,T61,I20),/)") &
2673 "HFX_INFO| Number of periodic cells considered: ", SIZE(x_data%neighbor_cells)
2674 ELSE
2675 WRITE (unit=iw, fmt="((T3,A,T77,A))") &
2676 "HFX_INFO| Number of periodic shells considered: ", "NONE"
2677 WRITE (unit=iw, fmt="((T3,A,T77,A),/)") &
2678 "HFX_INFO| Number of periodic cells considered: ", "NONE"
2679 END IF
2680 END IF
2681 END SUBROUTINE hfx_print_std_info
2682
2683! **************************************************************************************************
2684!> \brief ...
2685!> \param ri_data ...
2686!> \param hfx_section ...
2687! **************************************************************************************************
2688 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 NULLIFY (logger, ri_section)
2698 logger => cp_get_default_logger()
2699
2700 ri_section => ri_data%ri_section
2701
2702 iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2703 extension=".scfLog")
2704
2705 IF (iw > 0) THEN
2706
2707 associate(ri_metric => ri_data%ri_metric, hfx_pot => ri_data%hfx_pot)
2708 SELECT CASE (ri_metric%potential_type)
2709 CASE (do_potential_coulomb)
2710 WRITE (unit=iw, fmt="(/T3,A,T74,A)") &
2711 "HFX_RI_INFO| RI metric: ", "COULOMB"
2712 CASE (do_potential_short)
2713 WRITE (unit=iw, fmt="(T3,A,T71,A)") &
2714 "HFX_RI_INFO| RI metric: ", "SHORTRANGE"
2715 WRITE (iw, '(T3,A,T61,F20.10)') &
2716 "HFX_RI_INFO| Omega: ", ri_metric%omega
2717 rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
2718 WRITE (iw, '(T3,A,T61,F20.10)') &
2719 "HFX_RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
2720 CASE (do_potential_long)
2721 WRITE (unit=iw, fmt="(T3,A,T72,A)") &
2722 "HFX_RI_INFO| RI metric: ", "LONGRANGE"
2723 WRITE (iw, '(T3,A,T61,F20.10)') &
2724 "HFX_RI_INFO| Omega: ", ri_metric%omega
2725 CASE (do_potential_id)
2726 WRITE (unit=iw, fmt="(T3,A,T73,A)") &
2727 "HFX_RI_INFO| RI metric: ", "OVERLAP"
2728 CASE (do_potential_truncated)
2729 WRITE (unit=iw, fmt="(T3,A,T72,A)") &
2730 "HFX_RI_INFO| RI metric: ", "TRUNCATED COULOMB"
2731 rc_ang = cp_unit_from_cp2k(ri_metric%cutoff_radius, "angstrom")
2732 WRITE (iw, '(T3,A,T61,F20.10)') &
2733 "HFX_RI_INFO| Cutoff Radius [angstrom]: ", rc_ang
2734 END SELECT
2735
2736 END associate
2737 SELECT CASE (ri_data%flavor)
2738 CASE (ri_mo)
2739 WRITE (unit=iw, fmt="(T3, A, T79, A)") &
2740 "HFX_RI_INFO| RI flavor: ", "MO"
2741 CASE (ri_pmat)
2742 WRITE (unit=iw, fmt="(T3, A, T78, A)") &
2743 "HFX_RI_INFO| RI flavor: ", "RHO"
2744 END SELECT
2745 SELECT CASE (ri_data%t2c_method)
2746 CASE (hfx_ri_do_2c_iter)
2747 WRITE (unit=iw, fmt="(T3, A, T69, A)") &
2748 "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 "HFX_RI_INFO| Matrix SQRT/INV", "Dense / diag"
2752 END SELECT
2753 WRITE (unit=iw, fmt="(T3, A, T73, ES8.1)") &
2754 "HFX_RI_INFO| EPS_FILTER", ri_data%filter_eps
2755 WRITE (unit=iw, fmt="(T3, A, T73, ES8.1)") &
2756 "HFX_RI_INFO| EPS_FILTER 2-center", ri_data%filter_eps_2c
2757 WRITE (unit=iw, fmt="(T3, A, T73, ES8.1)") &
2758 "HFX_RI_INFO| EPS_FILTER storage", ri_data%filter_eps_storage
2759 WRITE (unit=iw, fmt="(T3, A, T73, ES8.1)") &
2760 "HFX_RI_INFO| EPS_FILTER MO", ri_data%filter_eps_mo
2761 WRITE (unit=iw, fmt="(T3, A, T73, ES8.1)") &
2762 "HFX_RI_INFO| EPS_PGF_ORB", ri_data%eps_pgf_orb
2763 WRITE (unit=iw, fmt="((T3, A, T73, ES8.1))") &
2764 "HFX_RI_INFO| EPS_SCHWARZ: ", ri_data%eps_schwarz
2765 WRITE (unit=iw, fmt="((T3, A, T73, ES8.1))") &
2766 "HFX_RI_INFO| EPS_SCHWARZ_FORCES: ", ri_data%eps_schwarz_forces
2767 WRITE (unit=iw, fmt="(T3, A, T78, I3)") &
2768 "HFX_RI_INFO| Minimum block size", ri_data%min_bsize
2769 WRITE (unit=iw, fmt="(T3, A, T78, I3)") &
2770 "HFX_RI_INFO| MO block size", ri_data%max_bsize_MO
2771 WRITE (unit=iw, fmt="(T3, A, T79, I2)") &
2772 "HFX_RI_INFO| Memory reduction factor", ri_data%n_mem_input
2773 END IF
2774
2775 END SUBROUTINE
2776
2777! **************************************************************************************************
2778!> \brief ...
2779!> \param x_data ...
2780!> \param hfx_section ...
2781!> \param i_rep ...
2782! **************************************************************************************************
2783 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 NULLIFY (logger)
2793 logger => cp_get_default_logger()
2794
2795 iw = cp_print_key_unit_nr(logger, hfx_section, "HF_INFO", &
2796 extension=".scfLog")
2797
2798 IF (iw > 0) THEN
2799 WRITE (unit=iw, fmt="(/,(T3,A,T61,I20))") &
2800 "HFX_INFO| Replica ID: ", i_rep
2801
2802 WRITE (iw, '(T3,A,T61,F20.10)') &
2803 "HFX_INFO| FRACTION: ", x_data%general_parameter%fraction
2804 SELECT CASE (x_data%potential_parameter%potential_type)
2805 CASE (do_potential_coulomb)
2806 WRITE (unit=iw, fmt="((T3,A,T74,A))") &
2807 "HFX_INFO| Interaction Potential: ", "COULOMB"
2808 CASE (do_potential_short)
2809 WRITE (unit=iw, fmt="((T3,A,T71,A))") &
2810 "HFX_INFO| Interaction Potential: ", "SHORTRANGE"
2811 WRITE (iw, '(T3,A,T61,F20.10)') &
2812 "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2813 rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2814 WRITE (iw, '(T3,A,T61,F20.10)') &
2815 "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2816 CASE (do_potential_long)
2817 WRITE (unit=iw, fmt="((T3,A,T72,A))") &
2818 "HFX_INFO| Interaction Potential: ", "LONGRANGE"
2819 WRITE (iw, '(T3,A,T61,F20.10)') &
2820 "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2821 CASE (do_potential_mix_cl)
2822 WRITE (unit=iw, fmt="((T3,A,T75,A))") &
2823 "HFX_INFO| Interaction Potential: ", "MIX_CL"
2824 WRITE (iw, '(T3,A,T61,F20.10)') &
2825 "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2826 WRITE (iw, '(T3,A,T61,F20.10)') &
2827 "HFX_INFO| SCALE_COULOMB: ", x_data%potential_parameter%scale_coulomb
2828 WRITE (iw, '(T3,A,T61,F20.10)') &
2829 "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 "HFX_INFO| Interaction Potential: ", "GAUSSIAN"
2833 WRITE (iw, '(T3,A,T61,F20.10)') &
2834 "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2835 CASE (do_potential_mix_lg)
2836 WRITE (unit=iw, fmt="((T3,A,T75,A))") &
2837 "HFX_INFO| Interaction Potential: ", "MIX_LG"
2838 WRITE (iw, '(T3,A,T61,F20.10)') &
2839 "HFX_INFO| Omega: ", x_data%potential_parameter%omega
2840 WRITE (iw, '(T3,A,T61,F20.10)') &
2841 "HFX_INFO| SCALE_LONGRANGE: ", x_data%potential_parameter%scale_longrange
2842 WRITE (iw, '(T3,A,T61,F20.10)') &
2843 "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 "HFX_INFO| Interaction Potential: ", "IDENTITY"
2847 CASE (do_potential_truncated)
2848 WRITE (unit=iw, fmt="((T3,A,T72,A))") &
2849 "HFX_INFO| Interaction Potential: ", "TRUNCATED"
2850 rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2851 WRITE (iw, '(T3,A,T61,F20.10)') &
2852 "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2853 CASE (do_potential_mix_cl_trunc)
2854 WRITE (unit=iw, fmt="((T3,A,T65,A))") &
2855 "HFX_INFO| Interaction Potential: ", "TRUNCATED MIX_CL"
2856 rc_ang = cp_unit_from_cp2k(x_data%potential_parameter%cutoff_radius, "angstrom")
2857 WRITE (iw, '(T3,A,T61,F20.10)') &
2858 "HFX_INFO| Cutoff Radius [angstrom]: ", rc_ang
2859 END SELECT
2860
2861 END IF
2862 IF (x_data%do_hfx_ri) THEN
2863 CALL hfx_print_ri_info(x_data%ri_data, hfx_section)
2864 ELSE
2865 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 "HF_INFO")
2870 END SUBROUTINE
2871
2872! **************************************************************************************************
2873!> \brief ...
2874!> \param DATA ...
2875!> \param memory_usage ...
2876! **************************************************************************************************
2877 SUBROUTINE dealloc_containers(DATA, memory_usage)
2878 TYPE(hfx_compression_type) :: data
2879 INTEGER :: memory_usage
2880
2881 INTEGER :: bin, i
2882
2883 DO bin = 1, SIZE(data%maxval_container)
2884 CALL hfx_init_container(data%maxval_container(bin), memory_usage, &
2885 .false.)
2886 DEALLOCATE (data%maxval_container(bin)%first)
2887 END DO
2888 DEALLOCATE (data%maxval_container)
2889 DEALLOCATE (data%maxval_cache)
2890
2891 DO bin = 1, SIZE(data%integral_containers, 2)
2892 DO i = 1, 64
2893 CALL hfx_init_container(data%integral_containers(i, bin), memory_usage, &
2894 .false.)
2895 DEALLOCATE (data%integral_containers(i, bin)%first)
2896 END DO
2897 END DO
2898 DEALLOCATE (data%integral_containers)
2899
2900 DEALLOCATE (data%integral_caches)
2901
2902 END SUBROUTINE dealloc_containers
2903
2904! **************************************************************************************************
2905!> \brief ...
2906!> \param DATA ...
2907!> \param bin_size ...
2908! **************************************************************************************************
2909 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 ALLOCATE (data%maxval_cache(bin_size))
2916 DO bin = 1, bin_size
2917 data%maxval_cache(bin)%element_counter = 1
2918 END DO
2919 ALLOCATE (data%maxval_container(bin_size))
2920 DO bin = 1, bin_size
2921 ALLOCATE (data%maxval_container(bin)%first)
2922 data%maxval_container(bin)%first%prev => null()
2923 data%maxval_container(bin)%first%next => null()
2924 data%maxval_container(bin)%current => data%maxval_container(bin)%first
2925 data%maxval_container(bin)%current%data = 0
2926 data%maxval_container(bin)%element_counter = 1
2927 END DO
2928
2929 ALLOCATE (data%integral_containers(64, bin_size))
2930 ALLOCATE (data%integral_caches(64, bin_size))
2931
2932 DO bin = 1, bin_size
2933 DO i = 1, 64
2934 data%integral_caches(i, bin)%element_counter = 1
2935 data%integral_caches(i, bin)%data = 0
2936 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 data%integral_containers(i, bin)%current => data%integral_containers(i, bin)%first
2940 data%integral_containers(i, bin)%current%data = 0
2941 data%integral_containers(i, bin)%element_counter = 1
2942 END DO
2943 END DO
2944
2945 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 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 is_identical = .true.
2969 IF (PRESENT(same_except_frac)) same_except_frac = .false.
2970
2971 CALL section_vals_get(hfx_section1, n_repetition=n_rep_hf1)
2972 CALL section_vals_get(hfx_section2, n_repetition=n_rep_hf2)
2973 is_identical = n_rep_hf1 == n_rep_hf2
2974 IF (.NOT. is_identical) RETURN
2975
2976 DO irep = 1, n_rep_hf1
2977 CALL section_vals_val_get(hfx_section1, "PW_HFX", l_val=lval1, i_rep_section=irep)
2978 CALL section_vals_val_get(hfx_section2, "PW_HFX", l_val=lval2, i_rep_section=irep)
2979 IF (lval1 .NEQV. lval2) is_identical = .false.
2980
2981 CALL section_vals_val_get(hfx_section1, "PW_HFX_BLOCKSIZE", i_val=ival1, i_rep_section=irep)
2982 CALL section_vals_val_get(hfx_section2, "PW_HFX_BLOCKSIZE", i_val=ival2, i_rep_section=irep)
2983 IF (ival1 .NE. ival2) is_identical = .false.
2984
2985 CALL section_vals_val_get(hfx_section1, "TREAT_LSD_IN_CORE", l_val=lval1, i_rep_section=irep)
2986 CALL section_vals_val_get(hfx_section2, "TREAT_LSD_IN_CORE", l_val=lval2, i_rep_section=irep)
2987 IF (lval1 .NEQV. lval2) is_identical = .false.
2988
2989 hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "INTERACTION_POTENTIAL", i_rep_section=irep)
2990 hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "INTERACTION_POTENTIAL", i_rep_section=irep)
2991
2992 CALL section_vals_val_get(hfx_sub_section1, "OMEGA", r_val=rval1, i_rep_section=irep)
2993 CALL section_vals_val_get(hfx_sub_section2, "OMEGA", r_val=rval2, i_rep_section=irep)
2994 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
2995
2996 CALL section_vals_val_get(hfx_sub_section1, "POTENTIAL_TYPE", i_val=ival1, i_rep_section=irep)
2997 CALL section_vals_val_get(hfx_sub_section2, "POTENTIAL_TYPE", i_val=ival2, i_rep_section=irep)
2998 IF (ival1 .NE. ival2) is_identical = .false.
2999 IF (.NOT. is_identical) RETURN
3000
3001 IF (ival1 == do_potential_truncated .OR. ival1 == do_potential_mix_cl_trunc) THEN
3002 CALL section_vals_val_get(hfx_sub_section1, "CUTOFF_RADIUS", r_val=rval1, i_rep_section=irep)
3003 CALL section_vals_val_get(hfx_sub_section2, "CUTOFF_RADIUS", r_val=rval2, i_rep_section=irep)
3004 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3005
3006 CALL section_vals_val_get(hfx_sub_section1, "T_C_G_DATA", c_val=cval1, i_rep_section=irep)
3007 CALL section_vals_val_get(hfx_sub_section2, "T_C_G_DATA", c_val=cval2, i_rep_section=irep)
3008 IF (cval1 .NE. cval2) is_identical = .false.
3009 END IF
3010
3011 CALL section_vals_val_get(hfx_sub_section1, "SCALE_COULOMB", r_val=rval1, i_rep_section=irep)
3012 CALL section_vals_val_get(hfx_sub_section2, "SCALE_COULOMB", r_val=rval2, i_rep_section=irep)
3013 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3014
3015 CALL section_vals_val_get(hfx_sub_section1, "SCALE_GAUSSIAN", r_val=rval1, i_rep_section=irep)
3016 CALL section_vals_val_get(hfx_sub_section2, "SCALE_GAUSSIAN", r_val=rval2, i_rep_section=irep)
3017 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3018
3019 CALL section_vals_val_get(hfx_sub_section1, "SCALE_LONGRANGE", r_val=rval1, i_rep_section=irep)
3020 CALL section_vals_val_get(hfx_sub_section2, "SCALE_LONGRANGE", r_val=rval2, i_rep_section=irep)
3021 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3022
3023 hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "PERIODIC", i_rep_section=irep)
3024 hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "PERIODIC", i_rep_section=irep)
3025
3026 CALL section_vals_val_get(hfx_sub_section1, "NUMBER_OF_SHELLS", i_val=ival1, i_rep_section=irep)
3027 CALL section_vals_val_get(hfx_sub_section2, "NUMBER_OF_SHELLS", i_val=ival2, i_rep_section=irep)
3028 IF (ival1 .NE. ival2) is_identical = .false.
3029
3030 hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "RI", i_rep_section=irep)
3031 hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "RI", i_rep_section=irep)
3032
3033 CALL section_vals_val_get(hfx_sub_section1, "_SECTION_PARAMETERS_", l_val=lval1, i_rep_section=irep)
3034 CALL section_vals_val_get(hfx_sub_section2, "_SECTION_PARAMETERS_", l_val=lval2, i_rep_section=irep)
3035 IF (lval1 .NEQV. lval2) is_identical = .false.
3036
3037 CALL section_vals_val_get(hfx_sub_section1, "CUTOFF_RADIUS", r_val=rval1, i_rep_section=irep)
3038 CALL section_vals_val_get(hfx_sub_section2, "CUTOFF_RADIUS", r_val=rval2, i_rep_section=irep)
3039 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3040
3041 CALL section_vals_val_get(hfx_sub_section1, "EPS_EIGVAL", r_val=rval1, i_rep_section=irep)
3042 CALL section_vals_val_get(hfx_sub_section2, "EPS_EIGVAL", r_val=rval2, i_rep_section=irep)
3043 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3044
3045 CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER", r_val=rval1, i_rep_section=irep)
3046 CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER", r_val=rval2, i_rep_section=irep)
3047 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3048
3049 CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER_2C", r_val=rval1, i_rep_section=irep)
3050 CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER_2C", r_val=rval2, i_rep_section=irep)
3051 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3052
3053 CALL section_vals_val_get(hfx_sub_section1, "EPS_FILTER_MO", r_val=rval1, i_rep_section=irep)
3054 CALL section_vals_val_get(hfx_sub_section2, "EPS_FILTER_MO", r_val=rval2, i_rep_section=irep)
3055 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3056
3057 CALL section_vals_val_get(hfx_sub_section1, "EPS_PGF_ORB", r_val=rval1, i_rep_section=irep)
3058 CALL section_vals_val_get(hfx_sub_section2, "EPS_PGF_ORB", r_val=rval2, i_rep_section=irep)
3059 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3060
3061 CALL section_vals_val_get(hfx_sub_section1, "MAX_BLOCK_SIZE_MO", i_val=ival1, i_rep_section=irep)
3062 CALL section_vals_val_get(hfx_sub_section2, "MAX_BLOCK_SIZE_MO", i_val=ival2, i_rep_section=irep)
3063 IF (ival1 .NE. ival2) is_identical = .false.
3064
3065 CALL section_vals_val_get(hfx_sub_section1, "MIN_BLOCK_SIZE", i_val=ival1, i_rep_section=irep)
3066 CALL section_vals_val_get(hfx_sub_section2, "MIN_BLOCK_SIZE", i_val=ival2, i_rep_section=irep)
3067 IF (ival1 .NE. ival2) is_identical = .false.
3068
3069 CALL section_vals_val_get(hfx_sub_section1, "OMEGA", r_val=rval1, i_rep_section=irep)
3070 CALL section_vals_val_get(hfx_sub_section2, "OMEGA", r_val=rval2, i_rep_section=irep)
3071 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3072
3073 CALL section_vals_val_get(hfx_sub_section1, "RI_FLAVOR", i_val=ival1, i_rep_section=irep)
3074 CALL section_vals_val_get(hfx_sub_section2, "RI_FLAVOR", i_val=ival2, i_rep_section=irep)
3075 IF (ival1 .NE. ival2) is_identical = .false.
3076
3077 CALL section_vals_val_get(hfx_sub_section1, "RI_METRIC", i_val=ival1, i_rep_section=irep)
3078 CALL section_vals_val_get(hfx_sub_section2, "RI_METRIC", i_val=ival2, i_rep_section=irep)
3079 IF (ival1 .NE. ival2) is_identical = .false.
3080
3081 hfx_sub_section1 => section_vals_get_subs_vals(hfx_section1, "SCREENING", i_rep_section=irep)
3082 hfx_sub_section2 => section_vals_get_subs_vals(hfx_section2, "SCREENING", i_rep_section=irep)
3083
3084 CALL section_vals_val_get(hfx_sub_section1, "EPS_SCHWARZ", r_val=rval1, i_rep_section=irep)
3085 CALL section_vals_val_get(hfx_sub_section2, "EPS_SCHWARZ", r_val=rval2, i_rep_section=irep)
3086 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3087
3088 CALL section_vals_val_get(hfx_sub_section1, "EPS_SCHWARZ_FORCES", r_val=rval1, i_rep_section=irep)
3089 CALL section_vals_val_get(hfx_sub_section2, "EPS_SCHWARZ_FORCES", r_val=rval2, i_rep_section=irep)
3090 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3091
3092 CALL section_vals_val_get(hfx_sub_section1, "P_SCREEN_CORRECTION_FACTOR", r_val=rval1, i_rep_section=irep)
3093 CALL section_vals_val_get(hfx_sub_section2, "P_SCREEN_CORRECTION_FACTOR", r_val=rval2, i_rep_section=irep)
3094 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3095
3096 CALL section_vals_val_get(hfx_sub_section1, "SCREEN_ON_INITIAL_P", l_val=lval1, i_rep_section=irep)
3097 CALL section_vals_val_get(hfx_sub_section2, "SCREEN_ON_INITIAL_P", l_val=lval2, i_rep_section=irep)
3098 IF (lval1 .NEQV. lval2) is_identical = .false.
3099
3100 CALL section_vals_val_get(hfx_sub_section1, "SCREEN_P_FORCES", l_val=lval1, i_rep_section=irep)
3101 CALL section_vals_val_get(hfx_sub_section2, "SCREEN_P_FORCES", l_val=lval2, i_rep_section=irep)
3102 IF (lval1 .NEQV. lval2) is_identical = .false.
3103
3104 END DO
3105
3106 !Test of the fraction
3107 IF (is_identical) THEN
3108 DO irep = 1, n_rep_hf1
3109 CALL section_vals_val_get(hfx_section1, "FRACTION", r_val=rval1, i_rep_section=irep)
3110 CALL section_vals_val_get(hfx_section2, "FRACTION", r_val=rval2, i_rep_section=irep)
3111 IF (abs(rval1 - rval2) > epsilon(1.0_dp)) is_identical = .false.
3112 END DO
3113
3114 IF (PRESENT(same_except_frac)) THEN
3115 IF (.NOT. is_identical) same_except_frac = .true.
3116 END IF
3117 END IF
3118
3119 END SUBROUTINE compare_hfx_sections
3120
3121END MODULE hfx_types
3122
ncoset
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition grid_common.h:76
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:156
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:630
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::guidon2008
integer, save, public guidon2008
Definition bibliography.F:36
bibliography::guidon2009
integer, save, public guidon2009
Definition bibliography.F:36
bibliography::bussy2023
integer, save, public bussy2023
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:516
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:195
cell_types::plane_distance
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:252
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1133
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
cp_files::file_exists
logical function, public file_exists(file_name)
Checks if file exists, considering also the file discovery mechanism.
Definition cp_files.F:501
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1178
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
hfx_helpers
Some auxiliary functions and subroutines needed for HFX calculations.
Definition hfx_helpers.F:14
hfx_helpers::count_cells_perd
integer function, public count_cells_perd(shell, perd)
Auxiliary function for creating periodic neighbor cells
Definition hfx_helpers.F:38
hfx_helpers::next_image_cell_perd
subroutine, public next_image_cell_perd(m, perd)
Auxiliary function for creating periodic neighbor cells
Definition hfx_helpers.F:62
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
hfx_types::hfx_create
subroutine, public hfx_create(x_data, para_env, hfx_section, atomic_kind_set, qs_kind_set, particle_set, dft_control, cell, orb_basis, ri_basis, nelectron_total, nkp_grid)
This routine allocates and initializes all types in hfx_data
Definition hfx_types.F:595
hfx_types::hfx_init_container
subroutine, public hfx_init_container(container, memory_usage, do_disk_storage)
This routine deletes all list entries in a container in order to deallocate the memory.
Definition hfx_types.F:2527
hfx_types::hfx_set_distr_energy
subroutine, public hfx_set_distr_energy(ptr_to_distr, x_data)
This routine stores the data obtained from the load balance routine for the energy
Definition hfx_types.F:2572
hfx_types::hfx_set_distr_forces
subroutine, public hfx_set_distr_forces(ptr_to_distr, x_data)
This routine stores the data obtained from the load balance routine for the forces
Definition hfx_types.F:2592
hfx_types::max_atom_block
integer, parameter, public max_atom_block
Definition hfx_types.F:117
hfx_types::parse_memory_section
subroutine, public parse_memory_section(memory_parameter, hf_sub_section, storage_id, i_thread, n_threads, para_env, irep, skip_disk, skip_in_core_forces)
Parses the memory section
Definition hfx_types.F:1818
hfx_types::hfx_release_basis_types
subroutine, public hfx_release_basis_types(basis_parameter)
...
Definition hfx_types.F:1786
hfx_types::init_t_c_g0_lmax
integer, save, public init_t_c_g0_lmax
Definition hfx_types.F:134
hfx_types::log_zero
real(dp), parameter, public log_zero
Definition hfx_types.F:119
hfx_types::max_images
integer, parameter, public max_images
Definition hfx_types.F:118
hfx_types::hfx_release
subroutine, public hfx_release(x_data)
This routine deallocates all data structures
Definition hfx_types.F:1910
hfx_types::alloc_containers
subroutine, public alloc_containers(data, bin_size)
...
Definition hfx_types.F:2910
hfx_types::hfx_create_neighbor_cells
subroutine, public hfx_create_neighbor_cells(x_data, pbc_shells, cell, i_thread, nkp_grid)
This routine computes the neighbor cells that are taken into account in periodic runs
Definition hfx_types.F:2054
hfx_types::dealloc_containers
subroutine, public dealloc_containers(data, memory_usage)
...
Definition hfx_types.F:2878
hfx_types::hfx_create_basis_types
subroutine, public hfx_create_basis_types(basis_parameter, basis_info, qs_kind_set, basis_type)
This routine allocates and initializes the basis_info and basis_parameter types
Definition hfx_types.F:1660
hfx_types::hfx_ri_init
subroutine, public hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
...
Definition hfx_types.F:1212
hfx_types::compare_hfx_sections
subroutine, public compare_hfx_sections(hfx_section1, hfx_section2, is_identical, same_except_frac)
Compares the non-technical parts of two HFX input section and check whether they are the same Ignore ...
Definition hfx_types.F:2957
hfx_types::mul_fact
real(kind=dp), dimension(0:10), parameter, public mul_fact
Definition hfx_types.F:121
hfx_types::powell_min_log
real(dp), parameter, public powell_min_log
Definition hfx_types.F:120
hfx_types::hfx_reset_memory_usage_counter
subroutine, public hfx_reset_memory_usage_counter(memory_parameter, subtr_size_mb)
resets the maximum memory usage for a HFX calculation subtracting all relevant buffers from the input...
Definition hfx_types.F:2613
hfx_types::hfx_ri_release
subroutine, public hfx_ri_release(ri_data, write_stats)
...
Definition hfx_types.F:1469
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::hfx_ri_do_2c_diag
integer, parameter, public hfx_ri_do_2c_diag
Definition input_constants.F:834
input_constants::do_potential_mix_cl
integer, parameter, public do_potential_mix_cl
Definition input_constants.F:808
input_constants::do_potential_gaussian
integer, parameter, public do_potential_gaussian
Definition input_constants.F:808
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition input_constants.F:808
input_constants::do_potential_mix_lg
integer, parameter, public do_potential_mix_lg
Definition input_constants.F:808
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:808
input_constants::hfx_ri_do_2c_iter
integer, parameter, public hfx_ri_do_2c_iter
Definition input_constants.F:834
input_constants::do_hfx_auto_shells
integer, parameter, public do_hfx_auto_shells
Definition input_constants.F:821
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition input_constants.F:808
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition input_constants.F:808
input_constants::do_potential_mix_cl_trunc
integer, parameter, public do_potential_mix_cl_trunc
Definition input_constants.F:808
input_constants::do_potential_long
integer, parameter, public do_potential_long
Definition input_constants.F:808
input_cp2k_hfx
function that builds the hartree fock exchange section of the input
Definition input_cp2k_hfx.F:14
input_cp2k_hfx::ri_pmat
integer, parameter, public ri_pmat
Definition input_cp2k_hfx.F:44
input_cp2k_hfx::ri_mo
integer, parameter, public ri_mo
Definition input_cp2k_hfx.F:44
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::compare_potential_types
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare libint_potential_types.
Definition libint_2c_3c.F:1062
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition libint_wrapper.F:15
libint_wrapper::cp_libint_init_eri1
subroutine, public cp_libint_init_eri1(lib, max_am)
Definition libint_wrapper.F:1116
libint_wrapper::prim_data_f_size
integer, parameter, public prim_data_f_size
Definition libint_wrapper.F:63
libint_wrapper::cp_libint_cleanup_eri1
subroutine, public cp_libint_cleanup_eri1(lib)
Definition libint_wrapper.F:1186
libint_wrapper::cp_libint_static_cleanup
subroutine, public cp_libint_static_cleanup()
Definition libint_wrapper.F:1244
libint_wrapper::cp_libint_init_eri
subroutine, public cp_libint_init_eri(lib, max_am)
Definition libint_wrapper.F:1086
libint_wrapper::cp_libint_static_init
subroutine, public cp_libint_static_init()
Definition libint_wrapper.F:1236
libint_wrapper::cp_libint_cleanup_eri
subroutine, public cp_libint_cleanup_eri(lib)
Definition libint_wrapper.F:1176
libint_wrapper::cp_libint_set_contrdepth
subroutine, public cp_libint_set_contrdepth(lib, contrdepth)
Definition libint_wrapper.F:1098
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_getcwd
subroutine, public m_getcwd(curdir)
...
Definition machine.F:616
machine::m_chdir
subroutine, public m_chdir(dir, ierror)
...
Definition machine.F:645
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::erfc_cutoff
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Definition mathlib.F:1686
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:454
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:879
qs_tensors_types
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors_types.F:12
qs_tensors_types::distribution_3d_create
subroutine, public distribution_3d_create(dist_3d, dist1, dist2, dist3, nkind, particle_set, mp_comm_3d, own_comm)
Create a 3d distribution.
Definition qs_tensors_types.F:91
qs_tensors_types::default_block_size
integer, parameter, public default_block_size
Definition qs_tensors_types.F:54
qs_tensors_types::create_2c_tensor
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Definition qs_tensors_types.F:382
qs_tensors_types::split_block_sizes
subroutine, public split_block_sizes(blk_sizes, blk_sizes_split, max_size)
...
Definition qs_tensors_types.F:428
qs_tensors_types::pgf_block_sizes
subroutine, public pgf_block_sizes(atomic_kind_set, basis, min_blk_size, pgf_blk_sizes)
...
Definition qs_tensors_types.F:464
qs_tensors_types::distribution_3d_destroy
subroutine, public distribution_3d_destroy(dist)
Destroy a 3d distribution.
Definition qs_tensors_types.F:135
qs_tensors_types::create_tensor_batches
subroutine, public create_tensor_batches(sizes, nbatches, starts_array, ends_array, starts_array_block, ends_array_block)
...
Definition qs_tensors_types.F:522
qs_tensors_types::create_3c_tensor
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Definition qs_tensors_types.F:335
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::compress
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Definition string_utilities.F:3191
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition t_c_g0.F:57
t_c_g0::free_c0
subroutine, public free_c0()
...
Definition t_c_g0.F:1388
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_array_utils::cp_1d_logical_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:282
cp_control_types::dft_control_type
Definition cp_control_types.F:599
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
hfx_types::block_ind_type
Definition hfx_types.F:368
hfx_types::hfx_2d_map
Definition hfx_types.F:302
hfx_types::hfx_basis_info_type
Definition hfx_types.F:285
hfx_types::hfx_basis_type
Definition hfx_types.F:264
hfx_types::hfx_block_range_type
Definition hfx_types.F:340
hfx_types::hfx_cache_type
Definition hfx_types.F:242
hfx_types::hfx_cell_type
Definition hfx_types.F:205
hfx_types::hfx_compression_type
Definition hfx_types.F:352
hfx_types::hfx_container_type
Definition hfx_types.F:254
hfx_types::hfx_distribution
Definition hfx_types.F:211
hfx_types::hfx_general_type
Definition hfx_types.F:199
hfx_types::hfx_load_balance_type
Definition hfx_types.F:189
hfx_types::hfx_memory_type
Definition hfx_types.F:158
hfx_types::hfx_p_kind
Definition hfx_types.F:297
hfx_types::hfx_pgf_list
Definition hfx_types.F:319
hfx_types::hfx_pgf_product_list
Definition hfx_types.F:330
hfx_types::hfx_potential_type
Definition hfx_types.F:139
hfx_types::hfx_ri_type
Definition hfx_types.F:372
hfx_types::hfx_screen_coeff_type
Definition hfx_types.F:292
hfx_types::hfx_screening_type
Definition hfx_types.F:150
hfx_types::hfx_task_list_type
Definition hfx_types.F:346
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:510
hfx_types::pair_list_element_type
Definition hfx_types.F:222
hfx_types::pair_list_type
Definition hfx_types.F:236
hfx_types::pair_set_list_type
Definition hfx_types.F:231
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
libint_2c_3c::libint_potential_type
Definition libint_2c_3c.F:67
libint_wrapper::cp_libint_t
Definition libint_wrapper.F:69
message_passing::mp_cart_type
Definition message_passing.F:740
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_tensors_types::distribution_3d_type
Definition qs_tensors_types.F:57