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