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