Line data Source code
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 : abc_contract_xsmm
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 10304 : 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 10304 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
156 : REAL(kind=dp) :: subcells
157 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: i_radius, j_radius
158 10304 : 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 10304 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
163 10304 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
164 10304 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
165 :
166 10304 : NULLIFY (atomic_kind_set, cell, local_particles, molecule_set, &
167 10304 : particle_set, dist_2d_prv)
168 :
169 10304 : IF (PRESENT(pot_to_rad)) THEN
170 1296 : 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 10304 : molecule_set=molecule_set)
183 :
184 10304 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
185 :
186 51520 : ALLOCATE (i_present(nkind), j_present(nkind))
187 51520 : ALLOCATE (i_radius(nkind), j_radius(nkind))
188 :
189 26848 : i_present = .FALSE.
190 26848 : j_present = .FALSE.
191 26848 : i_radius = 0.0_dp
192 26848 : j_radius = 0.0_dp
193 :
194 10304 : IF (PRESENT(dist_2d)) dist_2d_prv => dist_2d
195 :
196 : ! Set up the radii, depending on the operator type
197 10304 : IF (potential_parameter%potential_type == do_potential_id) THEN
198 :
199 : !overlap => use the kind radius for both i and j
200 22276 : DO ikind = 1, nkind
201 13720 : IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
202 13720 : i_present(ikind) = .TRUE.
203 13720 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
204 : END IF
205 22276 : IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
206 13720 : j_present(ikind) = .TRUE.
207 13720 : CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
208 : END IF
209 : END DO
210 :
211 1748 : 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 1086 : DO ikind = 1, nkind
215 724 : IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
216 724 : i_present(ikind) = .TRUE.
217 724 : IF (pot_to_rad_prv == 1) i_radius(ikind) = 1000000.0_dp
218 : END IF
219 1086 : IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
220 724 : j_present(ikind) = .TRUE.
221 724 : IF (pot_to_rad_prv == 2) j_radius(ikind) = 1000000.0_dp
222 : END IF
223 : END DO !ikind
224 :
225 1386 : 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 3486 : DO ikind = 1, nkind
230 2100 : IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
231 2100 : i_present(ikind) = .TRUE.
232 2100 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
233 2100 : IF (pot_to_rad_prv == 1) i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
234 : END IF
235 3486 : IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
236 2100 : j_present(ikind) = .TRUE.
237 2100 : CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
238 2100 : 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 0 : CPABORT("Operator not implemented.")
244 : END IF
245 :
246 41216 : ALLOCATE (pair_radius(nkind, nkind))
247 55908 : pair_radius = 0.0_dp
248 10304 : CALL pair_radius_setup(i_present, j_present, i_radius, j_radius, pair_radius)
249 :
250 30912 : ALLOCATE (atom2d(nkind))
251 :
252 : CALL atom2d_build(atom2d, local_particles, dist_2d_prv, atomic_kind_set, &
253 10304 : molecule_set, molecule_only=.FALSE., particle_set=particle_set)
254 : CALL build_neighbor_lists(ij_list, particle_set, atom2d, cell, pair_radius, subcells, &
255 10304 : symmetric=sym_ij, molecular=molecular, nlname=TRIM(name))
256 :
257 10304 : CALL atom2d_cleanup(atom2d)
258 :
259 30912 : 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 648 : 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 648 : CALL timeset(routineN, handle)
298 :
299 648 : IF (PRESENT(op_pos)) THEN
300 478 : op_pos_prv = op_pos
301 : ELSE
302 : op_pos_prv = 1
303 : END IF
304 :
305 466 : SELECT CASE (op_pos_prv)
306 : CASE (1)
307 466 : pot_par_1 = potential_parameter
308 466 : pot_par_2%potential_type = do_potential_id
309 : CASE (2)
310 182 : pot_par_2 = potential_parameter
311 478 : 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 648 : 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 648 : dist_2d=dist_3d%dist_2d_2, pot_to_rad=2)
321 :
322 648 : ijk_list%sym = symmetric_none
323 :
324 648 : sym_level = 0
325 648 : IF (PRESENT(sym_ij)) THEN
326 140 : IF (sym_ij) THEN
327 0 : ijk_list%sym = symmetric_ij
328 0 : sym_level = sym_level + 1
329 : END IF
330 : END IF
331 :
332 648 : IF (PRESENT(sym_jk)) THEN
333 472 : IF (sym_jk) THEN
334 418 : ijk_list%sym = symmetric_jk
335 418 : sym_level = sym_level + 1
336 : END IF
337 : END IF
338 :
339 648 : IF (PRESENT(sym_ik)) THEN
340 0 : IF (sym_ik) THEN
341 0 : ijk_list%sym = symmetrik_ik
342 0 : sym_level = sym_level + 1
343 : END IF
344 : END IF
345 :
346 648 : IF (sym_level >= 2) THEN
347 0 : ijk_list%sym = symmetric_ijk
348 : END IF
349 :
350 648 : ijk_list%dist_3d = dist_3d
351 648 : IF (PRESENT(own_dist)) THEN
352 648 : ijk_list%owns_dist = own_dist
353 : ELSE
354 0 : ijk_list%owns_dist = .FALSE.
355 : END IF
356 :
357 648 : CALL timestop(handle)
358 648 : END SUBROUTINE
359 :
360 : ! **************************************************************************************************
361 : !> \brief Symmetry criterion
362 : !> \param a ...
363 : !> \param b ...
364 : !> \return ...
365 : ! **************************************************************************************************
366 3048714 : PURE FUNCTION include_symmetric(a, b)
367 : INTEGER, INTENT(IN) :: a, b
368 : LOGICAL :: include_symmetric
369 :
370 3048714 : IF (a > b) THEN
371 1171848 : include_symmetric = (MODULO(a + b, 2) /= 0)
372 : ELSE
373 1876866 : include_symmetric = (MODULO(a + b, 2) == 0)
374 : END IF
375 :
376 3048714 : END FUNCTION
377 :
378 : ! **************************************************************************************************
379 : !> \brief Destroy 3c neighborlist
380 : !> \param ijk_list ...
381 : ! **************************************************************************************************
382 648 : SUBROUTINE neighbor_list_3c_destroy(ijk_list)
383 : TYPE(neighbor_list_3c_type), INTENT(INOUT) :: ijk_list
384 :
385 648 : CALL release_neighbor_list_sets(ijk_list%ij_list)
386 648 : CALL release_neighbor_list_sets(ijk_list%jk_list)
387 :
388 648 : IF (ijk_list%owns_dist) THEN
389 648 : CALL distribution_3d_destroy(ijk_list%dist_3d)
390 : END IF
391 :
392 648 : END SUBROUTINE
393 :
394 : ! **************************************************************************************************
395 : !> \brief Create a 3-center neighborlist iterator
396 : !> \param iterator ...
397 : !> \param ijk_nl ...
398 : ! **************************************************************************************************
399 56826 : 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 6314 : CALL timeset(routineN, handle)
408 :
409 6314 : CALL neighbor_list_iterator_create(iterator%iter_ij, ijk_nl%ij_list)
410 6314 : CALL neighbor_list_iterator_create(iterator%iter_jk, ijk_nl%jk_list, search=.TRUE.)
411 :
412 6314 : iterator%iter_level = 0
413 6314 : iterator%ijk_nl = ijk_nl
414 :
415 18942 : iterator%bounds_i = 0
416 18942 : iterator%bounds_j = 0
417 18942 : iterator%bounds_k = 0
418 :
419 6314 : CALL timestop(handle)
420 6314 : 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 6244 : 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 20476 : IF (PRESENT(bounds_i)) iterator%bounds_i = bounds_i
435 10492 : IF (PRESENT(bounds_j)) iterator%bounds_j = bounds_j
436 13624 : IF (PRESENT(bounds_k)) iterator%bounds_k = bounds_k
437 :
438 6244 : END SUBROUTINE
439 :
440 : ! **************************************************************************************************
441 : !> \brief Destroy 3c-nl iterator
442 : !> \param iterator ...
443 : ! **************************************************************************************************
444 6314 : 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 6314 : CALL timeset(routineN, handle)
453 6314 : CALL neighbor_list_iterator_release(iterator%iter_ij)
454 6314 : CALL neighbor_list_iterator_release(iterator%iter_jk)
455 6314 : NULLIFY (iterator%iter_ij)
456 6314 : NULLIFY (iterator%iter_jk)
457 :
458 6314 : CALL timestop(handle)
459 6314 : 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 6348959 : 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 6348959 : iter_level = iterator%iter_level
476 :
477 6348959 : IF (iter_level == 0) THEN
478 236346 : iter_stat = neighbor_list_iterate(iterator%iter_ij)
479 :
480 236346 : IF (iter_stat /= 0) THEN
481 : RETURN
482 : END IF
483 :
484 230032 : CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom)
485 230032 : skip_this = .FALSE.
486 : IF ((iterator%bounds_i(1) > 0 .AND. iatom < iterator%bounds_i(1)) &
487 230032 : .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 230032 : .OR. (iterator%bounds_j(2) > 0 .AND. jatom > iterator%bounds_j(2))) skip_this = .TRUE.
490 :
491 230032 : IF (skip_this) THEN
492 72744 : iter_stat = neighbor_list_3c_iterate(iterator)
493 72744 : RETURN
494 : END IF
495 :
496 : END IF
497 6269901 : iter_stat = nl_sub_iterate(iterator%iter_jk, iterator%iter_ij)
498 6269901 : IF (iter_stat /= 0) THEN
499 157288 : iterator%iter_level = 0
500 157288 : iter_stat = neighbor_list_3c_iterate(iterator)
501 157288 : RETURN
502 : ELSE
503 6112613 : iterator%iter_level = 1
504 : END IF
505 :
506 : CPASSERT(iter_stat == 0)
507 : CPASSERT(iterator%iter_level == 1)
508 6112613 : CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom_1)
509 6112613 : CALL get_iterator_info(iterator%iter_jk, iatom=jatom_2, jatom=katom)
510 :
511 6112613 : CPASSERT(jatom_1 == jatom_2)
512 6112613 : jatom = jatom_1
513 :
514 6112613 : skip_this = .FALSE.
515 : IF ((iterator%bounds_k(1) > 0 .AND. katom < iterator%bounds_k(1)) &
516 6112613 : .OR. (iterator%bounds_k(2) > 0 .AND. katom > iterator%bounds_k(2))) skip_this = .TRUE.
517 :
518 : IF (skip_this) THEN
519 241296 : iter_stat = neighbor_list_3c_iterate(iterator)
520 241296 : RETURN
521 : END IF
522 :
523 5871317 : SELECT CASE (iterator%ijk_nl%sym)
524 : CASE (symmetric_none)
525 0 : skip_this = .FALSE.
526 : CASE (symmetric_ij)
527 0 : skip_this = .NOT. include_symmetric(iatom, jatom)
528 : CASE (symmetric_jk)
529 3048714 : skip_this = .NOT. include_symmetric(jatom, katom)
530 : CASE (symmetrik_ik)
531 0 : skip_this = .NOT. include_symmetric(iatom, katom)
532 : CASE (symmetric_ijk)
533 0 : skip_this = .NOT. include_symmetric(iatom, jatom) .OR. .NOT. include_symmetric(jatom, katom)
534 : CASE DEFAULT
535 5871317 : CPABORT("should not happen")
536 : END SELECT
537 :
538 3048714 : IF (skip_this) THEN
539 1173987 : iter_stat = neighbor_list_3c_iterate(iterator)
540 1173987 : 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 4697330 : 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 4697330 : 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 4697330 : 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 4697330 : cell=cell_2)
586 :
587 4697330 : IF (PRESENT(ikind)) ikind = kinds_1(1)
588 4697330 : IF (PRESENT(jkind)) jkind = kinds_1(2)
589 4697330 : IF (PRESENT(kkind)) kkind = kinds_2(2)
590 4697330 : IF (PRESENT(iatom)) iatom = atoms_1(1)
591 4697330 : IF (PRESENT(jatom)) jatom = atoms_1(2)
592 4697330 : IF (PRESENT(katom)) katom = atoms_2(2)
593 :
594 4697330 : IF (PRESENT(rij)) rij = r_1
595 4697330 : IF (PRESENT(rjk)) rjk = r_2
596 23486650 : IF (PRESENT(rik)) rik = r_1 + r_2
597 :
598 4697330 : IF (PRESENT(cell_j)) cell_j = cell_1
599 22807450 : IF (PRESENT(cell_k)) cell_k = cell_1 + cell_2
600 :
601 4697330 : 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 2076 : 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 2076 : 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 2076 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:, :) :: nblk
645 2076 : 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 2076 : 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 2076 : 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 2076 : DIMENSION(:, :) :: alloc_i, alloc_j, alloc_k
662 2076 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
663 :
664 2076 : CALL timeset(routineN, handle)
665 2076 : NULLIFY (qs_kind_set, atomic_kind_set, cell)
666 :
667 2076 : IF (PRESENT(do_kpoints)) THEN
668 436 : do_kpoints_prv = do_kpoints
669 : ELSE
670 : do_kpoints_prv = .FALSE.
671 : END IF
672 2076 : IF (PRESENT(do_hfx_kpoints)) THEN
673 196 : do_hfx_kpoints_prv = do_hfx_kpoints
674 : ELSE
675 : do_hfx_kpoints_prv = .FALSE.
676 : END IF
677 196 : IF (do_hfx_kpoints_prv) THEN
678 196 : CPASSERT(do_kpoints_prv)
679 196 : CPASSERT(PRESENT(RI_range))
680 196 : CPASSERT(PRESENT(img_to_RI_cell))
681 : END IF
682 :
683 1880 : IF (PRESENT(img_to_RI_cell)) THEN
684 588 : ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
685 9564 : img_to_RI_cell_prv(:) = img_to_RI_cell
686 : END IF
687 :
688 2076 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
689 :
690 2076 : op_ij = do_potential_id; op_jk = do_potential_id
691 :
692 2076 : IF (PRESENT(op_pos)) THEN
693 2076 : op_pos_prv = op_pos
694 : ELSE
695 : op_pos_prv = 1
696 : END IF
697 :
698 1880 : SELECT CASE (op_pos_prv)
699 : CASE (1)
700 1880 : op_ij = potential_parameter%potential_type
701 : CASE (2)
702 2076 : op_jk = potential_parameter%potential_type
703 : END SELECT
704 :
705 2076 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
706 1054 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
707 1054 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
708 1022 : ELSEIF (op_ij == do_potential_coulomb) THEN
709 436 : dr_ij = 1000000.0_dp
710 436 : dr_ik = 1000000.0_dp
711 : END IF
712 :
713 2076 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
714 20 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
715 20 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
716 2056 : ELSEIF (op_jk == do_potential_coulomb) THEN
717 0 : dr_jk = 1000000.0_dp
718 0 : 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 2076 : dft_control=dft_control, kpoints=kpoints, para_env=para_env, cell=cell)
723 :
724 2076 : IF (do_kpoints_prv) THEN
725 202 : nimg = dft_control%nimages
726 202 : ncell_RI = nimg
727 202 : IF (do_hfx_kpoints_prv) THEN
728 196 : nimg = SIZE(t3c, 1)
729 196 : ncell_RI = SIZE(t3c, 2)
730 : END IF
731 202 : 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 2076 : IF (do_kpoints_prv) THEN
738 808 : kp_index_lbounds = LBOUND(cell_to_index)
739 808 : 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 8304 : ALLOCATE (nblk(nimg, ncell_RI))
744 11080 : nblk(:, :) = 0
745 :
746 2076 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
747 :
748 2076 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
749 :
750 1491748 : 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 1489672 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
753 :
754 5958688 : djk = NORM2(rjk)
755 5958688 : dij = NORM2(rij)
756 5958688 : dik = NORM2(rik)
757 :
758 1489672 : IF (do_kpoints_prv) THEN
759 :
760 603372 : 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 86196 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
764 86196 : IF (jcell > nimg .OR. jcell < 1) CYCLE
765 :
766 507674 : 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 68679 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
770 68679 : IF (kcell > nimg .OR. kcell < 1) CYCLE
771 48533 : IF (do_hfx_kpoints_prv) THEN
772 48071 : IF (dik > RI_range) CYCLE
773 : kcell = 1
774 : END IF
775 : ELSE
776 : jcell = 1; kcell = 1
777 : END IF
778 :
779 1412052 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
780 1412052 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
781 1412052 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
782 :
783 1412052 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
784 1412052 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
785 1412052 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
786 :
787 1484008 : nblk(jcell, kcell) = nblk(jcell, kcell) + 1
788 : END DO
789 2076 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
790 :
791 : !Do a second loop over the nl to give block indices
792 17308 : ALLOCATE (alloc_i(nimg, ncell_RI))
793 15232 : ALLOCATE (alloc_j(nimg, ncell_RI))
794 15232 : ALLOCATE (alloc_k(nimg, ncell_RI))
795 4176 : DO k_img = 1, ncell_RI
796 11080 : DO j_img = 1, nimg
797 18325 : ALLOCATE (alloc_i(j_img, k_img)%array(nblk(j_img, k_img)))
798 18325 : ALLOCATE (alloc_j(j_img, k_img)%array(nblk(j_img, k_img)))
799 20425 : ALLOCATE (alloc_k(j_img, k_img)%array(nblk(j_img, k_img)))
800 : END DO
801 : END DO
802 11080 : nblk(:, :) = 0
803 :
804 2076 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
805 :
806 2076 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
807 :
808 1491748 : 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 1489672 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
812 :
813 5958688 : djk = NORM2(rjk)
814 5958688 : dij = NORM2(rij)
815 5958688 : dik = NORM2(rik)
816 :
817 1489672 : IF (do_kpoints_prv) THEN
818 :
819 603372 : 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 86196 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
823 86196 : IF (jcell > nimg .OR. jcell < 1) CYCLE
824 :
825 507674 : 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 68679 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
829 68679 : IF (kcell > nimg .OR. kcell < 1) CYCLE
830 48533 : IF (do_hfx_kpoints_prv) THEN
831 48071 : IF (dik > RI_range) CYCLE
832 8114 : kcell = img_to_RI_cell_prv(kcell)
833 : END IF
834 : ELSE
835 : jcell = 1; kcell = 1
836 : END IF
837 :
838 5648208 : blk_idx = [iatom, jatom, katom]
839 1412052 : IF (do_hfx_kpoints_prv) THEN
840 8114 : blk_idx(3) = (kcell - 1)*natom + katom
841 8114 : kcell = 1
842 : END IF
843 :
844 1412052 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
845 1412052 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
846 1412052 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
847 :
848 1412052 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
849 1412052 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
850 1412052 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
851 :
852 693109 : 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 693109 : alloc_i(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(1)
856 693109 : alloc_j(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(2)
857 1484008 : alloc_k(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(3)
858 :
859 : END DO
860 2076 : 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 2076 : !$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 2076 : CALL timestop(handle)
881 :
882 41472 : 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 0 : 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 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: tmp
915 : INTEGER, DIMENSION(3) :: cell_j, cell_k, kp_index_lbounds, &
916 : kp_index_ubounds
917 0 : 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 0 : 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 0 : DIMENSION(:, :) :: alloc_i, alloc_j, alloc_k
930 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
931 :
932 0 : CALL timeset(routineN, handle)
933 0 : NULLIFY (qs_kind_set, atomic_kind_set)
934 :
935 0 : IF (PRESENT(do_kpoints)) THEN
936 0 : do_kpoints_prv = do_kpoints
937 : ELSE
938 : do_kpoints_prv = .FALSE.
939 : END IF
940 :
941 0 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
942 :
943 0 : op_ij = do_potential_id; op_jk = do_potential_id
944 :
945 0 : IF (PRESENT(op_pos)) THEN
946 0 : op_pos_prv = op_pos
947 : ELSE
948 : op_pos_prv = 1
949 : END IF
950 :
951 0 : SELECT CASE (op_pos_prv)
952 : CASE (1)
953 0 : op_ij = potential_parameter%potential_type
954 : CASE (2)
955 0 : op_jk = potential_parameter%potential_type
956 : END SELECT
957 :
958 0 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
959 0 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
960 0 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
961 0 : ELSEIF (op_ij == do_potential_coulomb) THEN
962 0 : dr_ij = 1000000.0_dp
963 0 : dr_ik = 1000000.0_dp
964 : END IF
965 :
966 0 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
967 0 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
968 0 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
969 0 : ELSEIF (op_jk == do_potential_coulomb) THEN
970 0 : dr_jk = 1000000.0_dp
971 0 : 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 0 : dft_control=dft_control, kpoints=kpoints, para_env=para_env)
976 :
977 0 : IF (do_kpoints_prv) THEN
978 0 : nimg = dft_control%nimages
979 0 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
980 : ELSE
981 : nimg = 1
982 : END IF
983 :
984 0 : ALLOCATE (alloc_i(nimg, nimg))
985 0 : ALLOCATE (alloc_j(nimg, nimg))
986 0 : ALLOCATE (alloc_k(nimg, nimg))
987 :
988 0 : IF (do_kpoints_prv) THEN
989 0 : kp_index_lbounds = LBOUND(cell_to_index)
990 0 : kp_index_ubounds = UBOUND(cell_to_index)
991 : END IF
992 :
993 0 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
994 0 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
995 0 : 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 0 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
999 :
1000 0 : IF (do_kpoints_prv) THEN
1001 :
1002 0 : 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 0 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1006 0 : IF (jcell > nimg .OR. jcell < 1) CYCLE
1007 :
1008 0 : 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 0 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
1012 0 : IF (kcell > nimg .OR. kcell < 1) CYCLE
1013 : ELSE
1014 : jcell = 1; kcell = 1
1015 : END IF
1016 :
1017 0 : djk = NORM2(rjk)
1018 0 : dij = NORM2(rij)
1019 0 : dik = NORM2(rik)
1020 :
1021 0 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
1022 0 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
1023 0 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
1024 :
1025 0 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
1026 0 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
1027 0 : 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 0 : ASSOCIATE (ai => alloc_i(jcell, kcell))
1035 0 : ASSOCIATE (aj => alloc_j(jcell, kcell))
1036 0 : ASSOCIATE (ak => alloc_k(jcell, kcell))
1037 :
1038 0 : new_block = .TRUE.
1039 0 : IF (ALLOCATED(aj%array)) THEN
1040 0 : DO iblk = 1, SIZE(aj%array)
1041 : IF (ai%array(iblk) == iatom .AND. &
1042 0 : aj%array(iblk) == jatom .AND. &
1043 0 : ak%array(iblk) == katom) THEN
1044 : new_block = .FALSE.
1045 : EXIT
1046 : END IF
1047 : END DO
1048 : END IF
1049 0 : IF (.NOT. new_block) CYCLE
1050 :
1051 0 : IF (ALLOCATED(ai%array)) THEN
1052 0 : blk_cnt = SIZE(ai%array)
1053 0 : ALLOCATE (tmp(blk_cnt))
1054 0 : tmp(:) = ai%array(:)
1055 0 : DEALLOCATE (ai%array)
1056 0 : ALLOCATE (ai%array(blk_cnt + 1))
1057 0 : ai%array(1:blk_cnt) = tmp(:)
1058 0 : ai%array(blk_cnt + 1) = iatom
1059 : ELSE
1060 0 : ALLOCATE (ai%array(1))
1061 0 : ai%array(1) = iatom
1062 : END IF
1063 :
1064 0 : IF (ALLOCATED(aj%array)) THEN
1065 0 : tmp(:) = aj%array(:)
1066 0 : DEALLOCATE (aj%array)
1067 0 : ALLOCATE (aj%array(blk_cnt + 1))
1068 0 : aj%array(1:blk_cnt) = tmp(:)
1069 0 : aj%array(blk_cnt + 1) = jatom
1070 : ELSE
1071 0 : ALLOCATE (aj%array(1))
1072 0 : aj%array(1) = jatom
1073 : END IF
1074 :
1075 0 : IF (ALLOCATED(ak%array)) THEN
1076 0 : tmp(:) = ak%array(:)
1077 0 : DEALLOCATE (ak%array)
1078 0 : ALLOCATE (ak%array(blk_cnt + 1))
1079 0 : ak%array(1:blk_cnt) = tmp(:)
1080 0 : ak%array(blk_cnt + 1) = katom
1081 : ELSE
1082 0 : ALLOCATE (ak%array(1))
1083 0 : ak%array(1) = katom
1084 : END IF
1085 :
1086 0 : IF (ALLOCATED(tmp)) DEALLOCATE (tmp)
1087 : END ASSOCIATE
1088 : END ASSOCIATE
1089 : END ASSOCIATE
1090 : END DO
1091 :
1092 0 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
1093 :
1094 0 : DO i_img = 1, nimg
1095 0 : DO j_img = 1, nimg
1096 0 : IF (ALLOCATED(alloc_i(i_img, j_img)%array)) THEN
1097 0 : 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 0 : iproc)
1102 0 : 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 0 : alloc_k(i_img, j_img)%array)
1109 : END IF
1110 : END DO
1111 : END DO
1112 :
1113 0 : CALL timestop(handle)
1114 :
1115 0 : 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 562 : SUBROUTINE build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, &
1142 562 : 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 562 : 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 562 : 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 562 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
1174 1124 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
1175 562 : nsgfj, nsgfk
1176 562 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
1177 562 : 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 562 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
1187 562 : max_contraction_i, max_contraction_j, &
1188 562 : max_contraction_k
1189 562 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dijk_contr, dummy_block_t
1190 562 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: block_t_i, block_t_j, block_t_k, dijk_j, &
1191 1124 : dijk_k, tmp_ijk_i, tmp_ijk_j
1192 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
1193 562 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
1194 562 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
1195 1124 : sphi_k, zeti, zetj, zetk
1196 562 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1197 1124 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
1198 : TYPE(cp_libint_t) :: lib
1199 3934 : TYPE(dbt_type) :: t3c_tmp
1200 562 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t3c_template
1201 562 : 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 562 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1208 :
1209 562 : CALL timeset(routineN, handle)
1210 :
1211 562 : IF (PRESENT(do_kpoints)) THEN
1212 74 : do_kpoints_prv = do_kpoints
1213 : ELSE
1214 : do_kpoints_prv = .FALSE.
1215 : END IF
1216 :
1217 562 : IF (PRESENT(do_hfx_kpoints)) THEN
1218 74 : do_hfx_kpoints_prv = do_hfx_kpoints
1219 : ELSE
1220 : do_hfx_kpoints_prv = .FALSE.
1221 : END IF
1222 74 : IF (do_hfx_kpoints_prv) THEN
1223 74 : CPASSERT(do_kpoints_prv)
1224 74 : CPASSERT(PRESENT(RI_range))
1225 74 : CPASSERT(PRESENT(img_to_RI_cell))
1226 : END IF
1227 :
1228 488 : IF (PRESENT(img_to_RI_cell)) THEN
1229 222 : ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
1230 3586 : img_to_RI_cell_prv(:) = img_to_RI_cell
1231 : END IF
1232 :
1233 562 : op_ij = do_potential_id; op_jk = do_potential_id
1234 :
1235 562 : IF (PRESENT(op_pos)) THEN
1236 562 : op_pos_prv = op_pos
1237 : ELSE
1238 0 : op_pos_prv = 1
1239 : END IF
1240 :
1241 488 : SELECT CASE (op_pos_prv)
1242 : CASE (1)
1243 488 : op_ij = potential_parameter%potential_type
1244 : CASE (2)
1245 562 : op_jk = potential_parameter%potential_type
1246 : END SELECT
1247 :
1248 562 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
1249 :
1250 562 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
1251 50 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
1252 50 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1253 512 : ELSEIF (op_ij == do_potential_coulomb) THEN
1254 244 : dr_ij = 1000000.0_dp
1255 244 : dr_ik = 1000000.0_dp
1256 : END IF
1257 :
1258 562 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
1259 8 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
1260 8 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1261 554 : ELSEIF (op_jk == do_potential_coulomb) THEN
1262 0 : dr_jk = 1000000.0_dp
1263 0 : dr_ik = 1000000.0_dp
1264 : END IF
1265 :
1266 562 : 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 562 : natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
1271 :
1272 562 : IF (do_kpoints_prv) THEN
1273 74 : nimg = dft_control%nimages
1274 74 : ncell_RI = nimg
1275 74 : IF (do_hfx_kpoints_prv) THEN
1276 74 : nimg = SIZE(t3c_der_k, 1)
1277 74 : ncell_RI = SIZE(t3c_der_k, 2)
1278 : END IF
1279 74 : 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 562 : IF (do_hfx_kpoints_prv) THEN
1286 74 : CPASSERT(op_pos_prv == 2)
1287 : ELSE
1288 1952 : CPASSERT(ALL(SHAPE(t3c_der_k) == [nimg, ncell_RI, 3]))
1289 : END IF
1290 :
1291 8446 : ALLOCATE (t3c_template(nimg, ncell_RI))
1292 1124 : DO j_img = 1, ncell_RI
1293 3388 : DO i_img = 1, nimg
1294 2826 : 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 1050 : RI_range=RI_range, img_to_RI_cell=img_to_RI_cell)
1302 2248 : DO i_xyz = 1, 3
1303 3934 : DO j_img = 1, ncell_RI
1304 10164 : DO i_img = 1, nimg
1305 6792 : CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_i(i_img, j_img, i_xyz))
1306 8478 : 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 1124 : DO j_img = 1, ncell_RI
1312 3388 : DO i_img = 1, nimg
1313 2826 : CALL dbt_destroy(t3c_template(i_img, j_img))
1314 : END DO
1315 : END DO
1316 2826 : DEALLOCATE (t3c_template)
1317 :
1318 562 : IF (nl_3c%sym == symmetric_jk) THEN
1319 9760 : ALLOCATE (t3c_der_j(nimg, ncell_RI, 3))
1320 1952 : DO i_xyz = 1, 3
1321 3416 : DO j_img = 1, ncell_RI
1322 4392 : DO i_img = 1, nimg
1323 1464 : CALL dbt_create(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
1324 2928 : 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 562 : nbasis = SIZE(basis_i)
1332 562 : max_nsgfi = 0
1333 562 : max_ncoi = 0
1334 562 : max_nset = 0
1335 562 : maxli = 0
1336 1614 : DO ibasis = 1, nbasis
1337 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
1338 1052 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
1339 1052 : maxli = MAX(maxli, imax)
1340 1052 : max_nset = MAX(max_nset, iset)
1341 5178 : max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
1342 6792 : max_ncoi = MAX(max_ncoi, MAXVAL(npgfi)*ncoset(maxli))
1343 : END DO
1344 : max_nsgfj = 0
1345 : max_ncoj = 0
1346 : maxlj = 0
1347 1614 : DO ibasis = 1, nbasis
1348 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
1349 1052 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
1350 1052 : maxlj = MAX(maxlj, imax)
1351 1052 : max_nset = MAX(max_nset, jset)
1352 3870 : max_nsgfj = MAX(max_nsgfj, MAXVAL(nsgfj))
1353 5484 : max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
1354 : END DO
1355 : max_nsgfk = 0
1356 : max_ncok = 0
1357 : maxlk = 0
1358 1614 : DO ibasis = 1, nbasis
1359 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
1360 1052 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
1361 1052 : maxlk = MAX(maxlk, imax)
1362 1052 : max_nset = MAX(max_nset, kset)
1363 3670 : max_nsgfk = MAX(max_nsgfk, MAXVAL(nsgfk))
1364 5284 : max_ncok = MAX(max_ncok, MAXVAL(npgfk)*ncoset(maxlk))
1365 : END DO
1366 562 : 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 562 : NULLIFY (tspj, spi, spk)
1372 26396 : ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
1373 :
1374 1614 : DO ibasis = 1, nbasis
1375 7862 : DO iset = 1, max_nset
1376 6248 : NULLIFY (spi(iset, ibasis)%array)
1377 6248 : NULLIFY (tspj(iset, ibasis)%array)
1378 :
1379 7300 : NULLIFY (spk(iset, ibasis)%array)
1380 : END DO
1381 : END DO
1382 :
1383 2248 : DO ilist = 1, 3
1384 5404 : DO ibasis = 1, nbasis
1385 3156 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
1386 3156 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
1387 3156 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
1388 :
1389 14404 : DO iset = 1, basis_set%nset
1390 :
1391 9562 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
1392 9562 : sgfi = basis_set%first_sgf(1, iset)
1393 9562 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
1394 :
1395 12718 : IF (ilist == 1) THEN
1396 16504 : ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1397 3073846 : spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1398 :
1399 5436 : ELSE IF (ilist == 2) THEN
1400 11272 : ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
1401 176530 : tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
1402 :
1403 : ELSE
1404 10472 : ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1405 833426 : 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 562 : IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN
1414 :
1415 58 : IF (m_max > get_lmax_init()) THEN
1416 8 : IF (para_env%mepos == 0) THEN
1417 4 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
1418 : END IF
1419 8 : CALL init(m_max, unit_id, para_env%mepos, para_env)
1420 8 : IF (para_env%mepos == 0) THEN
1421 4 : CALL close_file(unit_id)
1422 : END IF
1423 : END IF
1424 : END IF
1425 :
1426 562 : CALL init_md_ftable(nmax=m_max)
1427 :
1428 562 : IF (do_kpoints_prv) THEN
1429 296 : kp_index_lbounds = LBOUND(cell_to_index)
1430 296 : kp_index_ubounds = UBOUND(cell_to_index)
1431 : END IF
1432 :
1433 562 : nthread = 1
1434 562 : !$ 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 562 : !$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 abc_contract_xsmm(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, prefac)
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 : 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 abc_contract_xsmm(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, prefac)
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 : 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 562 : IF (do_kpoints_prv .AND. .NOT. do_hfx_kpoints_prv) THEN
1782 0 : DO i_xyz = 1, 3
1783 0 : DO kcell = 1, nimg
1784 0 : DO jcell = 1, nimg
1785 : ! need half of filter eps because afterwards we add transposed tensor
1786 0 : CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps/2)
1787 0 : CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps/2)
1788 : END DO
1789 : END DO
1790 : END DO
1791 :
1792 562 : 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 488 : CALL dbt_create(t3c_der_k(1, 1, 1), t3c_tmp)
1795 1952 : DO i_xyz = 1, 3
1796 3416 : DO kcell = 1, nimg
1797 4392 : DO jcell = 1, nimg
1798 : CALL dbt_copy(t3c_der_j(jcell, kcell, i_xyz), t3c_der_k(jcell, kcell, i_xyz), &
1799 1464 : order=[1, 3, 2], move_data=.TRUE., summation=.TRUE.)
1800 1464 : CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
1801 :
1802 1464 : 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 1464 : order=[1, 3, 2], move_data=.TRUE., summation=.TRUE.)
1805 2928 : CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
1806 : END DO
1807 : END DO
1808 : END DO
1809 488 : CALL dbt_destroy(t3c_tmp)
1810 :
1811 74 : ELSEIF (nl_3c%sym == symmetric_none) THEN
1812 296 : DO i_xyz = 1, 3
1813 518 : DO kcell = 1, ncell_RI
1814 5772 : DO jcell = 1, nimg
1815 5328 : CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
1816 5550 : CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
1817 : END DO
1818 : END DO
1819 : END DO
1820 : ELSE
1821 0 : CPABORT("requested symmetric case not implemented")
1822 : END IF
1823 :
1824 562 : IF (nl_3c%sym == symmetric_jk) THEN
1825 1952 : DO i_xyz = 1, 3
1826 3416 : DO j_img = 1, nimg
1827 4392 : DO i_img = 1, nimg
1828 2928 : 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 3846 : DO iset = 1, max_nset
1835 10094 : DO ibasis = 1, nbasis
1836 6248 : IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
1837 6248 : IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
1838 :
1839 9532 : IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
1840 : END DO
1841 : END DO
1842 :
1843 562 : DEALLOCATE (spi, tspj, spk)
1844 :
1845 562 : CALL timestop(handle)
1846 6522 : 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 16 : SUBROUTINE calc_3c_virial(work_virial, t3c_trace, pref, qs_env, &
1866 16 : 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 16 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
1889 32 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
1890 16 : nsgfj, nsgfk
1891 16 : 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 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
1901 16 : max_contraction_i, max_contraction_j, &
1902 16 : max_contraction_k
1903 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ablock, dijk_contr, tmp_block
1904 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: block_t_i, block_t_j, block_t_k, dijk_j, &
1905 16 : dijk_k
1906 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk, scoord
1907 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
1908 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
1909 16 : sphi_k, zeti, zetj, zetk
1910 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1911 : TYPE(cell_type), POINTER :: cell
1912 16 : 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 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1919 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1920 :
1921 16 : CALL timeset(routineN, handle)
1922 :
1923 16 : op_ij = do_potential_id; op_jk = do_potential_id
1924 :
1925 16 : IF (PRESENT(op_pos)) THEN
1926 16 : op_pos_prv = op_pos
1927 : ELSE
1928 : op_pos_prv = 1
1929 : END IF
1930 16 : CPASSERT(op_pos == 1)
1931 16 : CPASSERT(.NOT. nl_3c%sym == symmetric_jk)
1932 :
1933 16 : SELECT CASE (op_pos_prv)
1934 : CASE (1)
1935 16 : op_ij = potential_parameter%potential_type
1936 : CASE (2)
1937 16 : op_jk = potential_parameter%potential_type
1938 : END SELECT
1939 :
1940 16 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
1941 :
1942 16 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
1943 6 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
1944 6 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1945 10 : ELSEIF (op_ij == do_potential_coulomb) THEN
1946 0 : dr_ij = 1000000.0_dp
1947 0 : dr_ik = 1000000.0_dp
1948 : END IF
1949 :
1950 16 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
1951 0 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
1952 0 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1953 16 : ELSEIF (op_jk == do_potential_coulomb) THEN
1954 0 : dr_jk = 1000000.0_dp
1955 0 : dr_ik = 1000000.0_dp
1956 : END IF
1957 :
1958 16 : 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 16 : 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 16 : nbasis = SIZE(basis_i)
1967 16 : max_nsgfi = 0
1968 16 : max_ncoi = 0
1969 16 : max_nset = 0
1970 16 : maxli = 0
1971 48 : DO ibasis = 1, nbasis
1972 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
1973 32 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
1974 32 : maxli = MAX(maxli, imax)
1975 32 : max_nset = MAX(max_nset, iset)
1976 166 : max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
1977 214 : max_ncoi = MAX(max_ncoi, MAXVAL(npgfi)*ncoset(maxli))
1978 : END DO
1979 : max_nsgfj = 0
1980 : max_ncoj = 0
1981 : maxlj = 0
1982 48 : DO ibasis = 1, nbasis
1983 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
1984 32 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
1985 32 : maxlj = MAX(maxlj, imax)
1986 32 : max_nset = MAX(max_nset, jset)
1987 122 : max_nsgfj = MAX(max_nsgfj, MAXVAL(nsgfj))
1988 170 : max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
1989 : END DO
1990 : max_nsgfk = 0
1991 : max_ncok = 0
1992 : maxlk = 0
1993 48 : DO ibasis = 1, nbasis
1994 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
1995 32 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
1996 32 : maxlk = MAX(maxlk, imax)
1997 32 : max_nset = MAX(max_nset, kset)
1998 122 : max_nsgfk = MAX(max_nsgfk, MAXVAL(nsgfk))
1999 170 : max_ncok = MAX(max_ncok, MAXVAL(npgfk)*ncoset(maxlk))
2000 : END DO
2001 16 : 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 16 : NULLIFY (tspj, spi, spk)
2007 920 : ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
2008 :
2009 48 : DO ibasis = 1, nbasis
2010 280 : DO iset = 1, max_nset
2011 232 : NULLIFY (spi(iset, ibasis)%array)
2012 232 : NULLIFY (tspj(iset, ibasis)%array)
2013 :
2014 264 : NULLIFY (spk(iset, ibasis)%array)
2015 : END DO
2016 : END DO
2017 :
2018 64 : DO ilist = 1, 3
2019 160 : DO ibasis = 1, nbasis
2020 96 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
2021 96 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
2022 96 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
2023 :
2024 458 : DO iset = 1, basis_set%nset
2025 :
2026 314 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
2027 314 : sgfi = basis_set%first_sgf(1, iset)
2028 314 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
2029 :
2030 410 : IF (ilist == 1) THEN
2031 536 : ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2032 177990 : spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2033 :
2034 180 : ELSE IF (ilist == 2) THEN
2035 360 : ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
2036 5194 : tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
2037 :
2038 : ELSE
2039 360 : ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2040 4442 : 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 16 : IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN
2049 :
2050 6 : IF (m_max > get_lmax_init()) THEN
2051 0 : IF (para_env%mepos == 0) THEN
2052 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
2053 : END IF
2054 0 : CALL init(m_max, unit_id, para_env%mepos, para_env)
2055 0 : IF (para_env%mepos == 0) THEN
2056 0 : CALL close_file(unit_id)
2057 : END IF
2058 : END IF
2059 : END IF
2060 :
2061 16 : CALL init_md_ftable(nmax=m_max)
2062 :
2063 16 : nthread = 1
2064 16 : !$ 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 16 : !$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 abc_contract_xsmm(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 abc_contract_xsmm(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 132 : DO iset = 1, max_nset
2331 364 : DO ibasis = 1, nbasis
2332 232 : IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
2333 232 : IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
2334 :
2335 348 : IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
2336 : END DO
2337 : END DO
2338 :
2339 16 : DEALLOCATE (spi, tspj, spk)
2340 :
2341 16 : CALL timestop(handle)
2342 128 : 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 1514 : SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
2371 1514 : 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 1514 : 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 1514 : 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 1514 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
2403 3028 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
2404 1514 : nsgfj, nsgfk
2405 1514 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
2406 1514 : 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 1514 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
2415 1514 : max_contraction_j, max_contraction_k
2416 3028 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_t, dummy_block_t, sijk, &
2417 3028 : 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 1514 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
2421 1514 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
2422 3028 : sphi_k, zeti, zetj, zetk
2423 1514 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2424 : TYPE(cell_type), POINTER :: cell
2425 3028 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
2426 : TYPE(cp_libint_t) :: lib
2427 10598 : 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 1514 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2434 :
2435 1514 : CALL timeset(routineN, handle)
2436 :
2437 1514 : debug = .FALSE.
2438 :
2439 1514 : IF (PRESENT(do_kpoints)) THEN
2440 362 : do_kpoints_prv = do_kpoints
2441 : ELSE
2442 : do_kpoints_prv = .FALSE.
2443 : END IF
2444 :
2445 1514 : IF (PRESENT(do_hfx_kpoints)) THEN
2446 122 : do_hfx_kpoints_prv = do_hfx_kpoints
2447 : ELSE
2448 : do_hfx_kpoints_prv = .FALSE.
2449 : END IF
2450 122 : IF (do_hfx_kpoints_prv) THEN
2451 122 : CPASSERT(do_kpoints_prv)
2452 122 : CPASSERT(PRESENT(RI_range))
2453 122 : CPASSERT(PRESENT(img_to_RI_cell))
2454 : END IF
2455 :
2456 1392 : IF (PRESENT(img_to_RI_cell)) THEN
2457 366 : ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
2458 5978 : img_to_RI_cell_prv(:) = img_to_RI_cell
2459 : END IF
2460 :
2461 1514 : IF (PRESENT(desymmetrize)) THEN
2462 1514 : desymmetrize_prv = desymmetrize
2463 : ELSE
2464 : desymmetrize_prv = .TRUE.
2465 : END IF
2466 :
2467 1514 : op_ij = do_potential_id; op_jk = do_potential_id
2468 :
2469 1514 : IF (PRESENT(op_pos)) THEN
2470 438 : op_pos_prv = op_pos
2471 : ELSE
2472 1076 : op_pos_prv = 1
2473 : END IF
2474 :
2475 1392 : SELECT CASE (op_pos_prv)
2476 : CASE (1)
2477 1392 : op_ij = potential_parameter%potential_type
2478 : CASE (2)
2479 1514 : op_jk = potential_parameter%potential_type
2480 : END SELECT
2481 :
2482 1514 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
2483 :
2484 1514 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN
2485 1004 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
2486 1004 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
2487 510 : ELSEIF (op_ij == do_potential_coulomb) THEN
2488 192 : dr_ij = 1000000.0_dp
2489 192 : dr_ik = 1000000.0_dp
2490 : END IF
2491 :
2492 1514 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short) THEN
2493 12 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
2494 12 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
2495 1502 : ELSEIF (op_jk == do_potential_coulomb) THEN
2496 0 : dr_jk = 1000000.0_dp
2497 0 : dr_ik = 1000000.0_dp
2498 : END IF
2499 :
2500 1514 : 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 1514 : 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 2906 : RI_range=RI_range, img_to_RI_cell=img_to_RI_cell)
2510 :
2511 1514 : IF (do_kpoints_prv) THEN
2512 128 : nimg = dft_control%nimages
2513 128 : ncell_RI = nimg
2514 128 : IF (do_hfx_kpoints_prv) THEN
2515 122 : nimg = SIZE(t3c, 1)
2516 122 : ncell_RI = SIZE(t3c, 2)
2517 : END IF
2518 128 : 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 1514 : IF (do_hfx_kpoints_prv) THEN
2525 122 : CPASSERT(op_pos_prv == 2)
2526 122 : CPASSERT(.NOT. desymmetrize_prv)
2527 1392 : ELSE IF (do_kpoints_prv) THEN
2528 18 : 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 1514 : nbasis = SIZE(basis_i)
2533 1514 : max_nsgfi = 0
2534 1514 : max_ncoi = 0
2535 1514 : max_nset = 0
2536 1514 : maxli = 0
2537 4414 : DO ibasis = 1, nbasis
2538 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
2539 2900 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
2540 2900 : maxli = MAX(maxli, imax)
2541 2900 : max_nset = MAX(max_nset, iset)
2542 12366 : max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
2543 16780 : max_ncoi = MAX(max_ncoi, MAXVAL(npgfi)*ncoset(maxli))
2544 : END DO
2545 : max_nsgfj = 0
2546 : max_ncoj = 0
2547 : maxlj = 0
2548 4414 : DO ibasis = 1, nbasis
2549 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
2550 2900 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
2551 2900 : maxlj = MAX(maxlj, imax)
2552 2900 : max_nset = MAX(max_nset, jset)
2553 9868 : max_nsgfj = MAX(max_nsgfj, MAXVAL(nsgfj))
2554 14282 : max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
2555 : END DO
2556 : max_nsgfk = 0
2557 : max_ncok = 0
2558 : maxlk = 0
2559 4414 : DO ibasis = 1, nbasis
2560 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
2561 2900 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
2562 2900 : maxlk = MAX(maxlk, imax)
2563 2900 : max_nset = MAX(max_nset, kset)
2564 9568 : max_nsgfk = MAX(max_nsgfk, MAXVAL(nsgfk))
2565 13982 : max_ncok = MAX(max_ncok, MAXVAL(npgfk)*ncoset(maxlk))
2566 : END DO
2567 1514 : 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 1514 : NULLIFY (tspj, spi, spk)
2573 66184 : ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
2574 :
2575 4414 : DO ibasis = 1, nbasis
2576 19538 : DO iset = 1, max_nset
2577 15124 : NULLIFY (spi(iset, ibasis)%array)
2578 15124 : NULLIFY (tspj(iset, ibasis)%array)
2579 :
2580 18024 : NULLIFY (spk(iset, ibasis)%array)
2581 : END DO
2582 : END DO
2583 :
2584 6056 : DO ilist = 1, 3
2585 14756 : DO ibasis = 1, nbasis
2586 8700 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
2587 8700 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
2588 8700 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
2589 :
2590 36344 : DO iset = 1, basis_set%nset
2591 :
2592 23102 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
2593 23102 : sgfi = basis_set%first_sgf(1, iset)
2594 23102 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
2595 :
2596 31802 : IF (ilist == 1) THEN
2597 37864 : ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2598 3217162 : spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2599 :
2600 13636 : ELSE IF (ilist == 2) THEN
2601 27872 : ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
2602 428368 : tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
2603 :
2604 : ELSE
2605 26672 : ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2606 1491100 : 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 1514 : IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN
2615 :
2616 1010 : IF (m_max > get_lmax_init()) THEN
2617 76 : IF (para_env%mepos == 0) THEN
2618 38 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
2619 : END IF
2620 76 : CALL init(m_max, unit_id, para_env%mepos, para_env)
2621 76 : IF (para_env%mepos == 0) THEN
2622 38 : CALL close_file(unit_id)
2623 : END IF
2624 : END IF
2625 : END IF
2626 :
2627 1514 : CALL init_md_ftable(nmax=m_max)
2628 :
2629 1514 : IF (do_kpoints_prv) THEN
2630 512 : kp_index_lbounds = LBOUND(cell_to_index)
2631 512 : kp_index_ubounds = UBOUND(cell_to_index)
2632 : END IF
2633 :
2634 6056 : counter = 1.0_dp
2635 :
2636 1514 : nthread = 1
2637 1514 : !$ 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 1514 : !$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 abc_contract_xsmm(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, prefac)
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 : 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 1514 : IF (nl_3c%sym == symmetric_jk .OR. do_kpoints_prv) THEN
2904 :
2905 678 : IF (.NOT. do_hfx_kpoints_prv) THEN
2906 1136 : DO kcell = 1, nimg
2907 1836 : DO jcell = 1, nimg
2908 : ! need half of filter eps because afterwards we add transposed tensor
2909 1280 : CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
2910 : END DO
2911 : END DO
2912 : ELSE
2913 244 : DO kcell = 1, ncell_RI
2914 3348 : DO jcell = 1, nimg
2915 3226 : CALL dbt_filter(t3c(jcell, kcell), filter_eps)
2916 : END DO
2917 : END DO
2918 : END IF
2919 :
2920 678 : IF (desymmetrize_prv) THEN
2921 : ! add transposed of overlap integrals
2922 0 : CALL dbt_create(t3c(1, 1), t_3c_tmp)
2923 0 : DO kcell = 1, jcell
2924 0 : DO jcell = 1, nimg
2925 0 : CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
2926 0 : CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.)
2927 0 : CALL dbt_filter(t3c(kcell, jcell), filter_eps)
2928 : END DO
2929 : END DO
2930 0 : DO kcell = jcell + 1, nimg
2931 0 : DO jcell = 1, nimg
2932 0 : CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
2933 0 : CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.)
2934 0 : CALL dbt_filter(t3c(kcell, jcell), filter_eps)
2935 : END DO
2936 : END DO
2937 0 : CALL dbt_destroy(t_3c_tmp)
2938 : END IF
2939 836 : ELSEIF (nl_3c%sym == symmetric_ij) THEN
2940 0 : DO kcell = 1, nimg
2941 0 : DO jcell = 1, nimg
2942 0 : CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
2943 : END DO
2944 : END DO
2945 836 : ELSEIF (nl_3c%sym == symmetric_none) THEN
2946 1672 : DO kcell = 1, nimg
2947 2508 : DO jcell = 1, nimg
2948 1672 : CALL dbt_filter(t3c(jcell, kcell), filter_eps)
2949 : END DO
2950 : END DO
2951 : ELSE
2952 0 : CPABORT("requested symmetric case not implemented")
2953 : END IF
2954 :
2955 9364 : DO iset = 1, max_nset
2956 24488 : DO ibasis = 1, nbasis
2957 15124 : IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
2958 15124 : IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
2959 :
2960 22974 : IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
2961 : END DO
2962 : END DO
2963 1514 : DEALLOCATE (spi, tspj, spk)
2964 :
2965 1514 : CALL timestop(handle)
2966 13626 : 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 372 : SUBROUTINE build_2c_derivatives(t2c_der, filter_eps, qs_env, &
2981 372 : 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 372 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3001 372 : npgfj, nsgfi, nsgfj
3002 372 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3003 372 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3004 : LOGICAL :: do_kpoints_prv, do_symmetric, found, &
3005 : trans
3006 : REAL(KIND=dp) :: dab
3007 372 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij_contr, dij_rs
3008 372 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dij
3009 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj
3010 372 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3011 372 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3012 372 : zetj
3013 372 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3014 1488 : 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 372 : DIMENSION(:), POINTER :: nl_iterator
3021 372 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3022 :
3023 372 : CALL timeset(routineN, handle)
3024 :
3025 372 : IF (PRESENT(do_kpoints)) THEN
3026 84 : do_kpoints_prv = do_kpoints
3027 : ELSE
3028 : do_kpoints_prv = .FALSE.
3029 : END IF
3030 :
3031 372 : op_prv = potential_parameter%potential_type
3032 :
3033 372 : 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 372 : natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
3038 :
3039 372 : IF (do_kpoints_prv) THEN
3040 84 : nimg = SIZE(t2c_der, 1)
3041 84 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
3042 336 : kp_index_lbounds = LBOUND(cell_to_index)
3043 336 : kp_index_ubounds = UBOUND(cell_to_index)
3044 : ELSE
3045 : nimg = 1
3046 : END IF
3047 :
3048 : ! check for symmetry
3049 372 : CPASSERT(SIZE(nl_2c) > 0)
3050 372 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3051 :
3052 372 : IF (do_symmetric) THEN
3053 576 : DO img = 1, nimg
3054 : !Derivtive matrix is assymetric
3055 1440 : DO i_xyz = 1, 3
3056 1152 : CPASSERT(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_antisymmetric)
3057 : END DO
3058 : END DO
3059 : ELSE
3060 1960 : DO img = 1, nimg
3061 7588 : DO i_xyz = 1, 3
3062 7504 : 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 2536 : DO img = 1, nimg
3068 9028 : DO i_xyz = 1, 3
3069 8656 : CALL cp_dbcsr_alloc_block_from_nbl(t2c_der(img, i_xyz), nl_2c)
3070 : END DO
3071 : END DO
3072 :
3073 372 : maxli = 0
3074 1062 : DO ibasis = 1, SIZE(basis_i)
3075 690 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3076 1062 : maxli = MAX(maxli, imax)
3077 : END DO
3078 372 : maxlj = 0
3079 1062 : DO ibasis = 1, SIZE(basis_j)
3080 690 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3081 1062 : maxlj = MAX(maxlj, imax)
3082 : END DO
3083 :
3084 372 : m_max = maxli + maxlj + 1
3085 :
3086 : !Init the truncated Coulomb operator
3087 372 : IF (op_prv == do_potential_truncated) THEN
3088 :
3089 54 : IF (m_max > get_lmax_init()) THEN
3090 18 : IF (para_env%mepos == 0) THEN
3091 9 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3092 : END IF
3093 18 : CALL init(m_max, unit_id, para_env%mepos, para_env)
3094 18 : IF (para_env%mepos == 0) THEN
3095 9 : CALL close_file(unit_id)
3096 : END IF
3097 : END IF
3098 : END IF
3099 :
3100 372 : CALL init_md_ftable(nmax=m_max)
3101 :
3102 372 : CALL cp_libint_init_2eri1(lib, MAX(maxli, maxlj))
3103 372 : CALL cp_libint_set_contrdepth(lib, 1)
3104 :
3105 372 : CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3106 14603 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3107 :
3108 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3109 14231 : iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
3110 14231 : IF (do_kpoints_prv) THEN
3111 10010 : 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 1430 : img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3114 1430 : 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 14203 : 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 14203 : sphi=sphi_j, zet=zetj)
3126 :
3127 14203 : IF (do_symmetric) THEN
3128 12801 : IF (iatom <= jatom) THEN
3129 7824 : irow = iatom
3130 7824 : icol = jatom
3131 : ELSE
3132 4977 : irow = jatom
3133 4977 : icol = iatom
3134 : END IF
3135 : ELSE
3136 1402 : irow = iatom
3137 1402 : icol = jatom
3138 : END IF
3139 :
3140 56812 : dab = NORM2(rij)
3141 14203 : trans = do_symmetric .AND. (iatom > jatom)
3142 :
3143 56812 : DO i_xyz = 1, 3
3144 : CALL dbcsr_get_block_p(matrix=t2c_der(img, i_xyz), &
3145 42609 : row=irow, col=icol, BLOCK=block_t(i_xyz)%block, found=found)
3146 56812 : CPASSERT(found)
3147 : END DO
3148 :
3149 48405 : DO iset = 1, nseti
3150 :
3151 33830 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3152 33830 : sgfi = first_sgf_i(1, iset)
3153 :
3154 165921 : DO jset = 1, nsetj
3155 :
3156 117860 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3157 117860 : sgfj = first_sgf_j(1, jset)
3158 :
3159 151690 : IF (ncoi*ncoj > 0) THEN
3160 471440 : ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
3161 589300 : ALLOCATE (dij(ncoi, ncoj, 3))
3162 312725222 : dij(:, :, :) = 0.0_dp
3163 :
3164 117860 : ri = 0.0_dp
3165 117860 : 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 117860 : rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3170 :
3171 471440 : DO i_xyz = 1, 3
3172 :
3173 123133524 : dij_contr(:, :) = 0.0_dp
3174 : CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
3175 : sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3176 353580 : ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3177 :
3178 353580 : IF (trans) THEN
3179 : !if transpose, then -1 factor for antisymmetry
3180 448116 : ALLOCATE (dij_rs(nsgfj(jset), nsgfi(iset)))
3181 44937846 : dij_rs(:, :) = -1.0_dp*TRANSPOSE(dij_contr)
3182 : ELSE
3183 966204 : ALLOCATE (dij_rs(nsgfi(iset), nsgfj(jset)))
3184 78086691 : 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 353580 : sgfi, sgfj, trans=trans)
3190 471440 : DEALLOCATE (dij_rs)
3191 : END DO
3192 :
3193 117860 : DEALLOCATE (dij, dij_contr)
3194 : END IF
3195 : END DO
3196 : END DO
3197 : END DO
3198 :
3199 372 : CALL cp_libint_cleanup_2eri1(lib)
3200 :
3201 372 : CALL neighbor_list_iterator_release(nl_iterator)
3202 2536 : DO img = 1, nimg
3203 9028 : DO i_xyz = 1, 3
3204 6492 : CALL dbcsr_finalize(t2c_der(img, i_xyz))
3205 8656 : CALL dbcsr_filter(t2c_der(img, i_xyz), filter_eps)
3206 : END DO
3207 : END DO
3208 :
3209 372 : CALL timestop(handle)
3210 :
3211 1116 : 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 12 : 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 12 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3240 12 : npgfj, nsgfi, nsgfj
3241 12 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3242 : LOGICAL :: do_symmetric, found
3243 : REAL(dp) :: force
3244 12 : REAL(dp), DIMENSION(:, :), POINTER :: pblock
3245 : REAL(KIND=dp) :: dab
3246 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij_contr
3247 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dij
3248 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj, scoord
3249 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3250 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3251 12 : zetj
3252 12 : 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 12 : DIMENSION(:), POINTER :: nl_iterator
3259 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3260 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3261 :
3262 12 : CALL timeset(routineN, handle)
3263 :
3264 12 : op_prv = potential_parameter%potential_type
3265 :
3266 12 : 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 12 : particle_set=particle_set, cell=cell)
3272 :
3273 : ! check for symmetry
3274 12 : CPASSERT(SIZE(nl_2c) > 0)
3275 12 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3276 12 : CPASSERT(.NOT. do_symmetric)
3277 :
3278 12 : maxli = 0
3279 36 : DO ibasis = 1, SIZE(basis_i)
3280 24 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3281 36 : maxli = MAX(maxli, imax)
3282 : END DO
3283 12 : maxlj = 0
3284 36 : DO ibasis = 1, SIZE(basis_j)
3285 24 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3286 36 : maxlj = MAX(maxlj, imax)
3287 : END DO
3288 :
3289 12 : m_max = maxli + maxlj + 1
3290 :
3291 : !Init the truncated Coulomb operator
3292 12 : IF (op_prv == do_potential_truncated) THEN
3293 :
3294 2 : IF (m_max > get_lmax_init()) THEN
3295 0 : IF (para_env%mepos == 0) THEN
3296 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3297 : END IF
3298 0 : CALL init(m_max, unit_id, para_env%mepos, para_env)
3299 0 : IF (para_env%mepos == 0) THEN
3300 0 : CALL close_file(unit_id)
3301 : END IF
3302 : END IF
3303 : END IF
3304 :
3305 12 : CALL init_md_ftable(nmax=m_max)
3306 :
3307 12 : CALL cp_libint_init_2eri1(lib, MAX(maxli, maxlj))
3308 12 : CALL cp_libint_set_contrdepth(lib, 1)
3309 :
3310 12 : CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3311 1644 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3312 :
3313 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3314 1632 : 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 1632 : 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 1632 : sphi=sphi_j, zet=zetj)
3323 :
3324 6528 : dab = NORM2(rij)
3325 :
3326 1632 : CALL dbcsr_get_block_p(t2c_trace, iatom, jatom, pblock, found)
3327 1632 : IF (.NOT. found) CYCLE
3328 :
3329 6113 : DO iset = 1, nseti
3330 :
3331 4469 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3332 4469 : sgfi = first_sgf_i(1, iset)
3333 :
3334 25808 : DO jset = 1, nsetj
3335 :
3336 19707 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3337 19707 : sgfj = first_sgf_j(1, jset)
3338 :
3339 24176 : IF (ncoi*ncoj > 0) THEN
3340 78828 : ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
3341 98535 : ALLOCATE (dij(ncoi, ncoj, 3))
3342 5649186 : dij(:, :, :) = 0.0_dp
3343 :
3344 19707 : ri = 0.0_dp
3345 19707 : 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 19707 : rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3350 :
3351 78828 : DO i_xyz = 1, 3
3352 :
3353 2522610 : dij_contr(:, :) = 0.0_dp
3354 : CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
3355 : sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3356 59121 : ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3357 :
3358 2522610 : force = SUM(pblock(sgfi:sgfi + nsgfi(iset) - 1, sgfj:sgfj + nsgfj(jset) - 1)*dij_contr(:, :))
3359 59121 : force = pref*force
3360 :
3361 : !iatom virial
3362 59121 : CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
3363 236484 : DO j_xyz = 1, 3
3364 236484 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
3365 : END DO
3366 :
3367 : !jatom virial
3368 236484 : CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
3369 256191 : DO j_xyz = 1, 3
3370 236484 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - force*scoord(j_xyz)
3371 : END DO
3372 : END DO
3373 :
3374 19707 : DEALLOCATE (dij, dij_contr)
3375 : END IF
3376 : END DO
3377 : END DO
3378 : END DO
3379 12 : CALL neighbor_list_iterator_release(nl_iterator)
3380 :
3381 12 : CALL cp_libint_cleanup_2eri1(lib)
3382 :
3383 12 : CALL timestop(handle)
3384 :
3385 24 : 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 900 : SUBROUTINE build_2c_integrals(t2c, filter_eps, qs_env, &
3403 900 : 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 900 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3426 900 : npgfj, nsgfi, nsgfj
3427 900 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3428 900 : 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 900 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sij, sij_contr, sij_rs
3433 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj
3434 900 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3435 900 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3436 900 : zetj
3437 900 : 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 900 : DIMENSION(:), POINTER :: nl_iterator
3446 900 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3447 :
3448 900 : CALL timeset(routineN, handle)
3449 :
3450 900 : IF (PRESENT(do_kpoints)) THEN
3451 850 : do_kpoints_prv = do_kpoints
3452 : ELSE
3453 : do_kpoints_prv = .FALSE.
3454 : END IF
3455 :
3456 900 : IF (PRESENT(do_hfx_kpoints)) THEN
3457 266 : do_hfx_kpoints_prv = do_hfx_kpoints
3458 : ELSE
3459 : do_hfx_kpoints_prv = .FALSE.
3460 : END IF
3461 266 : IF (do_hfx_kpoints_prv) THEN
3462 140 : CPASSERT(do_kpoints_prv)
3463 : END IF
3464 :
3465 900 : my_regularization_RI = 0.0_dp
3466 900 : IF (PRESENT(regularization_RI)) my_regularization_RI = regularization_RI
3467 :
3468 900 : op_prv = potential_parameter%potential_type
3469 :
3470 900 : 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 900 : natom=natom, dft_control=dft_control, para_env=para_env, kpoints=kpoints)
3475 :
3476 900 : IF (PRESENT(ext_kpoints)) kpoints => ext_kpoints
3477 :
3478 900 : IF (do_kpoints_prv) THEN
3479 526 : nimg = SIZE(t2c)
3480 526 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
3481 2104 : kp_index_lbounds = LBOUND(cell_to_index)
3482 2104 : kp_index_ubounds = UBOUND(cell_to_index)
3483 : ELSE
3484 : nimg = 1
3485 : END IF
3486 :
3487 : ! check for symmetry
3488 900 : CPASSERT(SIZE(nl_2c) > 0)
3489 900 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3490 :
3491 900 : IF (do_symmetric) THEN
3492 11946 : DO img = 1, nimg
3493 11946 : CPASSERT(dbcsr_has_symmetry(t2c(img)))
3494 : END DO
3495 : ELSE
3496 3428 : DO img = 1, nimg
3497 3428 : CPASSERT(.NOT. dbcsr_has_symmetry(t2c(img)))
3498 : END DO
3499 : END IF
3500 :
3501 15374 : DO img = 1, nimg
3502 15374 : CALL cp_dbcsr_alloc_block_from_nbl(t2c(img), nl_2c)
3503 : END DO
3504 :
3505 900 : maxli = 0
3506 2602 : DO ibasis = 1, SIZE(basis_i)
3507 1702 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3508 2602 : maxli = MAX(maxli, imax)
3509 : END DO
3510 900 : maxlj = 0
3511 2602 : DO ibasis = 1, SIZE(basis_j)
3512 1702 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3513 2602 : maxlj = MAX(maxlj, imax)
3514 : END DO
3515 :
3516 900 : m_max = maxli + maxlj
3517 :
3518 : !Init the truncated Coulomb operator
3519 900 : IF (op_prv == do_potential_truncated) THEN
3520 :
3521 590 : IF (m_max > get_lmax_init()) THEN
3522 88 : IF (para_env%mepos == 0) THEN
3523 44 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3524 : END IF
3525 88 : CALL init(m_max, unit_id, para_env%mepos, para_env)
3526 88 : IF (para_env%mepos == 0) THEN
3527 44 : CALL close_file(unit_id)
3528 : END IF
3529 : END IF
3530 : END IF
3531 :
3532 900 : CALL init_md_ftable(nmax=m_max)
3533 :
3534 900 : CALL cp_libint_init_2eri(lib, MAX(maxli, maxlj))
3535 900 : CALL cp_libint_set_contrdepth(lib, 1)
3536 :
3537 900 : CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3538 28823 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3539 :
3540 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3541 27923 : iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
3542 27923 : IF (do_kpoints_prv) THEN
3543 47243 : 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 6749 : img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3546 6749 : 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 27873 : 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 27873 : sphi=sphi_j, zet=zetj)
3558 :
3559 27873 : IF (do_symmetric) THEN
3560 25253 : IF (iatom <= jatom) THEN
3561 15165 : irow = iatom
3562 15165 : icol = jatom
3563 : ELSE
3564 10088 : irow = jatom
3565 10088 : icol = iatom
3566 : END IF
3567 : ELSE
3568 2620 : irow = iatom
3569 2620 : icol = jatom
3570 : END IF
3571 :
3572 111492 : dab = NORM2(rij)
3573 :
3574 : CALL dbcsr_get_block_p(matrix=t2c(img), &
3575 27873 : row=irow, col=icol, BLOCK=block_t%block, found=found)
3576 27873 : CPASSERT(found)
3577 27873 : trans = do_symmetric .AND. (iatom > jatom)
3578 :
3579 105716 : DO iset = 1, nseti
3580 :
3581 76943 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3582 76943 : sgfi = first_sgf_i(1, iset)
3583 :
3584 492351 : DO jset = 1, nsetj
3585 :
3586 387485 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3587 387485 : sgfj = first_sgf_j(1, jset)
3588 :
3589 464428 : IF (ncoi*ncoj > 0) THEN
3590 1549940 : ALLOCATE (sij_contr(nsgfi(iset), nsgfj(jset)))
3591 76501075 : sij_contr(:, :) = 0.0_dp
3592 :
3593 1549940 : ALLOCATE (sij(ncoi, ncoj))
3594 199925913 : sij(:, :) = 0.0_dp
3595 :
3596 387485 : ri = 0.0_dp
3597 387485 : 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 387485 : rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3602 :
3603 : CALL ab_contract(sij_contr, sij, &
3604 : sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3605 387485 : ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3606 :
3607 387485 : DEALLOCATE (sij)
3608 387485 : IF (trans) THEN
3609 486560 : ALLOCATE (sij_rs(nsgfj(jset), nsgfi(iset)))
3610 29089897 : sij_rs(:, :) = TRANSPOSE(sij_contr)
3611 : ELSE
3612 1063380 : ALLOCATE (sij_rs(nsgfi(iset), nsgfj(jset)))
3613 47339704 : sij_rs(:, :) = sij_contr
3614 : END IF
3615 :
3616 387485 : DEALLOCATE (sij_contr)
3617 :
3618 : ! RI regularization
3619 387485 : 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 370316 : iatom == jatom .AND. iset == jset) THEN
3622 4058 : DO i_diag = 1, nsgfi(iset)
3623 : sij_rs(i_diag, i_diag) = sij_rs(i_diag, i_diag) + &
3624 9386 : 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 387485 : sgfi, sgfj, trans=trans)
3632 387485 : DEALLOCATE (sij_rs)
3633 :
3634 : END IF
3635 : END DO
3636 : END DO
3637 : END DO
3638 :
3639 900 : CALL cp_libint_cleanup_2eri(lib)
3640 :
3641 900 : CALL neighbor_list_iterator_release(nl_iterator)
3642 15374 : DO img = 1, nimg
3643 14474 : CALL dbcsr_finalize(t2c(img))
3644 15374 : CALL dbcsr_filter(t2c(img), filter_eps)
3645 : END DO
3646 :
3647 900 : CALL timestop(handle)
3648 :
3649 1800 : 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 20690 : 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 20690 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET :: blk_data
3676 20690 : REAL(dp), DIMENSION(:), POINTER :: blk_data_1d
3677 : TYPE(dbt_iterator_type) :: iter
3678 20690 : TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
3679 : TYPE(hfx_cache_type), POINTER :: maxval_cache
3680 20690 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
3681 : TYPE(hfx_container_type), POINTER :: maxval_container
3682 :
3683 20690 : CALL dealloc_containers(compressed, memory_usage)
3684 20690 : CALL alloc_containers(compressed, 1)
3685 :
3686 20690 : maxval_container => compressed%maxval_container(1)
3687 20690 : integral_containers => compressed%integral_containers(:, 1)
3688 :
3689 20690 : CALL hfx_init_container(maxval_container, memory_usage, .FALSE.)
3690 1344850 : DO i = 1, 64
3691 1344850 : CALL hfx_init_container(integral_containers(i), memory_usage, .FALSE.)
3692 : END DO
3693 :
3694 20690 : maxval_cache => compressed%maxval_cache(1)
3695 20690 : integral_caches => compressed%integral_caches(:, 1)
3696 :
3697 20690 : IF (ALLOCATED(blk_indices)) DEALLOCATE (blk_indices)
3698 55017 : ALLOCATE (blk_indices(dbt_get_num_blocks(tensor), 3))
3699 20690 : shared_offset = 0
3700 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices,shared_offset) &
3701 20690 : !$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 292942 : DO i = 1, SIZE(blk_indices, 1)
3717 1089008 : ind = blk_indices(i, :)
3718 272252 : CALL dbt_get_block(tensor, ind, blk_data, found)
3719 272252 : CPASSERT(found)
3720 1089008 : nints = SIZE(blk_data)
3721 272252 : blk_data_1d(1:nints) => blk_data
3722 141757668 : spherical_estimate = MAXVAL(ABS(blk_data_1d))
3723 272252 : IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate)
3724 272252 : estimate_to_store_int = EXPONENT(spherical_estimate)
3725 272252 : 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 272252 : .FALSE.)
3730 :
3731 272252 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
3732 :
3733 272252 : nbits = EXPONENT(ANINT(spherical_estimate/eps)) + 1
3734 272252 : IF (nbits > 64) THEN
3735 : CALL cp_abort(__LOCATION__, &
3736 0 : "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 601769 : DO WHILE (buffer_left > 0)
3743 329517 : 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 329517 : .FALSE.)
3751 329517 : buffer_left = buffer_left - buffer_size
3752 329517 : buffer_start = buffer_start + buffer_size
3753 : END DO
3754 :
3755 565194 : NULLIFY (blk_data_1d); DEALLOCATE (blk_data)
3756 : END DO
3757 :
3758 20690 : CALL dbt_clear(tensor)
3759 :
3760 20690 : storage_counter_integrals = memory_usage*cache_size
3761 20690 : 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 20690 : .FALSE.)
3767 1344850 : DO i = 1, 64
3768 : CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
3769 1344850 : memory_usage, .FALSE.)
3770 : END DO
3771 :
3772 20690 : CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .FALSE.)
3773 1344850 : DO i = 1, 64
3774 : CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
3775 1344850 : memory_usage, .FALSE.)
3776 : END DO
3777 :
3778 41380 : 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 60684 : 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 60684 : REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: blk_data
3801 60684 : REAL(dp), DIMENSION(:, :, :), POINTER :: blk_data_3d
3802 60684 : TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
3803 : TYPE(hfx_cache_type), POINTER :: maxval_cache
3804 60684 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
3805 : TYPE(hfx_container_type), POINTER :: maxval_container
3806 :
3807 60684 : maxval_cache => compressed%maxval_cache(1)
3808 60684 : maxval_container => compressed%maxval_container(1)
3809 60684 : integral_caches => compressed%integral_caches(:, 1)
3810 60684 : integral_containers => compressed%integral_containers(:, 1)
3811 :
3812 60684 : memory_usage = 0
3813 :
3814 60684 : CALL hfx_decompress_first_cache(6, maxval_cache, maxval_container, memory_usage, .FALSE.)
3815 :
3816 3944460 : DO i = 1, 64
3817 : CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), &
3818 3944460 : memory_usage, .FALSE.)
3819 : END DO
3820 :
3821 : !TODO: Parallelize creation of block list.
3822 60684 : !$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 1200555 : DO i = 1, SIZE(blk_indices, 1)
3831 4559484 : ind = blk_indices(i, :)
3832 1139871 : CALL dbt_blk_sizes(tensor, ind, blk_size)
3833 4559484 : 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 1139871 : .FALSE.)
3838 :
3839 1139871 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
3840 :
3841 1139871 : nbits = EXPONENT(ANINT(spherical_estimate/eps)) + 1
3842 :
3843 1139871 : buffer_left = nints
3844 1139871 : buffer_start = 1
3845 :
3846 3419613 : ALLOCATE (blk_data(nints))
3847 2394575 : DO WHILE (buffer_left > 0)
3848 1254704 : 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 1254704 : .FALSE.)
3856 1254704 : buffer_left = buffer_left - buffer_size
3857 1254704 : buffer_start = buffer_start + buffer_size
3858 : END DO
3859 :
3860 1139871 : blk_data_3d(1:blk_size(1), 1:blk_size(2), 1:blk_size(3)) => blk_data
3861 1139871 : CALL dbt_put_block(tensor, ind, blk_size, blk_data_3d)
3862 1200555 : NULLIFY (blk_data_3d); DEALLOCATE (blk_data)
3863 : END DO
3864 :
3865 60684 : CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .FALSE.)
3866 3944460 : DO i = 1, 64
3867 : CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
3868 3944460 : memory_usage, .FALSE.)
3869 : END DO
3870 121368 : END SUBROUTINE
3871 :
3872 : ! **************************************************************************************************
3873 : !> \brief ...
3874 : !> \param tensor ...
3875 : !> \param nze ...
3876 : !> \param occ ...
3877 : ! **************************************************************************************************
3878 111923 : 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 223846 : INTEGER, DIMENSION(dbt_ndims(tensor)) :: dims
3884 :
3885 111923 : nze = dbt_get_nze_total(tensor)
3886 111923 : CALL dbt_get_info(tensor, nfull_total=dims)
3887 430606 : occ = REAL(nze, dp)/PRODUCT(REAL(dims, dp))
3888 :
3889 111923 : END SUBROUTINE
3890 :
3891 0 : END MODULE
|
