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