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