(git:34ef472)
qs_tensors.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Utility methods to build 3-center integral tensors of various types.
10 ! **************************************************************************************************
11 MODULE qs_tensors
12  USE omp_lib, ONLY: omp_get_num_threads,&
13  omp_get_thread_num
14  USE ai_contraction, ONLY: block_add
15  USE ai_contraction_sphi, ONLY: ab_contract,&
16  libxsmm_abc_contract
17  USE atomic_kind_types, ONLY: atomic_kind_type
18  USE basis_set_types, ONLY: get_gto_basis_set,&
19  gto_basis_set_p_type,&
20  gto_basis_set_type
21  USE block_p_types, ONLY: block_p_type
22  USE cell_types, ONLY: cell_type,&
23  real_to_scaled
24  USE cp_array_utils, ONLY: cp_2d_r_p_type
25  USE cp_control_types, ONLY: dft_control_type
26  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
27  USE cp_files, ONLY: close_file,&
28  open_file
29  USE dbcsr_api, ONLY: dbcsr_filter,&
30  dbcsr_finalize,&
31  dbcsr_get_block_p,&
32  dbcsr_get_matrix_type,&
33  dbcsr_has_symmetry,&
34  dbcsr_type,&
35  dbcsr_type_antisymmetric,&
36  dbcsr_type_no_symmetry
37  USE dbt_api, ONLY: &
38  dbt_blk_sizes, dbt_clear, dbt_copy, dbt_create, dbt_destroy, dbt_filter, dbt_get_block, &
39  dbt_get_info, dbt_get_num_blocks, dbt_get_nze_total, dbt_get_stored_coordinates, &
40  dbt_iterator_next_block, dbt_iterator_num_blocks, dbt_iterator_start, dbt_iterator_stop, &
41  dbt_iterator_type, dbt_ndims, dbt_put_block, dbt_reserve_blocks, dbt_type
42  USE distribution_1d_types, ONLY: distribution_1d_type
43  USE distribution_2d_types, ONLY: distribution_2d_type
44  USE gamma, ONLY: init_md_ftable
45  USE hfx_compression_methods, ONLY: hfx_add_mult_cache_elements,&
46  hfx_add_single_cache_element,&
47  hfx_decompress_first_cache,&
48  hfx_flush_last_cache,&
49  hfx_get_mult_cache_elements,&
50  hfx_get_single_cache_element,&
51  hfx_reset_cache_and_container
52  USE hfx_types, ONLY: alloc_containers,&
53  dealloc_containers,&
54  hfx_cache_type,&
55  hfx_compression_type,&
56  hfx_container_type,&
57  hfx_init_container
58  USE input_constants, ONLY: do_potential_coulomb,&
59  do_potential_id,&
60  do_potential_short,&
61  do_potential_truncated
62  USE input_section_types, ONLY: section_vals_val_get
63  USE kinds, ONLY: dp,&
64  int_8
65  USE kpoint_types, ONLY: get_kpoint_info,&
66  kpoint_type
67  USE libint_2c_3c, ONLY: cutoff_screen_factor,&
68  eri_2center,&
69  eri_2center_derivs,&
70  eri_3center,&
71  eri_3center_derivs,&
72  libint_potential_type
73  USE libint_wrapper, ONLY: &
74  cp_libint_cleanup_2eri, cp_libint_cleanup_2eri1, cp_libint_cleanup_3eri, &
75  cp_libint_cleanup_3eri1, cp_libint_init_2eri, cp_libint_init_2eri1, cp_libint_init_3eri, &
76  cp_libint_init_3eri1, cp_libint_set_contrdepth, cp_libint_t
77  USE message_passing, ONLY: mp_para_env_type
78  USE molecule_types, ONLY: molecule_type
79  USE orbital_pointers, ONLY: ncoset
80  USE particle_types, ONLY: particle_type
81  USE qs_environment_types, ONLY: get_qs_env,&
82  qs_environment_type
83  USE qs_kind_types, ONLY: qs_kind_type
84  USE qs_neighbor_list_types, ONLY: &
85  get_iterator_info, get_neighbor_list_set_p, neighbor_list_iterate, &
86  neighbor_list_iterator_create, neighbor_list_iterator_p_type, &
87  neighbor_list_iterator_release, neighbor_list_set_p_type, nl_sub_iterate, &
88  release_neighbor_list_sets
89  USE qs_neighbor_lists, ONLY: atom2d_build,&
90  atom2d_cleanup,&
91  build_neighbor_lists,&
92  local_atoms_type,&
93  pair_radius_setup
94  USE qs_tensors_types, ONLY: &
95  distribution_3d_destroy, distribution_3d_type, neighbor_list_3c_iterator_type, &
96  neighbor_list_3c_type, symmetric_ij, symmetric_ijk, symmetric_jk, symmetric_none, &
97  symmetrik_ik
98  USE t_c_g0, ONLY: get_lmax_init,&
99  init
100  USE util, ONLY: get_limit
101 
102 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
103 #include "./base/base_uses.f90"
104 
105  IMPLICIT NONE
106 
107  PRIVATE
108 
109  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tensors'
110 
111  PUBLIC :: build_3c_neighbor_lists, &
112  neighbor_list_3c_destroy, neighbor_list_3c_iterate, neighbor_list_3c_iterator_create, &
113  neighbor_list_3c_iterator_destroy, get_3c_iterator_info, build_3c_integrals, &
114  build_2c_neighbor_lists, build_2c_integrals, cutoff_screen_factor, &
115  get_tensor_occupancy, compress_tensor, decompress_tensor, &
116  build_3c_derivatives, build_2c_derivatives, calc_2c_virial, calc_3c_virial
117 
118  TYPE one_dim_int_array
119  INTEGER, DIMENSION(:), ALLOCATABLE :: array
120  END TYPE
121 
122  ! cache size for integral compression
123  INTEGER, PARAMETER, PRIVATE :: cache_size = 1024
124 
125 CONTAINS
126 
127 ! **************************************************************************************************
128 !> \brief Build 2-center neighborlists adapted to different operators
129 !> This mainly wraps build_neighbor_lists for consistency with build_3c_neighbor_lists
130 !> \param ij_list 2c neighbor list for atom pairs i, j
131 !> \param basis_i basis object for atoms i
132 !> \param basis_j basis object for atoms j
133 !> \param potential_parameter ...
134 !> \param name name of 2c neighbor list
135 !> \param qs_env ...
136 !> \param sym_ij Symmetry in i, j (default .TRUE.)
137 !> \param molecular ...
138 !> \param dist_2d optionally a custom 2d distribution
139 !> \param pot_to_rad which radius (1 for i, 2 for j) should be adapted to incorporate potential
140 ! **************************************************************************************************
141  SUBROUTINE build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, &
142  sym_ij, molecular, dist_2d, pot_to_rad)
143  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
144  POINTER :: ij_list
145  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
146  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
147  CHARACTER(LEN=*), INTENT(IN) :: name
148  TYPE(qs_environment_type), POINTER :: qs_env
149  LOGICAL, INTENT(IN), OPTIONAL :: sym_ij, molecular
150  TYPE(distribution_2d_type), OPTIONAL, POINTER :: dist_2d
151  INTEGER, INTENT(IN), OPTIONAL :: pot_to_rad
152 
153  INTEGER :: ikind, nkind, pot_to_rad_prv
154  LOGICAL, ALLOCATABLE, DIMENSION(:) :: i_present, j_present
155  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
156  REAL(kind=dp) :: subcells
157  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: i_radius, j_radius
158  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
159  TYPE(cell_type), POINTER :: cell
160  TYPE(distribution_1d_type), POINTER :: local_particles
161  TYPE(distribution_2d_type), POINTER :: dist_2d_prv
162  TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
163  TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
164  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
165 
166  NULLIFY (atomic_kind_set, cell, local_particles, molecule_set, &
167  particle_set, dist_2d_prv)
168 
169  IF (PRESENT(pot_to_rad)) THEN
170  pot_to_rad_prv = pot_to_rad
171  ELSE
172  pot_to_rad_prv = 1
173  END IF
174 
175  CALL get_qs_env(qs_env, &
176  nkind=nkind, &
177  cell=cell, &
178  particle_set=particle_set, &
179  atomic_kind_set=atomic_kind_set, &
180  local_particles=local_particles, &
181  distribution_2d=dist_2d_prv, &
182  molecule_set=molecule_set)
183 
184  CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
185 
186  ALLOCATE (i_present(nkind), j_present(nkind))
187  ALLOCATE (i_radius(nkind), j_radius(nkind))
188 
189  i_present = .false.
190  j_present = .false.
191  i_radius = 0.0_dp
192  j_radius = 0.0_dp
193 
194  IF (PRESENT(dist_2d)) dist_2d_prv => dist_2d
195 
196  ! Set up the radii, depending on the operator type
197  IF (potential_parameter%potential_type == do_potential_id) THEN
198 
199  !overlap => use the kind radius for both i and j
200  DO ikind = 1, nkind
201  IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
202  i_present(ikind) = .true.
203  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
204  END IF
205  IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
206  j_present(ikind) = .true.
207  CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
208  END IF
209  END DO
210 
211  ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
212 
213  !Coulomb operator, virtually infinite range => set j_radius to arbitrarily large number
214  DO ikind = 1, nkind
215  IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
216  i_present(ikind) = .true.
217  IF (pot_to_rad_prv == 1) i_radius(ikind) = 1000000.0_dp
218  END IF
219  IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
220  j_present(ikind) = .true.
221  IF (pot_to_rad_prv == 2) j_radius(ikind) = 1000000.0_dp
222  END IF
223  END DO !ikind
224 
225  ELSE IF (potential_parameter%potential_type == do_potential_truncated .OR. &
226  potential_parameter%potential_type == do_potential_short) THEN
227 
228  !Truncated coulomb/short range: set j_radius to r_cutoff + the kind_radii
229  DO ikind = 1, nkind
230  IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
231  i_present(ikind) = .true.
232  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
233  IF (pot_to_rad_prv == 1) i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
234  END IF
235  IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
236  j_present(ikind) = .true.
237  CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
238  IF (pot_to_rad_prv == 2) j_radius(ikind) = j_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
239  END IF
240  END DO
241 
242  ELSE
243  cpabort("Operator not implemented.")
244  END IF
245 
246  ALLOCATE (pair_radius(nkind, nkind))
247  pair_radius = 0.0_dp
248  CALL pair_radius_setup(i_present, j_present, i_radius, j_radius, pair_radius)
249 
250  ALLOCATE (atom2d(nkind))
251 
252  CALL atom2d_build(atom2d, local_particles, dist_2d_prv, atomic_kind_set, &
253  molecule_set, molecule_only=.false., particle_set=particle_set)
254  CALL build_neighbor_lists(ij_list, particle_set, atom2d, cell, pair_radius, subcells, &
255  symmetric=sym_ij, molecular=molecular, nlname=trim(name))
256 
257  CALL atom2d_cleanup(atom2d)
258 
259  END SUBROUTINE
260 
261 ! **************************************************************************************************
262 !> \brief Build a 3-center neighbor list
263 !> \param ijk_list 3c neighbor list for atom triples i, j, k
264 !> \param basis_i basis object for atoms i
265 !> \param basis_j basis object for atoms j
266 !> \param basis_k basis object for atoms k
267 !> \param dist_3d 3d distribution object
268 !> \param potential_parameter ...
269 !> \param name name of 3c neighbor list
270 !> \param qs_env ...
271 !> \param sym_ij Symmetry in i, j (default .FALSE.)
272 !> \param sym_jk Symmetry in j, k (default .FALSE.)
273 !> \param sym_ik Symmetry in i, k (default .FALSE.)
274 !> \param molecular ??? not tested
275 !> \param op_pos ...
276 !> \param own_dist ...
277 ! **************************************************************************************************
278  SUBROUTINE build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, &
279  dist_3d, potential_parameter, name, qs_env, &
280  sym_ij, sym_jk, sym_ik, molecular, op_pos, &
281  own_dist)
282  TYPE(neighbor_list_3c_type), INTENT(OUT) :: ijk_list
283  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
284  TYPE(distribution_3d_type), INTENT(IN) :: dist_3d
285  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
286  CHARACTER(LEN=*), INTENT(IN) :: name
287  TYPE(qs_environment_type), POINTER :: qs_env
288  LOGICAL, INTENT(IN), OPTIONAL :: sym_ij, sym_jk, sym_ik, molecular
289  INTEGER, INTENT(IN), OPTIONAL :: op_pos
290  LOGICAL, INTENT(IN), OPTIONAL :: own_dist
291 
292  CHARACTER(len=*), PARAMETER :: routinen = 'build_3c_neighbor_lists'
293 
294  INTEGER :: handle, op_pos_prv, sym_level
295  TYPE(libint_potential_type) :: pot_par_1, pot_par_2
296 
297  CALL timeset(routinen, handle)
298 
299  IF (PRESENT(op_pos)) THEN
300  op_pos_prv = op_pos
301  ELSE
302  op_pos_prv = 1
303  END IF
304 
305  SELECT CASE (op_pos_prv)
306  CASE (1)
307  pot_par_1 = potential_parameter
308  pot_par_2%potential_type = do_potential_id
309  CASE (2)
310  pot_par_2 = potential_parameter
311  pot_par_1%potential_type = do_potential_id
312  END SELECT
313 
314  CALL build_2c_neighbor_lists(ijk_list%ij_list, basis_i, basis_j, pot_par_1, trim(name)//"_sub_1", &
315  qs_env, sym_ij=.false., molecular=molecular, &
316  dist_2d=dist_3d%dist_2d_1, pot_to_rad=1)
317 
318  CALL build_2c_neighbor_lists(ijk_list%jk_list, basis_j, basis_k, pot_par_2, trim(name)//"_sub_2", &
319  qs_env, sym_ij=.false., molecular=molecular, &
320  dist_2d=dist_3d%dist_2d_2, pot_to_rad=2)
321 
322  ijk_list%sym = symmetric_none
323 
324  sym_level = 0
325  IF (PRESENT(sym_ij)) THEN
326  IF (sym_ij) THEN
327  ijk_list%sym = symmetric_ij
328  sym_level = sym_level + 1
329  END IF
330  END IF
331 
332  IF (PRESENT(sym_jk)) THEN
333  IF (sym_jk) THEN
334  ijk_list%sym = symmetric_jk
335  sym_level = sym_level + 1
336  END IF
337  END IF
338 
339  IF (PRESENT(sym_ik)) THEN
340  IF (sym_ik) THEN
341  ijk_list%sym = symmetrik_ik
342  sym_level = sym_level + 1
343  END IF
344  END IF
345 
346  IF (sym_level >= 2) THEN
347  ijk_list%sym = symmetric_ijk
348  END IF
349 
350  ijk_list%dist_3d = dist_3d
351  IF (PRESENT(own_dist)) THEN
352  ijk_list%owns_dist = own_dist
353  ELSE
354  ijk_list%owns_dist = .false.
355  END IF
356 
357  CALL timestop(handle)
358  END SUBROUTINE
359 
360 ! **************************************************************************************************
361 !> \brief Symmetry criterion
362 !> \param a ...
363 !> \param b ...
364 !> \return ...
365 ! **************************************************************************************************
366  PURE FUNCTION include_symmetric(a, b)
367  INTEGER, INTENT(IN) :: a, b
368  LOGICAL :: include_symmetric
369 
370  IF (a > b) THEN
371  include_symmetric = (modulo(a + b, 2) /= 0)
372  ELSE
373  include_symmetric = (modulo(a + b, 2) == 0)
374  END IF
375 
376  END FUNCTION
377 
378 ! **************************************************************************************************
379 !> \brief Destroy 3c neighborlist
380 !> \param ijk_list ...
381 ! **************************************************************************************************
382  SUBROUTINE neighbor_list_3c_destroy(ijk_list)
383  TYPE(neighbor_list_3c_type), INTENT(INOUT) :: ijk_list
384 
385  CALL release_neighbor_list_sets(ijk_list%ij_list)
386  CALL release_neighbor_list_sets(ijk_list%jk_list)
387 
388  IF (ijk_list%owns_dist) THEN
389  CALL distribution_3d_destroy(ijk_list%dist_3d)
390  END IF
391 
392  END SUBROUTINE
393 
394 ! **************************************************************************************************
395 !> \brief Create a 3-center neighborlist iterator
396 !> \param iterator ...
397 !> \param ijk_nl ...
398 ! **************************************************************************************************
399  SUBROUTINE neighbor_list_3c_iterator_create(iterator, ijk_nl)
400  TYPE(neighbor_list_3c_iterator_type), INTENT(OUT) :: iterator
401  TYPE(neighbor_list_3c_type), INTENT(IN) :: ijk_nl
402 
403  CHARACTER(len=*), PARAMETER :: routinen = 'neighbor_list_3c_iterator_create'
404 
405  INTEGER :: handle
406 
407  CALL timeset(routinen, handle)
408 
409  CALL neighbor_list_iterator_create(iterator%iter_ij, ijk_nl%ij_list)
410  CALL neighbor_list_iterator_create(iterator%iter_jk, ijk_nl%jk_list, search=.true.)
411 
412  iterator%iter_level = 0
413  iterator%ijk_nl = ijk_nl
414 
415  iterator%bounds_i = 0
416  iterator%bounds_j = 0
417  iterator%bounds_k = 0
418 
419  CALL timestop(handle)
420  END SUBROUTINE
421 
422 ! **************************************************************************************************
423 !> \brief impose atomic upper and lower bounds
424 !> \param iterator ...
425 !> \param bounds_i ...
426 !> \param bounds_j ...
427 !> \param bounds_k ...
428 ! **************************************************************************************************
429  SUBROUTINE nl_3c_iter_set_bounds(iterator, bounds_i, bounds_j, bounds_k)
430  TYPE(neighbor_list_3c_iterator_type), &
431  INTENT(INOUT) :: iterator
432  INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
433 
434  IF (PRESENT(bounds_i)) iterator%bounds_i = bounds_i
435  IF (PRESENT(bounds_j)) iterator%bounds_j = bounds_j
436  IF (PRESENT(bounds_k)) iterator%bounds_k = bounds_k
437 
438  END SUBROUTINE
439 
440 ! **************************************************************************************************
441 !> \brief Destroy 3c-nl iterator
442 !> \param iterator ...
443 ! **************************************************************************************************
444  SUBROUTINE neighbor_list_3c_iterator_destroy(iterator)
445  TYPE(neighbor_list_3c_iterator_type), &
446  INTENT(INOUT) :: iterator
447 
448  CHARACTER(len=*), PARAMETER :: routinen = 'neighbor_list_3c_iterator_destroy'
449 
450  INTEGER :: handle
451 
452  CALL timeset(routinen, handle)
453  CALL neighbor_list_iterator_release(iterator%iter_ij)
454  CALL neighbor_list_iterator_release(iterator%iter_jk)
455  NULLIFY (iterator%iter_ij)
456  NULLIFY (iterator%iter_jk)
457 
458  CALL timestop(handle)
459  END SUBROUTINE
460 
461 ! **************************************************************************************************
462 !> \brief Iterate 3c-nl iterator
463 !> \param iterator ...
464 !> \return 0 if successful; 1 if end was reached
465 ! **************************************************************************************************
466  RECURSIVE FUNCTION neighbor_list_3c_iterate(iterator) RESULT(iter_stat)
467  TYPE(neighbor_list_3c_iterator_type), &
468  INTENT(INOUT) :: iterator
469  INTEGER :: iter_stat
470 
471  INTEGER :: iatom, iter_level, jatom, jatom_1, &
472  jatom_2, katom
473  LOGICAL :: skip_this
474 
475  iter_level = iterator%iter_level
476 
477  IF (iter_level == 0) THEN
478  iter_stat = neighbor_list_iterate(iterator%iter_ij)
479 
480  IF (iter_stat /= 0) THEN
481  RETURN
482  END IF
483 
484  CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom)
485  skip_this = .false.
486  IF ((iterator%bounds_i(1) > 0 .AND. iatom < iterator%bounds_i(1)) &
487  .OR. (iterator%bounds_i(2) > 0 .AND. iatom > iterator%bounds_i(2))) skip_this = .true.
488  IF ((iterator%bounds_j(1) > 0 .AND. jatom < iterator%bounds_j(1)) &
489  .OR. (iterator%bounds_j(2) > 0 .AND. jatom > iterator%bounds_j(2))) skip_this = .true.
490 
491  IF (skip_this) THEN
492  iter_stat = neighbor_list_3c_iterate(iterator)
493  RETURN
494  END IF
495 
496  END IF
497  iter_stat = nl_sub_iterate(iterator%iter_jk, iterator%iter_ij)
498  IF (iter_stat /= 0) THEN
499  iterator%iter_level = 0
500  iter_stat = neighbor_list_3c_iterate(iterator)
501  RETURN
502  ELSE
503  iterator%iter_level = 1
504  END IF
505 
506  cpassert(iter_stat == 0)
507  cpassert(iterator%iter_level == 1)
508  CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom_1)
509  CALL get_iterator_info(iterator%iter_jk, iatom=jatom_2, jatom=katom)
510 
511  cpassert(jatom_1 == jatom_2)
512  jatom = jatom_1
513 
514  skip_this = .false.
515  IF ((iterator%bounds_k(1) > 0 .AND. katom < iterator%bounds_k(1)) &
516  .OR. (iterator%bounds_k(2) > 0 .AND. katom > iterator%bounds_k(2))) skip_this = .true.
517 
518  IF (skip_this) THEN
519  iter_stat = neighbor_list_3c_iterate(iterator)
520  RETURN
521  END IF
522 
523  SELECT CASE (iterator%ijk_nl%sym)
524  CASE (symmetric_none)
525  skip_this = .false.
526  CASE (symmetric_ij)
527  skip_this = .NOT. include_symmetric(iatom, jatom)
528  CASE (symmetric_jk)
529  skip_this = .NOT. include_symmetric(jatom, katom)
530  CASE (symmetrik_ik)
531  skip_this = .NOT. include_symmetric(iatom, katom)
532  CASE (symmetric_ijk)
533  skip_this = .NOT. include_symmetric(iatom, jatom) .OR. .NOT. include_symmetric(jatom, katom)
534  CASE DEFAULT
535  cpabort("should not happen")
536  END SELECT
537 
538  IF (skip_this) THEN
539  iter_stat = neighbor_list_3c_iterate(iterator)
540  RETURN
541  END IF
542 
543  END FUNCTION
544 
545 ! **************************************************************************************************
546 !> \brief Get info of current iteration
547 !> \param iterator ...
548 !> \param ikind ...
549 !> \param jkind ...
550 !> \param kkind ...
551 !> \param nkind ...
552 !> \param iatom ...
553 !> \param jatom ...
554 !> \param katom ...
555 !> \param rij ...
556 !> \param rjk ...
557 !> \param rik ...
558 !> \param cell_j ...
559 !> \param cell_k ...
560 !> \return ...
561 ! **************************************************************************************************
562  SUBROUTINE get_3c_iterator_info(iterator, ikind, jkind, kkind, nkind, iatom, jatom, katom, &
563  rij, rjk, rik, cell_j, cell_k)
564  TYPE(neighbor_list_3c_iterator_type), &
565  INTENT(INOUT) :: iterator
566  INTEGER, INTENT(OUT), OPTIONAL :: ikind, jkind, kkind, nkind, iatom, &
567  jatom, katom
568  REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: rij, rjk, rik
569  INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: cell_j, cell_k
570 
571  INTEGER, DIMENSION(2) :: atoms_1, atoms_2, kinds_1, kinds_2
572  INTEGER, DIMENSION(3) :: cell_1, cell_2
573  REAL(kind=dp), DIMENSION(3) :: r_1, r_2
574 
575  cpassert(iterator%iter_level == 1)
576 
577  CALL get_iterator_info(iterator%iter_ij, &
578  ikind=kinds_1(1), jkind=kinds_1(2), nkind=nkind, &
579  iatom=atoms_1(1), jatom=atoms_1(2), r=r_1, &
580  cell=cell_1)
581 
582  CALL get_iterator_info(iterator%iter_jk, &
583  ikind=kinds_2(1), jkind=kinds_2(2), &
584  iatom=atoms_2(1), jatom=atoms_2(2), r=r_2, &
585  cell=cell_2)
586 
587  IF (PRESENT(ikind)) ikind = kinds_1(1)
588  IF (PRESENT(jkind)) jkind = kinds_1(2)
589  IF (PRESENT(kkind)) kkind = kinds_2(2)
590  IF (PRESENT(iatom)) iatom = atoms_1(1)
591  IF (PRESENT(jatom)) jatom = atoms_1(2)
592  IF (PRESENT(katom)) katom = atoms_2(2)
593 
594  IF (PRESENT(rij)) rij = r_1
595  IF (PRESENT(rjk)) rjk = r_2
596  IF (PRESENT(rik)) rik = r_1 + r_2
597 
598  IF (PRESENT(cell_j)) cell_j = cell_1
599  IF (PRESENT(cell_k)) cell_k = cell_1 + cell_2
600 
601  END SUBROUTINE
602 
603 ! **************************************************************************************************
604 !> \brief Allocate blocks of a 3-center tensor based on neighborlist
605 !> \param t3c empty DBCSR tensor
606 !> Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages)
607 !> if k-points are used
608 !> \param nl_3c 3-center neighborlist
609 !> \param basis_i ...
610 !> \param basis_j ...
611 !> \param basis_k ...
612 !> \param qs_env ...
613 !> \param potential_parameter ...
614 !> \param op_pos ...
615 !> \param do_kpoints ...
616 !> \param do_hfx_kpoints ...
617 !> \param bounds_i ...
618 !> \param bounds_j ...
619 !> \param bounds_k ...
620 !> \param RI_range ...
621 !> \param img_to_RI_cell ...
622 ! **************************************************************************************************
623  SUBROUTINE alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, op_pos, &
624  do_kpoints, do_hfx_kpoints, bounds_i, bounds_j, bounds_k, RI_range, &
625  img_to_RI_cell)
626  TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT) :: t3c
627  TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
628  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
629  TYPE(qs_environment_type), POINTER :: qs_env
630  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
631  INTEGER, INTENT(IN), OPTIONAL :: op_pos
632  LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints
633  INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
634  REAL(dp), INTENT(IN), OPTIONAL :: ri_range
635  INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: img_to_ri_cell
636 
637  CHARACTER(LEN=*), PARAMETER :: routinen = 'alloc_block_3c'
638 
639  INTEGER :: handle, iatom, ikind, j_img, jatom, &
640  jcell, jkind, k_img, katom, kcell, &
641  kkind, natom, ncell_ri, nimg, op_ij, &
642  op_jk, op_pos_prv
643  INTEGER(int_8) :: a, b, nblk_per_thread
644  INTEGER(int_8), ALLOCATABLE, DIMENSION(:, :) :: nblk
645  INTEGER, ALLOCATABLE, DIMENSION(:) :: img_to_ri_cell_prv
646  INTEGER, DIMENSION(3) :: blk_idx, cell_j, cell_k, &
647  kp_index_lbounds, kp_index_ubounds
648  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
649  LOGICAL :: do_hfx_kpoints_prv, do_kpoints_prv
650  REAL(kind=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
651  kind_radius_i, kind_radius_j, &
652  kind_radius_k
653  REAL(kind=dp), DIMENSION(3) :: rij, rik, rjk
654  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
655  TYPE(cell_type), POINTER :: cell
656  TYPE(dft_control_type), POINTER :: dft_control
657  TYPE(kpoint_type), POINTER :: kpoints
658  TYPE(mp_para_env_type), POINTER :: para_env
659  TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
660  TYPE(one_dim_int_array), ALLOCATABLE, &
661  DIMENSION(:, :) :: alloc_i, alloc_j, alloc_k
662  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
663 
664  CALL timeset(routinen, handle)
665  NULLIFY (qs_kind_set, atomic_kind_set, cell)
666 
667  IF (PRESENT(do_kpoints)) THEN
668  do_kpoints_prv = do_kpoints
669  ELSE
670  do_kpoints_prv = .false.
671  END IF
672  IF (PRESENT(do_hfx_kpoints)) THEN
673  do_hfx_kpoints_prv = do_hfx_kpoints
674  ELSE
675  do_hfx_kpoints_prv = .false.
676  END IF
677  IF (do_hfx_kpoints_prv) THEN
678  cpassert(do_kpoints_prv)
679  cpassert(PRESENT(ri_range))
680  cpassert(PRESENT(img_to_ri_cell))
681  END IF
682 
683  IF (PRESENT(img_to_ri_cell)) THEN
684  ALLOCATE (img_to_ri_cell_prv(SIZE(img_to_ri_cell)))
685  img_to_ri_cell_prv(:) = img_to_ri_cell
686  END IF
687 
688  dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
689 
690  op_ij = do_potential_id; op_jk = do_potential_id
691 
692  IF (PRESENT(op_pos)) THEN
693  op_pos_prv = op_pos
694  ELSE
695  op_pos_prv = 1
696  END IF
697 
698  SELECT CASE (op_pos_prv)
699  CASE (1)
700  op_ij = potential_parameter%potential_type
701  CASE (2)
702  op_jk = potential_parameter%potential_type
703  END SELECT
704 
705  IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
706  dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
707  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
708  ELSEIF (op_ij == do_potential_coulomb) THEN
709  dr_ij = 1000000.0_dp
710  dr_ik = 1000000.0_dp
711  END IF
712 
713  IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
714  dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
715  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
716  ELSEIF (op_jk == do_potential_coulomb) THEN
717  dr_jk = 1000000.0_dp
718  dr_ik = 1000000.0_dp
719  END IF
720 
721  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, natom=natom, &
722  dft_control=dft_control, kpoints=kpoints, para_env=para_env, cell=cell)
723 
724  IF (do_kpoints_prv) THEN
725  nimg = dft_control%nimages
726  ncell_ri = nimg
727  IF (do_hfx_kpoints_prv) THEN
728  nimg = SIZE(t3c, 1)
729  ncell_ri = SIZE(t3c, 2)
730  END IF
731  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
732  ELSE
733  nimg = 1
734  ncell_ri = 1
735  END IF
736 
737  IF (do_kpoints_prv) THEN
738  kp_index_lbounds = lbound(cell_to_index)
739  kp_index_ubounds = ubound(cell_to_index)
740  END IF
741 
742  !Do a first loop over the nl and count the blocks present
743  ALLOCATE (nblk(nimg, ncell_ri))
744  nblk(:, :) = 0
745 
746  CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
747 
748  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
749 
750  DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
751  CALL get_3c_iterator_info(nl_3c_iter, iatom=iatom, ikind=ikind, jkind=jkind, kkind=kkind, &
752  rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
753 
754  djk = norm2(rjk)
755  dij = norm2(rij)
756  dik = norm2(rik)
757 
758  IF (do_kpoints_prv) THEN
759 
760  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
761  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
762 
763  jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
764  IF (jcell > nimg .OR. jcell < 1) cycle
765 
766  IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
767  any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
768 
769  kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
770  IF (kcell > nimg .OR. kcell < 1) cycle
771  IF (do_hfx_kpoints_prv) THEN
772  IF (dik > ri_range) cycle
773  kcell = 1
774  END IF
775  ELSE
776  jcell = 1; kcell = 1
777  END IF
778 
779  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
780  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
781  CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
782 
783  IF (kind_radius_j + kind_radius_i + dr_ij < dij) cycle
784  IF (kind_radius_j + kind_radius_k + dr_jk < djk) cycle
785  IF (kind_radius_k + kind_radius_i + dr_ik < dik) cycle
786 
787  nblk(jcell, kcell) = nblk(jcell, kcell) + 1
788  END DO
789  CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
790 
791  !Do a second loop over the nl to give block indices
792  ALLOCATE (alloc_i(nimg, ncell_ri))
793  ALLOCATE (alloc_j(nimg, ncell_ri))
794  ALLOCATE (alloc_k(nimg, ncell_ri))
795  DO k_img = 1, ncell_ri
796  DO j_img = 1, nimg
797  ALLOCATE (alloc_i(j_img, k_img)%array(nblk(j_img, k_img)))
798  ALLOCATE (alloc_j(j_img, k_img)%array(nblk(j_img, k_img)))
799  ALLOCATE (alloc_k(j_img, k_img)%array(nblk(j_img, k_img)))
800  END DO
801  END DO
802  nblk(:, :) = 0
803 
804  CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
805 
806  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
807 
808  DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
809  CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
810  iatom=iatom, jatom=jatom, katom=katom, &
811  rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
812 
813  djk = norm2(rjk)
814  dij = norm2(rij)
815  dik = norm2(rik)
816 
817  IF (do_kpoints_prv) THEN
818 
819  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
820  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
821 
822  jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
823  IF (jcell > nimg .OR. jcell < 1) cycle
824 
825  IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
826  any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
827 
828  kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
829  IF (kcell > nimg .OR. kcell < 1) cycle
830  IF (do_hfx_kpoints_prv) THEN
831  IF (dik > ri_range) cycle
832  kcell = img_to_ri_cell_prv(kcell)
833  END IF
834  ELSE
835  jcell = 1; kcell = 1
836  END IF
837 
838  blk_idx = [iatom, jatom, katom]
839  IF (do_hfx_kpoints_prv) THEN
840  blk_idx(3) = (kcell - 1)*natom + katom
841  kcell = 1
842  END IF
843 
844  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
845  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
846  CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
847 
848  IF (kind_radius_j + kind_radius_i + dr_ij < dij) cycle
849  IF (kind_radius_j + kind_radius_k + dr_jk < djk) cycle
850  IF (kind_radius_k + kind_radius_i + dr_ik < dik) cycle
851 
852  nblk(jcell, kcell) = nblk(jcell, kcell) + 1
853 
854  !Note: there may be repeated indices due to periodic images => dbt_reserve_blocks takes care of it
855  alloc_i(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(1)
856  alloc_j(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(2)
857  alloc_k(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(3)
858 
859  END DO
860  CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
861 
862 !TODO: Parallelize creation of block list.
863 !$OMP PARALLEL DEFAULT(NONE) SHARED(t3c,nimg,nblk,alloc_i,alloc_j,alloc_k,ncell_RI) &
864 !$OMP PRIVATE(k_img,j_img,nblk_per_thread,A,b)
865  DO k_img = 1, ncell_ri
866  DO j_img = 1, nimg
867  IF (ALLOCATED(alloc_i(j_img, k_img)%array)) THEN
868  nblk_per_thread = nblk(j_img, k_img)/omp_get_num_threads() + 1
869  a = omp_get_thread_num()*nblk_per_thread + 1
870  b = min(a + nblk_per_thread, nblk(j_img, k_img))
871  CALL dbt_reserve_blocks(t3c(j_img, k_img), &
872  alloc_i(j_img, k_img)%array(a:b), &
873  alloc_j(j_img, k_img)%array(a:b), &
874  alloc_k(j_img, k_img)%array(a:b))
875  END IF
876  END DO
877  END DO
878 !$OMP END PARALLEL
879 
880  CALL timestop(handle)
881 
882  END SUBROUTINE
883 
884 ! **************************************************************************************************
885 !> \brief ...
886 !> \param t3c ...
887 !> \param nl_3c ...
888 !> \param basis_i ...
889 !> \param basis_j ...
890 !> \param basis_k ...
891 !> \param qs_env ...
892 !> \param potential_parameter ...
893 !> \param op_pos ...
894 !> \param do_kpoints ...
895 !> \param bounds_i ...
896 !> \param bounds_j ...
897 !> \param bounds_k ...
898 ! **************************************************************************************************
899  SUBROUTINE alloc_block_3c_old(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, op_pos, &
900  do_kpoints, bounds_i, bounds_j, bounds_k)
901  TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT) :: t3c
902  TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
903  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
904  TYPE(qs_environment_type), POINTER :: qs_env
905  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
906  INTEGER, INTENT(IN), OPTIONAL :: op_pos
907  LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints
908  INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
909 
910  CHARACTER(LEN=*), PARAMETER :: routinen = 'alloc_block_3c_old'
911 
912  INTEGER :: blk_cnt, handle, i, i_img, iatom, iblk, ikind, iproc, j_img, jatom, jcell, jkind, &
913  katom, kcell, kkind, natom, nimg, op_ij, op_jk, op_pos_prv
914  INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp
915  INTEGER, DIMENSION(3) :: cell_j, cell_k, kp_index_lbounds, &
916  kp_index_ubounds
917  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
918  LOGICAL :: do_kpoints_prv, new_block
919  REAL(kind=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
920  kind_radius_i, kind_radius_j, &
921  kind_radius_k
922  REAL(kind=dp), DIMENSION(3) :: rij, rik, rjk
923  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
924  TYPE(dft_control_type), POINTER :: dft_control
925  TYPE(kpoint_type), POINTER :: kpoints
926  TYPE(mp_para_env_type), POINTER :: para_env
927  TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
928  TYPE(one_dim_int_array), ALLOCATABLE, &
929  DIMENSION(:, :) :: alloc_i, alloc_j, alloc_k
930  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
931 
932  CALL timeset(routinen, handle)
933  NULLIFY (qs_kind_set, atomic_kind_set)
934 
935  IF (PRESENT(do_kpoints)) THEN
936  do_kpoints_prv = do_kpoints
937  ELSE
938  do_kpoints_prv = .false.
939  END IF
940 
941  dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
942 
943  op_ij = do_potential_id; op_jk = do_potential_id
944 
945  IF (PRESENT(op_pos)) THEN
946  op_pos_prv = op_pos
947  ELSE
948  op_pos_prv = 1
949  END IF
950 
951  SELECT CASE (op_pos_prv)
952  CASE (1)
953  op_ij = potential_parameter%potential_type
954  CASE (2)
955  op_jk = potential_parameter%potential_type
956  END SELECT
957 
958  IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
959  dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
960  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
961  ELSEIF (op_ij == do_potential_coulomb) THEN
962  dr_ij = 1000000.0_dp
963  dr_ik = 1000000.0_dp
964  END IF
965 
966  IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
967  dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
968  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
969  ELSEIF (op_jk == do_potential_coulomb) THEN
970  dr_jk = 1000000.0_dp
971  dr_ik = 1000000.0_dp
972  END IF
973 
974  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, natom=natom, &
975  dft_control=dft_control, kpoints=kpoints, para_env=para_env)
976 
977  IF (do_kpoints_prv) THEN
978  nimg = dft_control%nimages
979  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
980  ELSE
981  nimg = 1
982  END IF
983 
984  ALLOCATE (alloc_i(nimg, nimg))
985  ALLOCATE (alloc_j(nimg, nimg))
986  ALLOCATE (alloc_k(nimg, nimg))
987 
988  IF (do_kpoints_prv) THEN
989  kp_index_lbounds = lbound(cell_to_index)
990  kp_index_ubounds = ubound(cell_to_index)
991  END IF
992 
993  CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
994  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
995  DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
996  CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
997  iatom=iatom, jatom=jatom, katom=katom, &
998  rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
999 
1000  IF (do_kpoints_prv) THEN
1001 
1002  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
1003  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
1004 
1005  jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1006  IF (jcell > nimg .OR. jcell < 1) cycle
1007 
1008  IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
1009  any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
1010 
1011  kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
1012  IF (kcell > nimg .OR. kcell < 1) cycle
1013  ELSE
1014  jcell = 1; kcell = 1
1015  END IF
1016 
1017  djk = norm2(rjk)
1018  dij = norm2(rij)
1019  dik = norm2(rik)
1020 
1021  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
1022  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
1023  CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
1024 
1025  IF (kind_radius_j + kind_radius_i + dr_ij < dij) cycle
1026  IF (kind_radius_j + kind_radius_k + dr_jk < djk) cycle
1027  IF (kind_radius_k + kind_radius_i + dr_ik < dik) cycle
1028 
1029  ! tensor is not symmetric therefore need to allocate rows and columns in
1030  ! correspondence with neighborlist. Note that this only allocates half
1031  ! of the blocks (since neighborlist is symmetric). After filling the blocks,
1032  ! tensor will be added to its transposed
1033 
1034  associate(ai => alloc_i(jcell, kcell))
1035  associate(aj => alloc_j(jcell, kcell))
1036  associate(ak => alloc_k(jcell, kcell))
1037 
1038  new_block = .true.
1039  IF (ALLOCATED(aj%array)) THEN
1040  DO iblk = 1, SIZE(aj%array)
1041  IF (ai%array(iblk) == iatom .AND. &
1042  aj%array(iblk) == jatom .AND. &
1043  ak%array(iblk) == katom) THEN
1044  new_block = .false.
1045  EXIT
1046  END IF
1047  END DO
1048  END IF
1049  IF (.NOT. new_block) cycle
1050 
1051  IF (ALLOCATED(ai%array)) THEN
1052  blk_cnt = SIZE(ai%array)
1053  ALLOCATE (tmp(blk_cnt))
1054  tmp(:) = ai%array(:)
1055  DEALLOCATE (ai%array)
1056  ALLOCATE (ai%array(blk_cnt + 1))
1057  ai%array(1:blk_cnt) = tmp(:)
1058  ai%array(blk_cnt + 1) = iatom
1059  ELSE
1060  ALLOCATE (ai%array(1))
1061  ai%array(1) = iatom
1062  END IF
1063 
1064  IF (ALLOCATED(aj%array)) THEN
1065  tmp(:) = aj%array(:)
1066  DEALLOCATE (aj%array)
1067  ALLOCATE (aj%array(blk_cnt + 1))
1068  aj%array(1:blk_cnt) = tmp(:)
1069  aj%array(blk_cnt + 1) = jatom
1070  ELSE
1071  ALLOCATE (aj%array(1))
1072  aj%array(1) = jatom
1073  END IF
1074 
1075  IF (ALLOCATED(ak%array)) THEN
1076  tmp(:) = ak%array(:)
1077  DEALLOCATE (ak%array)
1078  ALLOCATE (ak%array(blk_cnt + 1))
1079  ak%array(1:blk_cnt) = tmp(:)
1080  ak%array(blk_cnt + 1) = katom
1081  ELSE
1082  ALLOCATE (ak%array(1))
1083  ak%array(1) = katom
1084  END IF
1085 
1086  IF (ALLOCATED(tmp)) DEALLOCATE (tmp)
1087  END associate
1088  END associate
1089  END associate
1090  END DO
1091 
1092  CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
1093 
1094  DO i_img = 1, nimg
1095  DO j_img = 1, nimg
1096  IF (ALLOCATED(alloc_i(i_img, j_img)%array)) THEN
1097  DO i = 1, SIZE(alloc_i(i_img, j_img)%array)
1098  CALL dbt_get_stored_coordinates(t3c(i_img, j_img), &
1099  [alloc_i(i_img, j_img)%array(i), alloc_j(i_img, j_img)%array(i), &
1100  alloc_k(i_img, j_img)%array(i)], &
1101  iproc)
1102  cpassert(iproc .EQ. para_env%mepos)
1103  END DO
1104 
1105  CALL dbt_reserve_blocks(t3c(i_img, j_img), &
1106  alloc_i(i_img, j_img)%array, &
1107  alloc_j(i_img, j_img)%array, &
1108  alloc_k(i_img, j_img)%array)
1109  END IF
1110  END DO
1111  END DO
1112 
1113  CALL timestop(handle)
1114 
1115  END SUBROUTINE
1116 
1117 ! **************************************************************************************************
1118 !> \brief Build 3-center derivative tensors
1119 !> \param t3c_der_i empty DBCSR tensor which will contain the 1st center derivatives
1120 !> \param t3c_der_k empty DBCSR tensor which will contain the 3rd center derivatives
1121 !> \param filter_eps Filter threshold for tensor blocks
1122 !> \param qs_env ...
1123 !> \param nl_3c 3-center neighborlist
1124 !> \param basis_i ...
1125 !> \param basis_j ...
1126 !> \param basis_k ...
1127 !> \param potential_parameter ...
1128 !> \param der_eps neglect integrals smaller than der_eps
1129 !> \param op_pos operator position.
1130 !> 1: calculate (i|jk) integrals,
1131 !> 2: calculate (ij|k) integrals
1132 !> \param do_kpoints ...
1133 !> this routine requires that libint has been static initialised somewhere else
1134 !> \param do_hfx_kpoints ...
1135 !> \param bounds_i ...
1136 !> \param bounds_j ...
1137 !> \param bounds_k ...
1138 !> \param RI_range ...
1139 !> \param img_to_RI_cell ...
1140 ! **************************************************************************************************
1141  SUBROUTINE build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, &
1142  nl_3c, basis_i, basis_j, basis_k, &
1143  potential_parameter, der_eps, &
1144  op_pos, do_kpoints, do_hfx_kpoints, &
1145  bounds_i, bounds_j, bounds_k, &
1146  RI_range, img_to_RI_cell)
1147 
1148  TYPE(dbt_type), DIMENSION(:, :, :), INTENT(INOUT) :: t3c_der_i, t3c_der_k
1149  REAL(kind=dp), INTENT(IN) :: filter_eps
1150  TYPE(qs_environment_type), POINTER :: qs_env
1151  TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
1152  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
1153  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1154  REAL(kind=dp), INTENT(IN), OPTIONAL :: der_eps
1155  INTEGER, INTENT(IN), OPTIONAL :: op_pos
1156  LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints
1157  INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
1158  REAL(dp), INTENT(IN), OPTIONAL :: ri_range
1159  INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: img_to_ri_cell
1160 
1161  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_3c_derivatives'
1162 
1163  INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
1164  block_start_k, egfi, handle, handle2, i, i_img, i_xyz, iatom, ibasis, ikind, ilist, imax, &
1165  iset, j_img, jatom, jcell, jkind, jset, katom, kcell, kkind, kset, m_max, max_ncoi, &
1166  max_ncoj, max_ncok, max_nset, max_nsgfi, max_nsgfj, max_nsgfk, maxli, maxlj, maxlk, &
1167  mepos, natom, nbasis, ncell_ri, ncoi, ncoj, ncok, nimg, nseti, nsetj, nsetk, nthread, &
1168  op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id
1169  INTEGER, ALLOCATABLE, DIMENSION(:) :: img_to_ri_cell_prv
1170  INTEGER, DIMENSION(2) :: bo
1171  INTEGER, DIMENSION(3) :: blk_idx, blk_size, cell_j, cell_k, &
1172  kp_index_lbounds, kp_index_ubounds, sp
1173  INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
1174  lmin_k, npgfi, npgfj, npgfk, nsgfi, &
1175  nsgfj, nsgfk
1176  INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
1177  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1178  LOGICAL :: do_hfx_kpoints_prv, do_kpoints_prv, &
1179  found, skip
1180  LOGICAL, DIMENSION(3) :: block_j_not_zero, block_k_not_zero, &
1181  der_j_zero, der_k_zero
1182  REAL(dp), DIMENSION(3) :: der_ext_i, der_ext_j, der_ext_k
1183  REAL(kind=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
1184  kind_radius_i, kind_radius_j, &
1185  kind_radius_k, prefac
1186  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
1187  max_contraction_i, max_contraction_j, &
1188  max_contraction_k
1189  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dijk_contr, dummy_block_t
1190  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: block_t_i, block_t_j, block_t_k, dijk_j, &
1191  dijk_k, tmp_ijk_i, tmp_ijk_j
1192  REAL(kind=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
1193  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
1194  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
1195  sphi_k, zeti, zetj, zetk
1196  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1197  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
1198  TYPE(cp_libint_t) :: lib
1199  TYPE(dbt_type) :: t3c_tmp
1200  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t3c_template
1201  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :, :) :: t3c_der_j
1202  TYPE(dft_control_type), POINTER :: dft_control
1203  TYPE(gto_basis_set_type), POINTER :: basis_set
1204  TYPE(kpoint_type), POINTER :: kpoints
1205  TYPE(mp_para_env_type), POINTER :: para_env
1206  TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
1207  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1208 
1209  CALL timeset(routinen, handle)
1210 
1211  IF (PRESENT(do_kpoints)) THEN
1212  do_kpoints_prv = do_kpoints
1213  ELSE
1214  do_kpoints_prv = .false.
1215  END IF
1216 
1217  IF (PRESENT(do_hfx_kpoints)) THEN
1218  do_hfx_kpoints_prv = do_hfx_kpoints
1219  ELSE
1220  do_hfx_kpoints_prv = .false.
1221  END IF
1222  IF (do_hfx_kpoints_prv) THEN
1223  cpassert(do_kpoints_prv)
1224  cpassert(PRESENT(ri_range))
1225  cpassert(PRESENT(img_to_ri_cell))
1226  END IF
1227 
1228  IF (PRESENT(img_to_ri_cell)) THEN
1229  ALLOCATE (img_to_ri_cell_prv(SIZE(img_to_ri_cell)))
1230  img_to_ri_cell_prv(:) = img_to_ri_cell
1231  END IF
1232 
1233  op_ij = do_potential_id; op_jk = do_potential_id
1234 
1235  IF (PRESENT(op_pos)) THEN
1236  op_pos_prv = op_pos
1237  ELSE
1238  op_pos_prv = 1
1239  END IF
1240 
1241  SELECT CASE (op_pos_prv)
1242  CASE (1)
1243  op_ij = potential_parameter%potential_type
1244  CASE (2)
1245  op_jk = potential_parameter%potential_type
1246  END SELECT
1247 
1248  dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
1249 
1250  IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
1251  dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
1252  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1253  ELSEIF (op_ij == do_potential_coulomb) THEN
1254  dr_ij = 1000000.0_dp
1255  dr_ik = 1000000.0_dp
1256  END IF
1257 
1258  IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
1259  dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
1260  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1261  ELSEIF (op_jk == do_potential_coulomb) THEN
1262  dr_jk = 1000000.0_dp
1263  dr_ik = 1000000.0_dp
1264  END IF
1265 
1266  NULLIFY (qs_kind_set, atomic_kind_set)
1267 
1268  ! get stuff
1269  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1270  natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
1271 
1272  IF (do_kpoints_prv) THEN
1273  nimg = dft_control%nimages
1274  ncell_ri = nimg
1275  IF (do_hfx_kpoints_prv) THEN
1276  nimg = SIZE(t3c_der_k, 1)
1277  ncell_ri = SIZE(t3c_der_k, 2)
1278  END IF
1279  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
1280  ELSE
1281  nimg = 1
1282  ncell_ri = 1
1283  END IF
1284 
1285  IF (do_hfx_kpoints_prv) THEN
1286  cpassert(op_pos_prv == 2)
1287  ELSE
1288  cpassert(all(shape(t3c_der_k) == [nimg, ncell_ri, 3]))
1289  END IF
1290 
1291  ALLOCATE (t3c_template(nimg, ncell_ri))
1292  DO j_img = 1, ncell_ri
1293  DO i_img = 1, nimg
1294  CALL dbt_create(t3c_der_i(i_img, j_img, 1), t3c_template(i_img, j_img))
1295  END DO
1296  END DO
1297 
1298  CALL alloc_block_3c(t3c_template, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, &
1299  op_pos=op_pos_prv, do_kpoints=do_kpoints, do_hfx_kpoints=do_hfx_kpoints, &
1300  bounds_i=bounds_i, bounds_j=bounds_j, bounds_k=bounds_k, &
1301  ri_range=ri_range, img_to_ri_cell=img_to_ri_cell)
1302  DO i_xyz = 1, 3
1303  DO j_img = 1, ncell_ri
1304  DO i_img = 1, nimg
1305  CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_i(i_img, j_img, i_xyz))
1306  CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_k(i_img, j_img, i_xyz))
1307  END DO
1308  END DO
1309  END DO
1310 
1311  DO j_img = 1, ncell_ri
1312  DO i_img = 1, nimg
1313  CALL dbt_destroy(t3c_template(i_img, j_img))
1314  END DO
1315  END DO
1316  DEALLOCATE (t3c_template)
1317 
1318  IF (nl_3c%sym == symmetric_jk) THEN
1319  ALLOCATE (t3c_der_j(nimg, ncell_ri, 3))
1320  DO i_xyz = 1, 3
1321  DO j_img = 1, ncell_ri
1322  DO i_img = 1, nimg
1323  CALL dbt_create(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
1324  CALL dbt_copy(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
1325  END DO
1326  END DO
1327  END DO
1328  END IF
1329 
1330  !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
1331  nbasis = SIZE(basis_i)
1332  max_nsgfi = 0
1333  max_ncoi = 0
1334  max_nset = 0
1335  maxli = 0
1336  DO ibasis = 1, nbasis
1337  CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
1338  nset=iset, nsgf_set=nsgfi, npgf=npgfi)
1339  maxli = max(maxli, imax)
1340  max_nset = max(max_nset, iset)
1341  max_nsgfi = max(max_nsgfi, maxval(nsgfi))
1342  max_ncoi = max(max_ncoi, maxval(npgfi)*ncoset(maxli))
1343  END DO
1344  max_nsgfj = 0
1345  max_ncoj = 0
1346  maxlj = 0
1347  DO ibasis = 1, nbasis
1348  CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
1349  nset=jset, nsgf_set=nsgfj, npgf=npgfj)
1350  maxlj = max(maxlj, imax)
1351  max_nset = max(max_nset, jset)
1352  max_nsgfj = max(max_nsgfj, maxval(nsgfj))
1353  max_ncoj = max(max_ncoj, maxval(npgfj)*ncoset(maxlj))
1354  END DO
1355  max_nsgfk = 0
1356  max_ncok = 0
1357  maxlk = 0
1358  DO ibasis = 1, nbasis
1359  CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
1360  nset=kset, nsgf_set=nsgfk, npgf=npgfk)
1361  maxlk = max(maxlk, imax)
1362  max_nset = max(max_nset, kset)
1363  max_nsgfk = max(max_nsgfk, maxval(nsgfk))
1364  max_ncok = max(max_ncok, maxval(npgfk)*ncoset(maxlk))
1365  END DO
1366  m_max = maxli + maxlj + maxlk + 1
1367 
1368  !To minimize expensive memory opsand generally optimize contraction, pre-allocate
1369  !contiguous sphi arrays (and transposed in the cas of sphi_i)
1370 
1371  NULLIFY (tspj, spi, spk)
1372  ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
1373 
1374  DO ibasis = 1, nbasis
1375  DO iset = 1, max_nset
1376  NULLIFY (spi(iset, ibasis)%array)
1377  NULLIFY (tspj(iset, ibasis)%array)
1378 
1379  NULLIFY (spk(iset, ibasis)%array)
1380  END DO
1381  END DO
1382 
1383  DO ilist = 1, 3
1384  DO ibasis = 1, nbasis
1385  IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
1386  IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
1387  IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
1388 
1389  DO iset = 1, basis_set%nset
1390 
1391  ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
1392  sgfi = basis_set%first_sgf(1, iset)
1393  egfi = sgfi + basis_set%nsgf_set(iset) - 1
1394 
1395  IF (ilist == 1) THEN
1396  ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1397  spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1398 
1399  ELSE IF (ilist == 2) THEN
1400  ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
1401  tspj(iset, ibasis)%array(:, :) = transpose(basis_set%sphi(1:ncoi, sgfi:egfi))
1402 
1403  ELSE
1404  ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1405  spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1406  END IF
1407 
1408  END DO !iset
1409  END DO !ibasis
1410  END DO !ilist
1411 
1412  !Init the truncated Coulomb operator
1413  IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN
1414 
1415  IF (m_max > get_lmax_init()) THEN
1416  IF (para_env%mepos == 0) THEN
1417  CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
1418  END IF
1419  CALL init(m_max, unit_id, para_env%mepos, para_env)
1420  IF (para_env%mepos == 0) THEN
1421  CALL close_file(unit_id)
1422  END IF
1423  END IF
1424  END IF
1425 
1426  CALL init_md_ftable(nmax=m_max)
1427 
1428  IF (do_kpoints_prv) THEN
1429  kp_index_lbounds = lbound(cell_to_index)
1430  kp_index_ubounds = ubound(cell_to_index)
1431  END IF
1432 
1433  nthread = 1
1434 !$ nthread = omp_get_max_threads()
1435 
1436 !$OMP PARALLEL DEFAULT(NONE) &
1437 !$OMP SHARED (nthread,do_kpoints_prv,kp_index_lbounds,kp_index_ubounds,maxli,maxlk,maxlj,bounds_i,&
1438 !$OMP bounds_j,bounds_k,nimg,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,op_pos_prv,&
1439 !$OMP potential_parameter,der_eps,tspj,spi,spk,cell_to_index,max_ncoi,max_nsgfk,RI_range,&
1440 !$OMP max_nsgfj,max_ncok,natom,nl_3c,t3c_der_i,t3c_der_k,t3c_der_j,ncell_RI,&
1441 !$OMP img_to_RI_cell_prv,do_hfx_kpoints_prv) &
1442 !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,cell_j,cell_k,&
1443 !$OMP prefac,jcell,kcell,first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
1444 !$OMP sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
1445 !$OMP set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
1446 !$OMP rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
1447 !$OMP sgfk,dijk_j,dijk_k,ri,rj,rk,max_contraction_i,max_contraction_j,blk_idx,&
1448 !$OMP max_contraction_k,iset,jset,kset,block_t_i,blk_size,dijk_contr,cpp_buffer,ccp_buffer,&
1449 !$OMP block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
1450 !$OMP dummy_block_t,sp,handle2,mepos,bo,block_t_k,der_ext_i,der_ext_j,der_ext_k,tmp_ijk_i,&
1451 !$OMP block_k_not_zero,der_k_zero,skip,der_j_zero,block_t_j,block_j_not_zero,tmp_ijk_j,i)
1452 
1453  mepos = 0
1454 !$ mepos = omp_get_thread_num()
1455 
1456  CALL cp_libint_init_3eri1(lib, max(maxli, maxlj, maxlk))
1457  CALL cp_libint_set_contrdepth(lib, 1)
1458 
1459  !pre-allocate contraction buffers
1460  ALLOCATE (cpp_buffer(max_nsgfj*max_ncok), ccp_buffer(max_nsgfj*max_nsgfk*max_ncoi))
1461 
1462  CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
1463 
1464  !We split the provided bounds among the threads such that each threads works on a different set of atoms
1465  IF (PRESENT(bounds_i)) THEN
1466  bo = get_limit(bounds_i(2) - bounds_i(1) + 1, nthread, mepos)
1467  bo(:) = bo(:) + bounds_i(1) - 1
1468  CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
1469  ELSE IF (PRESENT(bounds_j)) THEN
1470  bo = get_limit(bounds_j(2) - bounds_j(1) + 1, nthread, mepos)
1471  bo(:) = bo(:) + bounds_j(1) - 1
1472  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bo, bounds_k)
1473  ELSE IF (PRESENT(bounds_k)) THEN
1474  bo = get_limit(bounds_k(2) - bounds_k(1) + 1, nthread, mepos)
1475  bo(:) = bo(:) + bounds_k(1) - 1
1476  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bo)
1477  ELSE
1478  bo = get_limit(natom, nthread, mepos)
1479  CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
1480  END IF
1481 
1482  skip = .false.
1483  IF (bo(1) > bo(2)) skip = .true.
1484 
1485  DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
1486  CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
1487  iatom=iatom, jatom=jatom, katom=katom, &
1488  rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
1489  IF (skip) EXIT
1490 
1491  djk = norm2(rjk)
1492  dij = norm2(rij)
1493  dik = norm2(rik)
1494 
1495  IF (do_kpoints_prv) THEN
1496  prefac = 0.5_dp
1497  ELSEIF (nl_3c%sym == symmetric_jk) THEN
1498  IF (jatom == katom) THEN
1499  prefac = 0.5_dp
1500  ELSE
1501  prefac = 1.0_dp
1502  END IF
1503  ELSE
1504  prefac = 1.0_dp
1505  END IF
1506  IF (do_hfx_kpoints_prv) prefac = 1.0_dp
1507 
1508  IF (do_kpoints_prv) THEN
1509 
1510  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
1511  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
1512 
1513  jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1514  IF (jcell > nimg .OR. jcell < 1) cycle
1515 
1516  IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
1517  any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
1518 
1519  kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
1520  IF (kcell > nimg .OR. kcell < 1) cycle
1521  !In case of HFX k-points, we only consider P^k if d_ik <= RI_range
1522  IF (do_hfx_kpoints_prv) THEN
1523  IF (dik > ri_range) cycle
1524  kcell = img_to_ri_cell_prv(kcell)
1525  END IF
1526  ELSE
1527  jcell = 1; kcell = 1
1528  END IF
1529 
1530  blk_idx = [iatom, jatom, katom]
1531  IF (do_hfx_kpoints_prv) THEN
1532  blk_idx(3) = (kcell - 1)*natom + katom
1533  kcell = 1
1534  END IF
1535 
1536  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
1537  npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
1538  sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
1539 
1540  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
1541  npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
1542  sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
1543 
1544  CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
1545  npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
1546  sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
1547 
1548  IF (kind_radius_j + kind_radius_i + dr_ij < dij) cycle
1549  IF (kind_radius_j + kind_radius_k + dr_jk < djk) cycle
1550  IF (kind_radius_k + kind_radius_i + dr_ik < dik) cycle
1551 
1552  ALLOCATE (max_contraction_i(nseti))
1553  max_contraction_i = 0.0_dp
1554  DO iset = 1, nseti
1555  sgfi = first_sgf_i(1, iset)
1556  max_contraction_i(iset) = maxval((/(sum(abs(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/))
1557  END DO
1558 
1559  ALLOCATE (max_contraction_j(nsetj))
1560  max_contraction_j = 0.0_dp
1561  DO jset = 1, nsetj
1562  sgfj = first_sgf_j(1, jset)
1563  max_contraction_j(jset) = maxval((/(sum(abs(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/))
1564  END DO
1565 
1566  ALLOCATE (max_contraction_k(nsetk))
1567  max_contraction_k = 0.0_dp
1568  DO kset = 1, nsetk
1569  sgfk = first_sgf_k(1, kset)
1570  max_contraction_k(kset) = maxval((/(sum(abs(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/))
1571  END DO
1572 
1573  CALL dbt_blk_sizes(t3c_der_i(jcell, kcell, 1), blk_idx, blk_size)
1574 
1575  ALLOCATE (block_t_i(blk_size(2), blk_size(3), blk_size(1), 3))
1576  ALLOCATE (block_t_j(blk_size(2), blk_size(3), blk_size(1), 3))
1577  ALLOCATE (block_t_k(blk_size(2), blk_size(3), blk_size(1), 3))
1578 
1579  block_t_i = 0.0_dp
1580  block_t_j = 0.0_dp
1581  block_t_k = 0.0_dp
1582  block_j_not_zero = .false.
1583  block_k_not_zero = .false.
1584 
1585  DO iset = 1, nseti
1586 
1587  DO jset = 1, nsetj
1588 
1589  IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) cycle
1590 
1591  DO kset = 1, nsetk
1592 
1593  IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) cycle
1594  IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) cycle
1595 
1596  ncoi = npgfi(iset)*ncoset(lmax_i(iset))
1597  ncoj = npgfj(jset)*ncoset(lmax_j(jset))
1598  ncok = npgfk(kset)*ncoset(lmax_k(kset))
1599 
1600  sgfi = first_sgf_i(1, iset)
1601  sgfj = first_sgf_j(1, jset)
1602  sgfk = first_sgf_k(1, kset)
1603 
1604  IF (ncoj*ncok*ncoi > 0) THEN
1605  ALLOCATE (dijk_j(ncoj, ncok, ncoi, 3))
1606  ALLOCATE (dijk_k(ncoj, ncok, ncoi, 3))
1607  dijk_j(:, :, :, :) = 0.0_dp
1608  dijk_k(:, :, :, :) = 0.0_dp
1609 
1610  der_j_zero = .false.
1611  der_k_zero = .false.
1612 
1613  !need positions for libint. Only relative positions are needed => set ri to 0.0
1614  ri = 0.0_dp
1615  rj = rij ! ri + rij
1616  rk = rik ! ri + rik
1617 
1618  IF (op_pos_prv == 1) THEN
1619  CALL eri_3center_derivs(dijk_j, dijk_k, &
1620  lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
1621  lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
1622  lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
1623  djk, dij, dik, lib, potential_parameter, &
1624  der_abc_1_ext=der_ext_j, der_abc_2_ext=der_ext_k)
1625  ELSE
1626  ALLOCATE (tmp_ijk_i(ncoi, ncoj, ncok, 3))
1627  ALLOCATE (tmp_ijk_j(ncoi, ncoj, ncok, 3))
1628  tmp_ijk_i(:, :, :, :) = 0.0_dp
1629  tmp_ijk_j(:, :, :, :) = 0.0_dp
1630 
1631  CALL eri_3center_derivs(tmp_ijk_i, tmp_ijk_j, &
1632  lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
1633  lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
1634  lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
1635  dij, dik, djk, lib, potential_parameter, &
1636  der_abc_1_ext=der_ext_i, der_abc_2_ext=der_ext_j)
1637 
1638  !TODO: is that inefficient?
1639  der_ext_k = 0
1640  DO i_xyz = 1, 3
1641  DO i = 1, ncoi
1642  dijk_j(:, :, i, i_xyz) = tmp_ijk_j(i, :, :, i_xyz)
1643  dijk_k(:, :, i, i_xyz) = -(dijk_j(:, :, i, i_xyz) + tmp_ijk_i(i, :, :, i_xyz))
1644  der_ext_k(i_xyz) = max(der_ext_k(i_xyz), maxval(abs(dijk_k(:, :, i, i_xyz))))
1645  END DO
1646  END DO
1647  DEALLOCATE (tmp_ijk_i, tmp_ijk_j)
1648 
1649  END IF
1650 
1651  IF (PRESENT(der_eps)) THEN
1652  DO i_xyz = 1, 3
1653  IF (der_eps > der_ext_j(i_xyz)*(max_contraction_i(iset)* &
1654  max_contraction_j(jset)* &
1655  max_contraction_k(kset))) THEN
1656  der_j_zero(i_xyz) = .true.
1657  END IF
1658  END DO
1659 
1660  DO i_xyz = 1, 3
1661  IF (der_eps > der_ext_k(i_xyz)*(max_contraction_i(iset)* &
1662  max_contraction_j(jset)* &
1663  max_contraction_k(kset))) THEN
1664  der_k_zero(i_xyz) = .true.
1665  END IF
1666  END DO
1667  IF (all(der_j_zero) .AND. all(der_k_zero)) THEN
1668  DEALLOCATE (dijk_j, dijk_k)
1669  cycle
1670  END IF
1671  END IF
1672 
1673  ALLOCATE (dijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
1674 
1675  block_start_j = sgfj
1676  block_end_j = sgfj + nsgfj(jset) - 1
1677  block_start_k = sgfk
1678  block_end_k = sgfk + nsgfk(kset) - 1
1679  block_start_i = sgfi
1680  block_end_i = sgfi + nsgfi(iset) - 1
1681 
1682  DO i_xyz = 1, 3
1683  IF (der_j_zero(i_xyz)) cycle
1684 
1685  block_j_not_zero(i_xyz) = .true.
1686  CALL libxsmm_abc_contract(dijk_contr, dijk_j(:, :, :, i_xyz), tspj(jset, jkind)%array, &
1687  spk(kset, kkind)%array, spi(iset, ikind)%array, &
1688  ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
1689  nsgfi(iset), cpp_buffer, ccp_buffer)
1690 
1691  block_t_j(block_start_j:block_end_j, &
1692  block_start_k:block_end_k, &
1693  block_start_i:block_end_i, i_xyz) = &
1694  block_t_j(block_start_j:block_end_j, &
1695  block_start_k:block_end_k, &
1696  block_start_i:block_end_i, i_xyz) + &
1697  prefac*dijk_contr(:, :, :)
1698 
1699  END DO
1700 
1701  DO i_xyz = 1, 3
1702  IF (der_k_zero(i_xyz)) cycle
1703 
1704  block_k_not_zero(i_xyz) = .true.
1705  CALL libxsmm_abc_contract(dijk_contr, dijk_k(:, :, :, i_xyz), tspj(jset, jkind)%array, &
1706  spk(kset, kkind)%array, spi(iset, ikind)%array, &
1707  ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
1708  nsgfi(iset), cpp_buffer, ccp_buffer)
1709 
1710  block_t_k(block_start_j:block_end_j, &
1711  block_start_k:block_end_k, &
1712  block_start_i:block_end_i, i_xyz) = &
1713  block_t_k(block_start_j:block_end_j, &
1714  block_start_k:block_end_k, &
1715  block_start_i:block_end_i, i_xyz) + &
1716  prefac*dijk_contr(:, :, :)
1717 
1718  END DO
1719 
1720  DEALLOCATE (dijk_j, dijk_k, dijk_contr)
1721  END IF ! number of triples > 0
1722  END DO
1723  END DO
1724  END DO
1725 
1726  CALL timeset(routinen//"_put_dbcsr", handle2)
1727 !$OMP CRITICAL
1728  sp = shape(block_t_i(:, :, :, 1))
1729  sp([2, 3, 1]) = sp
1730 
1731  DO i_xyz = 1, 3
1732  !Derivatives wrt to center i can be obtained by translational invariance
1733  IF ((.NOT. block_j_not_zero(i_xyz)) .AND. (.NOT. block_k_not_zero(i_xyz))) cycle
1734  block_t_i(:, :, :, i_xyz) = -(block_t_j(:, :, :, i_xyz) + block_t_k(:, :, :, i_xyz))
1735 
1736  CALL dbt_put_block(t3c_der_i(jcell, kcell, i_xyz), blk_idx, sp, &
1737  reshape(block_t_i(:, :, :, i_xyz), shape=sp, order=[2, 3, 1]), &
1738  summation=.true.)
1739  END DO
1740 !$OMP END CRITICAL
1741 !$OMP CRITICAL
1742  IF (nl_3c%sym == symmetric_jk) THEN
1743 
1744  sp = shape(block_t_j(:, :, :, 1))
1745  sp([2, 3, 1]) = sp
1746 
1747  DO i_xyz = 1, 3
1748  IF (.NOT. block_j_not_zero(i_xyz)) cycle
1749  CALL dbt_put_block(t3c_der_j(jcell, kcell, i_xyz), blk_idx, sp, &
1750  reshape(block_t_j(:, :, :, i_xyz), shape=sp, order=[2, 3, 1]), &
1751  summation=.true.)
1752  END DO
1753  END IF
1754 !$OMP END CRITICAL
1755 !$OMP CRITICAL
1756  sp = shape(block_t_k(:, :, :, 1))
1757  sp([2, 3, 1]) = sp
1758 
1759  DO i_xyz = 1, 3
1760  IF (.NOT. block_k_not_zero(i_xyz)) cycle
1761  CALL dbt_put_block(t3c_der_k(jcell, kcell, i_xyz), blk_idx, sp, &
1762  reshape(block_t_k(:, :, :, i_xyz), shape=sp, order=[2, 3, 1]), &
1763  summation=.true.)
1764  END DO
1765 !$OMP END CRITICAL
1766 
1767  CALL timestop(handle2)
1768 
1769  DEALLOCATE (block_t_i)
1770  DEALLOCATE (block_t_j)
1771  DEALLOCATE (block_t_k)
1772 
1773  DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k)
1774  END DO
1775 
1776  CALL cp_libint_cleanup_3eri1(lib)
1777 
1778  CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
1779 !$OMP END PARALLEL
1780 
1781  IF (do_kpoints_prv .AND. .NOT. do_hfx_kpoints_prv) THEN
1782  DO i_xyz = 1, 3
1783  DO kcell = 1, nimg
1784  DO jcell = 1, nimg
1785  ! need half of filter eps because afterwards we add transposed tensor
1786  CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps/2)
1787  CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps/2)
1788  END DO
1789  END DO
1790  END DO
1791 
1792  ELSEIF (nl_3c%sym == symmetric_jk) THEN
1793  !Add the transpose of t3c_der_j to t3c_der_k to get the fully populated tensor
1794  CALL dbt_create(t3c_der_k(1, 1, 1), t3c_tmp)
1795  DO i_xyz = 1, 3
1796  DO kcell = 1, nimg
1797  DO jcell = 1, nimg
1798  CALL dbt_copy(t3c_der_j(jcell, kcell, i_xyz), t3c_der_k(jcell, kcell, i_xyz), &
1799  order=[1, 3, 2], move_data=.true., summation=.true.)
1800  CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
1801 
1802  CALL dbt_copy(t3c_der_i(jcell, kcell, i_xyz), t3c_tmp)
1803  CALL dbt_copy(t3c_tmp, t3c_der_i(jcell, kcell, i_xyz), &
1804  order=[1, 3, 2], move_data=.true., summation=.true.)
1805  CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
1806  END DO
1807  END DO
1808  END DO
1809  CALL dbt_destroy(t3c_tmp)
1810 
1811  ELSEIF (nl_3c%sym == symmetric_none) THEN
1812  DO i_xyz = 1, 3
1813  DO kcell = 1, ncell_ri
1814  DO jcell = 1, nimg
1815  CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
1816  CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
1817  END DO
1818  END DO
1819  END DO
1820  ELSE
1821  cpabort("requested symmetric case not implemented")
1822  END IF
1823 
1824  IF (nl_3c%sym == symmetric_jk) THEN
1825  DO i_xyz = 1, 3
1826  DO j_img = 1, nimg
1827  DO i_img = 1, nimg
1828  CALL dbt_destroy(t3c_der_j(i_img, j_img, i_xyz))
1829  END DO
1830  END DO
1831  END DO
1832  END IF
1833 
1834  DO iset = 1, max_nset
1835  DO ibasis = 1, nbasis
1836  IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
1837  IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
1838 
1839  IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
1840  END DO
1841  END DO
1842 
1843  DEALLOCATE (spi, tspj, spk)
1844 
1845  CALL timestop(handle)
1846  END SUBROUTINE build_3c_derivatives
1847 
1848 ! **************************************************************************************************
1849 !> \brief Calculates the 3c virial contributions on the fly
1850 !> \param work_virial ...
1851 !> \param t3c_trace the tensor with which the trace should be taken
1852 !> \param pref ...
1853 !> \param qs_env ...
1854 !> \param nl_3c 3-center neighborlist, with a distribution matching that of t3c_trace
1855 !> \param basis_i ...
1856 !> \param basis_j ...
1857 !> \param basis_k ...
1858 !> \param potential_parameter ...
1859 !> \param der_eps neglect integrals smaller than der_eps
1860 !> \param op_pos operator position.
1861 !> 1: calculate (i|jk) integrals,
1862 !> 2: calculate (ij|k) integrals
1863 !> this routine requires that libint has been static initialised somewhere else
1864 ! **************************************************************************************************
1865  SUBROUTINE calc_3c_virial(work_virial, t3c_trace, pref, qs_env, &
1866  nl_3c, basis_i, basis_j, basis_k, &
1867  potential_parameter, der_eps, op_pos)
1868 
1869  REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: work_virial
1870  TYPE(dbt_type), INTENT(INOUT) :: t3c_trace
1871  REAL(kind=dp), INTENT(IN) :: pref
1872  TYPE(qs_environment_type), POINTER :: qs_env
1873  TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
1874  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
1875  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1876  REAL(kind=dp), INTENT(IN), OPTIONAL :: der_eps
1877  INTEGER, INTENT(IN), OPTIONAL :: op_pos
1878 
1879  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_3c_virial'
1880 
1881  INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
1882  block_start_k, egfi, handle, i, i_xyz, iatom, ibasis, ikind, ilist, imax, iset, j_xyz, &
1883  jatom, jkind, jset, katom, kkind, kset, m_max, max_ncoi, max_ncoj, max_ncok, max_nset, &
1884  max_nsgfi, max_nsgfj, max_nsgfk, maxli, maxlj, maxlk, mepos, natom, nbasis, ncoi, ncoj, &
1885  ncok, nseti, nsetj, nsetk, nthread, op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id
1886  INTEGER, DIMENSION(2) :: bo
1887  INTEGER, DIMENSION(3) :: blk_size, sp
1888  INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
1889  lmin_k, npgfi, npgfj, npgfk, nsgfi, &
1890  nsgfj, nsgfk
1891  INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
1892  LOGICAL :: found, skip
1893  LOGICAL, DIMENSION(3) :: block_j_not_zero, block_k_not_zero, &
1894  der_j_zero, der_k_zero
1895  REAL(dp) :: force
1896  REAL(dp), DIMENSION(3) :: der_ext_i, der_ext_j, der_ext_k
1897  REAL(kind=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
1898  kind_radius_i, kind_radius_j, &
1899  kind_radius_k
1900  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
1901  max_contraction_i, max_contraction_j, &
1902  max_contraction_k
1903  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ablock, dijk_contr, tmp_block
1904  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: block_t_i, block_t_j, block_t_k, dijk_j, &
1905  dijk_k
1906  REAL(kind=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk, scoord
1907  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
1908  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
1909  sphi_k, zeti, zetj, zetk
1910  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1911  TYPE(cell_type), POINTER :: cell
1912  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
1913  TYPE(cp_libint_t) :: lib
1914  TYPE(dft_control_type), POINTER :: dft_control
1915  TYPE(gto_basis_set_type), POINTER :: basis_set
1916  TYPE(mp_para_env_type), POINTER :: para_env
1917  TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
1918  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1919  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1920 
1921  CALL timeset(routinen, handle)
1922 
1923  op_ij = do_potential_id; op_jk = do_potential_id
1924 
1925  IF (PRESENT(op_pos)) THEN
1926  op_pos_prv = op_pos
1927  ELSE
1928  op_pos_prv = 1
1929  END IF
1930  cpassert(op_pos == 1)
1931  cpassert(.NOT. nl_3c%sym == symmetric_jk)
1932 
1933  SELECT CASE (op_pos_prv)
1934  CASE (1)
1935  op_ij = potential_parameter%potential_type
1936  CASE (2)
1937  op_jk = potential_parameter%potential_type
1938  END SELECT
1939 
1940  dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
1941 
1942  IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
1943  dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
1944  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1945  ELSEIF (op_ij == do_potential_coulomb) THEN
1946  dr_ij = 1000000.0_dp
1947  dr_ik = 1000000.0_dp
1948  END IF
1949 
1950  IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
1951  dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
1952  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1953  ELSEIF (op_jk == do_potential_coulomb) THEN
1954  dr_jk = 1000000.0_dp
1955  dr_ik = 1000000.0_dp
1956  END IF
1957 
1958  NULLIFY (qs_kind_set, atomic_kind_set)
1959 
1960  ! get stuff
1961  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1962  natom=natom, dft_control=dft_control, para_env=para_env, &
1963  particle_set=particle_set, cell=cell)
1964 
1965  !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
1966  nbasis = SIZE(basis_i)
1967  max_nsgfi = 0
1968  max_ncoi = 0
1969  max_nset = 0
1970  maxli = 0
1971  DO ibasis = 1, nbasis
1972  CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
1973  nset=iset, nsgf_set=nsgfi, npgf=npgfi)
1974  maxli = max(maxli, imax)
1975  max_nset = max(max_nset, iset)
1976  max_nsgfi = max(max_nsgfi, maxval(nsgfi))
1977  max_ncoi = max(max_ncoi, maxval(npgfi)*ncoset(maxli))
1978  END DO
1979  max_nsgfj = 0
1980  max_ncoj = 0
1981  maxlj = 0
1982  DO ibasis = 1, nbasis
1983  CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
1984  nset=jset, nsgf_set=nsgfj, npgf=npgfj)
1985  maxlj = max(maxlj, imax)
1986  max_nset = max(max_nset, jset)
1987  max_nsgfj = max(max_nsgfj, maxval(nsgfj))
1988  max_ncoj = max(max_ncoj, maxval(npgfj)*ncoset(maxlj))
1989  END DO
1990  max_nsgfk = 0
1991  max_ncok = 0
1992  maxlk = 0
1993  DO ibasis = 1, nbasis
1994  CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
1995  nset=kset, nsgf_set=nsgfk, npgf=npgfk)
1996  maxlk = max(maxlk, imax)
1997  max_nset = max(max_nset, kset)
1998  max_nsgfk = max(max_nsgfk, maxval(nsgfk))
1999  max_ncok = max(max_ncok, maxval(npgfk)*ncoset(maxlk))
2000  END DO
2001  m_max = maxli + maxlj + maxlk + 1
2002 
2003  !To minimize expensive memory opsand generally optimize contraction, pre-allocate
2004  !contiguous sphi arrays (and transposed in the cas of sphi_i)
2005 
2006  NULLIFY (tspj, spi, spk)
2007  ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
2008 
2009  DO ibasis = 1, nbasis
2010  DO iset = 1, max_nset
2011  NULLIFY (spi(iset, ibasis)%array)
2012  NULLIFY (tspj(iset, ibasis)%array)
2013 
2014  NULLIFY (spk(iset, ibasis)%array)
2015  END DO
2016  END DO
2017 
2018  DO ilist = 1, 3
2019  DO ibasis = 1, nbasis
2020  IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
2021  IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
2022  IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
2023 
2024  DO iset = 1, basis_set%nset
2025 
2026  ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
2027  sgfi = basis_set%first_sgf(1, iset)
2028  egfi = sgfi + basis_set%nsgf_set(iset) - 1
2029 
2030  IF (ilist == 1) THEN
2031  ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2032  spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2033 
2034  ELSE IF (ilist == 2) THEN
2035  ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
2036  tspj(iset, ibasis)%array(:, :) = transpose(basis_set%sphi(1:ncoi, sgfi:egfi))
2037 
2038  ELSE
2039  ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2040  spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2041  END IF
2042 
2043  END DO !iset
2044  END DO !ibasis
2045  END DO !ilist
2046 
2047  !Init the truncated Coulomb operator
2048  IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN
2049 
2050  IF (m_max > get_lmax_init()) THEN
2051  IF (para_env%mepos == 0) THEN
2052  CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
2053  END IF
2054  CALL init(m_max, unit_id, para_env%mepos, para_env)
2055  IF (para_env%mepos == 0) THEN
2056  CALL close_file(unit_id)
2057  END IF
2058  END IF
2059  END IF
2060 
2061  CALL init_md_ftable(nmax=m_max)
2062 
2063  nthread = 1
2064 !$ nthread = omp_get_max_threads()
2065 
2066 !$OMP PARALLEL DEFAULT(NONE) &
2067 !$OMP SHARED (nthread,maxli,maxlk,maxlj,i,work_virial,pref,&
2068 !$OMP basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,&
2069 !$OMP potential_parameter,der_eps,tspj,spi,spk,max_ncoi,max_nsgfk,&
2070 !$OMP max_nsgfj,max_ncok,natom,nl_3c,t3c_trace,cell,particle_set) &
2071 !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,i_xyz,j_xyz,&
2072 !$OMP first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
2073 !$OMP sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
2074 !$OMP set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
2075 !$OMP rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
2076 !$OMP sgfk,dijk_j,dijk_k,ri,rj,rk,max_contraction_i,max_contraction_j,tmp_block,&
2077 !$OMP max_contraction_k,iset,jset,kset,block_t_i,blk_size,dijk_contr,cpp_buffer,ccp_buffer,&
2078 !$OMP block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
2079 !$OMP sp,mepos,bo,block_t_k,der_ext_i,der_ext_j,der_ext_k,ablock,force,scoord,&
2080 !$OMP block_k_not_zero,der_k_zero,skip,der_j_zero,block_t_j,block_j_not_zero)
2081 
2082  mepos = 0
2083 !$ mepos = omp_get_thread_num()
2084 
2085  CALL cp_libint_init_3eri1(lib, max(maxli, maxlj, maxlk))
2086  CALL cp_libint_set_contrdepth(lib, 1)
2087 
2088  !pre-allocate contraction buffers
2089  ALLOCATE (cpp_buffer(max_nsgfj*max_ncok), ccp_buffer(max_nsgfj*max_nsgfk*max_ncoi))
2090 
2091  CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
2092 
2093  !We split the provided bounds among the threads such that each threads works on a different set of atoms
2094 
2095  bo = get_limit(natom, nthread, mepos)
2096  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i=bo)
2097 
2098  skip = .false.
2099  IF (bo(1) > bo(2)) skip = .true.
2100 
2101  DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
2102  CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
2103  iatom=iatom, jatom=jatom, katom=katom, &
2104  rij=rij, rjk=rjk, rik=rik)
2105  IF (skip) EXIT
2106 
2107  CALL dbt_get_block(t3c_trace, [iatom, jatom, katom], tmp_block, found)
2108  IF (.NOT. found) cycle
2109 
2110  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
2111  npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
2112  sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
2113 
2114  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
2115  npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
2116  sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
2117 
2118  CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
2119  npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
2120  sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
2121 
2122  djk = norm2(rjk)
2123  dij = norm2(rij)
2124  dik = norm2(rik)
2125 
2126  IF (kind_radius_j + kind_radius_i + dr_ij < dij) cycle
2127  IF (kind_radius_j + kind_radius_k + dr_jk < djk) cycle
2128  IF (kind_radius_k + kind_radius_i + dr_ik < dik) cycle
2129 
2130  ALLOCATE (max_contraction_i(nseti))
2131  max_contraction_i = 0.0_dp
2132  DO iset = 1, nseti
2133  sgfi = first_sgf_i(1, iset)
2134  max_contraction_i(iset) = maxval((/(sum(abs(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/))
2135  END DO
2136 
2137  ALLOCATE (max_contraction_j(nsetj))
2138  max_contraction_j = 0.0_dp
2139  DO jset = 1, nsetj
2140  sgfj = first_sgf_j(1, jset)
2141  max_contraction_j(jset) = maxval((/(sum(abs(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/))
2142  END DO
2143 
2144  ALLOCATE (max_contraction_k(nsetk))
2145  max_contraction_k = 0.0_dp
2146  DO kset = 1, nsetk
2147  sgfk = first_sgf_k(1, kset)
2148  max_contraction_k(kset) = maxval((/(sum(abs(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/))
2149  END DO
2150 
2151  CALL dbt_blk_sizes(t3c_trace, [iatom, jatom, katom], blk_size)
2152 
2153  ALLOCATE (block_t_i(blk_size(2), blk_size(3), blk_size(1), 3))
2154  ALLOCATE (block_t_j(blk_size(2), blk_size(3), blk_size(1), 3))
2155  ALLOCATE (block_t_k(blk_size(2), blk_size(3), blk_size(1), 3))
2156 
2157  ALLOCATE (ablock(blk_size(2), blk_size(3), blk_size(1)))
2158  DO i = 1, blk_size(1)
2159  ablock(:, :, i) = tmp_block(i, :, :)
2160  END DO
2161  DEALLOCATE (tmp_block)
2162 
2163  block_t_i = 0.0_dp
2164  block_t_j = 0.0_dp
2165  block_t_k = 0.0_dp
2166  block_j_not_zero = .false.
2167  block_k_not_zero = .false.
2168 
2169  DO iset = 1, nseti
2170 
2171  DO jset = 1, nsetj
2172 
2173  IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) cycle
2174 
2175  DO kset = 1, nsetk
2176 
2177  IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) cycle
2178  IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) cycle
2179 
2180  ncoi = npgfi(iset)*ncoset(lmax_i(iset))
2181  ncoj = npgfj(jset)*ncoset(lmax_j(jset))
2182  ncok = npgfk(kset)*ncoset(lmax_k(kset))
2183 
2184  sgfi = first_sgf_i(1, iset)
2185  sgfj = first_sgf_j(1, jset)
2186  sgfk = first_sgf_k(1, kset)
2187 
2188  IF (ncoj*ncok*ncoi > 0) THEN
2189  ALLOCATE (dijk_j(ncoj, ncok, ncoi, 3))
2190  ALLOCATE (dijk_k(ncoj, ncok, ncoi, 3))
2191  dijk_j(:, :, :, :) = 0.0_dp
2192  dijk_k(:, :, :, :) = 0.0_dp
2193 
2194  der_j_zero = .false.
2195  der_k_zero = .false.
2196 
2197  !need positions for libint. Only relative positions are needed => set ri to 0.0
2198  ri = 0.0_dp
2199  rj = rij ! ri + rij
2200  rk = rik ! ri + rik
2201 
2202  CALL eri_3center_derivs(dijk_j, dijk_k, &
2203  lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
2204  lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
2205  lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
2206  djk, dij, dik, lib, potential_parameter, &
2207  der_abc_1_ext=der_ext_j, der_abc_2_ext=der_ext_k)
2208 
2209  IF (PRESENT(der_eps)) THEN
2210  DO i_xyz = 1, 3
2211  IF (der_eps > der_ext_j(i_xyz)*(max_contraction_i(iset)* &
2212  max_contraction_j(jset)* &
2213  max_contraction_k(kset))) THEN
2214  der_j_zero(i_xyz) = .true.
2215  END IF
2216  END DO
2217 
2218  DO i_xyz = 1, 3
2219  IF (der_eps > der_ext_k(i_xyz)*(max_contraction_i(iset)* &
2220  max_contraction_j(jset)* &
2221  max_contraction_k(kset))) THEN
2222  der_k_zero(i_xyz) = .true.
2223  END IF
2224  END DO
2225  IF (all(der_j_zero) .AND. all(der_k_zero)) THEN
2226  DEALLOCATE (dijk_j, dijk_k)
2227  cycle
2228  END IF
2229  END IF
2230 
2231  ALLOCATE (dijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
2232 
2233  block_start_j = sgfj
2234  block_end_j = sgfj + nsgfj(jset) - 1
2235  block_start_k = sgfk
2236  block_end_k = sgfk + nsgfk(kset) - 1
2237  block_start_i = sgfi
2238  block_end_i = sgfi + nsgfi(iset) - 1
2239 
2240  DO i_xyz = 1, 3
2241  IF (der_j_zero(i_xyz)) cycle
2242 
2243  block_j_not_zero(i_xyz) = .true.
2244  CALL libxsmm_abc_contract(dijk_contr, dijk_j(:, :, :, i_xyz), tspj(jset, jkind)%array, &
2245  spk(kset, kkind)%array, spi(iset, ikind)%array, &
2246  ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
2247  nsgfi(iset), cpp_buffer, ccp_buffer)
2248 
2249  block_t_j(block_start_j:block_end_j, &
2250  block_start_k:block_end_k, &
2251  block_start_i:block_end_i, i_xyz) = &
2252  block_t_j(block_start_j:block_end_j, &
2253  block_start_k:block_end_k, &
2254  block_start_i:block_end_i, i_xyz) + &
2255  dijk_contr(:, :, :)
2256 
2257  END DO
2258 
2259  DO i_xyz = 1, 3
2260  IF (der_k_zero(i_xyz)) cycle
2261 
2262  block_k_not_zero(i_xyz) = .true.
2263  CALL libxsmm_abc_contract(dijk_contr, dijk_k(:, :, :, i_xyz), tspj(jset, jkind)%array, &
2264  spk(kset, kkind)%array, spi(iset, ikind)%array, &
2265  ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
2266  nsgfi(iset), cpp_buffer, ccp_buffer)
2267 
2268  block_t_k(block_start_j:block_end_j, &
2269  block_start_k:block_end_k, &
2270  block_start_i:block_end_i, i_xyz) = &
2271  block_t_k(block_start_j:block_end_j, &
2272  block_start_k:block_end_k, &
2273  block_start_i:block_end_i, i_xyz) + &
2274  dijk_contr(:, :, :)
2275 
2276  END DO
2277 
2278  DEALLOCATE (dijk_j, dijk_k, dijk_contr)
2279  END IF ! number of triples > 0
2280  END DO
2281  END DO
2282  END DO
2283 
2284  !We obtain the derivative wrt to first center using translational invariance
2285  DO i_xyz = 1, 3
2286  block_t_i(:, :, :, i_xyz) = -block_t_j(:, :, :, i_xyz) - block_t_k(:, :, :, i_xyz)
2287  END DO
2288 
2289  !virial contribution coming from deriv wrt to first center
2290  DO i_xyz = 1, 3
2291  force = pref*sum(ablock(:, :, :)*block_t_i(:, :, :, i_xyz))
2292  CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
2293  DO j_xyz = 1, 3
2294 !$OMP ATOMIC
2295  work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
2296  END DO
2297  END DO
2298 
2299  !second center
2300  DO i_xyz = 1, 3
2301  force = pref*sum(ablock(:, :, :)*block_t_j(:, :, :, i_xyz))
2302  CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
2303  DO j_xyz = 1, 3
2304 !$OMP ATOMIC
2305  work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
2306  END DO
2307  END DO
2308 
2309  !third center
2310  DO i_xyz = 1, 3
2311  force = pref*sum(ablock(:, :, :)*block_t_k(:, :, :, i_xyz))
2312  CALL real_to_scaled(scoord, particle_set(iatom)%r + rik, cell)
2313  DO j_xyz = 1, 3
2314 !$OMP ATOMIC
2315  work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
2316  END DO
2317  END DO
2318 
2319  DEALLOCATE (block_t_i)
2320  DEALLOCATE (block_t_j)
2321  DEALLOCATE (block_t_k)
2322 
2323  DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k, ablock)
2324  END DO
2325 
2326  CALL cp_libint_cleanup_3eri1(lib)
2327  CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
2328 !$OMP END PARALLEL
2329 
2330  DO iset = 1, max_nset
2331  DO ibasis = 1, nbasis
2332  IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
2333  IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
2334 
2335  IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
2336  END DO
2337  END DO
2338 
2339  DEALLOCATE (spi, tspj, spk)
2340 
2341  CALL timestop(handle)
2342  END SUBROUTINE calc_3c_virial
2343 
2344 ! **************************************************************************************************
2345 !> \brief Build 3-center integral tensor
2346 !> \param t3c empty DBCSR tensor
2347 !> Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages)
2348 !> if k-points are used
2349 !> \param filter_eps Filter threshold for tensor blocks
2350 !> \param qs_env ...
2351 !> \param nl_3c 3-center neighborlist
2352 !> \param basis_i ...
2353 !> \param basis_j ...
2354 !> \param basis_k ...
2355 !> \param potential_parameter ...
2356 !> \param int_eps neglect integrals smaller than int_eps
2357 !> \param op_pos operator position.
2358 !> 1: calculate (i|jk) integrals,
2359 !> 2: calculate (ij|k) integrals
2360 !> \param do_kpoints ...
2361 !> this routine requires that libint has been static initialised somewhere else
2362 !> \param do_hfx_kpoints ...
2363 !> \param desymmetrize ...
2364 !> \param bounds_i ...
2365 !> \param bounds_j ...
2366 !> \param bounds_k ...
2367 !> \param RI_range ...
2368 !> \param img_to_RI_cell ...
2369 ! **************************************************************************************************
2370  SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
2371  nl_3c, basis_i, basis_j, basis_k, &
2372  potential_parameter, int_eps, &
2373  op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, &
2374  bounds_i, bounds_j, bounds_k, &
2375  RI_range, img_to_RI_cell)
2376 
2377  TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT) :: t3c
2378  REAL(kind=dp), INTENT(IN) :: filter_eps
2379  TYPE(qs_environment_type), POINTER :: qs_env
2380  TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
2381  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
2382  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
2383  REAL(kind=dp), INTENT(IN), OPTIONAL :: int_eps
2384  INTEGER, INTENT(IN), OPTIONAL :: op_pos
2385  LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints, desymmetrize
2386  INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
2387  REAL(dp), INTENT(IN), OPTIONAL :: ri_range
2388  INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: img_to_ri_cell
2389 
2390  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_3c_integrals'
2391 
2392  INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
2393  block_start_k, egfi, handle, handle2, i, iatom, ibasis, ikind, ilist, imax, iset, jatom, &
2394  jcell, jkind, jset, katom, kcell, kkind, kset, m_max, max_ncoi, max_ncoj, max_ncok, &
2395  max_nset, max_nsgfi, max_nsgfj, max_nsgfk, maxli, maxlj, maxlk, mepos, natom, nbasis, &
2396  ncell_ri, ncoi, ncoj, ncok, nimg, nseti, nsetj, nsetk, nthread, op_ij, op_jk, op_pos_prv, &
2397  sgfi, sgfj, sgfk, unit_id
2398  INTEGER, ALLOCATABLE, DIMENSION(:) :: img_to_ri_cell_prv
2399  INTEGER, DIMENSION(2) :: bo
2400  INTEGER, DIMENSION(3) :: blk_idx, blk_size, cell_j, cell_k, &
2401  kp_index_lbounds, kp_index_ubounds, sp
2402  INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
2403  lmin_k, npgfi, npgfj, npgfk, nsgfi, &
2404  nsgfj, nsgfk
2405  INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
2406  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2407  LOGICAL :: block_not_zero, debug, desymmetrize_prv, &
2408  do_hfx_kpoints_prv, do_kpoints_prv, &
2409  found, skip
2410  REAL(kind=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
2411  kind_radius_i, kind_radius_j, &
2412  kind_radius_k, max_contraction_i, &
2413  prefac, sijk_ext
2414  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
2415  max_contraction_j, max_contraction_k
2416  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_t, dummy_block_t, sijk, &
2417  sijk_contr, tmp_ijk
2418  REAL(kind=dp), DIMENSION(1, 1, 1) :: counter
2419  REAL(kind=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
2420  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
2421  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
2422  sphi_k, zeti, zetj, zetk
2423  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2424  TYPE(cell_type), POINTER :: cell
2425  TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
2426  TYPE(cp_libint_t) :: lib
2427  TYPE(dbt_type) :: t_3c_tmp
2428  TYPE(dft_control_type), POINTER :: dft_control
2429  TYPE(gto_basis_set_type), POINTER :: basis_set
2430  TYPE(kpoint_type), POINTER :: kpoints
2431  TYPE(mp_para_env_type), POINTER :: para_env
2432  TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
2433  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2434 
2435  CALL timeset(routinen, handle)
2436 
2437  debug = .false.
2438 
2439  IF (PRESENT(do_kpoints)) THEN
2440  do_kpoints_prv = do_kpoints
2441  ELSE
2442  do_kpoints_prv = .false.
2443  END IF
2444 
2445  IF (PRESENT(do_hfx_kpoints)) THEN
2446  do_hfx_kpoints_prv = do_hfx_kpoints
2447  ELSE
2448  do_hfx_kpoints_prv = .false.
2449  END IF
2450  IF (do_hfx_kpoints_prv) THEN
2451  cpassert(do_kpoints_prv)
2452  cpassert(PRESENT(ri_range))
2453  cpassert(PRESENT(img_to_ri_cell))
2454  END IF
2455 
2456  IF (PRESENT(img_to_ri_cell)) THEN
2457  ALLOCATE (img_to_ri_cell_prv(SIZE(img_to_ri_cell)))
2458  img_to_ri_cell_prv(:) = img_to_ri_cell
2459  END IF
2460 
2461  IF (PRESENT(desymmetrize)) THEN
2462  desymmetrize_prv = desymmetrize
2463  ELSE
2464  desymmetrize_prv = .true.
2465  END IF
2466 
2467  op_ij = do_potential_id; op_jk = do_potential_id
2468 
2469  IF (PRESENT(op_pos)) THEN
2470  op_pos_prv = op_pos
2471  ELSE
2472  op_pos_prv = 1
2473  END IF
2474 
2475  SELECT CASE (op_pos_prv)
2476  CASE (1)
2477  op_ij = potential_parameter%potential_type
2478  CASE (2)
2479  op_jk = potential_parameter%potential_type
2480  END SELECT
2481 
2482  dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
2483 
2484  IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
2485  dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
2486  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
2487  ELSEIF (op_ij == do_potential_coulomb) THEN
2488  dr_ij = 1000000.0_dp
2489  dr_ik = 1000000.0_dp
2490  END IF
2491 
2492  IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
2493  dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
2494  dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
2495  ELSEIF (op_jk == do_potential_coulomb) THEN
2496  dr_jk = 1000000.0_dp
2497  dr_ik = 1000000.0_dp
2498  END IF
2499 
2500  NULLIFY (qs_kind_set, atomic_kind_set)
2501 
2502  ! get stuff
2503  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
2504  natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
2505 
2506  CALL alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, &
2507  op_pos=op_pos_prv, do_kpoints=do_kpoints, do_hfx_kpoints=do_hfx_kpoints, &
2508  bounds_i=bounds_i, bounds_j=bounds_j, bounds_k=bounds_k, &
2509  ri_range=ri_range, img_to_ri_cell=img_to_ri_cell)
2510 
2511  IF (do_kpoints_prv) THEN
2512  nimg = dft_control%nimages
2513  ncell_ri = nimg
2514  IF (do_hfx_kpoints_prv) THEN
2515  nimg = SIZE(t3c, 1)
2516  ncell_ri = SIZE(t3c, 2)
2517  END IF
2518  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
2519  ELSE
2520  nimg = 1
2521  ncell_ri = 1
2522  END IF
2523 
2524  IF (do_hfx_kpoints_prv) THEN
2525  cpassert(op_pos_prv == 2)
2526  cpassert(.NOT. desymmetrize_prv)
2527  ELSE IF (do_kpoints_prv) THEN
2528  cpassert(all(shape(t3c) == [nimg, ncell_ri]))
2529  END IF
2530 
2531  !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
2532  nbasis = SIZE(basis_i)
2533  max_nsgfi = 0
2534  max_ncoi = 0
2535  max_nset = 0
2536  maxli = 0
2537  DO ibasis = 1, nbasis
2538  CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
2539  nset=iset, nsgf_set=nsgfi, npgf=npgfi)
2540  maxli = max(maxli, imax)
2541  max_nset = max(max_nset, iset)
2542  max_nsgfi = max(max_nsgfi, maxval(nsgfi))
2543  max_ncoi = max(max_ncoi, maxval(npgfi)*ncoset(maxli))
2544  END DO
2545  max_nsgfj = 0
2546  max_ncoj = 0
2547  maxlj = 0
2548  DO ibasis = 1, nbasis
2549  CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
2550  nset=jset, nsgf_set=nsgfj, npgf=npgfj)
2551  maxlj = max(maxlj, imax)
2552  max_nset = max(max_nset, jset)
2553  max_nsgfj = max(max_nsgfj, maxval(nsgfj))
2554  max_ncoj = max(max_ncoj, maxval(npgfj)*ncoset(maxlj))
2555  END DO
2556  max_nsgfk = 0
2557  max_ncok = 0
2558  maxlk = 0
2559  DO ibasis = 1, nbasis
2560  CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
2561  nset=kset, nsgf_set=nsgfk, npgf=npgfk)
2562  maxlk = max(maxlk, imax)
2563  max_nset = max(max_nset, kset)
2564  max_nsgfk = max(max_nsgfk, maxval(nsgfk))
2565  max_ncok = max(max_ncok, maxval(npgfk)*ncoset(maxlk))
2566  END DO
2567  m_max = maxli + maxlj + maxlk
2568 
2569  !To minimize expensive memory ops and generally optimize contraction, pre-allocate
2570  !contiguous sphi arrays (and transposed in the case of sphi_i)
2571 
2572  NULLIFY (tspj, spi, spk)
2573  ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
2574 
2575  DO ibasis = 1, nbasis
2576  DO iset = 1, max_nset
2577  NULLIFY (spi(iset, ibasis)%array)
2578  NULLIFY (tspj(iset, ibasis)%array)
2579 
2580  NULLIFY (spk(iset, ibasis)%array)
2581  END DO
2582  END DO
2583 
2584  DO ilist = 1, 3
2585  DO ibasis = 1, nbasis
2586  IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
2587  IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
2588  IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
2589 
2590  DO iset = 1, basis_set%nset
2591 
2592  ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
2593  sgfi = basis_set%first_sgf(1, iset)
2594  egfi = sgfi + basis_set%nsgf_set(iset) - 1
2595 
2596  IF (ilist == 1) THEN
2597  ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2598  spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2599 
2600  ELSE IF (ilist == 2) THEN
2601  ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
2602  tspj(iset, ibasis)%array(:, :) = transpose(basis_set%sphi(1:ncoi, sgfi:egfi))
2603 
2604  ELSE
2605  ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2606  spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2607  END IF
2608 
2609  END DO !iset
2610  END DO !ibasis
2611  END DO !ilist
2612 
2613  !Init the truncated Coulomb operator
2614  IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN
2615 
2616  IF (m_max > get_lmax_init()) THEN
2617  IF (para_env%mepos == 0) THEN
2618  CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
2619  END IF
2620  CALL init(m_max, unit_id, para_env%mepos, para_env)
2621  IF (para_env%mepos == 0) THEN
2622  CALL close_file(unit_id)
2623  END IF
2624  END IF
2625  END IF
2626 
2627  CALL init_md_ftable(nmax=m_max)
2628 
2629  IF (do_kpoints_prv) THEN
2630  kp_index_lbounds = lbound(cell_to_index)
2631  kp_index_ubounds = ubound(cell_to_index)
2632  END IF
2633 
2634  counter = 1.0_dp
2635 
2636  nthread = 1
2637 !$ nthread = omp_get_max_threads()
2638 
2639 !$OMP PARALLEL DEFAULT(NONE) &
2640 !$OMP SHARED (nthread,do_kpoints_prv,kp_index_lbounds,kp_index_ubounds,maxli,maxlk,maxlj,bounds_i,&
2641 !$OMP bounds_j,bounds_k,nimg,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,&
2642 !$OMP potential_parameter,int_eps,t3c,tspj,spi,spk,debug,cell_to_index,max_ncoi,max_nsgfk,&
2643 !$OMP max_nsgfj,max_ncok,natom,nl_3c,cell,op_pos_prv,do_hfx_kpoints_prv,RI_range,ncell_RI, &
2644 !$OMP img_to_RI_cell_prv) &
2645 !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,cell_j,cell_k,&
2646 !$OMP prefac,jcell,kcell,first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
2647 !$OMP sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
2648 !$OMP set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
2649 !$OMP rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
2650 !$OMP sgfk,sijk,ri,rj,rk,sijk_ext,block_not_zero,max_contraction_i,max_contraction_j,&
2651 !$OMP max_contraction_k,iset,jset,kset,block_t,blk_size,sijk_contr,cpp_buffer,ccp_buffer,&
2652 !$OMP block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
2653 !$OMP dummy_block_t,sp,handle2,mepos,bo,skip,tmp_ijk,i,blk_idx)
2654 
2655  mepos = 0
2656 !$ mepos = omp_get_thread_num()
2657 
2658  CALL cp_libint_init_3eri(lib, max(maxli, maxlj, maxlk))
2659  CALL cp_libint_set_contrdepth(lib, 1)
2660 
2661  !pre-allocate contraction buffers
2662  ALLOCATE (cpp_buffer(max_nsgfj*max_ncok), ccp_buffer(max_nsgfj*max_nsgfk*max_ncoi))
2663 
2664  CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
2665 
2666  !We split the provided bounds among the threads such that each threads works on a different set of atoms
2667  IF (PRESENT(bounds_i)) THEN
2668  bo = get_limit(bounds_i(2) - bounds_i(1) + 1, nthread, mepos)
2669  bo(:) = bo(:) + bounds_i(1) - 1
2670  CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
2671  ELSE IF (PRESENT(bounds_j)) THEN
2672 
2673  bo = get_limit(bounds_j(2) - bounds_j(1) + 1, nthread, mepos)
2674  bo(:) = bo(:) + bounds_j(1) - 1
2675  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bo, bounds_k)
2676  ELSE IF (PRESENT(bounds_k)) THEN
2677  bo = get_limit(bounds_k(2) - bounds_k(1) + 1, nthread, mepos)
2678  bo(:) = bo(:) + bounds_k(1) - 1
2679  CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bo)
2680  ELSE
2681  bo = get_limit(natom, nthread, mepos)
2682  CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
2683  END IF
2684 
2685  skip = .false.
2686  IF (bo(1) > bo(2)) skip = .true.
2687 
2688  DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
2689  CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
2690  iatom=iatom, jatom=jatom, katom=katom, &
2691  rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
2692  IF (skip) EXIT
2693 
2694  djk = norm2(rjk)
2695  dij = norm2(rij)
2696  dik = norm2(rik)
2697 
2698  IF (do_kpoints_prv) THEN
2699  prefac = 0.5_dp
2700  ELSEIF (nl_3c%sym == symmetric_jk) THEN
2701  IF (jatom == katom) THEN
2702  ! factor 0.5 due to double-counting of diagonal blocks
2703  ! (we desymmetrize by adding transpose)
2704  prefac = 0.5_dp
2705  ELSE
2706  prefac = 1.0_dp
2707  END IF
2708  ELSEIF (nl_3c%sym == symmetric_ij) THEN
2709  IF (iatom == jatom) THEN
2710  ! factor 0.5 due to double-counting of diagonal blocks
2711  ! (we desymmetrize by adding transpose)
2712  prefac = 0.5_dp
2713  ELSE
2714  prefac = 1.0_dp
2715  END IF
2716  ELSE
2717  prefac = 1.0_dp
2718  END IF
2719  IF (do_hfx_kpoints_prv) prefac = 1.0_dp
2720 
2721  IF (do_kpoints_prv) THEN
2722 
2723  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
2724  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
2725 
2726  jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2727  IF (jcell > nimg .OR. jcell < 1) cycle
2728 
2729  IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
2730  any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
2731 
2732  kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
2733  IF (kcell > nimg .OR. kcell < 1) cycle
2734  !In case of HFX k-points, we only consider P^k if d_ik <= RI_range
2735  IF (do_hfx_kpoints_prv) THEN
2736  IF (dik > ri_range) cycle
2737  kcell = img_to_ri_cell_prv(kcell)
2738  END IF
2739  ELSE
2740  jcell = 1; kcell = 1
2741  END IF
2742 
2743  blk_idx = [iatom, jatom, katom]
2744  IF (do_hfx_kpoints_prv) THEN
2745  blk_idx(3) = (kcell - 1)*natom + katom
2746  kcell = 1
2747  END IF
2748 
2749  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
2750  npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
2751  sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
2752 
2753  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
2754  npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
2755  sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
2756 
2757  CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
2758  npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
2759  sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
2760 
2761  IF (kind_radius_j + kind_radius_i + dr_ij < dij) cycle
2762  IF (kind_radius_j + kind_radius_k + dr_jk < djk) cycle
2763  IF (kind_radius_k + kind_radius_i + dr_ik < dik) cycle
2764 
2765  ALLOCATE (max_contraction_j(nsetj))
2766  DO jset = 1, nsetj
2767  sgfj = first_sgf_j(1, jset)
2768  max_contraction_j(jset) = maxval((/(sum(abs(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)/))
2769  END DO
2770 
2771  ALLOCATE (max_contraction_k(nsetk))
2772  DO kset = 1, nsetk
2773  sgfk = first_sgf_k(1, kset)
2774  max_contraction_k(kset) = maxval((/(sum(abs(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)/))
2775  END DO
2776 
2777  CALL dbt_blk_sizes(t3c(jcell, kcell), blk_idx, blk_size)
2778 
2779  ALLOCATE (block_t(blk_size(2), blk_size(3), blk_size(1)))
2780 
2781  block_t = 0.0_dp
2782  block_not_zero = .false.
2783  DO iset = 1, nseti
2784 
2785  sgfi = first_sgf_i(1, iset)
2786  max_contraction_i = maxval((/(sum(abs(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)/))
2787 
2788  DO jset = 1, nsetj
2789 
2790  IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) cycle
2791 
2792  DO kset = 1, nsetk
2793 
2794  IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) cycle
2795  IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) cycle
2796 
2797  ncoi = npgfi(iset)*ncoset(lmax_i(iset))
2798  ncoj = npgfj(jset)*ncoset(lmax_j(jset))
2799  ncok = npgfk(kset)*ncoset(lmax_k(kset))
2800 
2801  !ensure non-zero number of triples below
2802  IF (ncoj*ncok*ncoi == 0) cycle
2803 
2804  !need positions for libint. Only relative positions are needed => set ri to 0.0
2805  ri = 0.0_dp
2806  rj = rij ! ri + rij
2807  rk = rik ! ri + rik
2808 
2809  ALLOCATE (sijk(ncoj, ncok, ncoi))
2810  IF (op_pos_prv == 1) THEN
2811  sijk(:, :, :) = 0.0_dp
2812  CALL eri_3center(sijk, &
2813  lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
2814  lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
2815  lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
2816  djk, dij, dik, lib, potential_parameter, int_abc_ext=sijk_ext)
2817  ELSE
2818  ALLOCATE (tmp_ijk(ncoi, ncoj, ncok))
2819  tmp_ijk(:, :, :) = 0.0_dp
2820  CALL eri_3center(tmp_ijk, &
2821  lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
2822  lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
2823  lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
2824  dij, dik, djk, lib, potential_parameter, int_abc_ext=sijk_ext)
2825 
2826  !F08: sijk = RESHAPE(tmp_ijk, [ncoj, ncok, ncoi], ORDER=[2, 3, 1]) with sijk not allocated
2827  DO i = 1, ncoi !TODO: revise/check for efficiency
2828  sijk(:, :, i) = tmp_ijk(i, :, :)
2829  END DO
2830  DEALLOCATE (tmp_ijk)
2831  END IF
2832 
2833  IF (PRESENT(int_eps)) THEN
2834  IF (int_eps > sijk_ext*(max_contraction_i* &
2835  max_contraction_j(jset)* &
2836  max_contraction_k(kset))) THEN
2837  DEALLOCATE (sijk)
2838  cycle
2839  END IF
2840  END IF
2841 
2842  block_not_zero = .true.
2843  ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
2844  CALL libxsmm_abc_contract(sijk_contr, sijk, tspj(jset, jkind)%array, &
2845  spk(kset, kkind)%array, spi(iset, ikind)%array, &
2846  ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
2847  nsgfi(iset), cpp_buffer, ccp_buffer)
2848  DEALLOCATE (sijk)
2849 
2850  sgfj = first_sgf_j(1, jset)
2851  sgfk = first_sgf_k(1, kset)
2852 
2853  block_start_j = sgfj
2854  block_end_j = sgfj + nsgfj(jset) - 1
2855  block_start_k = sgfk
2856  block_end_k = sgfk + nsgfk(kset) - 1
2857  block_start_i = sgfi
2858  block_end_i = sgfi + nsgfi(iset) - 1
2859 
2860  block_t(block_start_j:block_end_j, &
2861  block_start_k:block_end_k, &
2862  block_start_i:block_end_i) = &
2863  block_t(block_start_j:block_end_j, &
2864  block_start_k:block_end_k, &
2865  block_start_i:block_end_i) + &
2866  prefac*sijk_contr(:, :, :)
2867  DEALLOCATE (sijk_contr)
2868 
2869  END DO
2870 
2871  END DO
2872 
2873  END DO
2874 
2875  IF (block_not_zero) THEN
2876 !$OMP CRITICAL
2877  CALL timeset(routinen//"_put_dbcsr", handle2)
2878  IF (debug) THEN
2879  CALL dbt_get_block(t3c(jcell, kcell), blk_idx, dummy_block_t, found=found)
2880  cpassert(found)
2881  END IF
2882 
2883  sp = shape(block_t)
2884  sp([2, 3, 1]) = sp ! sp = sp([2, 3, 1]) performs worse
2885 
2886  CALL dbt_put_block(t3c(jcell, kcell), blk_idx, sp, &
2887  reshape(block_t, shape=sp, order=[2, 3, 1]), summation=.true.)
2888 
2889  CALL timestop(handle2)
2890 !$OMP END CRITICAL
2891  END IF
2892 
2893  DEALLOCATE (block_t)
2894  DEALLOCATE (max_contraction_j, max_contraction_k)
2895  END DO
2896 
2897  CALL cp_libint_cleanup_3eri(lib)
2898 
2899  CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
2900 !$OMP END PARALLEL
2901 
2902  !TODO: deal with hfx_kpoints, because should not filter by 1/2
2903  IF (nl_3c%sym == symmetric_jk .OR. do_kpoints_prv) THEN
2904 
2905  IF (.NOT. do_hfx_kpoints_prv) THEN
2906  DO kcell = 1, nimg
2907  DO jcell = 1, nimg
2908  ! need half of filter eps because afterwards we add transposed tensor
2909  CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
2910  END DO
2911  END DO
2912  ELSE
2913  DO kcell = 1, ncell_ri
2914  DO jcell = 1, nimg
2915  CALL dbt_filter(t3c(jcell, kcell), filter_eps)
2916  END DO
2917  END DO
2918  END IF
2919 
2920  IF (desymmetrize_prv) THEN
2921  ! add transposed of overlap integrals
2922  CALL dbt_create(t3c(1, 1), t_3c_tmp)
2923  DO kcell = 1, jcell
2924  DO jcell = 1, nimg
2925  CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
2926  CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.true., move_data=.true.)
2927  CALL dbt_filter(t3c(kcell, jcell), filter_eps)
2928  END DO
2929  END DO
2930  DO kcell = jcell + 1, nimg
2931  DO jcell = 1, nimg
2932  CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
2933  CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.false., move_data=.true.)
2934  CALL dbt_filter(t3c(kcell, jcell), filter_eps)
2935  END DO
2936  END DO
2937  CALL dbt_destroy(t_3c_tmp)
2938  END IF
2939  ELSEIF (nl_3c%sym == symmetric_ij) THEN
2940  DO kcell = 1, nimg
2941  DO jcell = 1, nimg
2942  CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
2943  END DO
2944  END DO
2945  ELSEIF (nl_3c%sym == symmetric_none) THEN
2946  DO kcell = 1, nimg
2947  DO jcell = 1, nimg
2948  CALL dbt_filter(t3c(jcell, kcell), filter_eps)
2949  END DO
2950  END DO
2951  ELSE
2952  cpabort("requested symmetric case not implemented")
2953  END IF
2954 
2955  DO iset = 1, max_nset
2956  DO ibasis = 1, nbasis
2957  IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
2958  IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
2959 
2960  IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
2961  END DO
2962  END DO
2963  DEALLOCATE (spi, tspj, spk)
2964 
2965  CALL timestop(handle)
2966  END SUBROUTINE build_3c_integrals
2967 
2968 ! **************************************************************************************************
2969 !> \brief Calculates the derivatives of 2-center integrals, wrt to the first center
2970 !> \param t2c_der ...
2971 !> this routine requires that libint has been static initialised somewhere else
2972 !> \param filter_eps Filter threshold for matrix blocks
2973 !> \param qs_env ...
2974 !> \param nl_2c 2-center neighborlist
2975 !> \param basis_i ...
2976 !> \param basis_j ...
2977 !> \param potential_parameter ...
2978 !> \param do_kpoints ...
2979 ! **************************************************************************************************
2980  SUBROUTINE build_2c_derivatives(t2c_der, filter_eps, qs_env, &
2981  nl_2c, basis_i, basis_j, &
2982  potential_parameter, do_kpoints)
2983 
2984  TYPE(dbcsr_type), DIMENSION(:, :), INTENT(INOUT) :: t2c_der
2985  REAL(kind=dp), INTENT(IN) :: filter_eps
2986  TYPE(qs_environment_type), POINTER :: qs_env
2987  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2988  POINTER :: nl_2c
2989  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
2990  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
2991  LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints
2992 
2993  CHARACTER(len=*), PARAMETER :: routinen = 'build_2c_derivatives'
2994 
2995  INTEGER :: handle, i_xyz, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, &
2996  jset, m_max, maxli, maxlj, natom, ncoi, ncoj, nimg, nseti, nsetj, op_prv, sgfi, sgfj, &
2997  unit_id
2998  INTEGER, DIMENSION(3) :: cell_j, kp_index_lbounds, &
2999  kp_index_ubounds
3000  INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3001  npgfj, nsgfi, nsgfj
3002  INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3003  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3004  LOGICAL :: do_kpoints_prv, do_symmetric, found, &
3005  trans
3006  REAL(kind=dp) :: dab
3007  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dij_contr, dij_rs
3008  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dij
3009  REAL(kind=dp), DIMENSION(3) :: ri, rij, rj
3010  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3011  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3012  zetj
3013  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3014  TYPE(block_p_type), DIMENSION(3) :: block_t
3015  TYPE(cp_libint_t) :: lib
3016  TYPE(dft_control_type), POINTER :: dft_control
3017  TYPE(kpoint_type), POINTER :: kpoints
3018  TYPE(mp_para_env_type), POINTER :: para_env
3019  TYPE(neighbor_list_iterator_p_type), &
3020  DIMENSION(:), POINTER :: nl_iterator
3021  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3022 
3023  CALL timeset(routinen, handle)
3024 
3025  IF (PRESENT(do_kpoints)) THEN
3026  do_kpoints_prv = do_kpoints
3027  ELSE
3028  do_kpoints_prv = .false.
3029  END IF
3030 
3031  op_prv = potential_parameter%potential_type
3032 
3033  NULLIFY (qs_kind_set, atomic_kind_set, block_t(1)%block, block_t(2)%block, block_t(3)%block, cell_to_index)
3034 
3035  ! get stuff
3036  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
3037  natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
3038 
3039  IF (do_kpoints_prv) THEN
3040  nimg = SIZE(t2c_der, 1)
3041  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
3042  kp_index_lbounds = lbound(cell_to_index)
3043  kp_index_ubounds = ubound(cell_to_index)
3044  ELSE
3045  nimg = 1
3046  END IF
3047 
3048  ! check for symmetry
3049  cpassert(SIZE(nl_2c) > 0)
3050  CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3051 
3052  IF (do_symmetric) THEN
3053  DO img = 1, nimg
3054  !Derivtive matrix is assymetric
3055  DO i_xyz = 1, 3
3056  cpassert(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_antisymmetric)
3057  END DO
3058  END DO
3059  ELSE
3060  DO img = 1, nimg
3061  DO i_xyz = 1, 3
3062  cpassert(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_no_symmetry)
3063  END DO
3064  END DO
3065  END IF
3066 
3067  DO img = 1, nimg
3068  DO i_xyz = 1, 3
3069  CALL cp_dbcsr_alloc_block_from_nbl(t2c_der(img, i_xyz), nl_2c)
3070  END DO
3071  END DO
3072 
3073  maxli = 0
3074  DO ibasis = 1, SIZE(basis_i)
3075  CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3076  maxli = max(maxli, imax)
3077  END DO
3078  maxlj = 0
3079  DO ibasis = 1, SIZE(basis_j)
3080  CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3081  maxlj = max(maxlj, imax)
3082  END DO
3083 
3084  m_max = maxli + maxlj + 1
3085 
3086  !Init the truncated Coulomb operator
3087  IF (op_prv == do_potential_truncated) THEN
3088 
3089  IF (m_max > get_lmax_init()) THEN
3090  IF (para_env%mepos == 0) THEN
3091  CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3092  END IF
3093  CALL init(m_max, unit_id, para_env%mepos, para_env)
3094  IF (para_env%mepos == 0) THEN
3095  CALL close_file(unit_id)
3096  END IF
3097  END IF
3098  END IF
3099 
3100  CALL init_md_ftable(nmax=m_max)
3101 
3102  CALL cp_libint_init_2eri1(lib, max(maxli, maxlj))
3103  CALL cp_libint_set_contrdepth(lib, 1)
3104 
3105  CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3106  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3107 
3108  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3109  iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
3110  IF (do_kpoints_prv) THEN
3111  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
3112  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
3113  img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3114  IF (img > nimg .OR. img < 1) cycle
3115  ELSE
3116  img = 1
3117  END IF
3118 
3119  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
3120  npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
3121  sphi=sphi_i, zet=zeti)
3122 
3123  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
3124  npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
3125  sphi=sphi_j, zet=zetj)
3126 
3127  IF (do_symmetric) THEN
3128  IF (iatom <= jatom) THEN
3129  irow = iatom
3130  icol = jatom
3131  ELSE
3132  irow = jatom
3133  icol = iatom
3134  END IF
3135  ELSE
3136  irow = iatom
3137  icol = jatom
3138  END IF
3139 
3140  dab = norm2(rij)
3141  trans = do_symmetric .AND. (iatom > jatom)
3142 
3143  DO i_xyz = 1, 3
3144  CALL dbcsr_get_block_p(matrix=t2c_der(img, i_xyz), &
3145  row=irow, col=icol, block=block_t(i_xyz)%block, found=found)
3146  cpassert(found)
3147  END DO
3148 
3149  DO iset = 1, nseti
3150 
3151  ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3152  sgfi = first_sgf_i(1, iset)
3153 
3154  DO jset = 1, nsetj
3155 
3156  ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3157  sgfj = first_sgf_j(1, jset)
3158 
3159  IF (ncoi*ncoj > 0) THEN
3160  ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
3161  ALLOCATE (dij(ncoi, ncoj, 3))
3162  dij(:, :, :) = 0.0_dp
3163 
3164  ri = 0.0_dp
3165  rj = rij
3166 
3167  CALL eri_2center_derivs(dij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
3168  rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
3169  rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3170 
3171  DO i_xyz = 1, 3
3172 
3173  dij_contr(:, :) = 0.0_dp
3174  CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
3175  sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3176  ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3177 
3178  IF (trans) THEN
3179  !if transpose, then -1 factor for antisymmetry
3180  ALLOCATE (dij_rs(nsgfj(jset), nsgfi(iset)))
3181  dij_rs(:, :) = -1.0_dp*transpose(dij_contr)
3182  ELSE
3183  ALLOCATE (dij_rs(nsgfi(iset), nsgfj(jset)))
3184  dij_rs(:, :) = dij_contr
3185  END IF
3186 
3187  CALL block_add("IN", dij_rs, &
3188  nsgfi(iset), nsgfj(jset), block_t(i_xyz)%block, &
3189  sgfi, sgfj, trans=trans)
3190  DEALLOCATE (dij_rs)
3191  END DO
3192 
3193  DEALLOCATE (dij, dij_contr)
3194  END IF
3195  END DO
3196  END DO
3197  END DO
3198 
3199  CALL cp_libint_cleanup_2eri1(lib)
3200 
3201  CALL neighbor_list_iterator_release(nl_iterator)
3202  DO img = 1, nimg
3203  DO i_xyz = 1, 3
3204  CALL dbcsr_finalize(t2c_der(img, i_xyz))
3205  CALL dbcsr_filter(t2c_der(img, i_xyz), filter_eps)
3206  END DO
3207  END DO
3208 
3209  CALL timestop(handle)
3210 
3211  END SUBROUTINE build_2c_derivatives
3212 
3213 ! **************************************************************************************************
3214 !> \brief Calculates the virial coming from 2c derivatives on the fly
3215 !> \param work_virial ...
3216 !> \param t2c_trace the 2c tensor that we should trace with the derivatives
3217 !> \param pref ...
3218 !> \param qs_env ...
3219 !> \param nl_2c 2-center neighborlist. Assumed to have compatible distribution with t2c_trace,
3220 !> and to be non-symmetric
3221 !> \param basis_i ...
3222 !> \param basis_j ...
3223 !> \param potential_parameter ...
3224 ! **************************************************************************************************
3225  SUBROUTINE calc_2c_virial(work_virial, t2c_trace, pref, qs_env, nl_2c, basis_i, basis_j, potential_parameter)
3226  REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: work_virial
3227  TYPE(dbcsr_type), INTENT(INOUT) :: t2c_trace
3228  REAL(kind=dp), INTENT(IN) :: pref
3229  TYPE(qs_environment_type), POINTER :: qs_env
3230  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3231  POINTER :: nl_2c
3232  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
3233  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
3234 
3235  CHARACTER(len=*), PARAMETER :: routinen = 'calc_2c_virial'
3236 
3237  INTEGER :: handle, i_xyz, iatom, ibasis, ikind, imax, iset, j_xyz, jatom, jkind, jset, &
3238  m_max, maxli, maxlj, natom, ncoi, ncoj, nseti, nsetj, op_prv, sgfi, sgfj, unit_id
3239  INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3240  npgfj, nsgfi, nsgfj
3241  INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3242  LOGICAL :: do_symmetric, found
3243  REAL(dp) :: force
3244  REAL(dp), DIMENSION(:, :), POINTER :: pblock
3245  REAL(kind=dp) :: dab
3246  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dij_contr
3247  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dij
3248  REAL(kind=dp), DIMENSION(3) :: ri, rij, rj, scoord
3249  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3250  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3251  zetj
3252  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3253  TYPE(cell_type), POINTER :: cell
3254  TYPE(cp_libint_t) :: lib
3255  TYPE(dft_control_type), POINTER :: dft_control
3256  TYPE(mp_para_env_type), POINTER :: para_env
3257  TYPE(neighbor_list_iterator_p_type), &
3258  DIMENSION(:), POINTER :: nl_iterator
3259  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3260  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3261 
3262  CALL timeset(routinen, handle)
3263 
3264  op_prv = potential_parameter%potential_type
3265 
3266  NULLIFY (qs_kind_set, atomic_kind_set, pblock, particle_set, cell)
3267 
3268  ! get stuff
3269  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
3270  natom=natom, dft_control=dft_control, para_env=para_env, &
3271  particle_set=particle_set, cell=cell)
3272 
3273  ! check for symmetry
3274  cpassert(SIZE(nl_2c) > 0)
3275  CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3276  cpassert(.NOT. do_symmetric)
3277 
3278  maxli = 0
3279  DO ibasis = 1, SIZE(basis_i)
3280  CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3281  maxli = max(maxli, imax)
3282  END DO
3283  maxlj = 0
3284  DO ibasis = 1, SIZE(basis_j)
3285  CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3286  maxlj = max(maxlj, imax)
3287  END DO
3288 
3289  m_max = maxli + maxlj + 1
3290 
3291  !Init the truncated Coulomb operator
3292  IF (op_prv == do_potential_truncated) THEN
3293 
3294  IF (m_max > get_lmax_init()) THEN
3295  IF (para_env%mepos == 0) THEN
3296  CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3297  END IF
3298  CALL init(m_max, unit_id, para_env%mepos, para_env)
3299  IF (para_env%mepos == 0) THEN
3300  CALL close_file(unit_id)
3301  END IF
3302  END IF
3303  END IF
3304 
3305  CALL init_md_ftable(nmax=m_max)
3306 
3307  CALL cp_libint_init_2eri1(lib, max(maxli, maxlj))
3308  CALL cp_libint_set_contrdepth(lib, 1)
3309 
3310  CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3311  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3312 
3313  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3314  iatom=iatom, jatom=jatom, r=rij)
3315 
3316  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
3317  npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
3318  sphi=sphi_i, zet=zeti)
3319 
3320  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
3321  npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
3322  sphi=sphi_j, zet=zetj)
3323 
3324  dab = norm2(rij)
3325 
3326  CALL dbcsr_get_block_p(t2c_trace, iatom, jatom, pblock, found)
3327  IF (.NOT. found) cycle
3328 
3329  DO iset = 1, nseti
3330 
3331  ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3332  sgfi = first_sgf_i(1, iset)
3333 
3334  DO jset = 1, nsetj
3335 
3336  ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3337  sgfj = first_sgf_j(1, jset)
3338 
3339  IF (ncoi*ncoj > 0) THEN
3340  ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
3341  ALLOCATE (dij(ncoi, ncoj, 3))
3342  dij(:, :, :) = 0.0_dp
3343 
3344  ri = 0.0_dp
3345  rj = rij
3346 
3347  CALL eri_2center_derivs(dij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
3348  rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
3349  rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3350 
3351  DO i_xyz = 1, 3
3352 
3353  dij_contr(:, :) = 0.0_dp
3354  CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
3355  sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3356  ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3357 
3358  force = sum(pblock(sgfi:sgfi + nsgfi(iset) - 1, sgfj:sgfj + nsgfj(jset) - 1)*dij_contr(:, :))
3359  force = pref*force
3360 
3361  !iatom virial
3362  CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
3363  DO j_xyz = 1, 3
3364  work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
3365  END DO
3366 
3367  !jatom virial
3368  CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
3369  DO j_xyz = 1, 3
3370  work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - force*scoord(j_xyz)
3371  END DO
3372  END DO
3373 
3374  DEALLOCATE (dij, dij_contr)
3375  END IF
3376  END DO
3377  END DO
3378  END DO
3379  CALL neighbor_list_iterator_release(nl_iterator)
3380 
3381  CALL cp_libint_cleanup_2eri1(lib)
3382 
3383  CALL timestop(handle)
3384 
3385  END SUBROUTINE calc_2c_virial
3386 
3387 ! **************************************************************************************************
3388 !> \brief ...
3389 !> \param t2c empty DBCSR matrix
3390 !> \param filter_eps Filter threshold for matrix blocks
3391 !> \param qs_env ...
3392 !> \param nl_2c 2-center neighborlist
3393 !> \param basis_i ...
3394 !> \param basis_j ...
3395 !> \param potential_parameter ...
3396 !> \param do_kpoints ...
3397 !> \param do_hfx_kpoints ...
3398 !> this routine requires that libint has been static initialised somewhere else
3399 !> \param ext_kpoints ...
3400 !> \param regularization_RI ...
3401 ! **************************************************************************************************
3402  SUBROUTINE build_2c_integrals(t2c, filter_eps, qs_env, &
3403  nl_2c, basis_i, basis_j, &
3404  potential_parameter, do_kpoints, &
3405  do_hfx_kpoints, ext_kpoints, regularization_RI)
3406 
3407  TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: t2c
3408  REAL(kind=dp), INTENT(IN) :: filter_eps
3409  TYPE(qs_environment_type), POINTER :: qs_env
3410  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3411  POINTER :: nl_2c
3412  TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
3413  TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
3414  LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints
3415  TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
3416  REAL(kind=dp), OPTIONAL :: regularization_ri
3417 
3418  CHARACTER(len=*), PARAMETER :: routinen = 'build_2c_integrals'
3419 
3420  INTEGER :: handle, i_diag, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, &
3421  jset, m_max, maxli, maxlj, natom, ncoi, ncoj, nimg, nseti, nsetj, op_prv, sgfi, sgfj, &
3422  unit_id
3423  INTEGER, DIMENSION(3) :: cell_j, kp_index_lbounds, &
3424  kp_index_ubounds
3425  INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3426  npgfj, nsgfi, nsgfj
3427  INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3428  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3429  LOGICAL :: do_hfx_kpoints_prv, do_kpoints_prv, &
3430  do_symmetric, found, trans
3431  REAL(kind=dp) :: dab, my_regularization_ri
3432  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sij, sij_contr, sij_rs
3433  REAL(kind=dp), DIMENSION(3) :: ri, rij, rj
3434  REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3435  REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3436  zetj
3437  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3438  TYPE(block_p_type) :: block_t
3439  TYPE(cell_type), POINTER :: cell
3440  TYPE(cp_libint_t) :: lib
3441  TYPE(dft_control_type), POINTER :: dft_control
3442  TYPE(kpoint_type), POINTER :: kpoints
3443  TYPE(mp_para_env_type), POINTER :: para_env
3444  TYPE(neighbor_list_iterator_p_type), &
3445  DIMENSION(:), POINTER :: nl_iterator
3446  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3447 
3448  CALL timeset(routinen, handle)
3449 
3450  IF (PRESENT(do_kpoints)) THEN
3451  do_kpoints_prv = do_kpoints
3452  ELSE
3453  do_kpoints_prv = .false.
3454  END IF
3455 
3456  IF (PRESENT(do_hfx_kpoints)) THEN
3457  do_hfx_kpoints_prv = do_hfx_kpoints
3458  ELSE
3459  do_hfx_kpoints_prv = .false.
3460  END IF
3461  IF (do_hfx_kpoints_prv) THEN
3462  cpassert(do_kpoints_prv)
3463  END IF
3464 
3465  my_regularization_ri = 0.0_dp
3466  IF (PRESENT(regularization_ri)) my_regularization_ri = regularization_ri
3467 
3468  op_prv = potential_parameter%potential_type
3469 
3470  NULLIFY (qs_kind_set, atomic_kind_set, block_t%block, cell_to_index)
3471 
3472  ! get stuff
3473  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
3474  natom=natom, dft_control=dft_control, para_env=para_env, kpoints=kpoints)
3475 
3476  IF (PRESENT(ext_kpoints)) kpoints => ext_kpoints
3477 
3478  IF (do_kpoints_prv) THEN
3479  nimg = SIZE(t2c)
3480  CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
3481  kp_index_lbounds = lbound(cell_to_index)
3482  kp_index_ubounds = ubound(cell_to_index)
3483  ELSE
3484  nimg = 1
3485  END IF
3486 
3487  ! check for symmetry
3488  cpassert(SIZE(nl_2c) > 0)
3489  CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3490 
3491  IF (do_symmetric) THEN
3492  DO img = 1, nimg
3493  cpassert(dbcsr_has_symmetry(t2c(img)))
3494  END DO
3495  ELSE
3496  DO img = 1, nimg
3497  cpassert(.NOT. dbcsr_has_symmetry(t2c(img)))
3498  END DO
3499  END IF
3500 
3501  DO img = 1, nimg
3502  CALL cp_dbcsr_alloc_block_from_nbl(t2c(img), nl_2c)
3503  END DO
3504 
3505  maxli = 0
3506  DO ibasis = 1, SIZE(basis_i)
3507  CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3508  maxli = max(maxli, imax)
3509  END DO
3510  maxlj = 0
3511  DO ibasis = 1, SIZE(basis_j)
3512  CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3513  maxlj = max(maxlj, imax)
3514  END DO
3515 
3516  m_max = maxli + maxlj
3517 
3518  !Init the truncated Coulomb operator
3519  IF (op_prv == do_potential_truncated) THEN
3520 
3521  IF (m_max > get_lmax_init()) THEN
3522  IF (para_env%mepos == 0) THEN
3523  CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3524  END IF
3525  CALL init(m_max, unit_id, para_env%mepos, para_env)
3526  IF (para_env%mepos == 0) THEN
3527  CALL close_file(unit_id)
3528  END IF
3529  END IF
3530  END IF
3531 
3532  CALL init_md_ftable(nmax=m_max)
3533 
3534  CALL cp_libint_init_2eri(lib, max(maxli, maxlj))
3535  CALL cp_libint_set_contrdepth(lib, 1)
3536 
3537  CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3538  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3539 
3540  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3541  iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
3542  IF (do_kpoints_prv) THEN
3543  IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
3544  any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
3545  img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3546  IF (img > nimg .OR. img < 1) cycle
3547  ELSE
3548  img = 1
3549  END IF
3550 
3551  CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
3552  npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
3553  sphi=sphi_i, zet=zeti)
3554 
3555  CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
3556  npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
3557  sphi=sphi_j, zet=zetj)
3558 
3559  IF (do_symmetric) THEN
3560  IF (iatom <= jatom) THEN
3561  irow = iatom
3562  icol = jatom
3563  ELSE
3564  irow = jatom
3565  icol = iatom
3566  END IF
3567  ELSE
3568  irow = iatom
3569  icol = jatom
3570  END IF
3571 
3572  dab = norm2(rij)
3573 
3574  CALL dbcsr_get_block_p(matrix=t2c(img), &
3575  row=irow, col=icol, block=block_t%block, found=found)
3576  cpassert(found)
3577  trans = do_symmetric .AND. (iatom > jatom)
3578 
3579  DO iset = 1, nseti
3580 
3581  ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3582  sgfi = first_sgf_i(1, iset)
3583 
3584  DO jset = 1, nsetj
3585 
3586  ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3587  sgfj = first_sgf_j(1, jset)
3588 
3589  IF (ncoi*ncoj > 0) THEN
3590  ALLOCATE (sij_contr(nsgfi(iset), nsgfj(jset)))
3591  sij_contr(:, :) = 0.0_dp
3592 
3593  ALLOCATE (sij(ncoi, ncoj))
3594  sij(:, :) = 0.0_dp
3595 
3596  ri = 0.0_dp
3597  rj = rij
3598 
3599  CALL eri_2center(sij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
3600  rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
3601  rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3602 
3603  CALL ab_contract(sij_contr, sij, &
3604  sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3605  ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3606 
3607  DEALLOCATE (sij)
3608  IF (trans) THEN
3609  ALLOCATE (sij_rs(nsgfj(jset), nsgfi(iset)))
3610  sij_rs(:, :) = transpose(sij_contr)
3611  ELSE
3612  ALLOCATE (sij_rs(nsgfi(iset), nsgfj(jset)))
3613  sij_rs(:, :) = sij_contr
3614  END IF
3615 
3616  DEALLOCATE (sij_contr)
3617 
3618  ! RI regularization
3619  IF (.NOT. do_hfx_kpoints_prv) THEN
3620  IF (do_kpoints_prv .AND. cell_j(1) == 0 .AND. cell_j(2) == 0 .AND. cell_j(3) == 0 .AND. &
3621  iatom == jatom .AND. iset == jset) THEN
3622  DO i_diag = 1, nsgfi(iset)
3623  sij_rs(i_diag, i_diag) = sij_rs(i_diag, i_diag) + &
3624  regularization_ri*max(1.0_dp, 1.0_dp/minval(zeti(:, iset)))
3625  END DO
3626  END IF
3627  END IF
3628 
3629  CALL block_add("IN", sij_rs, &
3630  nsgfi(iset), nsgfj(jset), block_t%block, &
3631  sgfi, sgfj, trans=trans)
3632  DEALLOCATE (sij_rs)
3633 
3634  END IF
3635  END DO
3636  END DO
3637  END DO
3638 
3639  CALL cp_libint_cleanup_2eri(lib)
3640 
3641  CALL neighbor_list_iterator_release(nl_iterator)
3642  DO img = 1, nimg
3643  CALL dbcsr_finalize(t2c(img))
3644  CALL dbcsr_filter(t2c(img), filter_eps)
3645  END DO
3646 
3647  CALL timestop(handle)
3648 
3649  END SUBROUTINE build_2c_integrals
3650 
3651 ! **************************************************************************************************
3652 !> \brief ...
3653 !> \param tensor tensor with data. Data is cleared after compression.
3654 !> \param blk_indices ...
3655 !> \param compressed compressed tensor data
3656 !> \param eps all entries < eps are discarded
3657 !> \param memory ...
3658 ! **************************************************************************************************
3659  SUBROUTINE compress_tensor(tensor, blk_indices, compressed, eps, memory)
3660  TYPE(dbt_type), INTENT(INOUT) :: tensor
3661  INTEGER, ALLOCATABLE, DIMENSION(:, :), &
3662  INTENT(INOUT) :: blk_indices
3663  TYPE(hfx_compression_type), INTENT(INOUT) :: compressed
3664  REAL(dp), INTENT(IN) :: eps
3665  REAL(dp), INTENT(INOUT) :: memory
3666 
3667  INTEGER :: buffer_left, buffer_size, buffer_start, &
3668  i, iblk, memory_usage, nbits, nblk, &
3669  nints, offset, shared_offset
3670  INTEGER(int_8) :: estimate_to_store_int, &
3671  storage_counter_integrals
3672  INTEGER, DIMENSION(3) :: ind
3673  LOGICAL :: found
3674  REAL(dp) :: spherical_estimate
3675  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET :: blk_data
3676  REAL(dp), DIMENSION(:), POINTER :: blk_data_1d
3677  TYPE(dbt_iterator_type) :: iter
3678  TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
3679  TYPE(hfx_cache_type), POINTER :: maxval_cache
3680  TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
3681  TYPE(hfx_container_type), POINTER :: maxval_container
3682 
3683  CALL dealloc_containers(compressed, memory_usage)
3684  CALL alloc_containers(compressed, 1)
3685 
3686  maxval_container => compressed%maxval_container(1)
3687  integral_containers => compressed%integral_containers(:, 1)
3688 
3689  CALL hfx_init_container(maxval_container, memory_usage, .false.)
3690  DO i = 1, 64
3691  CALL hfx_init_container(integral_containers(i), memory_usage, .false.)
3692  END DO
3693 
3694  maxval_cache => compressed%maxval_cache(1)
3695  integral_caches => compressed%integral_caches(:, 1)
3696 
3697  IF (ALLOCATED(blk_indices)) DEALLOCATE (blk_indices)
3698  ALLOCATE (blk_indices(dbt_get_num_blocks(tensor), 3))
3699  shared_offset = 0
3700 !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices,shared_offset) &
3701 !$OMP PRIVATE(iter,ind,offset,nblk,iblk)
3702  CALL dbt_iterator_start(iter, tensor)
3703  nblk = dbt_iterator_num_blocks(iter)
3704 !$OMP CRITICAL
3705  offset = shared_offset
3706  shared_offset = shared_offset + nblk
3707 !$OMP END CRITICAL
3708  DO iblk = 1, nblk
3709  CALL dbt_iterator_next_block(iter, ind)
3710  blk_indices(offset + iblk, :) = ind(:)
3711  END DO
3712  CALL dbt_iterator_stop(iter)
3713 !$OMP END PARALLEL
3714 
3715  ! Can not use the tensor iterator here because the order of the blocks is not guaranteed.
3716  DO i = 1, SIZE(blk_indices, 1)
3717  ind = blk_indices(i, :)
3718  CALL dbt_get_block(tensor, ind, blk_data, found)
3719  cpassert(found)
3720  nints = SIZE(blk_data)
3721  blk_data_1d(1:nints) => blk_data
3722  spherical_estimate = maxval(abs(blk_data_1d))
3723  IF (spherical_estimate == 0.0_dp) spherical_estimate = tiny(spherical_estimate)
3724  estimate_to_store_int = exponent(spherical_estimate)
3725  estimate_to_store_int = max(estimate_to_store_int, -15_int_8)
3726 
3727  CALL hfx_add_single_cache_element(estimate_to_store_int, 6, &
3728  maxval_cache, maxval_container, memory_usage, &
3729  .false.)
3730 
3731  spherical_estimate = set_exponent(1.0_dp, estimate_to_store_int + 1)
3732 
3733  nbits = exponent(anint(spherical_estimate/eps)) + 1
3734  IF (nbits > 64) THEN
3735  CALL cp_abort(__location__, &
3736  "Overflow during tensor compression. Please use a larger EPS_FILTER or EPS_STORAGE_SCALING")
3737  END IF
3738 
3739  buffer_left = nints
3740  buffer_start = 1
3741 
3742  DO WHILE (buffer_left > 0)
3743  buffer_size = min(buffer_left, cache_size)
3744  CALL hfx_add_mult_cache_elements(blk_data_1d(buffer_start:), &
3745  buffer_size, nbits, &
3746  integral_caches(nbits), &
3747  integral_containers(nbits), &
3748  eps, 1.0_dp, &
3749  memory_usage, &
3750  .false.)
3751  buffer_left = buffer_left - buffer_size
3752  buffer_start = buffer_start + buffer_size
3753  END DO
3754 
3755  NULLIFY (blk_data_1d); DEALLOCATE (blk_data)
3756  END DO
3757 
3758  CALL dbt_clear(tensor)
3759 
3760  storage_counter_integrals = memory_usage*cache_size
3761  memory = memory + real(storage_counter_integrals, dp)/1024/128
3762  !WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
3763  ! "HFX_MEM_INFO| Total memory consumption ERI's RAM [MiB]: ", memory
3764 
3765  CALL hfx_flush_last_cache(6, maxval_cache, maxval_container, memory_usage, &
3766  .false.)
3767  DO i = 1, 64
3768  CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
3769  memory_usage, .false.)
3770  END DO
3771 
3772  CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .false.)
3773  DO i = 1, 64
3774  CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
3775  memory_usage, .false.)
3776  END DO
3777 
3778  END SUBROUTINE
3779 
3780 ! **************************************************************************************************
3781 !> \brief ...
3782 !> \param tensor empty tensor which is filled by decompressed data
3783 !> \param blk_indices indices of blocks to be reserved
3784 !> \param compressed compressed data
3785 !> \param eps all entries < eps are discarded
3786 ! **************************************************************************************************
3787  SUBROUTINE decompress_tensor(tensor, blk_indices, compressed, eps)
3788 
3789  TYPE(dbt_type), INTENT(INOUT) :: tensor
3790  INTEGER, DIMENSION(:, :) :: blk_indices
3791  TYPE(hfx_compression_type), INTENT(INOUT) :: compressed
3792  REAL(dp), INTENT(IN) :: eps
3793 
3794  INTEGER :: a, b, buffer_left, buffer_size, &
3795  buffer_start, i, memory_usage, nbits, &
3796  nblk_per_thread, nints
3797  INTEGER(int_8) :: estimate_to_store_int
3798  INTEGER, DIMENSION(3) :: blk_size, ind
3799  REAL(dp) :: spherical_estimate
3800  REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: blk_data
3801  REAL(dp), DIMENSION(:, :, :), POINTER :: blk_data_3d
3802  TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
3803  TYPE(hfx_cache_type), POINTER :: maxval_cache
3804  TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
3805  TYPE(hfx_container_type), POINTER :: maxval_container
3806 
3807  maxval_cache => compressed%maxval_cache(1)
3808  maxval_container => compressed%maxval_container(1)
3809  integral_caches => compressed%integral_caches(:, 1)
3810  integral_containers => compressed%integral_containers(:, 1)
3811 
3812  memory_usage = 0
3813 
3814  CALL hfx_decompress_first_cache(6, maxval_cache, maxval_container, memory_usage, .false.)
3815 
3816  DO i = 1, 64
3817  CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), &
3818  memory_usage, .false.)
3819  END DO
3820 
3821 !TODO: Parallelize creation of block list.
3822 !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices) PRIVATE(nblk_per_thread,A,b)
3823  nblk_per_thread = SIZE(blk_indices, 1)/omp_get_num_threads() + 1
3824  a = omp_get_thread_num()*nblk_per_thread + 1
3825  b = min(a + nblk_per_thread, SIZE(blk_indices, 1))
3826  CALL dbt_reserve_blocks(tensor, blk_indices(a:b, :))
3827 !$OMP END PARALLEL
3828 
3829  ! Can not use the tensor iterator here because the order of the blocks is not guaranteed.
3830  DO i = 1, SIZE(blk_indices, 1)
3831  ind = blk_indices(i, :)
3832  CALL dbt_blk_sizes(tensor, ind, blk_size)
3833  nints = product(blk_size)
3834  CALL hfx_get_single_cache_element( &
3835  estimate_to_store_int, 6, &
3836  maxval_cache, maxval_container, memory_usage, &
3837  .false.)
3838 
3839  spherical_estimate = set_exponent(1.0_dp, estimate_to_store_int + 1)
3840 
3841  nbits = exponent(anint(spherical_estimate/eps)) + 1
3842 
3843  buffer_left = nints
3844  buffer_start = 1
3845 
3846  ALLOCATE (blk_data(nints))
3847  DO WHILE (buffer_left > 0)
3848  buffer_size = min(buffer_left, cache_size)
3849  CALL hfx_get_mult_cache_elements(blk_data(buffer_start), &
3850  buffer_size, nbits, &
3851  integral_caches(nbits), &
3852  integral_containers(nbits), &
3853  eps, 1.0_dp, &
3854  memory_usage, &
3855  .false.)
3856  buffer_left = buffer_left - buffer_size
3857  buffer_start = buffer_start + buffer_size
3858  END DO
3859 
3860  blk_data_3d(1:blk_size(1), 1:blk_size(2), 1:blk_size(3)) => blk_data
3861  CALL dbt_put_block(tensor, ind, blk_size, blk_data_3d)
3862  NULLIFY (blk_data_3d); DEALLOCATE (blk_data)
3863  END DO
3864 
3865  CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .false.)
3866  DO i = 1, 64
3867  CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
3868  memory_usage, .false.)
3869  END DO
3870  END SUBROUTINE
3871 
3872 ! **************************************************************************************************
3873 !> \brief ...
3874 !> \param tensor ...
3875 !> \param nze ...
3876 !> \param occ ...
3877 ! **************************************************************************************************
3878  SUBROUTINE get_tensor_occupancy(tensor, nze, occ)
3879  TYPE(dbt_type), INTENT(IN) :: tensor
3880  INTEGER(int_8), INTENT(OUT) :: nze
3881  REAL(dp), INTENT(OUT) :: occ
3882 
3883  INTEGER, DIMENSION(dbt_ndims(tensor)) :: dims
3884 
3885  nze = dbt_get_nze_total(tensor)
3886  CALL dbt_get_info(tensor, nfull_total=dims)
3887  occ = real(nze, dp)/product(real(dims, dp))
3888 
3889  END SUBROUTINE
3890 
3891 END MODULE
imax
static int imax(int x, int y)
Returns the larger of two given integer (missing from the C standard)
Definition: dbm_distribution.c:73
ncoset
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition: grid_common.h:73
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
tensor
struct tensor_ tensor
ai_contraction_sphi
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
Definition: ai_contraction_sphi.F:15
ai_contraction_sphi::libxsmm_abc_contract
subroutine, public libxsmm_abc_contract(abcint, sabc, tsphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer)
3-center contraction routine from primitive cartesain Gaussians to spherical Gaussian functions....
Definition: ai_contraction_sphi.F:247
ai_contraction_sphi::ab_contract
subroutine, public ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
contract overlap integrals (a,b) and transfer to spherical Gaussians
Definition: ai_contraction_sphi.F:52
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition: ai_contraction.F:18
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
basis_set_types
Definition: basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius)
...
Definition: basis_set_types.F:615
block_p_types
collect pointers to a block of reals
Definition: block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition: cell_types.F:486
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition: cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition: dbt_api.F:17
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition: distribution_2d_types.F:15
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
gamma::init_md_ftable
subroutine, public init_md_ftable(nmax)
Initialize a table of F_n(t) values in the range 0 <= t <= 12 with a stepsize of 0....
Definition: gamma.F:540
hfx_compression_methods
routines and types for Hartree-Fock-Exchange
Definition: hfx_compression_methods.F:14
hfx_compression_methods::hfx_add_single_cache_element
subroutine, public hfx_add_single_cache_element(value, nbits, cache, container, memory_usage, use_disk_storage, max_val_memory)
This routine adds an int_8 value to a cache. If the cache is full a compression routine is invoked an...
Definition: hfx_compression_methods.F:72
hfx_compression_methods::hfx_get_mult_cache_elements
subroutine, public hfx_get_mult_cache_elements(values, nints, nbits, cache, container, eps_schwarz, pmax_entry, memory_usage, use_disk_storage)
This routine returns a bunch real values from a cache. If the cache is empty a decompression routine ...
Definition: hfx_compression_methods.F:421
hfx_compression_methods::hfx_flush_last_cache
subroutine, public hfx_flush_last_cache(nbits, cache, container, memory_usage, use_disk_storage)
This routine compresses the last probably not yet compressed cache into a container
Definition: hfx_compression_methods.F:305
hfx_compression_methods::hfx_get_single_cache_element
subroutine, public hfx_get_single_cache_element(value, nbits, cache, container, memory_usage, use_disk_storage)
This routine returns an int_8 value from a cache. If the cache is empty a decompression routine is in...
Definition: hfx_compression_methods.F:168
hfx_compression_methods::hfx_decompress_first_cache
subroutine, public hfx_decompress_first_cache(nbits, cache, container, memory_usage, use_disk_storage)
This routine decompresses the first bunch of data in a container and copies them into a cache
Definition: hfx_compression_methods.F:282
hfx_compression_methods::hfx_add_mult_cache_elements
subroutine, public hfx_add_mult_cache_elements(values, nints, nbits, cache, container, eps_schwarz, pmax_entry, memory_usage, use_disk_storage)
This routine adds an a few real values to a cache. If the cache is full a compression routine is invo...
Definition: hfx_compression_methods.F:341
hfx_compression_methods::hfx_reset_cache_and_container
subroutine, public hfx_reset_cache_and_container(cache, container, memory_usage, do_disk_storage)
This routine resets the containers list pointer to the first element and moves the element counters o...
Definition: hfx_compression_methods.F:251
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
hfx_types::dealloc_containers
subroutine, public dealloc_containers(DATA, memory_usage)
...
Definition: hfx_types.F:2874
hfx_types::alloc_containers
subroutine, public alloc_containers(DATA, bin_size)
...
Definition: hfx_types.F:2906
hfx_types::hfx_init_container
subroutine, public hfx_init_container(container, memory_usage, do_disk_storage)
This routine deletes all list entries in a container in order to deallocate the memory.
Definition: hfx_types.F:2523
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_potential_truncated
integer, parameter, public do_potential_truncated
Definition: input_constants.F:801
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition: input_constants.F:801
input_constants::do_potential_coulomb
integer, parameter, public do_potential_coulomb
Definition: input_constants.F:801
input_constants::do_potential_short
integer, parameter, public do_potential_short
Definition: input_constants.F:801
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition: kinds.F:54
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition: libint_2c_3c.F:14
libint_2c_3c::eri_2center
subroutine, public eri_2center(int_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center electron repulsion integrals (a|b) for a given set of cartesian gaussian orbita...
Definition: libint_2c_3c.F:795
libint_2c_3c::eri_3center
subroutine, public eri_3center(int_abc, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, int_abc_ext)
Computes the 3-center electron repulsion integrals (ab|c) for a given set of cartesian gaussian orbit...
Definition: libint_2c_3c.F:114
libint_2c_3c::cutoff_screen_factor
real(kind=dp), parameter, public cutoff_screen_factor
Definition: libint_2c_3c.F:48
libint_2c_3c::eri_2center_derivs
subroutine, public eri_2center_derivs(der_ab, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, dab, lib, potential_parameter)
Computes the 2-center derivatives of the electron repulsion integrals (a|b) for a given set of cartes...
Definition: libint_2c_3c.F:1005
libint_2c_3c::eri_3center_derivs
subroutine, public eri_3center_derivs(der_abc_1, der_abc_2, la_min, la_max, npgfa, zeta, rpgfa, ra, lb_min, lb_max, npgfb, zetb, rpgfb, rb, lc_min, lc_max, npgfc, zetc, rpgfc, rc, dab, dac, dbc, lib, potential_parameter, der_abc_1_ext, der_abc_2_ext)
Computes the derivatives of the 3-center electron repulsion integrals (ab|c) for a given set of carte...
Definition: libint_2c_3c.F:443
libint_wrapper
Interface to the Libint-Library or a c++ wrapper.
Definition: libint_wrapper.F:15
libint_wrapper::cp_libint_cleanup_3eri1
subroutine, public cp_libint_cleanup_3eri1(lib)
Definition: libint_wrapper.F:1206
libint_wrapper::cp_libint_init_3eri1
subroutine, public cp_libint_init_3eri1(lib, max_am)
Definition: libint_wrapper.F:1140
libint_wrapper::cp_libint_cleanup_2eri1
subroutine, public cp_libint_cleanup_2eri1(lib)
Definition: libint_wrapper.F:1216
libint_wrapper::cp_libint_init_2eri1
subroutine, public cp_libint_init_2eri1(lib, max_am)
Definition: libint_wrapper.F:1152
libint_wrapper::cp_libint_init_2eri
subroutine, public cp_libint_init_2eri(lib, max_am)
Definition: libint_wrapper.F:1164
libint_wrapper::cp_libint_init_3eri
subroutine, public cp_libint_init_3eri(lib, max_am)
Definition: libint_wrapper.F:1128
libint_wrapper::cp_libint_cleanup_3eri
subroutine, public cp_libint_cleanup_3eri(lib)
Definition: libint_wrapper.F:1196
libint_wrapper::cp_libint_set_contrdepth
subroutine, public cp_libint_set_contrdepth(lib, contrdepth)
Definition: libint_wrapper.F:1098
libint_wrapper::cp_libint_cleanup_2eri
subroutine, public cp_libint_cleanup_2eri(lib)
Definition: libint_wrapper.F:1226
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition: molecule_types.F:16
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition: qs_neighbor_list_types.F:1097
qs_neighbor_list_types::nl_sub_iterate
integer function, public nl_sub_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:440
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::get_neighbor_list_set_p
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Definition: qs_neighbor_list_types.F:993
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
qs_neighbor_lists
Generate the atomic neighbor lists.
Definition: qs_neighbor_lists.F:19
qs_neighbor_lists::atom2d_cleanup
subroutine, public atom2d_cleanup(atom2d)
free the internals of atom2d
Definition: qs_neighbor_lists.F:131
qs_neighbor_lists::pair_radius_setup
subroutine, public pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
...
Definition: qs_neighbor_lists.F:1618
qs_neighbor_lists::build_neighbor_lists
subroutine, public build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, mic, symmetric, molecular, subset_of_mol, current_subset, operator_type, nlname, atomb_to_keep)
Build simple pair neighbor lists.
Definition: qs_neighbor_lists.F:1020
qs_neighbor_lists::atom2d_build
subroutine, public atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, molecule_set, molecule_only, particle_set)
Build some distribution structure of atoms, refactored from build_qs_neighbor_lists.
Definition: qs_neighbor_lists.F:173
qs_tensors_types
Utility methods to build 3-center integral tensors of various types.
Definition: qs_tensors_types.F:12
qs_tensors_types::symmetrik_ik
integer, parameter, public symmetrik_ik
Definition: qs_tensors_types.F:52
qs_tensors_types::symmetric_jk
integer, parameter, public symmetric_jk
Definition: qs_tensors_types.F:52
qs_tensors_types::symmetric_ijk
integer, parameter, public symmetric_ijk
Definition: qs_tensors_types.F:52
qs_tensors_types::symmetric_ij
integer, parameter, public symmetric_ij
Definition: qs_tensors_types.F:52
qs_tensors_types::symmetric_none
integer, parameter, public symmetric_none
Definition: qs_tensors_types.F:52
qs_tensors_types::distribution_3d_destroy
subroutine, public distribution_3d_destroy(dist)
Destroy a 3d distribution.
Definition: qs_tensors_types.F:135
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition: qs_tensors.F:11
qs_tensors::calc_3c_virial
subroutine, public calc_3c_virial(work_virial, t3c_trace, pref, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos)
Calculates the 3c virial contributions on the fly.
Definition: qs_tensors.F:1868
qs_tensors::build_2c_neighbor_lists
subroutine, public build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, sym_ij, molecular, dist_2d, pot_to_rad)
Build 2-center neighborlists adapted to different operators This mainly wraps build_neighbor_lists fo...
Definition: qs_tensors.F:143
qs_tensors::build_3c_derivatives
subroutine, public build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos, do_kpoints, do_hfx_kpoints, bounds_i, bounds_j, bounds_k, RI_range, img_to_RI_cell)
Build 3-center derivative tensors.
Definition: qs_tensors.F:1147
qs_tensors::neighbor_list_3c_iterate
recursive integer function, public neighbor_list_3c_iterate(iterator)
Iterate 3c-nl iterator.
Definition: qs_tensors.F:467
qs_tensors::compress_tensor
subroutine, public compress_tensor(tensor, blk_indices, compressed, eps, memory)
...
Definition: qs_tensors.F:3660
qs_tensors::neighbor_list_3c_iterator_destroy
subroutine, public neighbor_list_3c_iterator_destroy(iterator)
Destroy 3c-nl iterator.
Definition: qs_tensors.F:445
qs_tensors::neighbor_list_3c_destroy
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
Definition: qs_tensors.F:383
qs_tensors::calc_2c_virial
subroutine, public calc_2c_virial(work_virial, t2c_trace, pref, qs_env, nl_2c, basis_i, basis_j, potential_parameter)
Calculates the virial coming from 2c derivatives on the fly.
Definition: qs_tensors.F:3226
qs_tensors::decompress_tensor
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
Definition: qs_tensors.F:3788
qs_tensors::build_2c_derivatives
subroutine, public build_2c_derivatives(t2c_der, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints)
Calculates the derivatives of 2-center integrals, wrt to the first center.
Definition: qs_tensors.F:2983
qs_tensors::get_tensor_occupancy
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Definition: qs_tensors.F:3879
qs_tensors::build_3c_neighbor_lists
subroutine, public build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, dist_3d, potential_parameter, name, qs_env, sym_ij, sym_jk, sym_ik, molecular, op_pos, own_dist)
Build a 3-center neighbor list.
Definition: qs_tensors.F:282
qs_tensors::build_2c_integrals
subroutine, public build_2c_integrals(t2c, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints, do_hfx_kpoints, ext_kpoints, regularization_RI)
...
Definition: qs_tensors.F:3406
qs_tensors::build_3c_integrals
subroutine, public build_3c_integrals(t3c, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, int_eps, op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, bounds_i, bounds_j, bounds_k, RI_range, img_to_RI_cell)
Build 3-center integral tensor.
Definition: qs_tensors.F:2376
qs_tensors::neighbor_list_3c_iterator_create
subroutine, public neighbor_list_3c_iterator_create(iterator, ijk_nl)
Create a 3-center neighborlist iterator.
Definition: qs_tensors.F:400
qs_tensors::get_3c_iterator_info
subroutine, public get_3c_iterator_info(iterator, ikind, jkind, kkind, nkind, iatom, jatom, katom, rij, rjk, rik, cell_j, cell_k)
Get info of current iteration.
Definition: qs_tensors.F:564
t_c_g0
This module computes the basic integrals for the truncated coulomb operator.
Definition: t_c_g0.F:57
t_c_g0::init
subroutine, public init(Nder, iunit, mepos, group)
...
Definition: t_c_g0.F:1357
t_c_g0::get_lmax_init
integer function, public get_lmax_init()
Returns the value of nderiv_init so that one can check if opening the potential file is worhtwhile.
Definition: t_c_g0.F:1464
util
All kind of helpful little routines.
Definition: util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333