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