Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_api, ONLY: dbcsr_filter,&
27 : dbcsr_finalize,&
28 : dbcsr_get_block_p,&
29 : dbcsr_get_matrix_type,&
30 : dbcsr_has_symmetry,&
31 : dbcsr_type,&
32 : dbcsr_type_antisymmetric,&
33 : dbcsr_type_no_symmetry
34 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
35 : USE cp_files, ONLY: close_file,&
36 : open_file
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_iterator_next_block, &
40 : dbt_iterator_num_blocks, dbt_iterator_start, dbt_iterator_stop, dbt_iterator_type, &
41 : 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_mix_cl_trunc,&
61 : do_potential_short,&
62 : do_potential_truncated
63 : USE input_section_types, ONLY: section_vals_val_get
64 : USE kinds, ONLY: dp,&
65 : int_8
66 : USE kpoint_types, ONLY: get_kpoint_info,&
67 : kpoint_type
68 : USE libint_2c_3c, ONLY: cutoff_screen_factor,&
69 : eri_2center,&
70 : eri_2center_derivs,&
71 : eri_3center,&
72 : eri_3center_derivs,&
73 : libint_potential_type
74 : USE libint_wrapper, ONLY: &
75 : cp_libint_cleanup_2eri, cp_libint_cleanup_2eri1, cp_libint_cleanup_3eri, &
76 : cp_libint_cleanup_3eri1, cp_libint_init_2eri, cp_libint_init_2eri1, cp_libint_init_3eri, &
77 : cp_libint_init_3eri1, cp_libint_set_contrdepth, cp_libint_t
78 : USE message_passing, ONLY: mp_para_env_type
79 : USE molecule_types, ONLY: molecule_type
80 : USE orbital_pointers, ONLY: ncoset
81 : USE particle_types, ONLY: particle_type
82 : USE qs_environment_types, ONLY: get_qs_env,&
83 : qs_environment_type
84 : USE qs_kind_types, ONLY: qs_kind_type
85 : USE qs_neighbor_list_types, ONLY: &
86 : get_iterator_info, get_neighbor_list_set_p, neighbor_list_iterate, &
87 : neighbor_list_iterator_create, neighbor_list_iterator_p_type, &
88 : neighbor_list_iterator_release, neighbor_list_set_p_type, nl_sub_iterate, &
89 : release_neighbor_list_sets
90 : USE qs_neighbor_lists, ONLY: atom2d_build,&
91 : atom2d_cleanup,&
92 : build_neighbor_lists,&
93 : local_atoms_type,&
94 : pair_radius_setup
95 : USE qs_tensors_types, ONLY: &
96 : distribution_3d_destroy, distribution_3d_type, neighbor_list_3c_iterator_type, &
97 : neighbor_list_3c_type, symmetric_ij, symmetric_ijk, symmetric_jk, symmetric_none, &
98 : symmetrik_ik
99 : USE t_c_g0, ONLY: get_lmax_init,&
100 : init
101 : USE util, ONLY: get_limit
102 :
103 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
104 : #include "./base/base_uses.f90"
105 :
106 : IMPLICIT NONE
107 :
108 : PRIVATE
109 :
110 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tensors'
111 :
112 : PUBLIC :: build_3c_neighbor_lists, &
113 : neighbor_list_3c_destroy, neighbor_list_3c_iterate, neighbor_list_3c_iterator_create, &
114 : neighbor_list_3c_iterator_destroy, get_3c_iterator_info, build_3c_integrals, &
115 : build_2c_neighbor_lists, build_2c_integrals, cutoff_screen_factor, &
116 : get_tensor_occupancy, compress_tensor, decompress_tensor, &
117 : build_3c_derivatives, build_2c_derivatives, calc_2c_virial, calc_3c_virial
118 :
119 : TYPE one_dim_int_array
120 : INTEGER, DIMENSION(:), ALLOCATABLE :: array
121 : END TYPE one_dim_int_array
122 :
123 : ! cache size for integral compression
124 : INTEGER, PARAMETER, PRIVATE :: cache_size = 1024
125 :
126 : CONTAINS
127 :
128 : ! **************************************************************************************************
129 : !> \brief Build 2-center neighborlists adapted to different operators
130 : !> This mainly wraps build_neighbor_lists for consistency with build_3c_neighbor_lists
131 : !> \param ij_list 2c neighbor list for atom pairs i, j
132 : !> \param basis_i basis object for atoms i
133 : !> \param basis_j basis object for atoms j
134 : !> \param potential_parameter ...
135 : !> \param name name of 2c neighbor list
136 : !> \param qs_env ...
137 : !> \param sym_ij Symmetry in i, j (default .TRUE.)
138 : !> \param molecular ...
139 : !> \param dist_2d optionally a custom 2d distribution
140 : !> \param pot_to_rad which radius (1 for i, 2 for j) should be adapted to incorporate potential
141 : ! **************************************************************************************************
142 11908 : SUBROUTINE build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, &
143 : sym_ij, molecular, dist_2d, pot_to_rad)
144 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
145 : POINTER :: ij_list
146 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
147 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
148 : CHARACTER(LEN=*), INTENT(IN) :: name
149 : TYPE(qs_environment_type), POINTER :: qs_env
150 : LOGICAL, INTENT(IN), OPTIONAL :: sym_ij, molecular
151 : TYPE(distribution_2d_type), OPTIONAL, POINTER :: dist_2d
152 : INTEGER, INTENT(IN), OPTIONAL :: pot_to_rad
153 :
154 : INTEGER :: ikind, nkind, pot_to_rad_prv
155 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: i_present, j_present
156 11908 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
157 : REAL(kind=dp) :: subcells
158 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: i_radius, j_radius
159 11908 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
160 : TYPE(cell_type), POINTER :: cell
161 : TYPE(distribution_1d_type), POINTER :: local_particles
162 : TYPE(distribution_2d_type), POINTER :: dist_2d_prv
163 11908 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
164 11908 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
165 11908 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
166 :
167 11908 : NULLIFY (atomic_kind_set, cell, local_particles, molecule_set, &
168 11908 : particle_set, dist_2d_prv)
169 :
170 11908 : IF (PRESENT(pot_to_rad)) THEN
171 1392 : pot_to_rad_prv = pot_to_rad
172 : ELSE
173 : pot_to_rad_prv = 1
174 : END IF
175 :
176 : CALL get_qs_env(qs_env, &
177 : nkind=nkind, &
178 : cell=cell, &
179 : particle_set=particle_set, &
180 : atomic_kind_set=atomic_kind_set, &
181 : local_particles=local_particles, &
182 : distribution_2d=dist_2d_prv, &
183 11908 : molecule_set=molecule_set)
184 :
185 11908 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
186 :
187 53577 : ALLOCATE (i_present(nkind), SOURCE=.FALSE.)
188 41669 : ALLOCATE (j_present(nkind), SOURCE=.FALSE.)
189 53577 : ALLOCATE (i_radius(nkind), SOURCE=0.0_dp)
190 41669 : ALLOCATE (j_radius(nkind), SOURCE=0.0_dp)
191 :
192 11908 : IF (PRESENT(dist_2d)) dist_2d_prv => dist_2d
193 :
194 : ! Set up the radii, depending on the operator type
195 11908 : IF (potential_parameter%potential_type == do_potential_id) THEN
196 :
197 : !overlap => use the kind radius for both i and j
198 25295 : DO ikind = 1, nkind
199 15133 : IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
200 15133 : i_present(ikind) = .TRUE.
201 15133 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
202 : END IF
203 25295 : IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
204 15133 : j_present(ikind) = .TRUE.
205 15133 : CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
206 : END IF
207 : END DO
208 :
209 : ELSE IF (potential_parameter%potential_type == do_potential_coulomb) THEN
210 :
211 : !Coulomb operator, virtually infinite range => set j_radius to arbitrarily large number
212 1110 : DO ikind = 1, nkind
213 736 : IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
214 736 : i_present(ikind) = .TRUE.
215 736 : IF (pot_to_rad_prv == 1) i_radius(ikind) = 1000000.0_dp
216 : END IF
217 1110 : IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
218 736 : j_present(ikind) = .TRUE.
219 736 : IF (pot_to_rad_prv == 2) j_radius(ikind) = 1000000.0_dp
220 : END IF
221 : END DO !ikind
222 :
223 : ELSE IF (potential_parameter%potential_type == do_potential_truncated .OR. &
224 : potential_parameter%potential_type == do_potential_short .OR. &
225 : potential_parameter%potential_type == do_potential_mix_cl_trunc) THEN
226 :
227 : !Truncated coulomb/short range: set j_radius to r_cutoff + the kind_radii
228 3356 : DO ikind = 1, nkind
229 1984 : IF (ASSOCIATED(basis_i(ikind)%gto_basis_set)) THEN
230 1984 : i_present(ikind) = .TRUE.
231 1984 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=i_radius(ikind))
232 1984 : IF (pot_to_rad_prv == 1) i_radius(ikind) = i_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
233 : END IF
234 3356 : IF (ASSOCIATED(basis_j(ikind)%gto_basis_set)) THEN
235 1984 : j_present(ikind) = .TRUE.
236 1984 : CALL get_gto_basis_set(basis_j(ikind)%gto_basis_set, kind_radius=j_radius(ikind))
237 1984 : IF (pot_to_rad_prv == 2) j_radius(ikind) = j_radius(ikind) + cutoff_screen_factor*potential_parameter%cutoff_radius
238 : END IF
239 : END DO
240 :
241 : ELSE
242 0 : CPABORT("Operator not implemented.")
243 : END IF
244 :
245 95264 : ALLOCATE (pair_radius(nkind, nkind), SOURCE=0.0_dp)
246 11908 : CALL pair_radius_setup(i_present, j_present, i_radius, j_radius, pair_radius)
247 :
248 53577 : ALLOCATE (atom2d(nkind))
249 :
250 : CALL atom2d_build(atom2d, local_particles, dist_2d_prv, atomic_kind_set, &
251 11908 : molecule_set, molecule_only=.FALSE., particle_set=particle_set)
252 : CALL build_neighbor_lists(ij_list, particle_set, atom2d, cell, pair_radius, subcells, &
253 11908 : symmetric=sym_ij, molecular=molecular, nlname=TRIM(name))
254 :
255 11908 : CALL atom2d_cleanup(atom2d)
256 :
257 35724 : END SUBROUTINE build_2c_neighbor_lists
258 :
259 : ! **************************************************************************************************
260 : !> \brief Build a 3-center neighbor list
261 : !> \param ijk_list 3c neighbor list for atom triples i, j, k
262 : !> \param basis_i basis object for atoms i
263 : !> \param basis_j basis object for atoms j
264 : !> \param basis_k basis object for atoms k
265 : !> \param dist_3d 3d distribution object
266 : !> \param potential_parameter ...
267 : !> \param name name of 3c neighbor list
268 : !> \param qs_env ...
269 : !> \param sym_ij Symmetry in i, j (default .FALSE.)
270 : !> \param sym_jk Symmetry in j, k (default .FALSE.)
271 : !> \param sym_ik Symmetry in i, k (default .FALSE.)
272 : !> \param molecular ??? not tested
273 : !> \param op_pos ...
274 : !> \param own_dist ...
275 : ! **************************************************************************************************
276 696 : SUBROUTINE build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, &
277 : dist_3d, potential_parameter, name, qs_env, &
278 : sym_ij, sym_jk, sym_ik, molecular, op_pos, &
279 : own_dist)
280 : TYPE(neighbor_list_3c_type), INTENT(OUT) :: ijk_list
281 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
282 : TYPE(distribution_3d_type), INTENT(IN) :: dist_3d
283 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
284 : CHARACTER(LEN=*), INTENT(IN) :: name
285 : TYPE(qs_environment_type), POINTER :: qs_env
286 : LOGICAL, INTENT(IN), OPTIONAL :: sym_ij, sym_jk, sym_ik, molecular
287 : INTEGER, INTENT(IN), OPTIONAL :: op_pos
288 : LOGICAL, INTENT(IN), OPTIONAL :: own_dist
289 :
290 : CHARACTER(len=*), PARAMETER :: routineN = 'build_3c_neighbor_lists'
291 :
292 : INTEGER :: handle, op_pos_prv, sym_level
293 : TYPE(libint_potential_type) :: pot_par_1, pot_par_2
294 :
295 696 : CALL timeset(routineN, handle)
296 :
297 696 : IF (PRESENT(op_pos)) THEN
298 512 : op_pos_prv = op_pos
299 : ELSE
300 : op_pos_prv = 1
301 : END IF
302 :
303 490 : SELECT CASE (op_pos_prv)
304 : CASE (1)
305 490 : pot_par_1 = potential_parameter
306 490 : pot_par_2%potential_type = do_potential_id
307 : CASE (2)
308 206 : pot_par_2 = potential_parameter
309 512 : pot_par_1%potential_type = do_potential_id
310 : END SELECT
311 :
312 : CALL build_2c_neighbor_lists(ijk_list%ij_list, basis_i, basis_j, pot_par_1, TRIM(name)//"_sub_1", &
313 : qs_env, sym_ij=.FALSE., molecular=molecular, &
314 696 : dist_2d=dist_3d%dist_2d_1, pot_to_rad=1)
315 :
316 : CALL build_2c_neighbor_lists(ijk_list%jk_list, basis_j, basis_k, pot_par_2, TRIM(name)//"_sub_2", &
317 : qs_env, sym_ij=.FALSE., molecular=molecular, &
318 696 : dist_2d=dist_3d%dist_2d_2, pot_to_rad=2)
319 :
320 696 : ijk_list%sym = symmetric_none
321 :
322 696 : sym_level = 0
323 696 : IF (PRESENT(sym_ij)) THEN
324 160 : IF (sym_ij) THEN
325 0 : ijk_list%sym = symmetric_ij
326 0 : sym_level = sym_level + 1
327 : END IF
328 : END IF
329 :
330 696 : IF (PRESENT(sym_jk)) THEN
331 486 : IF (sym_jk) THEN
332 428 : ijk_list%sym = symmetric_jk
333 428 : sym_level = sym_level + 1
334 : END IF
335 : END IF
336 :
337 696 : IF (PRESENT(sym_ik)) THEN
338 0 : IF (sym_ik) THEN
339 0 : ijk_list%sym = symmetrik_ik
340 0 : sym_level = sym_level + 1
341 : END IF
342 : END IF
343 :
344 696 : IF (sym_level >= 2) THEN
345 0 : ijk_list%sym = symmetric_ijk
346 : END IF
347 :
348 696 : ijk_list%dist_3d = dist_3d
349 696 : IF (PRESENT(own_dist)) THEN
350 696 : ijk_list%owns_dist = own_dist
351 : ELSE
352 0 : ijk_list%owns_dist = .FALSE.
353 : END IF
354 :
355 696 : CALL timestop(handle)
356 696 : END SUBROUTINE build_3c_neighbor_lists
357 :
358 : ! **************************************************************************************************
359 : !> \brief Symmetry criterion
360 : !> \param a ...
361 : !> \param b ...
362 : !> \return ...
363 : ! **************************************************************************************************
364 4207758 : PURE FUNCTION include_symmetric(a, b)
365 : INTEGER, INTENT(IN) :: a, b
366 : LOGICAL :: include_symmetric
367 :
368 4207758 : IF (a > b) THEN
369 1691358 : include_symmetric = (MODULO(a + b, 2) /= 0)
370 : ELSE
371 2516400 : include_symmetric = (MODULO(a + b, 2) == 0)
372 : END IF
373 :
374 4207758 : END FUNCTION include_symmetric
375 :
376 : ! **************************************************************************************************
377 : !> \brief Destroy 3c neighborlist
378 : !> \param ijk_list ...
379 : ! **************************************************************************************************
380 696 : SUBROUTINE neighbor_list_3c_destroy(ijk_list)
381 : TYPE(neighbor_list_3c_type), INTENT(INOUT) :: ijk_list
382 :
383 696 : CALL release_neighbor_list_sets(ijk_list%ij_list)
384 696 : CALL release_neighbor_list_sets(ijk_list%jk_list)
385 :
386 696 : IF (ijk_list%owns_dist) THEN
387 696 : CALL distribution_3d_destroy(ijk_list%dist_3d)
388 : END IF
389 :
390 696 : END SUBROUTINE neighbor_list_3c_destroy
391 :
392 : ! **************************************************************************************************
393 : !> \brief Create a 3-center neighborlist iterator
394 : !> \param iterator ...
395 : !> \param ijk_nl ...
396 : ! **************************************************************************************************
397 69768 : SUBROUTINE neighbor_list_3c_iterator_create(iterator, ijk_nl)
398 : TYPE(neighbor_list_3c_iterator_type), INTENT(OUT) :: iterator
399 : TYPE(neighbor_list_3c_type), INTENT(IN) :: ijk_nl
400 :
401 : CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterator_create'
402 :
403 : INTEGER :: handle
404 :
405 7752 : CALL timeset(routineN, handle)
406 :
407 7752 : CALL neighbor_list_iterator_create(iterator%iter_ij, ijk_nl%ij_list)
408 7752 : CALL neighbor_list_iterator_create(iterator%iter_jk, ijk_nl%jk_list, search=.TRUE.)
409 :
410 7752 : iterator%iter_level = 0
411 7752 : iterator%ijk_nl = ijk_nl
412 :
413 23256 : iterator%bounds_i = 0
414 23256 : iterator%bounds_j = 0
415 23256 : iterator%bounds_k = 0
416 :
417 7752 : CALL timestop(handle)
418 7752 : END SUBROUTINE neighbor_list_3c_iterator_create
419 :
420 : ! **************************************************************************************************
421 : !> \brief impose atomic upper and lower bounds
422 : !> \param iterator ...
423 : !> \param bounds_i ...
424 : !> \param bounds_j ...
425 : !> \param bounds_k ...
426 : ! **************************************************************************************************
427 7672 : SUBROUTINE nl_3c_iter_set_bounds(iterator, bounds_i, bounds_j, bounds_k)
428 : TYPE(neighbor_list_3c_iterator_type), &
429 : INTENT(INOUT) :: iterator
430 : INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
431 :
432 24190 : IF (PRESENT(bounds_i)) iterator%bounds_i = bounds_i
433 14260 : IF (PRESENT(bounds_j)) iterator%bounds_j = bounds_j
434 18292 : IF (PRESENT(bounds_k)) iterator%bounds_k = bounds_k
435 :
436 7672 : END SUBROUTINE nl_3c_iter_set_bounds
437 :
438 : ! **************************************************************************************************
439 : !> \brief Destroy 3c-nl iterator
440 : !> \param iterator ...
441 : ! **************************************************************************************************
442 7752 : SUBROUTINE neighbor_list_3c_iterator_destroy(iterator)
443 : TYPE(neighbor_list_3c_iterator_type), &
444 : INTENT(INOUT) :: iterator
445 :
446 : CHARACTER(len=*), PARAMETER :: routineN = 'neighbor_list_3c_iterator_destroy'
447 :
448 : INTEGER :: handle
449 :
450 7752 : CALL timeset(routineN, handle)
451 7752 : CALL neighbor_list_iterator_release(iterator%iter_ij)
452 7752 : CALL neighbor_list_iterator_release(iterator%iter_jk)
453 7752 : NULLIFY (iterator%iter_ij)
454 7752 : NULLIFY (iterator%iter_jk)
455 :
456 7752 : CALL timestop(handle)
457 7752 : END SUBROUTINE neighbor_list_3c_iterator_destroy
458 :
459 : ! **************************************************************************************************
460 : !> \brief Iterate 3c-nl iterator
461 : !> \param iterator ...
462 : !> \return 0 if successful; 1 if end was reached
463 : ! **************************************************************************************************
464 5706732 : RECURSIVE FUNCTION neighbor_list_3c_iterate(iterator) RESULT(iter_stat)
465 : TYPE(neighbor_list_3c_iterator_type), &
466 : INTENT(INOUT) :: iterator
467 : INTEGER :: iter_stat
468 :
469 : INTEGER :: iatom, iter_level, jatom, jatom_1, &
470 : jatom_2, katom
471 : LOGICAL :: skip_this
472 :
473 5706732 : iter_level = iterator%iter_level
474 :
475 5706732 : IF (iter_level == 0) THEN
476 195772 : iter_stat = neighbor_list_iterate(iterator%iter_ij)
477 :
478 195772 : IF (iter_stat /= 0) THEN
479 : RETURN
480 : END IF
481 :
482 188020 : CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom)
483 188020 : skip_this = .FALSE.
484 : IF ((iterator%bounds_i(1) > 0 .AND. iatom < iterator%bounds_i(1)) &
485 188020 : .OR. (iterator%bounds_i(2) > 0 .AND. iatom > iterator%bounds_i(2))) skip_this = .TRUE.
486 : IF ((iterator%bounds_j(1) > 0 .AND. jatom < iterator%bounds_j(1)) &
487 188020 : .OR. (iterator%bounds_j(2) > 0 .AND. jatom > iterator%bounds_j(2))) skip_this = .TRUE.
488 :
489 188020 : IF (skip_this) THEN
490 85323 : iter_stat = neighbor_list_3c_iterate(iterator)
491 85323 : RETURN
492 : END IF
493 : END IF
494 :
495 5613657 : iter_stat = nl_sub_iterate(iterator%iter_jk, iterator%iter_ij)
496 5613657 : IF (iter_stat /= 0) THEN
497 102697 : iterator%iter_level = 0
498 102697 : iter_stat = neighbor_list_3c_iterate(iterator)
499 102697 : RETURN
500 : ELSE
501 5510960 : iterator%iter_level = 1
502 : END IF
503 :
504 : CPASSERT(iter_stat == 0)
505 : CPASSERT(iterator%iter_level == 1)
506 5510960 : CALL get_iterator_info(iterator%iter_ij, iatom=iatom, jatom=jatom_1)
507 5510960 : CALL get_iterator_info(iterator%iter_jk, iatom=jatom_2, jatom=katom)
508 :
509 5510960 : CPASSERT(jatom_1 == jatom_2)
510 5510960 : jatom = jatom_1
511 :
512 5510960 : skip_this = .FALSE.
513 : IF ((iterator%bounds_k(1) > 0 .AND. katom < iterator%bounds_k(1)) &
514 5510960 : .OR. (iterator%bounds_k(2) > 0 .AND. katom > iterator%bounds_k(2))) skip_this = .TRUE.
515 :
516 : IF (skip_this) THEN
517 312756 : iter_stat = neighbor_list_3c_iterate(iterator)
518 312756 : RETURN
519 : END IF
520 :
521 5198204 : SELECT CASE (iterator%ijk_nl%sym)
522 : CASE (symmetric_none)
523 0 : skip_this = .FALSE.
524 : CASE (symmetric_ij)
525 0 : skip_this = .NOT. include_symmetric(iatom, jatom)
526 : CASE (symmetric_jk)
527 4207758 : skip_this = .NOT. include_symmetric(jatom, katom)
528 : CASE (symmetrik_ik)
529 0 : skip_this = .NOT. include_symmetric(iatom, katom)
530 : CASE (symmetric_ijk)
531 0 : skip_this = .NOT. include_symmetric(iatom, jatom) .OR. .NOT. include_symmetric(jatom, katom)
532 : CASE DEFAULT
533 5198204 : CPABORT("should not happen")
534 : END SELECT
535 :
536 4207758 : IF (skip_this) THEN
537 1692660 : iter_stat = neighbor_list_3c_iterate(iterator)
538 1692660 : RETURN
539 : END IF
540 :
541 : END FUNCTION neighbor_list_3c_iterate
542 :
543 : ! **************************************************************************************************
544 : !> \brief Get info of current iteration
545 : !> \param iterator ...
546 : !> \param ikind ...
547 : !> \param jkind ...
548 : !> \param kkind ...
549 : !> \param nkind ...
550 : !> \param iatom ...
551 : !> \param jatom ...
552 : !> \param katom ...
553 : !> \param rij ...
554 : !> \param rjk ...
555 : !> \param rik ...
556 : !> \param cell_j ...
557 : !> \param cell_k ...
558 : !> \return ...
559 : ! **************************************************************************************************
560 3505544 : SUBROUTINE get_3c_iterator_info(iterator, ikind, jkind, kkind, nkind, iatom, jatom, katom, &
561 : rij, rjk, rik, cell_j, cell_k)
562 : TYPE(neighbor_list_3c_iterator_type), &
563 : INTENT(INOUT) :: iterator
564 : INTEGER, INTENT(OUT), OPTIONAL :: ikind, jkind, kkind, nkind, iatom, &
565 : jatom, katom
566 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: rij, rjk, rik
567 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: cell_j, cell_k
568 :
569 : INTEGER, DIMENSION(2) :: atoms_1, atoms_2, kinds_1, kinds_2
570 : INTEGER, DIMENSION(3) :: cell_1, cell_2
571 : REAL(KIND=dp), DIMENSION(3) :: r_1, r_2
572 :
573 3505544 : CPASSERT(iterator%iter_level == 1)
574 :
575 : CALL get_iterator_info(iterator%iter_ij, &
576 : ikind=kinds_1(1), jkind=kinds_1(2), nkind=nkind, &
577 : iatom=atoms_1(1), jatom=atoms_1(2), r=r_1, &
578 3505544 : cell=cell_1)
579 :
580 : CALL get_iterator_info(iterator%iter_jk, &
581 : ikind=kinds_2(1), jkind=kinds_2(2), &
582 : iatom=atoms_2(1), jatom=atoms_2(2), r=r_2, &
583 3505544 : cell=cell_2)
584 :
585 3505544 : IF (PRESENT(ikind)) ikind = kinds_1(1)
586 3505544 : IF (PRESENT(jkind)) jkind = kinds_1(2)
587 3505544 : IF (PRESENT(kkind)) kkind = kinds_2(2)
588 3505544 : IF (PRESENT(iatom)) iatom = atoms_1(1)
589 3505544 : IF (PRESENT(jatom)) jatom = atoms_1(2)
590 3505544 : IF (PRESENT(katom)) katom = atoms_2(2)
591 :
592 3505544 : IF (PRESENT(rij)) rij = r_1
593 3505544 : IF (PRESENT(rjk)) rjk = r_2
594 17527720 : IF (PRESENT(rik)) rik = r_1 + r_2
595 :
596 3505544 : IF (PRESENT(cell_j)) cell_j = cell_1
597 16848520 : IF (PRESENT(cell_k)) cell_k = cell_1 + cell_2
598 :
599 3505544 : END SUBROUTINE get_3c_iterator_info
600 :
601 : ! **************************************************************************************************
602 : !> \brief Allocate blocks of a 3-center tensor based on neighborlist
603 : !> \param t3c empty DBCSR tensor
604 : !> Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages)
605 : !> if k-points are used
606 : !> \param nl_3c 3-center neighborlist
607 : !> \param basis_i ...
608 : !> \param basis_j ...
609 : !> \param basis_k ...
610 : !> \param qs_env ...
611 : !> \param potential_parameter ...
612 : !> \param op_pos ...
613 : !> \param do_kpoints ...
614 : !> \param do_hfx_kpoints ...
615 : !> \param bounds_i ...
616 : !> \param bounds_j ...
617 : !> \param bounds_k ...
618 : !> \param RI_range ...
619 : !> \param img_to_RI_cell ...
620 : !> \param cell_to_index ...
621 : !> \param cell_sym ...
622 : ! **************************************************************************************************
623 2552 : 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 2552 : img_to_RI_cell, cell_to_index, cell_sym)
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 : INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: cell_to_index
637 : LOGICAL, INTENT(IN), OPTIONAL :: cell_sym
638 :
639 : CHARACTER(LEN=*), PARAMETER :: routineN = 'alloc_block_3c'
640 :
641 : INTEGER :: handle, iatom, ikind, j_img, jatom, &
642 : jcell, jkind, k_img, katom, kcell, &
643 : kkind, natom, ncell_RI, nimg, op_ij, &
644 : op_jk, op_pos_prv
645 : INTEGER(int_8) :: a, b, nblk_per_thread
646 2552 : INTEGER(int_8), ALLOCATABLE, DIMENSION(:, :) :: nblk
647 2552 : INTEGER, ALLOCATABLE, DIMENSION(:) :: img_to_RI_cell_prv
648 : INTEGER, DIMENSION(3) :: blk_idx, cell_j, cell_k, &
649 : kp_index_lbounds, kp_index_ubounds
650 : LOGICAL :: cell_sym_prv, do_hfx_kpoints_prv, &
651 : do_kpoints_prv
652 : REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
653 : kind_radius_i, kind_radius_j, &
654 : kind_radius_k
655 : REAL(KIND=dp), DIMENSION(3) :: rij, rik, rjk
656 2552 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
657 : TYPE(cell_type), POINTER :: cell
658 : TYPE(dft_control_type), POINTER :: dft_control
659 : TYPE(mp_para_env_type), POINTER :: para_env
660 : TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
661 : TYPE(one_dim_int_array), ALLOCATABLE, &
662 2552 : DIMENSION(:, :) :: alloc_i, alloc_j, alloc_k
663 2552 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
664 :
665 2552 : CALL timeset(routineN, handle)
666 2552 : NULLIFY (qs_kind_set, atomic_kind_set, cell)
667 :
668 2552 : IF (PRESENT(do_kpoints)) THEN
669 464 : do_kpoints_prv = do_kpoints
670 : ELSE
671 : do_kpoints_prv = .FALSE.
672 : END IF
673 2552 : IF (PRESENT(do_hfx_kpoints)) THEN
674 222 : do_hfx_kpoints_prv = do_hfx_kpoints
675 : ELSE
676 : do_hfx_kpoints_prv = .FALSE.
677 : END IF
678 222 : IF (do_hfx_kpoints_prv) THEN
679 222 : CPASSERT(do_kpoints_prv)
680 222 : CPASSERT(PRESENT(RI_range))
681 222 : CPASSERT(PRESENT(img_to_RI_cell))
682 : END IF
683 :
684 2330 : IF (PRESENT(img_to_RI_cell)) THEN
685 666 : ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
686 10432 : img_to_RI_cell_prv(:) = img_to_RI_cell
687 : END IF
688 :
689 2552 : IF (PRESENT(cell_sym)) THEN
690 1978 : cell_sym_prv = cell_sym
691 : ELSE
692 : cell_sym_prv = .FALSE.
693 : END IF
694 :
695 2552 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
696 :
697 2552 : op_ij = do_potential_id; op_jk = do_potential_id
698 :
699 2552 : IF (PRESENT(op_pos)) THEN
700 2552 : op_pos_prv = op_pos
701 : ELSE
702 : op_pos_prv = 1
703 : END IF
704 :
705 2330 : SELECT CASE (op_pos_prv)
706 : CASE (1)
707 2330 : op_ij = potential_parameter%potential_type
708 : CASE (2)
709 2552 : op_jk = potential_parameter%potential_type
710 : END SELECT
711 :
712 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
713 2552 : .OR. op_ij == do_potential_mix_cl_trunc) THEN
714 1488 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
715 1488 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
716 1064 : ELSEIF (op_ij == do_potential_coulomb) THEN
717 436 : dr_ij = 1000000.0_dp
718 436 : dr_ik = 1000000.0_dp
719 : END IF
720 :
721 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
722 2552 : .OR. op_jk == do_potential_mix_cl_trunc) THEN
723 20 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
724 20 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
725 2532 : ELSEIF (op_jk == do_potential_coulomb) THEN
726 0 : dr_jk = 1000000.0_dp
727 0 : dr_ik = 1000000.0_dp
728 : END IF
729 :
730 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, natom=natom, &
731 2552 : dft_control=dft_control, para_env=para_env, cell=cell)
732 :
733 2552 : IF (do_kpoints_prv) THEN
734 234 : CPASSERT(PRESENT(cell_to_index))
735 234 : CPASSERT(ASSOCIATED(cell_to_index))
736 : ! nimg = dft_control%nimages
737 33696 : nimg = MAXVAL(cell_to_index)
738 234 : ncell_RI = nimg
739 234 : IF (do_hfx_kpoints_prv) THEN
740 222 : nimg = SIZE(t3c, 1)
741 222 : ncell_RI = SIZE(t3c, 2)
742 : END IF
743 : ELSE
744 : nimg = 1
745 : ncell_RI = 1
746 : END IF
747 :
748 2552 : IF (do_kpoints_prv) THEN
749 936 : kp_index_lbounds = LBOUND(cell_to_index)
750 936 : kp_index_ubounds = UBOUND(cell_to_index)
751 : END IF
752 :
753 : !Do a first loop over the nl and count the blocks present
754 10208 : ALLOCATE (nblk(nimg, ncell_RI))
755 14156 : nblk(:, :) = 0
756 :
757 2552 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
758 :
759 2552 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
760 :
761 1089046 : DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
762 : CALL get_3c_iterator_info(nl_3c_iter, iatom=iatom, ikind=ikind, jkind=jkind, kkind=kkind, &
763 1086494 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
764 :
765 4345976 : djk = NORM2(rjk)
766 4345976 : dij = NORM2(rij)
767 4345976 : dik = NORM2(rik)
768 :
769 1086494 : IF (do_kpoints_prv) THEN
770 :
771 801066 : IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
772 : ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
773 :
774 114438 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
775 114438 : IF (jcell > nimg .OR. jcell < 1) CYCLE
776 :
777 702612 : IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
778 : ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
779 :
780 96315 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
781 96315 : IF (kcell > nimg .OR. kcell < 1) CYCLE
782 66939 : IF (do_hfx_kpoints_prv) THEN
783 63309 : IF (dik > RI_range) CYCLE
784 : kcell = 1
785 : END IF
786 : ELSE
787 : jcell = 1; kcell = 1
788 : END IF
789 :
790 985236 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
791 985236 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
792 985236 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
793 :
794 985236 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
795 985236 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
796 985236 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
797 :
798 1080522 : nblk(jcell, kcell) = nblk(jcell, kcell) + 1
799 : END DO
800 2552 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
801 :
802 : !Do a second loop over the nl to give block indices
803 21812 : ALLOCATE (alloc_i(nimg, ncell_RI))
804 19260 : ALLOCATE (alloc_j(nimg, ncell_RI))
805 19260 : ALLOCATE (alloc_k(nimg, ncell_RI))
806 5180 : DO k_img = 1, ncell_RI
807 14156 : DO j_img = 1, nimg
808 23711 : ALLOCATE (alloc_i(j_img, k_img)%array(nblk(j_img, k_img)))
809 23711 : ALLOCATE (alloc_j(j_img, k_img)%array(nblk(j_img, k_img)))
810 26339 : ALLOCATE (alloc_k(j_img, k_img)%array(nblk(j_img, k_img)))
811 : END DO
812 : END DO
813 14156 : nblk(:, :) = 0
814 :
815 2552 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
816 :
817 2552 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)
818 :
819 1089046 : DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
820 : CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
821 : iatom=iatom, jatom=jatom, katom=katom, &
822 1086494 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
823 :
824 4345976 : djk = NORM2(rjk)
825 4345976 : dij = NORM2(rij)
826 4345976 : dik = NORM2(rik)
827 :
828 1086494 : IF (do_kpoints_prv) THEN
829 :
830 801066 : IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
831 : ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
832 :
833 114438 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
834 114438 : IF (jcell > nimg .OR. jcell < 1) CYCLE
835 :
836 702612 : IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
837 : ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
838 :
839 96315 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
840 96315 : IF (kcell > nimg .OR. kcell < 1) CYCLE
841 66939 : IF (do_hfx_kpoints_prv) THEN
842 63309 : IF (dik > RI_range) CYCLE
843 9550 : kcell = img_to_RI_cell_prv(kcell)
844 : END IF
845 : ELSE
846 : jcell = 1; kcell = 1
847 : END IF
848 :
849 3940944 : blk_idx = [iatom, jatom, katom]
850 985236 : IF (do_hfx_kpoints_prv) THEN
851 9550 : blk_idx(3) = (kcell - 1)*natom + katom
852 9550 : kcell = 1
853 : END IF
854 :
855 985236 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, kind_radius=kind_radius_i)
856 985236 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, kind_radius=kind_radius_j)
857 985236 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, kind_radius=kind_radius_k)
858 :
859 985236 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
860 985236 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
861 985236 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
862 :
863 497535 : nblk(jcell, kcell) = nblk(jcell, kcell) + 1
864 :
865 : !Note: there may be repeated indices due to periodic images => dbt_reserve_blocks takes care of it
866 497535 : alloc_i(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(1)
867 497535 : alloc_j(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(2)
868 1080522 : alloc_k(jcell, kcell)%array(nblk(jcell, kcell)) = blk_idx(3)
869 :
870 : END DO
871 2552 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
872 :
873 : !TODO: Parallelize creation of block list.
874 : !$OMP PARALLEL DEFAULT(NONE) SHARED(t3c,nimg,nblk,alloc_i,alloc_j,alloc_k,ncell_RI,cell_sym_prv) &
875 2552 : !$OMP PRIVATE(k_img,j_img,nblk_per_thread,A,b)
876 : DO k_img = 1, ncell_RI
877 : DO j_img = 1, nimg
878 : IF (cell_sym_prv .AND. j_img < k_img) CYCLE
879 : IF (ALLOCATED(alloc_i(j_img, k_img)%array)) THEN
880 : nblk_per_thread = nblk(j_img, k_img)/omp_get_num_threads() + 1
881 : a = omp_get_thread_num()*nblk_per_thread + 1
882 : b = MIN(a + nblk_per_thread, nblk(j_img, k_img))
883 : CALL dbt_reserve_blocks(t3c(j_img, k_img), &
884 : alloc_i(j_img, k_img)%array(a:b), &
885 : alloc_j(j_img, k_img)%array(a:b), &
886 : alloc_k(j_img, k_img)%array(a:b))
887 : END IF
888 : END DO
889 : END DO
890 : !$OMP END PARALLEL
891 :
892 2552 : CALL timestop(handle)
893 :
894 52448 : END SUBROUTINE alloc_block_3c
895 :
896 : ! **************************************************************************************************
897 : !> \brief Build 3-center derivative tensors
898 : !> \param t3c_der_i empty DBCSR tensor which will contain the 1st center derivatives
899 : !> \param t3c_der_k empty DBCSR tensor which will contain the 3rd center derivatives
900 : !> \param filter_eps Filter threshold for tensor blocks
901 : !> \param qs_env ...
902 : !> \param nl_3c 3-center neighborlist
903 : !> \param basis_i ...
904 : !> \param basis_j ...
905 : !> \param basis_k ...
906 : !> \param potential_parameter ...
907 : !> \param der_eps neglect integrals smaller than der_eps
908 : !> \param op_pos operator position.
909 : !> 1: calculate (i|jk) integrals,
910 : !> 2: calculate (ij|k) integrals
911 : !> \param do_kpoints ...
912 : !> this routine requires that libint has been static initialised somewhere else
913 : !> \param do_hfx_kpoints ...
914 : !> \param bounds_i ...
915 : !> \param bounds_j ...
916 : !> \param bounds_k ...
917 : !> \param RI_range ...
918 : !> \param img_to_RI_cell ...
919 : ! **************************************************************************************************
920 574 : SUBROUTINE build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, &
921 574 : nl_3c, basis_i, basis_j, basis_k, &
922 : potential_parameter, der_eps, &
923 : op_pos, do_kpoints, do_hfx_kpoints, &
924 : bounds_i, bounds_j, bounds_k, &
925 574 : RI_range, img_to_RI_cell)
926 :
927 : TYPE(dbt_type), DIMENSION(:, :, :), INTENT(INOUT) :: t3c_der_i, t3c_der_k
928 : REAL(KIND=dp), INTENT(IN) :: filter_eps
929 : TYPE(qs_environment_type), POINTER :: qs_env
930 : TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
931 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
932 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
933 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: der_eps
934 : INTEGER, INTENT(IN), OPTIONAL :: op_pos
935 : LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints
936 : INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
937 : REAL(dp), INTENT(IN), OPTIONAL :: RI_range
938 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: img_to_RI_cell
939 :
940 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_derivatives'
941 :
942 : INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
943 : block_start_k, egfi, handle, handle2, i, i_img, i_xyz, iatom, ibasis, ikind, ilist, imax, &
944 : iset, j_img, jatom, jcell, jkind, jset, katom, kcell, kkind, kset, m_max, max_ncoj, &
945 : max_nset, max_nsgfi, maxli, maxlj, maxlk, mepos, natom, nbasis, ncell_RI, ncoi, ncoj, &
946 : ncok, nimg, nseti, nsetj, nsetk, nthread, op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, &
947 : unit_id
948 574 : INTEGER, ALLOCATABLE, DIMENSION(:) :: img_to_RI_cell_prv
949 : INTEGER, DIMENSION(2) :: bo
950 : INTEGER, DIMENSION(3) :: blk_idx, blk_size, cell_j, cell_k, &
951 : kp_index_lbounds, kp_index_ubounds, sp
952 574 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
953 1148 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
954 574 : nsgfj, nsgfk
955 574 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
956 574 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
957 : LOGICAL :: do_hfx_kpoints_prv, do_kpoints_prv, &
958 : found, skip
959 : LOGICAL, DIMENSION(3) :: block_j_not_zero, block_k_not_zero, &
960 : der_j_zero, der_k_zero
961 : REAL(dp), DIMENSION(3) :: der_ext_i, der_ext_j, der_ext_k
962 : REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
963 : kind_radius_i, kind_radius_j, &
964 : kind_radius_k, prefac
965 574 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
966 574 : max_contraction_i, max_contraction_j, &
967 574 : max_contraction_k
968 574 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dijk_contr, dummy_block_t
969 574 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: block_t_i, block_t_j, block_t_k, dijk_j, &
970 1148 : dijk_k, tmp_ijk_i, tmp_ijk_j
971 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
972 574 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
973 574 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
974 1148 : sphi_k, zeti, zetj, zetk
975 574 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
976 1148 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
977 : TYPE(cp_libint_t) :: lib
978 4018 : TYPE(dbt_type) :: t3c_tmp
979 574 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t3c_template
980 574 : TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :, :) :: t3c_der_j
981 : TYPE(dft_control_type), POINTER :: dft_control
982 : TYPE(gto_basis_set_type), POINTER :: basis_set
983 : TYPE(kpoint_type), POINTER :: kpoints
984 : TYPE(mp_para_env_type), POINTER :: para_env
985 : TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
986 574 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
987 :
988 574 : CALL timeset(routineN, handle)
989 :
990 574 : IF (PRESENT(do_kpoints)) THEN
991 82 : do_kpoints_prv = do_kpoints
992 : ELSE
993 : do_kpoints_prv = .FALSE.
994 : END IF
995 :
996 574 : IF (PRESENT(do_hfx_kpoints)) THEN
997 82 : do_hfx_kpoints_prv = do_hfx_kpoints
998 : ELSE
999 : do_hfx_kpoints_prv = .FALSE.
1000 : END IF
1001 82 : IF (do_hfx_kpoints_prv) THEN
1002 82 : CPASSERT(do_kpoints_prv)
1003 82 : CPASSERT(PRESENT(RI_range))
1004 82 : CPASSERT(PRESENT(img_to_RI_cell))
1005 : END IF
1006 :
1007 492 : IF (PRESENT(img_to_RI_cell)) THEN
1008 246 : ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
1009 3834 : img_to_RI_cell_prv(:) = img_to_RI_cell
1010 : END IF
1011 :
1012 574 : op_ij = do_potential_id; op_jk = do_potential_id
1013 :
1014 574 : IF (PRESENT(op_pos)) THEN
1015 574 : op_pos_prv = op_pos
1016 : ELSE
1017 0 : op_pos_prv = 1
1018 : END IF
1019 :
1020 492 : SELECT CASE (op_pos_prv)
1021 : CASE (1)
1022 492 : op_ij = potential_parameter%potential_type
1023 : CASE (2)
1024 574 : op_jk = potential_parameter%potential_type
1025 : END SELECT
1026 :
1027 574 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
1028 :
1029 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
1030 574 : .OR. op_ij == do_potential_mix_cl_trunc) THEN
1031 50 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
1032 50 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1033 524 : ELSEIF (op_ij == do_potential_coulomb) THEN
1034 244 : dr_ij = 1000000.0_dp
1035 244 : dr_ik = 1000000.0_dp
1036 : END IF
1037 :
1038 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
1039 574 : .OR. op_jk == do_potential_mix_cl_trunc) THEN
1040 8 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
1041 8 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1042 566 : ELSEIF (op_jk == do_potential_coulomb) THEN
1043 0 : dr_jk = 1000000.0_dp
1044 0 : dr_ik = 1000000.0_dp
1045 : END IF
1046 :
1047 574 : NULLIFY (qs_kind_set, atomic_kind_set)
1048 :
1049 : ! get stuff
1050 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1051 574 : natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
1052 :
1053 574 : IF (do_kpoints_prv) THEN
1054 82 : nimg = dft_control%nimages
1055 82 : ncell_RI = nimg
1056 82 : IF (do_hfx_kpoints_prv) THEN
1057 82 : nimg = SIZE(t3c_der_k, 1)
1058 82 : ncell_RI = SIZE(t3c_der_k, 2)
1059 : END IF
1060 82 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
1061 : ELSE
1062 : nimg = 1
1063 : ncell_RI = 1
1064 : END IF
1065 :
1066 574 : IF (do_hfx_kpoints_prv) THEN
1067 82 : CPASSERT(op_pos_prv == 2)
1068 : ELSE
1069 1968 : CPASSERT(ALL(SHAPE(t3c_der_k) == [nimg, ncell_RI, 3]))
1070 : END IF
1071 :
1072 8910 : ALLOCATE (t3c_template(nimg, ncell_RI))
1073 1148 : DO j_img = 1, ncell_RI
1074 3744 : DO i_img = 1, nimg
1075 3170 : CALL dbt_create(t3c_der_i(i_img, j_img, 1), t3c_template(i_img, j_img))
1076 : END DO
1077 : END DO
1078 :
1079 : CALL alloc_block_3c(t3c_template, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, &
1080 : op_pos=op_pos_prv, do_kpoints=do_kpoints, do_hfx_kpoints=do_hfx_kpoints, &
1081 : bounds_i=bounds_i, bounds_j=bounds_j, bounds_k=bounds_k, &
1082 1066 : RI_range=RI_range, img_to_RI_cell=img_to_RI_cell, cell_to_index=cell_to_index)
1083 2296 : DO i_xyz = 1, 3
1084 4018 : DO j_img = 1, ncell_RI
1085 11232 : DO i_img = 1, nimg
1086 7788 : CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_i(i_img, j_img, i_xyz))
1087 9510 : CALL dbt_copy(t3c_template(i_img, j_img), t3c_der_k(i_img, j_img, i_xyz))
1088 : END DO
1089 : END DO
1090 : END DO
1091 :
1092 1148 : DO j_img = 1, ncell_RI
1093 3744 : DO i_img = 1, nimg
1094 3170 : CALL dbt_destroy(t3c_template(i_img, j_img))
1095 : END DO
1096 : END DO
1097 3170 : DEALLOCATE (t3c_template)
1098 :
1099 574 : IF (nl_3c%sym == symmetric_jk) THEN
1100 9840 : ALLOCATE (t3c_der_j(nimg, ncell_RI, 3))
1101 1968 : DO i_xyz = 1, 3
1102 3444 : DO j_img = 1, ncell_RI
1103 4428 : DO i_img = 1, nimg
1104 1476 : CALL dbt_create(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
1105 2952 : CALL dbt_copy(t3c_der_k(i_img, j_img, i_xyz), t3c_der_j(i_img, j_img, i_xyz))
1106 : END DO
1107 : END DO
1108 : END DO
1109 : END IF
1110 :
1111 : !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
1112 574 : nbasis = SIZE(basis_i)
1113 574 : max_nsgfi = 0
1114 574 : max_nset = 0
1115 574 : maxli = 0
1116 1638 : DO ibasis = 1, nbasis
1117 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
1118 1064 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
1119 1064 : maxli = MAX(maxli, imax)
1120 1064 : max_nset = MAX(max_nset, iset)
1121 5792 : max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
1122 : END DO
1123 : max_ncoj = 0
1124 : maxlj = 0
1125 1638 : DO ibasis = 1, nbasis
1126 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
1127 1064 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
1128 1064 : maxlj = MAX(maxlj, imax)
1129 1064 : max_nset = MAX(max_nset, jset)
1130 4484 : max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
1131 : END DO
1132 : maxlk = 0
1133 1638 : DO ibasis = 1, nbasis
1134 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
1135 1064 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
1136 1064 : maxlk = MAX(maxlk, imax)
1137 1638 : max_nset = MAX(max_nset, kset)
1138 : END DO
1139 574 : m_max = maxli + maxlj + maxlk + 1
1140 :
1141 : !To minimize expensive memory ops and generally optimize contraction, pre-allocate
1142 : !contiguous sphi arrays (and transposed in the case of sphi_i)
1143 :
1144 574 : NULLIFY (tspj, spi, spk)
1145 26756 : ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
1146 :
1147 1638 : DO ibasis = 1, nbasis
1148 7962 : DO iset = 1, max_nset
1149 6324 : NULLIFY (spi(iset, ibasis)%array)
1150 6324 : NULLIFY (tspj(iset, ibasis)%array)
1151 7388 : NULLIFY (spk(iset, ibasis)%array)
1152 : END DO
1153 : END DO
1154 :
1155 2296 : DO ilist = 1, 3
1156 5488 : DO ibasis = 1, nbasis
1157 3192 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
1158 3192 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
1159 3192 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
1160 :
1161 14608 : DO iset = 1, basis_set%nset
1162 :
1163 9694 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
1164 9694 : sgfi = basis_set%first_sgf(1, iset)
1165 9694 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
1166 :
1167 12886 : IF (ilist == 1) THEN
1168 16616 : ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1169 3075498 : spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1170 :
1171 5540 : ELSE IF (ilist == 2) THEN
1172 11384 : ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
1173 177774 : tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
1174 :
1175 : ELSE
1176 10776 : ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1177 838246 : spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1178 : END IF
1179 :
1180 : END DO !iset
1181 : END DO !ibasis
1182 : END DO !ilist
1183 :
1184 : !Init the truncated Coulomb operator
1185 : IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated .OR. &
1186 574 : op_ij == do_potential_mix_cl_trunc .OR. op_jk == do_potential_mix_cl_trunc) THEN
1187 :
1188 58 : IF (m_max > get_lmax_init()) THEN
1189 8 : IF (para_env%mepos == 0) THEN
1190 4 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
1191 : END IF
1192 8 : CALL init(m_max, unit_id, para_env%mepos, para_env)
1193 8 : IF (para_env%mepos == 0) THEN
1194 4 : CALL close_file(unit_id)
1195 : END IF
1196 : END IF
1197 : END IF
1198 :
1199 574 : CALL init_md_ftable(nmax=m_max)
1200 :
1201 574 : IF (do_kpoints_prv) THEN
1202 328 : kp_index_lbounds = LBOUND(cell_to_index)
1203 328 : kp_index_ubounds = UBOUND(cell_to_index)
1204 : END IF
1205 :
1206 574 : nthread = 1
1207 574 : !$ nthread = omp_get_max_threads()
1208 :
1209 : !$OMP PARALLEL DEFAULT(NONE) &
1210 : !$OMP SHARED (nthread,do_kpoints_prv,kp_index_lbounds,kp_index_ubounds,maxli,maxlk,maxlj,bounds_i,&
1211 : !$OMP bounds_j,bounds_k,nimg,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,op_pos_prv,&
1212 : !$OMP potential_parameter,der_eps,tspj,spi,spk,cell_to_index,RI_range,natom,nl_3c,&
1213 : !$OMP t3c_der_i,t3c_der_k,t3c_der_j,ncell_RI,img_to_RI_cell_prv,do_hfx_kpoints_prv) &
1214 : !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,cell_j,cell_k,&
1215 : !$OMP prefac,jcell,kcell,first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
1216 : !$OMP sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
1217 : !$OMP set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
1218 : !$OMP rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
1219 : !$OMP sgfk,dijk_j,dijk_k,ri,rj,rk,max_contraction_i,max_contraction_j,blk_idx,&
1220 : !$OMP max_contraction_k,iset,jset,kset,block_t_i,blk_size,dijk_contr,cpp_buffer,ccp_buffer,&
1221 : !$OMP block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
1222 : !$OMP dummy_block_t,sp,handle2,mepos,bo,block_t_k,der_ext_i,der_ext_j,der_ext_k,tmp_ijk_i,&
1223 574 : !$OMP block_k_not_zero,der_k_zero,skip,der_j_zero,block_t_j,block_j_not_zero,tmp_ijk_j,i)
1224 :
1225 : mepos = 0
1226 : !$ mepos = omp_get_thread_num()
1227 :
1228 : CALL cp_libint_init_3eri1(lib, MAX(maxli, maxlj, maxlk))
1229 : CALL cp_libint_set_contrdepth(lib, 1)
1230 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
1231 :
1232 : !We split the provided bounds among the threads such that each threads works on a different set of atoms
1233 : IF (PRESENT(bounds_i)) THEN
1234 : bo = get_limit(bounds_i(2) - bounds_i(1) + 1, nthread, mepos)
1235 : bo(:) = bo(:) + bounds_i(1) - 1
1236 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
1237 : ELSE IF (PRESENT(bounds_j)) THEN
1238 : bo = get_limit(bounds_j(2) - bounds_j(1) + 1, nthread, mepos)
1239 : bo(:) = bo(:) + bounds_j(1) - 1
1240 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bo, bounds_k)
1241 : ELSE IF (PRESENT(bounds_k)) THEN
1242 : bo = get_limit(bounds_k(2) - bounds_k(1) + 1, nthread, mepos)
1243 : bo(:) = bo(:) + bounds_k(1) - 1
1244 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bo)
1245 : ELSE
1246 : bo = get_limit(natom, nthread, mepos)
1247 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
1248 : END IF
1249 :
1250 : skip = .FALSE.
1251 : IF (bo(1) > bo(2)) skip = .TRUE.
1252 :
1253 : DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
1254 : CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
1255 : iatom=iatom, jatom=jatom, katom=katom, &
1256 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
1257 : IF (skip) EXIT
1258 :
1259 : djk = NORM2(rjk)
1260 : dij = NORM2(rij)
1261 : dik = NORM2(rik)
1262 :
1263 : IF (do_kpoints_prv) THEN
1264 : prefac = 0.5_dp
1265 : ELSEIF (nl_3c%sym == symmetric_jk) THEN
1266 : IF (jatom == katom) THEN
1267 : prefac = 0.5_dp
1268 : ELSE
1269 : prefac = 1.0_dp
1270 : END IF
1271 : ELSE
1272 : prefac = 1.0_dp
1273 : END IF
1274 : IF (do_hfx_kpoints_prv) prefac = 1.0_dp
1275 :
1276 : IF (do_kpoints_prv) THEN
1277 :
1278 : IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
1279 : ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
1280 :
1281 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1282 : IF (jcell > nimg .OR. jcell < 1) CYCLE
1283 :
1284 : IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
1285 : ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
1286 :
1287 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
1288 : IF (kcell > nimg .OR. kcell < 1) CYCLE
1289 : !In case of HFX k-points, we only consider P^k if d_ik <= RI_range
1290 : IF (do_hfx_kpoints_prv) THEN
1291 : IF (dik > RI_range) CYCLE
1292 : kcell = img_to_RI_cell_prv(kcell)
1293 : END IF
1294 : ELSE
1295 : jcell = 1; kcell = 1
1296 : END IF
1297 :
1298 : blk_idx = [iatom, jatom, katom]
1299 : IF (do_hfx_kpoints_prv) THEN
1300 : blk_idx(3) = (kcell - 1)*natom + katom
1301 : kcell = 1
1302 : END IF
1303 :
1304 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
1305 : npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
1306 : sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
1307 :
1308 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
1309 : npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
1310 : sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
1311 :
1312 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
1313 : npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
1314 : sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
1315 :
1316 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
1317 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
1318 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
1319 :
1320 : IF (PRESENT(der_eps)) THEN
1321 : ALLOCATE (max_contraction_i(nseti))
1322 : DO iset = 1, nseti
1323 : sgfi = first_sgf_i(1, iset)
1324 : max_contraction_i(iset) = MAXVAL([(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)])
1325 : END DO
1326 :
1327 : ALLOCATE (max_contraction_j(nsetj))
1328 : DO jset = 1, nsetj
1329 : sgfj = first_sgf_j(1, jset)
1330 : max_contraction_j(jset) = MAXVAL([(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)])
1331 : END DO
1332 :
1333 : ALLOCATE (max_contraction_k(nsetk))
1334 : DO kset = 1, nsetk
1335 : sgfk = first_sgf_k(1, kset)
1336 : max_contraction_k(kset) = MAXVAL([(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)])
1337 : END DO
1338 : END IF
1339 :
1340 : CALL dbt_blk_sizes(t3c_der_i(jcell, kcell, 1), blk_idx, blk_size)
1341 :
1342 : ALLOCATE (block_t_i(blk_size(2), blk_size(3), blk_size(1), 3), SOURCE=0.0_dp)
1343 : ALLOCATE (block_t_j(blk_size(2), blk_size(3), blk_size(1), 3), SOURCE=0.0_dp)
1344 : ALLOCATE (block_t_k(blk_size(2), blk_size(3), blk_size(1), 3), SOURCE=0.0_dp)
1345 :
1346 : block_j_not_zero = .FALSE.
1347 : block_k_not_zero = .FALSE.
1348 :
1349 : DO iset = 1, nseti
1350 :
1351 : DO jset = 1, nsetj
1352 :
1353 : IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE
1354 :
1355 : DO kset = 1, nsetk
1356 :
1357 : IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE
1358 : IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE
1359 :
1360 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
1361 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
1362 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
1363 :
1364 : sgfi = first_sgf_i(1, iset)
1365 : sgfj = first_sgf_j(1, jset)
1366 : sgfk = first_sgf_k(1, kset)
1367 :
1368 : IF (ncoj*ncok*ncoi > 0) THEN
1369 : ALLOCATE (dijk_j(ncoj, ncok, ncoi, 3), SOURCE=0.0_dp)
1370 : ALLOCATE (dijk_k(ncoj, ncok, ncoi, 3), SOURCE=0.0_dp)
1371 :
1372 : der_j_zero = .FALSE.
1373 : der_k_zero = .FALSE.
1374 :
1375 : !need positions for libint. Only relative positions are needed => set ri to 0.0
1376 : ri = 0.0_dp
1377 : rj = rij ! ri + rij
1378 : rk = rik ! ri + rik
1379 :
1380 : IF (op_pos_prv == 1) THEN
1381 : CALL eri_3center_derivs(dijk_j, dijk_k, &
1382 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
1383 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
1384 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
1385 : djk, dij, dik, lib, potential_parameter, &
1386 : der_abc_1_ext=der_ext_j, der_abc_2_ext=der_ext_k)
1387 : ELSE
1388 : ALLOCATE (tmp_ijk_i(ncoi, ncoj, ncok, 3), SOURCE=0.0_dp)
1389 : ALLOCATE (tmp_ijk_j(ncoi, ncoj, ncok, 3), SOURCE=0.0_dp)
1390 :
1391 : CALL eri_3center_derivs(tmp_ijk_i, tmp_ijk_j, &
1392 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
1393 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
1394 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
1395 : dij, dik, djk, lib, potential_parameter, &
1396 : der_abc_1_ext=der_ext_i, der_abc_2_ext=der_ext_j)
1397 :
1398 : !TODO: is that inefficient?
1399 : der_ext_k = 0
1400 : DO i_xyz = 1, 3
1401 : DO i = 1, ncoi
1402 : dijk_j(:, :, i, i_xyz) = tmp_ijk_j(i, :, :, i_xyz)
1403 : dijk_k(:, :, i, i_xyz) = -(dijk_j(:, :, i, i_xyz) + tmp_ijk_i(i, :, :, i_xyz))
1404 : der_ext_k(i_xyz) = MAX(der_ext_k(i_xyz), MAXVAL(ABS(dijk_k(:, :, i, i_xyz))))
1405 : END DO
1406 : END DO
1407 : DEALLOCATE (tmp_ijk_i, tmp_ijk_j)
1408 :
1409 : END IF
1410 :
1411 : IF (PRESENT(der_eps)) THEN
1412 : DO i_xyz = 1, 3
1413 : IF (der_eps > der_ext_j(i_xyz)*(max_contraction_i(iset)* &
1414 : max_contraction_j(jset)* &
1415 : max_contraction_k(kset))) THEN
1416 : der_j_zero(i_xyz) = .TRUE.
1417 : END IF
1418 : END DO
1419 :
1420 : DO i_xyz = 1, 3
1421 : IF (der_eps > der_ext_k(i_xyz)*(max_contraction_i(iset)* &
1422 : max_contraction_j(jset)* &
1423 : max_contraction_k(kset))) THEN
1424 : der_k_zero(i_xyz) = .TRUE.
1425 : END IF
1426 : END DO
1427 : IF (ALL(der_j_zero) .AND. ALL(der_k_zero)) THEN
1428 : DEALLOCATE (dijk_j, dijk_k)
1429 : CYCLE
1430 : END IF
1431 : END IF
1432 :
1433 : ALLOCATE (dijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
1434 :
1435 : block_start_j = sgfj
1436 : block_end_j = sgfj + nsgfj(jset) - 1
1437 : block_start_k = sgfk
1438 : block_end_k = sgfk + nsgfk(kset) - 1
1439 : block_start_i = sgfi
1440 : block_end_i = sgfi + nsgfi(iset) - 1
1441 :
1442 : DO i_xyz = 1, 3
1443 : IF (der_j_zero(i_xyz)) CYCLE
1444 :
1445 : block_j_not_zero(i_xyz) = .TRUE.
1446 : CALL abc_contract_xsmm(dijk_contr, dijk_j(:, :, :, i_xyz), tspj(jset, jkind)%array, &
1447 : spk(kset, kkind)%array, spi(iset, ikind)%array, &
1448 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
1449 : nsgfi(iset), cpp_buffer, ccp_buffer, prefac)
1450 :
1451 : block_t_j(block_start_j:block_end_j, &
1452 : block_start_k:block_end_k, &
1453 : block_start_i:block_end_i, i_xyz) = &
1454 : block_t_j(block_start_j:block_end_j, &
1455 : block_start_k:block_end_k, &
1456 : block_start_i:block_end_i, i_xyz) + &
1457 : dijk_contr(:, :, :)
1458 :
1459 : END DO
1460 :
1461 : DO i_xyz = 1, 3
1462 : IF (der_k_zero(i_xyz)) CYCLE
1463 :
1464 : block_k_not_zero(i_xyz) = .TRUE.
1465 : CALL abc_contract_xsmm(dijk_contr, dijk_k(:, :, :, i_xyz), tspj(jset, jkind)%array, &
1466 : spk(kset, kkind)%array, spi(iset, ikind)%array, &
1467 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
1468 : nsgfi(iset), cpp_buffer, ccp_buffer, prefac)
1469 :
1470 : block_t_k(block_start_j:block_end_j, &
1471 : block_start_k:block_end_k, &
1472 : block_start_i:block_end_i, i_xyz) = &
1473 : block_t_k(block_start_j:block_end_j, &
1474 : block_start_k:block_end_k, &
1475 : block_start_i:block_end_i, i_xyz) + &
1476 : dijk_contr(:, :, :)
1477 :
1478 : END DO
1479 :
1480 : DEALLOCATE (dijk_j, dijk_k, dijk_contr)
1481 : END IF ! number of triples > 0
1482 : END DO
1483 : END DO
1484 : END DO
1485 :
1486 : CALL timeset(routineN//"_put_dbcsr", handle2)
1487 : !$OMP CRITICAL(qs_tensors_t3c_der_i)
1488 : sp = SHAPE(block_t_i(:, :, :, 1))
1489 : sp([2, 3, 1]) = sp
1490 :
1491 : DO i_xyz = 1, 3
1492 : !Derivatives wrt to center i can be obtained by translational invariance
1493 : IF ((.NOT. block_j_not_zero(i_xyz)) .AND. (.NOT. block_k_not_zero(i_xyz))) CYCLE
1494 : block_t_i(:, :, :, i_xyz) = -(block_t_j(:, :, :, i_xyz) + block_t_k(:, :, :, i_xyz))
1495 :
1496 : CALL dbt_put_block(t3c_der_i(jcell, kcell, i_xyz), blk_idx, sp, &
1497 : RESHAPE(block_t_i(:, :, :, i_xyz), SHAPE=sp, ORDER=[2, 3, 1]), &
1498 : summation=.TRUE.)
1499 : END DO
1500 :
1501 : sp = SHAPE(block_t_k(:, :, :, 1))
1502 : sp([2, 3, 1]) = sp
1503 :
1504 : DO i_xyz = 1, 3
1505 : IF (.NOT. block_k_not_zero(i_xyz)) CYCLE
1506 : CALL dbt_put_block(t3c_der_k(jcell, kcell, i_xyz), blk_idx, sp, &
1507 : RESHAPE(block_t_k(:, :, :, i_xyz), SHAPE=sp, ORDER=[2, 3, 1]), &
1508 : summation=.TRUE.)
1509 : END DO
1510 : !$OMP END CRITICAL(qs_tensors_t3c_der_i)
1511 :
1512 : IF (nl_3c%sym == symmetric_jk) THEN
1513 : sp = SHAPE(block_t_j(:, :, :, 1))
1514 : sp([2, 3, 1]) = sp
1515 : !$OMP CRITICAL(qs_tensors_t3c_der_j)
1516 : DO i_xyz = 1, 3
1517 : IF (.NOT. block_j_not_zero(i_xyz)) CYCLE
1518 : CALL dbt_put_block(t3c_der_j(jcell, kcell, i_xyz), blk_idx, sp, &
1519 : RESHAPE(block_t_j(:, :, :, i_xyz), SHAPE=sp, ORDER=[2, 3, 1]), &
1520 : summation=.TRUE.)
1521 : END DO
1522 : !$OMP END CRITICAL(qs_tensors_t3c_der_j)
1523 : END IF
1524 :
1525 : CALL timestop(handle2)
1526 :
1527 : DEALLOCATE (block_t_i, block_t_j, block_t_k)
1528 : IF (PRESENT(der_eps)) THEN
1529 : DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k)
1530 : END IF
1531 : END DO
1532 :
1533 : IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
1534 : IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
1535 :
1536 : CALL cp_libint_cleanup_3eri1(lib)
1537 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
1538 : !$OMP END PARALLEL
1539 :
1540 574 : IF (do_kpoints_prv .AND. .NOT. do_hfx_kpoints_prv) THEN
1541 0 : DO i_xyz = 1, 3
1542 0 : DO kcell = 1, nimg
1543 0 : DO jcell = 1, nimg
1544 : ! need half of filter eps because afterwards we add transposed tensor
1545 0 : CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps/2)
1546 0 : CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps/2)
1547 : END DO
1548 : END DO
1549 : END DO
1550 :
1551 574 : ELSEIF (nl_3c%sym == symmetric_jk) THEN
1552 : !Add the transpose of t3c_der_j to t3c_der_k to get the fully populated tensor
1553 492 : CALL dbt_create(t3c_der_k(1, 1, 1), t3c_tmp)
1554 1968 : DO i_xyz = 1, 3
1555 3444 : DO kcell = 1, nimg
1556 4428 : DO jcell = 1, nimg
1557 : CALL dbt_copy(t3c_der_j(jcell, kcell, i_xyz), t3c_der_k(jcell, kcell, i_xyz), &
1558 1476 : order=[1, 3, 2], move_data=.TRUE., summation=.TRUE.)
1559 1476 : CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
1560 :
1561 1476 : CALL dbt_copy(t3c_der_i(jcell, kcell, i_xyz), t3c_tmp)
1562 : CALL dbt_copy(t3c_tmp, t3c_der_i(jcell, kcell, i_xyz), &
1563 1476 : order=[1, 3, 2], move_data=.TRUE., summation=.TRUE.)
1564 2952 : CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
1565 : END DO
1566 : END DO
1567 : END DO
1568 492 : CALL dbt_destroy(t3c_tmp)
1569 :
1570 82 : ELSEIF (nl_3c%sym == symmetric_none) THEN
1571 328 : DO i_xyz = 1, 3
1572 574 : DO kcell = 1, ncell_RI
1573 6804 : DO jcell = 1, nimg
1574 6312 : CALL dbt_filter(t3c_der_i(jcell, kcell, i_xyz), filter_eps)
1575 6558 : CALL dbt_filter(t3c_der_k(jcell, kcell, i_xyz), filter_eps)
1576 : END DO
1577 : END DO
1578 : END DO
1579 : ELSE
1580 0 : CPABORT("requested symmetric case not implemented")
1581 : END IF
1582 :
1583 574 : IF (nl_3c%sym == symmetric_jk) THEN
1584 1968 : DO i_xyz = 1, 3
1585 3444 : DO j_img = 1, nimg
1586 4428 : DO i_img = 1, nimg
1587 2952 : CALL dbt_destroy(t3c_der_j(i_img, j_img, i_xyz))
1588 : END DO
1589 : END DO
1590 : END DO
1591 : END IF
1592 :
1593 3934 : DO iset = 1, max_nset
1594 10258 : DO ibasis = 1, nbasis
1595 6324 : IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
1596 6324 : IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
1597 9684 : IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
1598 : END DO
1599 : END DO
1600 :
1601 574 : DEALLOCATE (spi, tspj, spk)
1602 :
1603 574 : CALL timestop(handle)
1604 6642 : END SUBROUTINE build_3c_derivatives
1605 :
1606 : ! **************************************************************************************************
1607 : !> \brief Calculates the 3c virial contributions on the fly
1608 : !> \param work_virial ...
1609 : !> \param t3c_trace the tensor with which the trace should be taken
1610 : !> \param pref ...
1611 : !> \param qs_env ...
1612 : !> \param nl_3c 3-center neighborlist, with a distribution matching that of t3c_trace
1613 : !> \param basis_i ...
1614 : !> \param basis_j ...
1615 : !> \param basis_k ...
1616 : !> \param potential_parameter ...
1617 : !> \param der_eps neglect integrals smaller than der_eps
1618 : !> \param op_pos operator position.
1619 : !> 1: calculate (i|jk) integrals,
1620 : !> 2: calculate (ij|k) integrals
1621 : !> this routine requires that libint has been static initialised somewhere else
1622 : ! **************************************************************************************************
1623 16 : SUBROUTINE calc_3c_virial(work_virial, t3c_trace, pref, qs_env, &
1624 16 : nl_3c, basis_i, basis_j, basis_k, &
1625 : potential_parameter, der_eps, op_pos)
1626 :
1627 : REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: work_virial
1628 : TYPE(dbt_type), INTENT(INOUT) :: t3c_trace
1629 : REAL(KIND=dp), INTENT(IN) :: pref
1630 : TYPE(qs_environment_type), POINTER :: qs_env
1631 : TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
1632 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
1633 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
1634 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: der_eps
1635 : INTEGER, INTENT(IN), OPTIONAL :: op_pos
1636 :
1637 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_3c_virial'
1638 :
1639 : INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
1640 : block_start_k, egfi, handle, i, i_xyz, iatom, ibasis, ikind, ilist, imax, iset, j_xyz, &
1641 : jatom, jkind, jset, katom, kkind, kset, m_max, max_ncoj, max_nset, max_nsgfi, maxli, &
1642 : maxlj, maxlk, mepos, natom, nbasis, ncoi, ncoj, ncok, nseti, nsetj, nsetk, nthread, &
1643 : op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id
1644 : INTEGER, DIMENSION(2) :: bo
1645 : INTEGER, DIMENSION(3) :: blk_size, sp
1646 16 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
1647 32 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
1648 16 : nsgfj, nsgfk
1649 16 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
1650 : LOGICAL :: found, skip
1651 : LOGICAL, DIMENSION(3) :: block_j_not_zero, block_k_not_zero, &
1652 : der_j_zero, der_k_zero
1653 : REAL(dp) :: force
1654 : REAL(dp), DIMENSION(3) :: der_ext_i, der_ext_j, der_ext_k
1655 : REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
1656 : kind_radius_i, kind_radius_j, &
1657 : kind_radius_k
1658 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
1659 16 : max_contraction_i, max_contraction_j, &
1660 16 : max_contraction_k
1661 32 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ablock, dijk_contr, tmp_block
1662 16 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: block_t_i, block_t_j, block_t_k, dijk_j, &
1663 16 : dijk_k
1664 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk, scoord
1665 16 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
1666 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
1667 16 : sphi_k, zeti, zetj, zetk
1668 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1669 : TYPE(cell_type), POINTER :: cell
1670 16 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
1671 : TYPE(cp_libint_t) :: lib
1672 : TYPE(dft_control_type), POINTER :: dft_control
1673 : TYPE(gto_basis_set_type), POINTER :: basis_set
1674 : TYPE(mp_para_env_type), POINTER :: para_env
1675 : TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
1676 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1677 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1678 :
1679 16 : CALL timeset(routineN, handle)
1680 :
1681 16 : op_ij = do_potential_id; op_jk = do_potential_id
1682 :
1683 16 : IF (PRESENT(op_pos)) THEN
1684 16 : op_pos_prv = op_pos
1685 : ELSE
1686 : op_pos_prv = 1
1687 : END IF
1688 16 : CPASSERT(op_pos == 1)
1689 16 : CPASSERT(.NOT. nl_3c%sym == symmetric_jk)
1690 :
1691 16 : SELECT CASE (op_pos_prv)
1692 : CASE (1)
1693 16 : op_ij = potential_parameter%potential_type
1694 : CASE (2)
1695 16 : op_jk = potential_parameter%potential_type
1696 : END SELECT
1697 :
1698 16 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
1699 :
1700 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
1701 16 : .OR. op_ij == do_potential_mix_cl_trunc) THEN
1702 6 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
1703 6 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1704 10 : ELSEIF (op_ij == do_potential_coulomb) THEN
1705 0 : dr_ij = 1000000.0_dp
1706 0 : dr_ik = 1000000.0_dp
1707 : END IF
1708 :
1709 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
1710 16 : .OR. op_jk == do_potential_mix_cl_trunc) THEN
1711 0 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
1712 0 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
1713 16 : ELSEIF (op_jk == do_potential_coulomb) THEN
1714 0 : dr_jk = 1000000.0_dp
1715 0 : dr_ik = 1000000.0_dp
1716 : END IF
1717 :
1718 16 : NULLIFY (qs_kind_set, atomic_kind_set)
1719 :
1720 : ! get stuff
1721 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1722 : natom=natom, dft_control=dft_control, para_env=para_env, &
1723 16 : particle_set=particle_set, cell=cell)
1724 :
1725 : !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
1726 16 : nbasis = SIZE(basis_i)
1727 16 : max_nsgfi = 0
1728 16 : max_nset = 0
1729 16 : maxli = 0
1730 48 : DO ibasis = 1, nbasis
1731 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
1732 32 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
1733 32 : maxli = MAX(maxli, imax)
1734 32 : max_nset = MAX(max_nset, iset)
1735 182 : max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
1736 : END DO
1737 : max_ncoj = 0
1738 : maxlj = 0
1739 48 : DO ibasis = 1, nbasis
1740 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
1741 32 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
1742 32 : maxlj = MAX(maxlj, imax)
1743 32 : max_nset = MAX(max_nset, jset)
1744 138 : max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
1745 : END DO
1746 : maxlk = 0
1747 48 : DO ibasis = 1, nbasis
1748 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
1749 32 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
1750 32 : maxlk = MAX(maxlk, imax)
1751 48 : max_nset = MAX(max_nset, kset)
1752 : END DO
1753 16 : m_max = maxli + maxlj + maxlk + 1
1754 :
1755 : !To minimize expensive memory ops and generally optimize contraction, pre-allocate
1756 : !contiguous sphi arrays (and transposed in the case of sphi_i)
1757 :
1758 16 : NULLIFY (tspj, spi, spk)
1759 920 : ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
1760 :
1761 48 : DO ibasis = 1, nbasis
1762 280 : DO iset = 1, max_nset
1763 232 : NULLIFY (spi(iset, ibasis)%array)
1764 232 : NULLIFY (tspj(iset, ibasis)%array)
1765 :
1766 264 : NULLIFY (spk(iset, ibasis)%array)
1767 : END DO
1768 : END DO
1769 :
1770 64 : DO ilist = 1, 3
1771 160 : DO ibasis = 1, nbasis
1772 96 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
1773 96 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
1774 96 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
1775 :
1776 458 : DO iset = 1, basis_set%nset
1777 :
1778 314 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
1779 314 : sgfi = basis_set%first_sgf(1, iset)
1780 314 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
1781 :
1782 410 : IF (ilist == 1) THEN
1783 536 : ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1784 177990 : spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1785 :
1786 180 : ELSE IF (ilist == 2) THEN
1787 360 : ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
1788 5194 : tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
1789 :
1790 : ELSE
1791 360 : ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
1792 4442 : spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
1793 : END IF
1794 :
1795 : END DO !iset
1796 : END DO !ibasis
1797 : END DO !ilist
1798 :
1799 : !Init the truncated Coulomb operator
1800 : IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated .OR. &
1801 16 : op_ij == do_potential_mix_cl_trunc .OR. op_jk == do_potential_mix_cl_trunc) THEN
1802 :
1803 6 : IF (m_max > get_lmax_init()) THEN
1804 0 : IF (para_env%mepos == 0) THEN
1805 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
1806 : END IF
1807 0 : CALL init(m_max, unit_id, para_env%mepos, para_env)
1808 0 : IF (para_env%mepos == 0) THEN
1809 0 : CALL close_file(unit_id)
1810 : END IF
1811 : END IF
1812 : END IF
1813 :
1814 16 : CALL init_md_ftable(nmax=m_max)
1815 :
1816 16 : nthread = 1
1817 16 : !$ nthread = omp_get_max_threads()
1818 :
1819 : !$OMP PARALLEL DEFAULT(NONE) &
1820 : !$OMP SHARED (nthread,maxli,maxlk,maxlj,work_virial,pref,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,&
1821 : !$OMP ncoset,potential_parameter,der_eps,tspj,spi,spk,natom,nl_3c,t3c_trace,cell,particle_set) &
1822 : !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,i_xyz,j_xyz,first_sgf_i,&
1823 : !$OMP lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,sphi_i,zeti,kind_radius_i,first_sgf_j,&
1824 : !$OMP lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,&
1825 : !$OMP lmax_k,lmin_k,npgfk,nsetk,nsgfk,rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,&
1826 : !$OMP ncoi,ncoj,ncok,sgfi,sgfj,sgfk,dijk_j,dijk_k,ri,rj,rk,max_contraction_i,max_contraction_j,&
1827 : !$OMP tmp_block,max_contraction_k,iset,jset,kset,block_t_i,blk_size,dijk_contr,found,sp,mepos,&
1828 : !$OMP block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,block_t_k,&
1829 : !$OMP bo,der_ext_i,der_ext_j,der_ext_k,ablock,force,scoord,skip,cpp_buffer,ccp_buffer,&
1830 16 : !$OMP block_k_not_zero,der_k_zero,der_j_zero,block_t_j,block_j_not_zero,i)
1831 :
1832 : mepos = 0
1833 : !$ mepos = omp_get_thread_num()
1834 :
1835 : CALL cp_libint_init_3eri1(lib, MAX(maxli, maxlj, maxlk))
1836 : CALL cp_libint_set_contrdepth(lib, 1)
1837 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
1838 :
1839 : !We split the provided bounds among the threads such that each threads works on a different set of atoms
1840 :
1841 : bo = get_limit(natom, nthread, mepos)
1842 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i=bo)
1843 :
1844 : skip = .FALSE.
1845 : IF (bo(1) > bo(2)) skip = .TRUE.
1846 :
1847 : DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
1848 : CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
1849 : iatom=iatom, jatom=jatom, katom=katom, &
1850 : rij=rij, rjk=rjk, rik=rik)
1851 : IF (skip) EXIT
1852 :
1853 : CALL dbt_get_block(t3c_trace, [iatom, jatom, katom], tmp_block, found)
1854 : IF (.NOT. found) CYCLE
1855 :
1856 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
1857 : npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
1858 : sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
1859 :
1860 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
1861 : npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
1862 : sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
1863 :
1864 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
1865 : npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
1866 : sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
1867 :
1868 : djk = NORM2(rjk)
1869 : dij = NORM2(rij)
1870 : dik = NORM2(rik)
1871 :
1872 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
1873 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
1874 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
1875 :
1876 : IF (PRESENT(der_eps)) THEN
1877 : ALLOCATE (max_contraction_i(nseti))
1878 : DO iset = 1, nseti
1879 : sgfi = first_sgf_i(1, iset)
1880 : max_contraction_i(iset) = MAXVAL([(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)])
1881 : END DO
1882 :
1883 : ALLOCATE (max_contraction_j(nsetj))
1884 : DO jset = 1, nsetj
1885 : sgfj = first_sgf_j(1, jset)
1886 : max_contraction_j(jset) = MAXVAL([(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)])
1887 : END DO
1888 :
1889 : ALLOCATE (max_contraction_k(nsetk))
1890 : DO kset = 1, nsetk
1891 : sgfk = first_sgf_k(1, kset)
1892 : max_contraction_k(kset) = MAXVAL([(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)])
1893 : END DO
1894 : END IF
1895 :
1896 : CALL dbt_blk_sizes(t3c_trace, [iatom, jatom, katom], blk_size)
1897 :
1898 : ALLOCATE (block_t_i(blk_size(2), blk_size(3), blk_size(1), 3))
1899 : ALLOCATE (block_t_j(blk_size(2), blk_size(3), blk_size(1), 3))
1900 : ALLOCATE (block_t_k(blk_size(2), blk_size(3), blk_size(1), 3))
1901 :
1902 : ALLOCATE (ablock(blk_size(2), blk_size(3), blk_size(1)))
1903 : DO i = 1, blk_size(1)
1904 : ablock(:, :, i) = tmp_block(i, :, :)
1905 : END DO
1906 : DEALLOCATE (tmp_block)
1907 :
1908 : block_t_i = 0.0_dp
1909 : block_t_j = 0.0_dp
1910 : block_t_k = 0.0_dp
1911 : block_j_not_zero = .FALSE.
1912 : block_k_not_zero = .FALSE.
1913 :
1914 : DO iset = 1, nseti
1915 :
1916 : DO jset = 1, nsetj
1917 :
1918 : IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE
1919 :
1920 : DO kset = 1, nsetk
1921 :
1922 : IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE
1923 : IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE
1924 :
1925 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
1926 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
1927 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
1928 :
1929 : sgfi = first_sgf_i(1, iset)
1930 : sgfj = first_sgf_j(1, jset)
1931 : sgfk = first_sgf_k(1, kset)
1932 :
1933 : IF (ncoj*ncok*ncoi > 0) THEN
1934 : ALLOCATE (dijk_j(ncoj, ncok, ncoi, 3), SOURCE=0.0_dp)
1935 : ALLOCATE (dijk_k(ncoj, ncok, ncoi, 3), SOURCE=0.0_dp)
1936 :
1937 : der_j_zero = .FALSE.
1938 : der_k_zero = .FALSE.
1939 :
1940 : !need positions for libint. Only relative positions are needed => set ri to 0.0
1941 : ri = 0.0_dp
1942 : rj = rij ! ri + rij
1943 : rk = rik ! ri + rik
1944 :
1945 : CALL eri_3center_derivs(dijk_j, dijk_k, &
1946 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
1947 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
1948 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
1949 : djk, dij, dik, lib, potential_parameter, &
1950 : der_abc_1_ext=der_ext_j, der_abc_2_ext=der_ext_k)
1951 :
1952 : IF (PRESENT(der_eps)) THEN
1953 : DO i_xyz = 1, 3
1954 : IF (der_eps > der_ext_j(i_xyz)*(max_contraction_i(iset)* &
1955 : max_contraction_j(jset)* &
1956 : max_contraction_k(kset))) THEN
1957 : der_j_zero(i_xyz) = .TRUE.
1958 : END IF
1959 : END DO
1960 :
1961 : DO i_xyz = 1, 3
1962 : IF (der_eps > der_ext_k(i_xyz)*(max_contraction_i(iset)* &
1963 : max_contraction_j(jset)* &
1964 : max_contraction_k(kset))) THEN
1965 : der_k_zero(i_xyz) = .TRUE.
1966 : END IF
1967 : END DO
1968 : IF (ALL(der_j_zero) .AND. ALL(der_k_zero)) THEN
1969 : DEALLOCATE (dijk_j, dijk_k)
1970 : CYCLE
1971 : END IF
1972 : END IF
1973 :
1974 : ALLOCATE (dijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
1975 :
1976 : block_start_j = sgfj
1977 : block_end_j = sgfj + nsgfj(jset) - 1
1978 : block_start_k = sgfk
1979 : block_end_k = sgfk + nsgfk(kset) - 1
1980 : block_start_i = sgfi
1981 : block_end_i = sgfi + nsgfi(iset) - 1
1982 :
1983 : DO i_xyz = 1, 3
1984 : IF (der_j_zero(i_xyz)) CYCLE
1985 :
1986 : block_j_not_zero(i_xyz) = .TRUE.
1987 : CALL abc_contract_xsmm(dijk_contr, dijk_j(:, :, :, i_xyz), tspj(jset, jkind)%array, &
1988 : spk(kset, kkind)%array, spi(iset, ikind)%array, &
1989 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
1990 : nsgfi(iset), cpp_buffer, ccp_buffer)
1991 :
1992 : block_t_j(block_start_j:block_end_j, &
1993 : block_start_k:block_end_k, &
1994 : block_start_i:block_end_i, i_xyz) = &
1995 : block_t_j(block_start_j:block_end_j, &
1996 : block_start_k:block_end_k, &
1997 : block_start_i:block_end_i, i_xyz) + &
1998 : dijk_contr(:, :, :)
1999 :
2000 : END DO
2001 :
2002 : DO i_xyz = 1, 3
2003 : IF (der_k_zero(i_xyz)) CYCLE
2004 :
2005 : block_k_not_zero(i_xyz) = .TRUE.
2006 : CALL abc_contract_xsmm(dijk_contr, dijk_k(:, :, :, i_xyz), tspj(jset, jkind)%array, &
2007 : spk(kset, kkind)%array, spi(iset, ikind)%array, &
2008 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
2009 : nsgfi(iset), cpp_buffer, ccp_buffer)
2010 :
2011 : block_t_k(block_start_j:block_end_j, &
2012 : block_start_k:block_end_k, &
2013 : block_start_i:block_end_i, i_xyz) = &
2014 : block_t_k(block_start_j:block_end_j, &
2015 : block_start_k:block_end_k, &
2016 : block_start_i:block_end_i, i_xyz) + &
2017 : dijk_contr(:, :, :)
2018 :
2019 : END DO
2020 :
2021 : DEALLOCATE (dijk_j, dijk_k, dijk_contr)
2022 : END IF ! number of triples > 0
2023 : END DO
2024 : END DO
2025 : END DO
2026 :
2027 : !We obtain the derivative wrt to first center using translational invariance
2028 : DO i_xyz = 1, 3
2029 : block_t_i(:, :, :, i_xyz) = -block_t_j(:, :, :, i_xyz) - block_t_k(:, :, :, i_xyz)
2030 : END DO
2031 :
2032 : !virial contribution coming from deriv wrt to first center
2033 : DO i_xyz = 1, 3
2034 : force = pref*SUM(ablock(:, :, :)*block_t_i(:, :, :, i_xyz))
2035 : CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
2036 : DO j_xyz = 1, 3
2037 : !$OMP ATOMIC
2038 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
2039 : END DO
2040 : END DO
2041 :
2042 : !second center
2043 : DO i_xyz = 1, 3
2044 : force = pref*SUM(ablock(:, :, :)*block_t_j(:, :, :, i_xyz))
2045 : CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
2046 : DO j_xyz = 1, 3
2047 : !$OMP ATOMIC
2048 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
2049 : END DO
2050 : END DO
2051 :
2052 : !third center
2053 : DO i_xyz = 1, 3
2054 : force = pref*SUM(ablock(:, :, :)*block_t_k(:, :, :, i_xyz))
2055 : CALL real_to_scaled(scoord, particle_set(iatom)%r + rik, cell)
2056 : DO j_xyz = 1, 3
2057 : !$OMP ATOMIC
2058 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
2059 : END DO
2060 : END DO
2061 :
2062 : DEALLOCATE (block_t_i, block_t_j, block_t_k)
2063 : IF (PRESENT(der_eps)) THEN
2064 : DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k, ablock)
2065 : END IF
2066 : END DO
2067 :
2068 : IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
2069 : IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
2070 :
2071 : CALL cp_libint_cleanup_3eri1(lib)
2072 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
2073 : !$OMP END PARALLEL
2074 :
2075 132 : DO iset = 1, max_nset
2076 364 : DO ibasis = 1, nbasis
2077 232 : IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
2078 232 : IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
2079 348 : IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
2080 : END DO
2081 : END DO
2082 :
2083 16 : DEALLOCATE (spi, tspj, spk)
2084 :
2085 16 : CALL timestop(handle)
2086 128 : END SUBROUTINE calc_3c_virial
2087 :
2088 : ! **************************************************************************************************
2089 : !> \brief Build 3-center integral tensor
2090 : !> \param t3c empty DBCSR tensor
2091 : !> Should be of shape (1,1) if no kpoints are used and of shape (nimages, nimages)
2092 : !> if k-points are used
2093 : !> \param filter_eps Filter threshold for tensor blocks
2094 : !> \param qs_env ...
2095 : !> \param nl_3c 3-center neighborlist
2096 : !> \param basis_i ...
2097 : !> \param basis_j ...
2098 : !> \param basis_k ...
2099 : !> \param potential_parameter ...
2100 : !> \param int_eps neglect integrals smaller than int_eps
2101 : !> \param op_pos operator position.
2102 : !> 1: calculate (i|jk) integrals,
2103 : !> 2: calculate (ij|k) integrals
2104 : !> \param do_kpoints ...
2105 : !> this routine requires that libint has been static initialised somewhere else
2106 : !> \param do_hfx_kpoints ...
2107 : !> \param desymmetrize ...
2108 : !> \param cell_sym ...
2109 : !> \param bounds_i ...
2110 : !> \param bounds_j ...
2111 : !> \param bounds_k ...
2112 : !> \param RI_range ...
2113 : !> \param img_to_RI_cell ...
2114 : !> \param cell_to_index_ext ...
2115 : ! **************************************************************************************************
2116 1978 : SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
2117 1978 : nl_3c, basis_i, basis_j, basis_k, &
2118 : potential_parameter, int_eps, &
2119 : op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, cell_sym, &
2120 : bounds_i, bounds_j, bounds_k, &
2121 1978 : RI_range, img_to_RI_cell, cell_to_index_ext)
2122 :
2123 : TYPE(dbt_type), DIMENSION(:, :), INTENT(INOUT) :: t3c
2124 : REAL(KIND=dp), INTENT(IN) :: filter_eps
2125 : TYPE(qs_environment_type), POINTER :: qs_env
2126 : TYPE(neighbor_list_3c_type), INTENT(INOUT) :: nl_3c
2127 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j, basis_k
2128 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
2129 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: int_eps
2130 : INTEGER, INTENT(IN), OPTIONAL :: op_pos
2131 : LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints, &
2132 : desymmetrize, cell_sym
2133 : INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k
2134 : REAL(dp), INTENT(IN), OPTIONAL :: RI_range
2135 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: img_to_RI_cell
2136 : INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: cell_to_index_ext
2137 :
2138 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integrals'
2139 :
2140 : INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, &
2141 : block_start_k, egfi, handle, handle2, i, iatom, ibasis, ikind, ilist, imax, iset, jatom, &
2142 : jcell, jkind, jset, katom, kcell, kkind, kset, m_max, max_ncoj, max_nset, max_nsgfi, &
2143 : maxli, maxlj, maxlk, mepos, natom, nbasis, ncell_RI, ncoi, ncoj, ncok, nimg, nseti, &
2144 : nsetj, nsetk, nthread, op_ij, op_jk, op_pos_prv, sgfi, sgfj, sgfk, unit_id
2145 1978 : INTEGER, ALLOCATABLE, DIMENSION(:) :: img_to_RI_cell_prv
2146 : INTEGER, DIMENSION(2) :: bo
2147 : INTEGER, DIMENSION(3) :: blk_idx, blk_size, cell_j, cell_k, &
2148 : kp_index_lbounds, kp_index_ubounds, sp
2149 1978 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, &
2150 3956 : lmin_k, npgfi, npgfj, npgfk, nsgfi, &
2151 1978 : nsgfj, nsgfk
2152 1978 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k
2153 1978 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2154 : LOGICAL :: block_not_zero, cell_sym_prv, debug, &
2155 : desymmetrize_prv, do_hfx_kpoints_prv, &
2156 : do_kpoints_prv, found, skip
2157 : REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, &
2158 : kind_radius_i, kind_radius_j, &
2159 : kind_radius_k, max_contraction_i, &
2160 : prefac, sijk_ext
2161 1978 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, &
2162 1978 : max_contraction_j, max_contraction_k
2163 3956 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: block_t, dummy_block_t, sijk, &
2164 3956 : sijk_contr, tmp_ijk
2165 : REAL(KIND=dp), DIMENSION(1, 1, 1) :: counter
2166 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk
2167 1978 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k
2168 1978 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, &
2169 3956 : sphi_k, zeti, zetj, zetk
2170 1978 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2171 : TYPE(cell_type), POINTER :: cell
2172 3956 : TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj
2173 : TYPE(cp_libint_t) :: lib
2174 13846 : TYPE(dbt_type) :: t_3c_tmp
2175 : TYPE(dft_control_type), POINTER :: dft_control
2176 : TYPE(gto_basis_set_type), POINTER :: basis_set
2177 : TYPE(kpoint_type), POINTER :: kpoints
2178 : TYPE(mp_para_env_type), POINTER :: para_env
2179 : TYPE(neighbor_list_3c_iterator_type) :: nl_3c_iter
2180 1978 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2181 :
2182 1978 : CALL timeset(routineN, handle)
2183 :
2184 1978 : debug = .FALSE.
2185 :
2186 1978 : IF (PRESENT(do_kpoints)) THEN
2187 382 : do_kpoints_prv = do_kpoints
2188 : ELSE
2189 : do_kpoints_prv = .FALSE.
2190 : END IF
2191 :
2192 1978 : IF (PRESENT(do_hfx_kpoints)) THEN
2193 140 : do_hfx_kpoints_prv = do_hfx_kpoints
2194 : ELSE
2195 : do_hfx_kpoints_prv = .FALSE.
2196 : END IF
2197 140 : IF (do_hfx_kpoints_prv) THEN
2198 140 : CPASSERT(do_kpoints_prv)
2199 140 : CPASSERT(PRESENT(RI_range))
2200 140 : CPASSERT(PRESENT(img_to_RI_cell))
2201 : END IF
2202 :
2203 1838 : IF (PRESENT(img_to_RI_cell)) THEN
2204 420 : ALLOCATE (img_to_RI_cell_prv(SIZE(img_to_RI_cell)))
2205 6598 : img_to_RI_cell_prv(:) = img_to_RI_cell
2206 : END IF
2207 :
2208 1978 : IF (PRESENT(desymmetrize)) THEN
2209 1978 : desymmetrize_prv = desymmetrize
2210 : ELSE
2211 : desymmetrize_prv = .TRUE.
2212 : END IF
2213 :
2214 1978 : IF (PRESENT(cell_sym)) THEN
2215 6 : cell_sym_prv = cell_sym
2216 : ELSE
2217 1972 : cell_sym_prv = .FALSE.
2218 : END IF
2219 :
2220 1978 : op_ij = do_potential_id; op_jk = do_potential_id
2221 :
2222 1978 : IF (PRESENT(op_pos)) THEN
2223 474 : op_pos_prv = op_pos
2224 : ELSE
2225 1504 : op_pos_prv = 1
2226 : END IF
2227 :
2228 1838 : SELECT CASE (op_pos_prv)
2229 : CASE (1)
2230 1838 : op_ij = potential_parameter%potential_type
2231 : CASE (2)
2232 1978 : op_jk = potential_parameter%potential_type
2233 : END SELECT
2234 :
2235 1978 : dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp
2236 :
2237 : IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short &
2238 1978 : .OR. op_ij == do_potential_mix_cl_trunc) THEN
2239 1438 : dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor
2240 1438 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
2241 540 : ELSEIF (op_ij == do_potential_coulomb) THEN
2242 192 : dr_ij = 1000000.0_dp
2243 192 : dr_ik = 1000000.0_dp
2244 : END IF
2245 :
2246 : IF (op_jk == do_potential_truncated .OR. op_jk == do_potential_short &
2247 1978 : .OR. op_jk == do_potential_mix_cl_trunc) THEN
2248 12 : dr_jk = potential_parameter%cutoff_radius*cutoff_screen_factor
2249 12 : dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor
2250 1966 : ELSEIF (op_jk == do_potential_coulomb) THEN
2251 0 : dr_jk = 1000000.0_dp
2252 0 : dr_ik = 1000000.0_dp
2253 : END IF
2254 :
2255 1978 : NULLIFY (qs_kind_set, atomic_kind_set)
2256 :
2257 : ! get stuff
2258 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
2259 1978 : natom=natom, dft_control=dft_control, para_env=para_env)
2260 :
2261 1978 : IF (do_kpoints_prv) THEN
2262 152 : IF (PRESENT(cell_to_index_ext)) THEN
2263 6 : cell_to_index => cell_to_index_ext
2264 288 : nimg = MAXVAL(cell_to_index)
2265 : ELSE
2266 146 : CALL get_qs_env(qs_env, kpoints=kpoints)
2267 146 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
2268 146 : nimg = dft_control%nimages
2269 : END IF
2270 152 : ncell_RI = nimg
2271 152 : IF (do_hfx_kpoints_prv) THEN
2272 140 : nimg = SIZE(t3c, 1)
2273 140 : ncell_RI = SIZE(t3c, 2)
2274 : END IF
2275 : ELSE
2276 : nimg = 1
2277 : ncell_RI = 1
2278 : END IF
2279 :
2280 : CALL alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potential_parameter, &
2281 : op_pos=op_pos_prv, do_kpoints=do_kpoints, do_hfx_kpoints=do_hfx_kpoints, &
2282 : bounds_i=bounds_i, bounds_j=bounds_j, bounds_k=bounds_k, &
2283 : RI_range=RI_range, img_to_RI_cell=img_to_RI_cell, cell_sym=cell_sym_prv, &
2284 3816 : cell_to_index=cell_to_index)
2285 :
2286 1978 : IF (do_hfx_kpoints_prv) THEN
2287 140 : CPASSERT(op_pos_prv == 2)
2288 140 : CPASSERT(.NOT. desymmetrize_prv)
2289 1838 : ELSE IF (do_kpoints_prv) THEN
2290 36 : CPASSERT(ALL(SHAPE(t3c) == [nimg, ncell_RI]))
2291 : END IF
2292 :
2293 : !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction
2294 1978 : nbasis = SIZE(basis_i)
2295 1978 : max_nsgfi = 0
2296 1978 : max_nset = 0
2297 1978 : maxli = 0
2298 4932 : DO ibasis = 1, nbasis
2299 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, &
2300 2954 : nset=iset, nsgf_set=nsgfi, npgf=npgfi)
2301 2954 : maxli = MAX(maxli, imax)
2302 2954 : max_nset = MAX(max_nset, iset)
2303 17516 : max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi))
2304 : END DO
2305 : max_ncoj = 0
2306 : maxlj = 0
2307 4932 : DO ibasis = 1, nbasis
2308 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, &
2309 2954 : nset=jset, nsgf_set=nsgfj, npgf=npgfj)
2310 2954 : maxlj = MAX(maxlj, imax)
2311 2954 : max_nset = MAX(max_nset, jset)
2312 11314 : max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj))
2313 : END DO
2314 : maxlk = 0
2315 4932 : DO ibasis = 1, nbasis
2316 : CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, &
2317 2954 : nset=kset, nsgf_set=nsgfk, npgf=npgfk)
2318 2954 : maxlk = MAX(maxlk, imax)
2319 4932 : max_nset = MAX(max_nset, kset)
2320 : END DO
2321 1978 : m_max = maxli + maxlj + maxlk
2322 :
2323 : !To minimize expensive memory ops and generally optimize contraction, pre-allocate
2324 : !contiguous sphi arrays (and transposed in the case of sphi_i)
2325 :
2326 1978 : NULLIFY (tspj, spi, spk)
2327 75842 : ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis))
2328 :
2329 4932 : DO ibasis = 1, nbasis
2330 21984 : DO iset = 1, max_nset
2331 17052 : NULLIFY (spi(iset, ibasis)%array)
2332 17052 : NULLIFY (tspj(iset, ibasis)%array)
2333 :
2334 20006 : NULLIFY (spk(iset, ibasis)%array)
2335 : END DO
2336 : END DO
2337 :
2338 7912 : DO ilist = 1, 3
2339 16774 : DO ibasis = 1, nbasis
2340 8862 : IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set
2341 8862 : IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set
2342 8862 : IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set
2343 :
2344 39876 : DO iset = 1, basis_set%nset
2345 :
2346 25080 : ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset))
2347 25080 : sgfi = basis_set%first_sgf(1, iset)
2348 25080 : egfi = sgfi + basis_set%nsgf_set(iset) - 1
2349 :
2350 33942 : IF (ilist == 1) THEN
2351 50336 : ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2352 3349236 : spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2353 :
2354 12496 : ELSE IF (ilist == 2) THEN
2355 25528 : ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi))
2356 358742 : tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi))
2357 :
2358 : ELSE
2359 24456 : ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset)))
2360 1495290 : spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi)
2361 : END IF
2362 :
2363 : END DO !iset
2364 : END DO !ibasis
2365 : END DO !ilist
2366 :
2367 : !Init the truncated Coulomb operator
2368 : IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated &
2369 1978 : .OR. op_ij == do_potential_mix_cl_trunc .OR. op_jk == do_potential_mix_cl_trunc) THEN
2370 :
2371 1444 : IF (m_max > get_lmax_init()) THEN
2372 84 : IF (para_env%mepos == 0) THEN
2373 42 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
2374 : END IF
2375 84 : CALL init(m_max, unit_id, para_env%mepos, para_env)
2376 84 : IF (para_env%mepos == 0) THEN
2377 42 : CALL close_file(unit_id)
2378 : END IF
2379 : END IF
2380 : END IF
2381 :
2382 1978 : CALL init_md_ftable(nmax=m_max)
2383 :
2384 1978 : IF (do_kpoints_prv) THEN
2385 608 : kp_index_lbounds = LBOUND(cell_to_index)
2386 608 : kp_index_ubounds = UBOUND(cell_to_index)
2387 : END IF
2388 :
2389 7912 : counter = 1.0_dp
2390 :
2391 1978 : nthread = 1
2392 1978 : !$ nthread = omp_get_max_threads()
2393 :
2394 : !$OMP PARALLEL DEFAULT(NONE) &
2395 : !$OMP SHARED (nthread,do_kpoints_prv,kp_index_lbounds,kp_index_ubounds,maxli,maxlk,maxlj,bounds_i,&
2396 : !$OMP bounds_j,bounds_k,nimg,basis_i,basis_j,basis_k,dr_ij,dr_jk,dr_ik,ncoset,&
2397 : !$OMP potential_parameter,int_eps,t3c,tspj,spi,spk,debug,cell_to_index,&
2398 : !$OMP natom,nl_3c,cell,op_pos_prv,do_hfx_kpoints_prv,RI_range,ncell_RI, &
2399 : !$OMP img_to_RI_cell_prv, cell_sym_prv) &
2400 : !$OMP PRIVATE (lib,nl_3c_iter,ikind,jkind,kkind,iatom,jatom,katom,rij,rjk,rik,cell_j,cell_k,&
2401 : !$OMP prefac,jcell,kcell,first_sgf_i,lmax_i,lmin_i,npgfi,nseti,nsgfi,rpgf_i,set_radius_i,&
2402 : !$OMP sphi_i,zeti,kind_radius_i,first_sgf_j,lmax_j,lmin_j,npgfj,nsetj,nsgfj,rpgf_j,&
2403 : !$OMP set_radius_j,sphi_j,zetj,kind_radius_j,first_sgf_k,lmax_k,lmin_k,npgfk,nsetk,nsgfk,&
2404 : !$OMP rpgf_k,set_radius_k,sphi_k,zetk,kind_radius_k,djk,dij,dik,ncoi,ncoj,ncok,sgfi,sgfj,&
2405 : !$OMP sgfk,sijk,ri,rj,rk,sijk_ext,block_not_zero,max_contraction_i,max_contraction_j,&
2406 : !$OMP max_contraction_k,iset,jset,kset,block_t,blk_size,sijk_contr,cpp_buffer,ccp_buffer,&
2407 : !$OMP block_start_j,block_end_j,block_start_k,block_end_k,block_start_i,block_end_i,found,&
2408 1978 : !$OMP dummy_block_t,sp,handle2,mepos,bo,skip,tmp_ijk,i,blk_idx)
2409 :
2410 : mepos = 0
2411 : !$ mepos = omp_get_thread_num()
2412 :
2413 : CALL cp_libint_init_3eri(lib, MAX(maxli, maxlj, maxlk))
2414 : CALL cp_libint_set_contrdepth(lib, 1)
2415 : CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
2416 :
2417 : !We split the provided bounds among the threads such that each threads works on a different set of atoms
2418 : IF (PRESENT(bounds_i)) THEN
2419 : bo = get_limit(bounds_i(2) - bounds_i(1) + 1, nthread, mepos)
2420 : bo(:) = bo(:) + bounds_i(1) - 1
2421 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
2422 : ELSE IF (PRESENT(bounds_j)) THEN
2423 :
2424 : bo = get_limit(bounds_j(2) - bounds_j(1) + 1, nthread, mepos)
2425 : bo(:) = bo(:) + bounds_j(1) - 1
2426 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bo, bounds_k)
2427 : ELSE IF (PRESENT(bounds_k)) THEN
2428 : bo = get_limit(bounds_k(2) - bounds_k(1) + 1, nthread, mepos)
2429 : bo(:) = bo(:) + bounds_k(1) - 1
2430 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bo)
2431 : ELSE
2432 : bo = get_limit(natom, nthread, mepos)
2433 : CALL nl_3c_iter_set_bounds(nl_3c_iter, bo, bounds_j, bounds_k)
2434 : END IF
2435 :
2436 : skip = .FALSE.
2437 : IF (bo(1) > bo(2)) skip = .TRUE.
2438 :
2439 : DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
2440 : CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
2441 : iatom=iatom, jatom=jatom, katom=katom, &
2442 : rij=rij, rjk=rjk, rik=rik, cell_j=cell_j, cell_k=cell_k)
2443 : IF (skip) EXIT
2444 :
2445 : djk = NORM2(rjk)
2446 : dij = NORM2(rij)
2447 : dik = NORM2(rik)
2448 :
2449 : IF (nl_3c%sym == symmetric_jk) THEN
2450 : IF (jatom == katom) THEN
2451 : ! factor 0.5 due to double-counting of diagonal blocks
2452 : ! (we desymmetrize by adding transpose)
2453 : prefac = 0.5_dp
2454 : ELSE
2455 : prefac = 1.0_dp
2456 : END IF
2457 : ELSEIF (nl_3c%sym == symmetric_ij) THEN
2458 : IF (iatom == jatom) THEN
2459 : ! factor 0.5 due to double-counting of diagonal blocks
2460 : ! (we desymmetrize by adding transpose)
2461 : prefac = 0.5_dp
2462 : ELSE
2463 : prefac = 1.0_dp
2464 : END IF
2465 : ELSE
2466 : prefac = 1.0_dp
2467 : END IF
2468 : IF (do_kpoints_prv) prefac = 1.0_dp
2469 :
2470 : IF (do_kpoints_prv) THEN
2471 :
2472 : IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
2473 : ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
2474 :
2475 : jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2476 : IF (jcell > nimg .OR. jcell < 1) CYCLE
2477 :
2478 : IF (ANY([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
2479 : ANY([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) CYCLE
2480 :
2481 : kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
2482 : IF (kcell > nimg .OR. kcell < 1) CYCLE
2483 : !In case of HFX k-points, we only consider P^k if d_ik <= RI_range
2484 : IF (do_hfx_kpoints_prv) THEN
2485 : IF (dik > RI_range) CYCLE
2486 : kcell = img_to_RI_cell_prv(kcell)
2487 : END IF
2488 : ELSE
2489 : jcell = 1; kcell = 1
2490 : END IF
2491 :
2492 : IF (cell_sym_prv .AND. jcell < kcell) CYCLE
2493 :
2494 : blk_idx = [iatom, jatom, katom]
2495 : IF (do_hfx_kpoints_prv) THEN
2496 : blk_idx(3) = (kcell - 1)*natom + katom
2497 : kcell = 1
2498 : END IF
2499 :
2500 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
2501 : npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
2502 : sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i)
2503 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
2504 : npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
2505 : sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j)
2506 : IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE
2507 :
2508 : CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, lmax=lmax_k, lmin=lmin_k, &
2509 : npgf=npgfk, nset=nsetk, nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, &
2510 : sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k)
2511 : IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE
2512 : IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE
2513 :
2514 : IF (PRESENT(int_eps)) THEN
2515 : ALLOCATE (max_contraction_j(nsetj))
2516 : DO jset = 1, nsetj
2517 : sgfj = first_sgf_j(1, jset)
2518 : max_contraction_j(jset) = MAXVAL([(SUM(ABS(sphi_j(:, i))), i=sgfj, sgfj + nsgfj(jset) - 1)])
2519 : END DO
2520 :
2521 : ALLOCATE (max_contraction_k(nsetk))
2522 : DO kset = 1, nsetk
2523 : sgfk = first_sgf_k(1, kset)
2524 : max_contraction_k(kset) = MAXVAL([(SUM(ABS(sphi_k(:, i))), i=sgfk, sgfk + nsgfk(kset) - 1)])
2525 : END DO
2526 : END IF
2527 :
2528 : CALL dbt_blk_sizes(t3c(jcell, kcell), blk_idx, blk_size)
2529 :
2530 : ALLOCATE (block_t(blk_size(2), blk_size(3), blk_size(1)), SOURCE=0.0_dp)
2531 :
2532 : block_not_zero = .FALSE.
2533 : DO iset = 1, nseti
2534 :
2535 : sgfi = first_sgf_i(1, iset)
2536 : max_contraction_i = MAXVAL([(SUM(ABS(sphi_i(:, i))), i=sgfi, sgfi + nsgfi(iset) - 1)])
2537 :
2538 : DO jset = 1, nsetj
2539 :
2540 : IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE
2541 :
2542 : DO kset = 1, nsetk
2543 :
2544 : IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE
2545 : IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE
2546 :
2547 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
2548 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
2549 : ncok = npgfk(kset)*ncoset(lmax_k(kset))
2550 :
2551 : !ensure non-zero number of triples below
2552 : IF (ncoj*ncok*ncoi == 0) CYCLE
2553 :
2554 : !need positions for libint. Only relative positions are needed => set ri to 0.0
2555 : ri = 0.0_dp
2556 : rj = rij ! ri + rij
2557 : rk = rik ! ri + rik
2558 :
2559 : ALLOCATE (sijk(ncoj, ncok, ncoi))
2560 : IF (op_pos_prv == 1) THEN
2561 : sijk(:, :, :) = 0.0_dp
2562 : CALL eri_3center(sijk, &
2563 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
2564 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
2565 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
2566 : djk, dij, dik, lib, potential_parameter, int_abc_ext=sijk_ext)
2567 : ELSE
2568 : ALLOCATE (tmp_ijk(ncoi, ncoj, ncok), SOURCE=0.0_dp)
2569 : CALL eri_3center(tmp_ijk, &
2570 : lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), rpgf_i(:, iset), ri, &
2571 : lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), rpgf_j(:, jset), rj, &
2572 : lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), rpgf_k(:, kset), rk, &
2573 : dij, dik, djk, lib, potential_parameter, int_abc_ext=sijk_ext)
2574 :
2575 : !F08: sijk = RESHAPE(tmp_ijk, [ncoj, ncok, ncoi], ORDER=[2, 3, 1]) with sijk not allocated
2576 : DO i = 1, ncoi !TODO: revise/check for efficiency
2577 : sijk(:, :, i) = tmp_ijk(i, :, :)
2578 : END DO
2579 : DEALLOCATE (tmp_ijk)
2580 : END IF
2581 :
2582 : IF (PRESENT(int_eps)) THEN
2583 : IF (int_eps > sijk_ext*(max_contraction_i* &
2584 : max_contraction_j(jset)* &
2585 : max_contraction_k(kset))) THEN
2586 : DEALLOCATE (sijk)
2587 : CYCLE
2588 : END IF
2589 : END IF
2590 :
2591 : block_not_zero = .TRUE.
2592 : ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset)))
2593 : CALL abc_contract_xsmm(sijk_contr, sijk, tspj(jset, jkind)%array, &
2594 : spk(kset, kkind)%array, spi(iset, ikind)%array, &
2595 : ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), &
2596 : nsgfi(iset), cpp_buffer, ccp_buffer, prefac)
2597 : DEALLOCATE (sijk)
2598 :
2599 : sgfj = first_sgf_j(1, jset)
2600 : sgfk = first_sgf_k(1, kset)
2601 :
2602 : block_start_j = sgfj
2603 : block_end_j = sgfj + nsgfj(jset) - 1
2604 : block_start_k = sgfk
2605 : block_end_k = sgfk + nsgfk(kset) - 1
2606 : block_start_i = sgfi
2607 : block_end_i = sgfi + nsgfi(iset) - 1
2608 :
2609 : block_t(block_start_j:block_end_j, &
2610 : block_start_k:block_end_k, &
2611 : block_start_i:block_end_i) = &
2612 : block_t(block_start_j:block_end_j, &
2613 : block_start_k:block_end_k, &
2614 : block_start_i:block_end_i) + &
2615 : sijk_contr(:, :, :)
2616 : DEALLOCATE (sijk_contr)
2617 :
2618 : END DO
2619 :
2620 : END DO
2621 :
2622 : END DO
2623 :
2624 : IF (block_not_zero) THEN
2625 : !$OMP CRITICAL(qs_tensors_t3c)
2626 : CALL timeset(routineN//"_put_dbcsr", handle2)
2627 : IF (debug) THEN
2628 : CALL dbt_get_block(t3c(jcell, kcell), blk_idx, dummy_block_t, found=found)
2629 : CPASSERT(found)
2630 : END IF
2631 :
2632 : sp = SHAPE(block_t)
2633 : sp([2, 3, 1]) = sp ! sp = sp([2, 3, 1]) performs worse
2634 :
2635 : CALL dbt_put_block(t3c(jcell, kcell), blk_idx, sp, &
2636 : RESHAPE(block_t, SHAPE=sp, ORDER=[2, 3, 1]), summation=.TRUE.)
2637 :
2638 : CALL timestop(handle2)
2639 : !$OMP END CRITICAL(qs_tensors_t3c)
2640 : END IF
2641 :
2642 : DEALLOCATE (block_t)
2643 : IF (PRESENT(int_eps)) THEN
2644 : DEALLOCATE (max_contraction_j, max_contraction_k)
2645 : END IF
2646 : END DO
2647 :
2648 : IF (ALLOCATED(ccp_buffer)) DEALLOCATE (ccp_buffer)
2649 : IF (ALLOCATED(cpp_buffer)) DEALLOCATE (cpp_buffer)
2650 :
2651 : CALL cp_libint_cleanup_3eri(lib)
2652 : CALL neighbor_list_3c_iterator_destroy(nl_3c_iter)
2653 : !$OMP END PARALLEL
2654 :
2655 : !TODO: deal with hfx_kpoints, because should not filter by 1/2
2656 1978 : IF (nl_3c%sym == symmetric_jk .OR. do_kpoints_prv) THEN
2657 :
2658 716 : IF (.NOT. do_hfx_kpoints_prv) THEN
2659 1228 : DO kcell = 1, nimg
2660 2568 : DO jcell = 1, nimg
2661 : ! need half of filter eps because afterwards we add transposed tensor
2662 1992 : CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
2663 : END DO
2664 : END DO
2665 : ELSE
2666 280 : DO kcell = 1, ncell_RI
2667 4058 : DO jcell = 1, nimg
2668 3918 : CALL dbt_filter(t3c(jcell, kcell), filter_eps)
2669 : END DO
2670 : END DO
2671 : END IF
2672 :
2673 716 : IF (desymmetrize_prv) THEN
2674 : ! add transposed of overlap integrals
2675 0 : CALL dbt_create(t3c(1, 1), t_3c_tmp)
2676 0 : DO kcell = 1, jcell
2677 0 : DO jcell = 1, nimg
2678 0 : CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
2679 0 : CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.)
2680 0 : CALL dbt_filter(t3c(kcell, jcell), filter_eps)
2681 : END DO
2682 : END DO
2683 0 : DO kcell = jcell + 1, nimg
2684 0 : DO jcell = 1, nimg
2685 0 : CALL dbt_copy(t3c(jcell, kcell), t_3c_tmp)
2686 0 : CALL dbt_copy(t_3c_tmp, t3c(kcell, jcell), order=[1, 3, 2], summation=.FALSE., move_data=.TRUE.)
2687 0 : CALL dbt_filter(t3c(kcell, jcell), filter_eps)
2688 : END DO
2689 : END DO
2690 0 : CALL dbt_destroy(t_3c_tmp)
2691 : END IF
2692 1262 : ELSEIF (nl_3c%sym == symmetric_ij) THEN
2693 0 : DO kcell = 1, nimg
2694 0 : DO jcell = 1, nimg
2695 0 : CALL dbt_filter(t3c(jcell, kcell), filter_eps/2)
2696 : END DO
2697 : END DO
2698 1262 : ELSEIF (nl_3c%sym == symmetric_none) THEN
2699 2524 : DO kcell = 1, nimg
2700 3786 : DO jcell = 1, nimg
2701 2524 : CALL dbt_filter(t3c(jcell, kcell), filter_eps)
2702 : END DO
2703 : END DO
2704 : ELSE
2705 0 : CPABORT("requested symmetric case not implemented")
2706 : END IF
2707 :
2708 4932 : DO ibasis = 1, nbasis
2709 21984 : DO iset = 1, max_nset
2710 17052 : IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array)
2711 17052 : IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array)
2712 20006 : IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array)
2713 : END DO
2714 : END DO
2715 1978 : DEALLOCATE (spi, tspj, spk)
2716 :
2717 1978 : CALL timestop(handle)
2718 17802 : END SUBROUTINE build_3c_integrals
2719 :
2720 : ! **************************************************************************************************
2721 : !> \brief Calculates the derivatives of 2-center integrals, wrt to the first center
2722 : !> \param t2c_der ...
2723 : !> this routine requires that libint has been static initialised somewhere else
2724 : !> \param filter_eps Filter threshold for matrix blocks
2725 : !> \param qs_env ...
2726 : !> \param nl_2c 2-center neighborlist
2727 : !> \param basis_i ...
2728 : !> \param basis_j ...
2729 : !> \param potential_parameter ...
2730 : !> \param do_kpoints ...
2731 : ! **************************************************************************************************
2732 384 : SUBROUTINE build_2c_derivatives(t2c_der, filter_eps, qs_env, &
2733 384 : nl_2c, basis_i, basis_j, &
2734 : potential_parameter, do_kpoints)
2735 :
2736 : TYPE(dbcsr_type), DIMENSION(:, :), INTENT(INOUT) :: t2c_der
2737 : REAL(KIND=dp), INTENT(IN) :: filter_eps
2738 : TYPE(qs_environment_type), POINTER :: qs_env
2739 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2740 : POINTER :: nl_2c
2741 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
2742 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
2743 : LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints
2744 :
2745 : CHARACTER(len=*), PARAMETER :: routineN = 'build_2c_derivatives'
2746 :
2747 : INTEGER :: handle, i_xyz, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, &
2748 : jset, m_max, maxli, maxlj, natom, ncoi, ncoj, nimg, nseti, nsetj, op_prv, sgfi, sgfj, &
2749 : unit_id
2750 : INTEGER, DIMENSION(3) :: cell_j, kp_index_lbounds, &
2751 : kp_index_ubounds
2752 384 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
2753 384 : npgfj, nsgfi, nsgfj
2754 384 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
2755 384 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
2756 : LOGICAL :: do_kpoints_prv, do_symmetric, found, &
2757 : trans
2758 : REAL(KIND=dp) :: dab
2759 384 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij_contr, dij_rs
2760 384 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dij
2761 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj
2762 384 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
2763 384 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
2764 384 : zetj
2765 384 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2766 1536 : TYPE(block_p_type), DIMENSION(3) :: block_t
2767 : TYPE(cp_libint_t) :: lib
2768 : TYPE(dft_control_type), POINTER :: dft_control
2769 : TYPE(kpoint_type), POINTER :: kpoints
2770 : TYPE(mp_para_env_type), POINTER :: para_env
2771 : TYPE(neighbor_list_iterator_p_type), &
2772 384 : DIMENSION(:), POINTER :: nl_iterator
2773 384 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2774 :
2775 384 : CALL timeset(routineN, handle)
2776 :
2777 384 : IF (PRESENT(do_kpoints)) THEN
2778 92 : do_kpoints_prv = do_kpoints
2779 : ELSE
2780 : do_kpoints_prv = .FALSE.
2781 : END IF
2782 :
2783 384 : op_prv = potential_parameter%potential_type
2784 :
2785 384 : NULLIFY (qs_kind_set, atomic_kind_set, block_t(1)%block, block_t(2)%block, block_t(3)%block, cell_to_index)
2786 :
2787 : ! get stuff
2788 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
2789 384 : natom=natom, kpoints=kpoints, dft_control=dft_control, para_env=para_env)
2790 :
2791 384 : IF (do_kpoints_prv) THEN
2792 92 : nimg = SIZE(t2c_der, 1)
2793 92 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
2794 368 : kp_index_lbounds = LBOUND(cell_to_index)
2795 368 : kp_index_ubounds = UBOUND(cell_to_index)
2796 : ELSE
2797 : nimg = 1
2798 : END IF
2799 :
2800 : ! check for symmetry
2801 384 : CPASSERT(SIZE(nl_2c) > 0)
2802 384 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
2803 :
2804 384 : IF (do_symmetric) THEN
2805 584 : DO img = 1, nimg
2806 : !Derivtive matrix is assymetric
2807 1460 : DO i_xyz = 1, 3
2808 1168 : CPASSERT(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_antisymmetric)
2809 : END DO
2810 : END DO
2811 : ELSE
2812 2296 : DO img = 1, nimg
2813 8908 : DO i_xyz = 1, 3
2814 8816 : CPASSERT(dbcsr_get_matrix_type(t2c_der(img, i_xyz)) == dbcsr_type_no_symmetry)
2815 : END DO
2816 : END DO
2817 : END IF
2818 :
2819 2880 : DO img = 1, nimg
2820 10368 : DO i_xyz = 1, 3
2821 9984 : CALL cp_dbcsr_alloc_block_from_nbl(t2c_der(img, i_xyz), nl_2c)
2822 : END DO
2823 : END DO
2824 :
2825 384 : maxli = 0
2826 1086 : DO ibasis = 1, SIZE(basis_i)
2827 702 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
2828 1086 : maxli = MAX(maxli, imax)
2829 : END DO
2830 384 : maxlj = 0
2831 1086 : DO ibasis = 1, SIZE(basis_j)
2832 702 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
2833 1086 : maxlj = MAX(maxlj, imax)
2834 : END DO
2835 :
2836 384 : m_max = maxli + maxlj + 1
2837 :
2838 : !Init the truncated Coulomb operator
2839 384 : IF (op_prv == do_potential_truncated .OR. op_prv == do_potential_mix_cl_trunc) THEN
2840 :
2841 60 : IF (m_max > get_lmax_init()) THEN
2842 22 : IF (para_env%mepos == 0) THEN
2843 11 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
2844 : END IF
2845 22 : CALL init(m_max, unit_id, para_env%mepos, para_env)
2846 22 : IF (para_env%mepos == 0) THEN
2847 11 : CALL close_file(unit_id)
2848 : END IF
2849 : END IF
2850 : END IF
2851 :
2852 384 : CALL init_md_ftable(nmax=m_max)
2853 :
2854 384 : CALL cp_libint_init_2eri1(lib, MAX(maxli, maxlj))
2855 384 : CALL cp_libint_set_contrdepth(lib, 1)
2856 :
2857 384 : CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
2858 14908 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
2859 :
2860 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
2861 14524 : iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
2862 14524 : IF (do_kpoints_prv) THEN
2863 11326 : IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
2864 : ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
2865 1618 : img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2866 1618 : IF (img > nimg .OR. img < 1) CYCLE
2867 : ELSE
2868 : img = 1
2869 : END IF
2870 :
2871 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
2872 : npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
2873 14496 : sphi=sphi_i, zet=zeti)
2874 :
2875 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
2876 : npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
2877 14496 : sphi=sphi_j, zet=zetj)
2878 :
2879 14496 : IF (do_symmetric) THEN
2880 12906 : IF (iatom <= jatom) THEN
2881 7900 : irow = iatom
2882 7900 : icol = jatom
2883 : ELSE
2884 5006 : irow = jatom
2885 5006 : icol = iatom
2886 : END IF
2887 : ELSE
2888 1590 : irow = iatom
2889 1590 : icol = jatom
2890 : END IF
2891 :
2892 57984 : dab = NORM2(rij)
2893 14496 : trans = do_symmetric .AND. (iatom > jatom)
2894 :
2895 57984 : DO i_xyz = 1, 3
2896 : CALL dbcsr_get_block_p(matrix=t2c_der(img, i_xyz), &
2897 43488 : row=irow, col=icol, BLOCK=block_t(i_xyz)%block, found=found)
2898 57984 : CPASSERT(found)
2899 : END DO
2900 :
2901 49556 : DO iset = 1, nseti
2902 :
2903 34676 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
2904 34676 : sgfi = first_sgf_i(1, iset)
2905 :
2906 172652 : DO jset = 1, nsetj
2907 :
2908 123452 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
2909 123452 : sgfj = first_sgf_j(1, jset)
2910 :
2911 158128 : IF (ncoi*ncoj > 0) THEN
2912 493808 : ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
2913 313937563 : ALLOCATE (dij(ncoi, ncoj, 3), SOURCE=0.0_dp)
2914 :
2915 123452 : ri = 0.0_dp
2916 123452 : rj = rij
2917 :
2918 : CALL eri_2center_derivs(dij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
2919 : rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
2920 123452 : rpgf_j(:, jset), rj, dab, lib, potential_parameter)
2921 :
2922 493808 : DO i_xyz = 1, 3
2923 :
2924 123412836 : dij_contr(:, :) = 0.0_dp
2925 : CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
2926 : sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
2927 370356 : ncoi, ncoj, nsgfi(iset), nsgfj(jset))
2928 :
2929 370356 : IF (trans) THEN
2930 : !if transpose, then -1 factor for antisymmetry
2931 449508 : ALLOCATE (dij_rs(nsgfj(jset), nsgfi(iset)))
2932 44948634 : dij_rs(:, :) = -1.0_dp*TRANSPOSE(dij_contr)
2933 : ELSE
2934 1031916 : ALLOCATE (dij_rs(nsgfi(iset), nsgfj(jset)))
2935 78355215 : dij_rs(:, :) = dij_contr
2936 : END IF
2937 :
2938 : CALL block_add("IN", dij_rs, &
2939 : nsgfi(iset), nsgfj(jset), block_t(i_xyz)%block, &
2940 370356 : sgfi, sgfj, trans=trans)
2941 493808 : DEALLOCATE (dij_rs)
2942 : END DO
2943 :
2944 123452 : DEALLOCATE (dij, dij_contr)
2945 : END IF
2946 : END DO
2947 : END DO
2948 : END DO
2949 :
2950 384 : CALL cp_libint_cleanup_2eri1(lib)
2951 384 : CALL neighbor_list_iterator_release(nl_iterator)
2952 :
2953 2880 : DO img = 1, nimg
2954 10368 : DO i_xyz = 1, 3
2955 7488 : CALL dbcsr_finalize(t2c_der(img, i_xyz))
2956 9984 : CALL dbcsr_filter(t2c_der(img, i_xyz), filter_eps)
2957 : END DO
2958 : END DO
2959 :
2960 384 : CALL timestop(handle)
2961 :
2962 1152 : END SUBROUTINE build_2c_derivatives
2963 :
2964 : ! **************************************************************************************************
2965 : !> \brief Calculates the virial coming from 2c derivatives on the fly
2966 : !> \param work_virial ...
2967 : !> \param t2c_trace the 2c tensor that we should trace with the derivatives
2968 : !> \param pref ...
2969 : !> \param qs_env ...
2970 : !> \param nl_2c 2-center neighborlist. Assumed to have compatible distribution with t2c_trace,
2971 : !> and to be non-symmetric
2972 : !> \param basis_i ...
2973 : !> \param basis_j ...
2974 : !> \param potential_parameter ...
2975 : ! **************************************************************************************************
2976 12 : SUBROUTINE calc_2c_virial(work_virial, t2c_trace, pref, qs_env, nl_2c, basis_i, basis_j, potential_parameter)
2977 : REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: work_virial
2978 : TYPE(dbcsr_type), INTENT(INOUT) :: t2c_trace
2979 : REAL(KIND=dp), INTENT(IN) :: pref
2980 : TYPE(qs_environment_type), POINTER :: qs_env
2981 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2982 : POINTER :: nl_2c
2983 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
2984 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
2985 :
2986 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_2c_virial'
2987 :
2988 : INTEGER :: handle, i_xyz, iatom, ibasis, ikind, imax, iset, j_xyz, jatom, jkind, jset, &
2989 : m_max, maxli, maxlj, natom, ncoi, ncoj, nseti, nsetj, op_prv, sgfi, sgfj, unit_id
2990 12 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
2991 12 : npgfj, nsgfi, nsgfj
2992 12 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
2993 : LOGICAL :: do_symmetric, found
2994 : REAL(dp) :: force
2995 12 : REAL(dp), DIMENSION(:, :), POINTER :: pblock
2996 : REAL(KIND=dp) :: dab
2997 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dij_contr
2998 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: dij
2999 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj, scoord
3000 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3001 12 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3002 12 : zetj
3003 12 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3004 : TYPE(cell_type), POINTER :: cell
3005 : TYPE(cp_libint_t) :: lib
3006 : TYPE(dft_control_type), POINTER :: dft_control
3007 : TYPE(mp_para_env_type), POINTER :: para_env
3008 : TYPE(neighbor_list_iterator_p_type), &
3009 12 : DIMENSION(:), POINTER :: nl_iterator
3010 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3011 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3012 :
3013 12 : CALL timeset(routineN, handle)
3014 :
3015 12 : op_prv = potential_parameter%potential_type
3016 :
3017 12 : NULLIFY (qs_kind_set, atomic_kind_set, pblock, particle_set, cell)
3018 :
3019 : ! get stuff
3020 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
3021 : natom=natom, dft_control=dft_control, para_env=para_env, &
3022 12 : particle_set=particle_set, cell=cell)
3023 :
3024 : ! check for symmetry
3025 12 : CPASSERT(SIZE(nl_2c) > 0)
3026 12 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3027 12 : CPASSERT(.NOT. do_symmetric)
3028 :
3029 12 : maxli = 0
3030 36 : DO ibasis = 1, SIZE(basis_i)
3031 24 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3032 36 : maxli = MAX(maxli, imax)
3033 : END DO
3034 12 : maxlj = 0
3035 36 : DO ibasis = 1, SIZE(basis_j)
3036 24 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3037 36 : maxlj = MAX(maxlj, imax)
3038 : END DO
3039 :
3040 12 : m_max = maxli + maxlj + 1
3041 :
3042 : !Init the truncated Coulomb operator
3043 12 : IF (op_prv == do_potential_truncated .OR. op_prv == do_potential_mix_cl_trunc) THEN
3044 :
3045 2 : IF (m_max > get_lmax_init()) THEN
3046 0 : IF (para_env%mepos == 0) THEN
3047 0 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3048 : END IF
3049 0 : CALL init(m_max, unit_id, para_env%mepos, para_env)
3050 0 : IF (para_env%mepos == 0) THEN
3051 0 : CALL close_file(unit_id)
3052 : END IF
3053 : END IF
3054 : END IF
3055 :
3056 12 : CALL init_md_ftable(nmax=m_max)
3057 :
3058 12 : CALL cp_libint_init_2eri1(lib, MAX(maxli, maxlj))
3059 12 : CALL cp_libint_set_contrdepth(lib, 1)
3060 :
3061 12 : CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3062 1644 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3063 :
3064 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3065 1632 : iatom=iatom, jatom=jatom, r=rij)
3066 :
3067 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
3068 : npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
3069 1632 : sphi=sphi_i, zet=zeti)
3070 :
3071 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
3072 : npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
3073 1632 : sphi=sphi_j, zet=zetj)
3074 :
3075 6528 : dab = NORM2(rij)
3076 :
3077 1632 : CALL dbcsr_get_block_p(t2c_trace, iatom, jatom, pblock, found)
3078 1632 : IF (.NOT. found) CYCLE
3079 :
3080 6113 : DO iset = 1, nseti
3081 :
3082 4469 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3083 4469 : sgfi = first_sgf_i(1, iset)
3084 :
3085 25808 : DO jset = 1, nsetj
3086 :
3087 19707 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3088 19707 : sgfj = first_sgf_j(1, jset)
3089 :
3090 24176 : IF (ncoi*ncoj > 0) THEN
3091 78828 : ALLOCATE (dij_contr(nsgfi(iset), nsgfj(jset)))
3092 5728014 : ALLOCATE (dij(ncoi, ncoj, 3), SOURCE=0.0_dp)
3093 :
3094 19707 : ri = 0.0_dp
3095 19707 : rj = rij
3096 :
3097 : CALL eri_2center_derivs(dij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
3098 : rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
3099 19707 : rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3100 :
3101 78828 : DO i_xyz = 1, 3
3102 :
3103 2522610 : dij_contr(:, :) = 0.0_dp
3104 : CALL ab_contract(dij_contr, dij(:, :, i_xyz), &
3105 : sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3106 59121 : ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3107 :
3108 2522610 : force = SUM(pblock(sgfi:sgfi + nsgfi(iset) - 1, sgfj:sgfj + nsgfj(jset) - 1)*dij_contr(:, :))
3109 59121 : force = pref*force
3110 :
3111 : !iatom virial
3112 59121 : CALL real_to_scaled(scoord, particle_set(iatom)%r, cell)
3113 236484 : DO j_xyz = 1, 3
3114 236484 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + force*scoord(j_xyz)
3115 : END DO
3116 :
3117 : !jatom virial
3118 236484 : CALL real_to_scaled(scoord, particle_set(iatom)%r + rij, cell)
3119 256191 : DO j_xyz = 1, 3
3120 236484 : work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - force*scoord(j_xyz)
3121 : END DO
3122 : END DO
3123 :
3124 19707 : DEALLOCATE (dij, dij_contr)
3125 : END IF
3126 : END DO
3127 : END DO
3128 : END DO
3129 :
3130 12 : CALL neighbor_list_iterator_release(nl_iterator)
3131 12 : CALL cp_libint_cleanup_2eri1(lib)
3132 :
3133 12 : CALL timestop(handle)
3134 :
3135 24 : END SUBROUTINE calc_2c_virial
3136 :
3137 : ! **************************************************************************************************
3138 : !> \brief ...
3139 : !> \param t2c empty DBCSR matrix
3140 : !> \param filter_eps Filter threshold for matrix blocks
3141 : !> \param qs_env ...
3142 : !> \param nl_2c 2-center neighborlist
3143 : !> \param basis_i ...
3144 : !> \param basis_j ...
3145 : !> \param potential_parameter ...
3146 : !> \param do_kpoints ...
3147 : !> \param do_hfx_kpoints ...
3148 : !> this routine requires that libint has been static initialised somewhere else
3149 : !> \param ext_kpoints ...
3150 : !> \param regularization_RI ...
3151 : ! **************************************************************************************************
3152 882 : SUBROUTINE build_2c_integrals(t2c, filter_eps, qs_env, &
3153 882 : nl_2c, basis_i, basis_j, &
3154 : potential_parameter, do_kpoints, &
3155 : do_hfx_kpoints, ext_kpoints, regularization_RI)
3156 :
3157 : TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: t2c
3158 : REAL(KIND=dp), INTENT(IN) :: filter_eps
3159 : TYPE(qs_environment_type), POINTER :: qs_env
3160 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3161 : POINTER :: nl_2c
3162 : TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_i, basis_j
3163 : TYPE(libint_potential_type), INTENT(IN) :: potential_parameter
3164 : LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, do_hfx_kpoints
3165 : TYPE(kpoint_type), OPTIONAL, POINTER :: ext_kpoints
3166 : REAL(KIND=dp), OPTIONAL :: regularization_RI
3167 :
3168 : CHARACTER(len=*), PARAMETER :: routineN = 'build_2c_integrals'
3169 :
3170 : INTEGER :: handle, i_diag, iatom, ibasis, icol, ikind, imax, img, irow, iset, jatom, jkind, &
3171 : jset, m_max, maxli, maxlj, natom, ncoi, ncoj, nimg, nseti, nsetj, op_prv, sgfi, sgfj, &
3172 : unit_id
3173 : INTEGER, DIMENSION(3) :: cell_j, kp_index_lbounds, &
3174 : kp_index_ubounds
3175 882 : INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmin_i, lmin_j, npgfi, &
3176 882 : npgfj, nsgfi, nsgfj
3177 882 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j
3178 882 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3179 : LOGICAL :: do_hfx_kpoints_prv, do_kpoints_prv, &
3180 : do_symmetric, found, trans
3181 : REAL(KIND=dp) :: dab, min_zet
3182 882 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sij, sij_contr, sij_rs
3183 : REAL(KIND=dp), DIMENSION(3) :: ri, rij, rj
3184 882 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j
3185 882 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, sphi_i, sphi_j, zeti, &
3186 882 : zetj
3187 882 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3188 : TYPE(block_p_type) :: block_t
3189 : TYPE(cell_type), POINTER :: cell
3190 : TYPE(cp_libint_t) :: lib
3191 : TYPE(dft_control_type), POINTER :: dft_control
3192 : TYPE(kpoint_type), POINTER :: kpoints
3193 : TYPE(mp_para_env_type), POINTER :: para_env
3194 : TYPE(neighbor_list_iterator_p_type), &
3195 882 : DIMENSION(:), POINTER :: nl_iterator
3196 882 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3197 :
3198 882 : CALL timeset(routineN, handle)
3199 :
3200 882 : IF (PRESENT(do_kpoints)) THEN
3201 832 : do_kpoints_prv = do_kpoints
3202 : ELSE
3203 : do_kpoints_prv = .FALSE.
3204 : END IF
3205 :
3206 882 : IF (PRESENT(do_hfx_kpoints)) THEN
3207 302 : do_hfx_kpoints_prv = do_hfx_kpoints
3208 : ELSE
3209 : do_hfx_kpoints_prv = .FALSE.
3210 : END IF
3211 302 : IF (do_hfx_kpoints_prv) THEN
3212 160 : CPASSERT(do_kpoints_prv)
3213 : END IF
3214 :
3215 882 : op_prv = potential_parameter%potential_type
3216 :
3217 882 : NULLIFY (qs_kind_set, atomic_kind_set, block_t%block, cell_to_index)
3218 :
3219 : ! get stuff
3220 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
3221 882 : natom=natom, dft_control=dft_control, para_env=para_env, kpoints=kpoints)
3222 :
3223 882 : IF (PRESENT(ext_kpoints)) kpoints => ext_kpoints
3224 :
3225 882 : IF (do_kpoints_prv) THEN
3226 464 : nimg = SIZE(t2c)
3227 464 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
3228 1856 : kp_index_lbounds = LBOUND(cell_to_index)
3229 1856 : kp_index_ubounds = UBOUND(cell_to_index)
3230 : ELSE
3231 : nimg = 1
3232 : END IF
3233 :
3234 : ! check for symmetry
3235 882 : CPASSERT(SIZE(nl_2c) > 0)
3236 882 : CALL get_neighbor_list_set_p(neighbor_list_sets=nl_2c, symmetric=do_symmetric)
3237 :
3238 882 : IF (do_symmetric) THEN
3239 2154 : DO img = 1, nimg
3240 2154 : CPASSERT(dbcsr_has_symmetry(t2c(img)))
3241 : END DO
3242 : ELSE
3243 4884 : DO img = 1, nimg
3244 4884 : CPASSERT(.NOT. dbcsr_has_symmetry(t2c(img)))
3245 : END DO
3246 : END IF
3247 :
3248 7038 : DO img = 1, nimg
3249 7038 : CALL cp_dbcsr_alloc_block_from_nbl(t2c(img), nl_2c)
3250 : END DO
3251 :
3252 882 : maxli = 0
3253 2470 : DO ibasis = 1, SIZE(basis_i)
3254 1588 : CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax)
3255 2470 : maxli = MAX(maxli, imax)
3256 : END DO
3257 882 : maxlj = 0
3258 2470 : DO ibasis = 1, SIZE(basis_j)
3259 1588 : CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax)
3260 2470 : maxlj = MAX(maxlj, imax)
3261 : END DO
3262 :
3263 882 : m_max = maxli + maxlj
3264 :
3265 : !Init the truncated Coulomb operator
3266 882 : IF (op_prv == do_potential_truncated .OR. op_prv == do_potential_mix_cl_trunc) THEN
3267 :
3268 544 : IF (m_max > get_lmax_init()) THEN
3269 106 : IF (para_env%mepos == 0) THEN
3270 53 : CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename)
3271 : END IF
3272 106 : CALL init(m_max, unit_id, para_env%mepos, para_env)
3273 106 : IF (para_env%mepos == 0) THEN
3274 53 : CALL close_file(unit_id)
3275 : END IF
3276 : END IF
3277 : END IF
3278 :
3279 882 : CALL init_md_ftable(nmax=m_max)
3280 :
3281 882 : CALL cp_libint_init_2eri(lib, MAX(maxli, maxlj))
3282 882 : CALL cp_libint_set_contrdepth(lib, 1)
3283 :
3284 882 : CALL neighbor_list_iterator_create(nl_iterator, nl_2c)
3285 29975 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3286 :
3287 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3288 29093 : iatom=iatom, jatom=jatom, r=rij, cell=cell_j)
3289 29093 : IF (do_kpoints_prv) THEN
3290 37009 : IF (ANY([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
3291 : ANY([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) CYCLE
3292 5247 : img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3293 5247 : IF (img > nimg .OR. img < 1) CYCLE
3294 : ELSE
3295 : img = 1
3296 : END IF
3297 :
3298 : CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, lmax=lmax_i, lmin=lmin_i, &
3299 : npgf=npgfi, nset=nseti, nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, &
3300 28985 : sphi=sphi_i, zet=zeti)
3301 :
3302 : CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, lmax=lmax_j, lmin=lmin_j, &
3303 : npgf=npgfj, nset=nsetj, nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, &
3304 28985 : sphi=sphi_j, zet=zetj)
3305 :
3306 28985 : IF (do_symmetric) THEN
3307 25105 : IF (iatom <= jatom) THEN
3308 14962 : irow = iatom
3309 14962 : icol = jatom
3310 : ELSE
3311 10143 : irow = jatom
3312 10143 : icol = iatom
3313 : END IF
3314 : ELSE
3315 3880 : irow = iatom
3316 3880 : icol = jatom
3317 : END IF
3318 :
3319 115940 : dab = NORM2(rij)
3320 :
3321 : CALL dbcsr_get_block_p(matrix=t2c(img), &
3322 28985 : row=irow, col=icol, BLOCK=block_t%block, found=found)
3323 28985 : CPASSERT(found)
3324 28985 : trans = do_symmetric .AND. (iatom > jatom)
3325 :
3326 110025 : DO iset = 1, nseti
3327 :
3328 80158 : ncoi = npgfi(iset)*ncoset(lmax_i(iset))
3329 80158 : sgfi = first_sgf_i(1, iset)
3330 :
3331 507394 : DO jset = 1, nsetj
3332 :
3333 398195 : ncoj = npgfj(jset)*ncoset(lmax_j(jset))
3334 398195 : sgfj = first_sgf_j(1, jset)
3335 :
3336 478353 : IF (ncoi*ncoj > 0) THEN
3337 80841058 : ALLOCATE (sij_contr(nsgfi(iset), nsgfj(jset)), SOURCE=0.0_dp)
3338 207189999 : ALLOCATE (sij(ncoi, ncoj), SOURCE=0.0_dp)
3339 :
3340 398195 : ri = 0.0_dp
3341 398195 : rj = rij
3342 :
3343 : CALL eri_2center(sij, lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), &
3344 : rpgf_i(:, iset), ri, lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), &
3345 398195 : rpgf_j(:, jset), rj, dab, lib, potential_parameter)
3346 :
3347 : CALL ab_contract(sij_contr, sij, &
3348 : sphi_i(:, sgfi:), sphi_j(:, sgfj:), &
3349 398195 : ncoi, ncoj, nsgfi(iset), nsgfj(jset))
3350 :
3351 398195 : DEALLOCATE (sij)
3352 398195 : IF (trans) THEN
3353 492840 : ALLOCATE (sij_rs(nsgfj(jset), nsgfi(iset)))
3354 30206681 : sij_rs(:, :) = TRANSPOSE(sij_contr)
3355 : ELSE
3356 1099940 : ALLOCATE (sij_rs(nsgfi(iset), nsgfj(jset)))
3357 49354022 : sij_rs(:, :) = sij_contr
3358 : END IF
3359 :
3360 398195 : DEALLOCATE (sij_contr)
3361 :
3362 : ! RI regularization
3363 : IF (.NOT. do_hfx_kpoints_prv .AND. PRESENT(regularization_RI) .AND. &
3364 : iatom == jatom .AND. iset == jset .AND. &
3365 398195 : cell_j(1) == 0 .AND. cell_j(2) == 0 .AND. cell_j(3) == 0) THEN
3366 11990 : DO i_diag = 1, nsgfi(iset)
3367 26016 : min_zet = MINVAL(zeti(:, iset))
3368 8672 : CPASSERT(min_zet > 1.0E-10_dp)
3369 : sij_rs(i_diag, i_diag) = sij_rs(i_diag, i_diag) + &
3370 11990 : regularization_RI*MAX(1.0_dp, 1.0_dp/min_zet)
3371 : END DO
3372 : END IF
3373 :
3374 : CALL block_add("IN", sij_rs, &
3375 : nsgfi(iset), nsgfj(jset), block_t%block, &
3376 398195 : sgfi, sgfj, trans=trans)
3377 398195 : DEALLOCATE (sij_rs)
3378 :
3379 : END IF
3380 : END DO
3381 : END DO
3382 : END DO
3383 :
3384 882 : CALL cp_libint_cleanup_2eri(lib)
3385 :
3386 882 : CALL neighbor_list_iterator_release(nl_iterator)
3387 7038 : DO img = 1, nimg
3388 6156 : CALL dbcsr_finalize(t2c(img))
3389 7038 : CALL dbcsr_filter(t2c(img), filter_eps)
3390 : END DO
3391 :
3392 882 : CALL timestop(handle)
3393 :
3394 1764 : END SUBROUTINE build_2c_integrals
3395 :
3396 : ! **************************************************************************************************
3397 : !> \brief ...
3398 : !> \param tensor tensor with data. Data is cleared after compression.
3399 : !> \param blk_indices ...
3400 : !> \param compressed compressed tensor data
3401 : !> \param eps all entries < eps are discarded
3402 : !> \param memory ...
3403 : ! **************************************************************************************************
3404 22202 : SUBROUTINE compress_tensor(tensor, blk_indices, compressed, eps, memory)
3405 : TYPE(dbt_type), INTENT(INOUT) :: tensor
3406 : INTEGER, ALLOCATABLE, DIMENSION(:, :), &
3407 : INTENT(INOUT) :: blk_indices
3408 : TYPE(hfx_compression_type), INTENT(INOUT) :: compressed
3409 : REAL(dp), INTENT(IN) :: eps
3410 : REAL(dp), INTENT(INOUT) :: memory
3411 :
3412 : INTEGER :: buffer_left, buffer_size, buffer_start, &
3413 : i, iblk, memory_usage, nbits, nblk, &
3414 : nints, offset, shared_offset
3415 : INTEGER(int_8) :: estimate_to_store_int, &
3416 : storage_counter_integrals
3417 : INTEGER, DIMENSION(3) :: ind
3418 : LOGICAL :: found
3419 : REAL(dp) :: spherical_estimate
3420 22202 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET :: blk_data
3421 22202 : REAL(dp), DIMENSION(:), POINTER :: blk_data_1d
3422 : TYPE(dbt_iterator_type) :: iter
3423 22202 : TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
3424 : TYPE(hfx_cache_type), POINTER :: maxval_cache
3425 22202 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
3426 : TYPE(hfx_container_type), POINTER :: maxval_container
3427 :
3428 22202 : CALL dealloc_containers(compressed, memory_usage)
3429 22202 : CALL alloc_containers(compressed, 1)
3430 :
3431 22202 : maxval_container => compressed%maxval_container(1)
3432 22202 : integral_containers => compressed%integral_containers(:, 1)
3433 :
3434 22202 : CALL hfx_init_container(maxval_container, memory_usage, .FALSE.)
3435 1443130 : DO i = 1, 64
3436 1443130 : CALL hfx_init_container(integral_containers(i), memory_usage, .FALSE.)
3437 : END DO
3438 :
3439 22202 : maxval_cache => compressed%maxval_cache(1)
3440 22202 : integral_caches => compressed%integral_caches(:, 1)
3441 :
3442 22202 : IF (ALLOCATED(blk_indices)) DEALLOCATE (blk_indices)
3443 59132 : ALLOCATE (blk_indices(dbt_get_num_blocks(tensor), 3))
3444 22202 : shared_offset = 0
3445 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices,shared_offset) &
3446 22202 : !$OMP PRIVATE(iter,ind,offset,nblk,iblk)
3447 : CALL dbt_iterator_start(iter, tensor)
3448 : nblk = dbt_iterator_num_blocks(iter)
3449 : !$OMP CRITICAL(qs_tensors_offset)
3450 : offset = shared_offset
3451 : shared_offset = shared_offset + nblk
3452 : !$OMP END CRITICAL(qs_tensors_offset)
3453 : DO iblk = 1, nblk
3454 : CALL dbt_iterator_next_block(iter, ind)
3455 : blk_indices(offset + iblk, :) = ind(:)
3456 : END DO
3457 : CALL dbt_iterator_stop(iter)
3458 : !$OMP END PARALLEL
3459 :
3460 : ! Can not use the tensor iterator here because the order of the blocks is not guaranteed.
3461 449765 : DO i = 1, SIZE(blk_indices, 1)
3462 1710252 : ind = blk_indices(i, :)
3463 427563 : CALL dbt_get_block(tensor, ind, blk_data, found)
3464 427563 : CPASSERT(found)
3465 1710252 : nints = SIZE(blk_data)
3466 427563 : blk_data_1d(1:nints) => blk_data
3467 225399711 : spherical_estimate = MAXVAL(ABS(blk_data_1d))
3468 427563 : IF (spherical_estimate == 0.0_dp) spherical_estimate = TINY(spherical_estimate)
3469 427563 : estimate_to_store_int = EXPONENT(spherical_estimate)
3470 427563 : estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8)
3471 :
3472 : CALL hfx_add_single_cache_element(estimate_to_store_int, 6, &
3473 : maxval_cache, maxval_container, memory_usage, &
3474 427563 : .FALSE.)
3475 :
3476 427563 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
3477 :
3478 427563 : nbits = EXPONENT(ANINT(spherical_estimate/eps)) + 1
3479 427563 : IF (nbits > 64) THEN
3480 : CALL cp_abort(__LOCATION__, &
3481 0 : "Overflow during tensor compression. Please use a larger EPS_FILTER or EPS_STORAGE_SCALING")
3482 : END IF
3483 :
3484 : buffer_left = nints
3485 : buffer_start = 1
3486 :
3487 953780 : DO WHILE (buffer_left > 0)
3488 526217 : buffer_size = MIN(buffer_left, cache_size)
3489 : CALL hfx_add_mult_cache_elements(blk_data_1d(buffer_start:), &
3490 : buffer_size, nbits, &
3491 : integral_caches(nbits), &
3492 : integral_containers(nbits), &
3493 : eps, 1.0_dp, &
3494 : memory_usage, &
3495 526217 : .FALSE.)
3496 526217 : buffer_left = buffer_left - buffer_size
3497 526217 : buffer_start = buffer_start + buffer_size
3498 : END DO
3499 :
3500 877328 : NULLIFY (blk_data_1d); DEALLOCATE (blk_data)
3501 : END DO
3502 :
3503 22202 : CALL dbt_clear(tensor)
3504 :
3505 22202 : storage_counter_integrals = memory_usage*cache_size
3506 22202 : memory = memory + REAL(storage_counter_integrals, dp)/1024/128
3507 : !WRITE (UNIT=iw, FMT="((T3,A,T60,I21))") &
3508 : ! "HFX_MEM_INFO| Total memory consumption ERI's RAM [MiB]: ", memory
3509 :
3510 : CALL hfx_flush_last_cache(6, maxval_cache, maxval_container, memory_usage, &
3511 22202 : .FALSE.)
3512 1443130 : DO i = 1, 64
3513 : CALL hfx_flush_last_cache(i, integral_caches(i), integral_containers(i), &
3514 1443130 : memory_usage, .FALSE.)
3515 : END DO
3516 :
3517 22202 : CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .FALSE.)
3518 1443130 : DO i = 1, 64
3519 : CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
3520 1443130 : memory_usage, .FALSE.)
3521 : END DO
3522 :
3523 44404 : END SUBROUTINE compress_tensor
3524 :
3525 : ! **************************************************************************************************
3526 : !> \brief ...
3527 : !> \param tensor empty tensor which is filled by decompressed data
3528 : !> \param blk_indices indices of blocks to be reserved
3529 : !> \param compressed compressed data
3530 : !> \param eps all entries < eps are discarded
3531 : ! **************************************************************************************************
3532 63866 : SUBROUTINE decompress_tensor(tensor, blk_indices, compressed, eps)
3533 :
3534 : TYPE(dbt_type), INTENT(INOUT) :: tensor
3535 : INTEGER, DIMENSION(:, :) :: blk_indices
3536 : TYPE(hfx_compression_type), INTENT(INOUT) :: compressed
3537 : REAL(dp), INTENT(IN) :: eps
3538 :
3539 : INTEGER :: A, b, buffer_left, buffer_size, &
3540 : buffer_start, i, memory_usage, nbits, &
3541 : nblk_per_thread, nints
3542 : INTEGER(int_8) :: estimate_to_store_int
3543 : INTEGER, DIMENSION(3) :: blk_size, ind
3544 : REAL(dp) :: spherical_estimate
3545 63866 : REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET :: blk_data
3546 63866 : REAL(dp), DIMENSION(:, :, :), POINTER :: blk_data_3d
3547 63866 : TYPE(hfx_cache_type), DIMENSION(:), POINTER :: integral_caches
3548 : TYPE(hfx_cache_type), POINTER :: maxval_cache
3549 63866 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
3550 : TYPE(hfx_container_type), POINTER :: maxval_container
3551 :
3552 63866 : maxval_cache => compressed%maxval_cache(1)
3553 63866 : maxval_container => compressed%maxval_container(1)
3554 63866 : integral_caches => compressed%integral_caches(:, 1)
3555 63866 : integral_containers => compressed%integral_containers(:, 1)
3556 :
3557 63866 : memory_usage = 0
3558 :
3559 63866 : CALL hfx_decompress_first_cache(6, maxval_cache, maxval_container, memory_usage, .FALSE.)
3560 :
3561 4151290 : DO i = 1, 64
3562 : CALL hfx_decompress_first_cache(i, integral_caches(i), integral_containers(i), &
3563 4151290 : memory_usage, .FALSE.)
3564 : END DO
3565 :
3566 : !TODO: Parallelize creation of block list.
3567 63866 : !$OMP PARALLEL DEFAULT(NONE) SHARED(tensor,blk_indices) PRIVATE(nblk_per_thread,A,b)
3568 : nblk_per_thread = SIZE(blk_indices, 1)/omp_get_num_threads() + 1
3569 : a = omp_get_thread_num()*nblk_per_thread + 1
3570 : b = MIN(a + nblk_per_thread, SIZE(blk_indices, 1))
3571 : CALL dbt_reserve_blocks(tensor, blk_indices(a:b, :))
3572 : !$OMP END PARALLEL
3573 :
3574 : ! Can not use the tensor iterator here because the order of the blocks is not guaranteed.
3575 1527616 : DO i = 1, SIZE(blk_indices, 1)
3576 5855000 : ind = blk_indices(i, :)
3577 1463750 : CALL dbt_blk_sizes(tensor, ind, blk_size)
3578 5855000 : nints = PRODUCT(blk_size)
3579 : CALL hfx_get_single_cache_element( &
3580 : estimate_to_store_int, 6, &
3581 : maxval_cache, maxval_container, memory_usage, &
3582 1463750 : .FALSE.)
3583 :
3584 1463750 : spherical_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1)
3585 :
3586 1463750 : nbits = EXPONENT(ANINT(spherical_estimate/eps)) + 1
3587 :
3588 1463750 : buffer_left = nints
3589 1463750 : buffer_start = 1
3590 :
3591 4391250 : ALLOCATE (blk_data(nints))
3592 3122039 : DO WHILE (buffer_left > 0)
3593 1658289 : buffer_size = MIN(buffer_left, cache_size)
3594 : CALL hfx_get_mult_cache_elements(blk_data(buffer_start), &
3595 : buffer_size, nbits, &
3596 : integral_caches(nbits), &
3597 : integral_containers(nbits), &
3598 : eps, 1.0_dp, &
3599 : memory_usage, &
3600 1658289 : .FALSE.)
3601 1658289 : buffer_left = buffer_left - buffer_size
3602 1658289 : buffer_start = buffer_start + buffer_size
3603 : END DO
3604 :
3605 1463750 : blk_data_3d(1:blk_size(1), 1:blk_size(2), 1:blk_size(3)) => blk_data
3606 1463750 : CALL dbt_put_block(tensor, ind, blk_size, blk_data_3d)
3607 1527616 : NULLIFY (blk_data_3d); DEALLOCATE (blk_data)
3608 : END DO
3609 :
3610 63866 : CALL hfx_reset_cache_and_container(maxval_cache, maxval_container, memory_usage, .FALSE.)
3611 4151290 : DO i = 1, 64
3612 : CALL hfx_reset_cache_and_container(integral_caches(i), integral_containers(i), &
3613 4151290 : memory_usage, .FALSE.)
3614 : END DO
3615 127732 : END SUBROUTINE decompress_tensor
3616 :
3617 : ! **************************************************************************************************
3618 : !> \brief ...
3619 : !> \param tensor ...
3620 : !> \param nze ...
3621 : !> \param occ ...
3622 : ! **************************************************************************************************
3623 139344 : SUBROUTINE get_tensor_occupancy(tensor, nze, occ)
3624 : TYPE(dbt_type), INTENT(IN) :: tensor
3625 : INTEGER(int_8), INTENT(OUT) :: nze
3626 : REAL(dp), INTENT(OUT) :: occ
3627 :
3628 278688 : INTEGER, DIMENSION(dbt_ndims(tensor)) :: dims
3629 :
3630 139344 : nze = dbt_get_nze_total(tensor)
3631 139344 : CALL dbt_get_info(tensor, nfull_total=dims)
3632 529578 : occ = REAL(nze, dp)/PRODUCT(REAL(dims, dp))
3633 :
3634 139344 : END SUBROUTINE get_tensor_occupancy
3635 :
3636 0 : END MODULE qs_tensors
|
