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