46 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
47 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
48 dbt_distribution_destroy, dbt_distribution_new, dbt_distribution_type, dbt_filter, &
49 dbt_finalize, dbt_get_block, dbt_get_info, dbt_get_stored_coordinates, &
50 dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, &
51 dbt_iterator_type, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
52 dbt_pgrid_type, dbt_put_block, dbt_scale, dbt_type
114#include "./base/base_uses.f90"
123 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'hfx_ri_kp'
133 SUBROUTINE adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
134 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
135 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
136 TYPE(qs_environment_type),
POINTER :: qs_env
138 INTEGER :: i_img, i_RI, i_spin, iatom, natom, &
139 nblks_RI, nimg, nkind, nspins
140 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_RI_ext, dist1, dist2, dist3
141 TYPE(dft_control_type),
POINTER :: dft_control
142 TYPE(mp_para_env_type),
POINTER :: para_env
144 NULLIFY (dft_control, para_env)
151 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, para_env=para_env, nkind=nkind)
155 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
156 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
157 DO i_ri = 1, ri_data%ncell_RI
158 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
161 ALLOCATE (ri_data%t_3c_int_ctr_1(1, nimg))
163 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
164 ri_data%bsizes_AO_split, [1, 2], [3], name=
"(AO RI | AO)")
167 CALL dbt_create(ri_data%t_3c_int_ctr_1(1, 1), ri_data%t_3c_int_ctr_1(1, i_img))
169 DEALLOCATE (dist1, dist2, dist3)
171 ALLOCATE (ri_data%t_3c_int_ctr_2(1, 1))
173 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
174 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
175 DEALLOCATE (dist1, dist2, dist3)
178 DEALLOCATE (bsizes_ri_ext)
179 nblks_ri =
SIZE(ri_data%bsizes_RI)
180 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
181 DO i_ri = 1, ri_data%ncell_RI
182 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
185 ALLOCATE (ri_data%t_2c_inv(1, natom), ri_data%t_2c_int(1, natom), ri_data%t_2c_pot(1, natom))
186 CALL create_2c_tensor(ri_data%t_2c_inv(1, 1), dist1, dist2, ri_data%pgrid_2d, &
187 bsizes_ri_ext, bsizes_ri_ext, &
189 DEALLOCATE (dist1, dist2)
190 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, 1))
191 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, 1))
193 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_inv(1, iatom))
194 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, iatom))
195 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, iatom))
198 ALLOCATE (ri_data%kp_cost(natom, natom, nimg))
199 ri_data%kp_cost = 0.0_dp
202 nspins = dft_control%nspins
203 ALLOCATE (ri_data%rho_ao_t(nspins, nimg), ri_data%ks_t(nspins, nimg))
204 CALL create_2c_tensor(ri_data%rho_ao_t(1, 1), dist1, dist2, ri_data%pgrid_2d, &
205 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
207 DEALLOCATE (dist1, dist2)
209 CALL dbt_create(dbcsr_template, ri_data%ks_t(1, 1))
211 IF (nspins == 2)
THEN
212 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
213 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
216 DO i_spin = 1, nspins
217 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(i_spin, i_img))
218 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(i_spin, i_img))
222 END SUBROUTINE adapt_ri_data_to_kp
230 SUBROUTINE hfx_ri_pre_scf_kp(dbcsr_template, ri_data, qs_env)
231 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
232 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
233 TYPE(qs_environment_type),
POINTER :: qs_env
235 CHARACTER(LEN=*),
PARAMETER :: routineN =
'hfx_ri_pre_scf_kp'
237 INTEGER :: handle, i_img, iatom, natom, nimg, nkind
238 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_op_pot, t_2c_op_RI
239 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_int
240 TYPE(dft_control_type),
POINTER :: dft_control
242 NULLIFY (dft_control)
244 CALL timeset(routinen, handle)
246 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, nkind=nkind)
248 CALL cleanup_kp(ri_data)
251 IF (ri_data%flavor /=
ri_pmat) cpabort(
"K-points RI-HFX only with RHO flavor")
252 IF (ri_data%same_op) ri_data%same_op = .false.
253 IF (abs(ri_data%eps_pgf_orb - dft_control%qs_control%eps_pgf_orb) > 1.0e-16_dp) &
254 cpabort(
"RI%EPS_PGF_ORB and QS%EPS_PGF_ORB must be identical for RI-HFX k-points")
256 CALL get_kp_and_ri_images(ri_data, qs_env)
260 ALLOCATE (t_2c_op_pot(nimg), t_2c_op_ri(nimg))
261 ALLOCATE (t_3c_int(1, nimg))
265 CALL adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
270 CALL get_ext_2c_int(ri_data%t_2c_inv(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
274 CALL get_ext_2c_int(ri_data%t_2c_int(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, &
275 qs_env, off_diagonal=.true.)
276 CALL apply_bump(ri_data%t_2c_int(1, iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
279 CALL get_ext_2c_int(ri_data%t_2c_pot(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
280 do_inverse=.true., skip_inverse=.true.)
281 CALL apply_bump(ri_data%t_2c_pot(1, iatom), iatom, ri_data, qs_env, from_left=.true., &
282 from_right=.true., debump=.true.)
290 ALLOCATE (ri_data%kp_mat_2c_pot(1, nimg))
292 CALL dbcsr_create(ri_data%kp_mat_2c_pot(1, i_img), template=t_2c_op_pot(i_img))
293 CALL dbcsr_copy(ri_data%kp_mat_2c_pot(1, i_img), t_2c_op_pot(i_img))
298 CALL reorder_3c_ints(t_3c_int(1, :), ri_data)
302 CALL precontract_3c_ints(t_3c_int, ri_data, qs_env)
304 CALL timestop(handle)
306 END SUBROUTINE hfx_ri_pre_scf_kp
312 SUBROUTINE cleanup_kp(ri_data)
313 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
317 IF (
ALLOCATED(ri_data%kp_cost))
DEALLOCATE (ri_data%kp_cost)
318 IF (
ALLOCATED(ri_data%idx_to_img))
DEALLOCATE (ri_data%idx_to_img)
319 IF (
ALLOCATED(ri_data%img_to_idx))
DEALLOCATE (ri_data%img_to_idx)
320 IF (
ALLOCATED(ri_data%present_images))
DEALLOCATE (ri_data%present_images)
321 IF (
ALLOCATED(ri_data%img_to_RI_cell))
DEALLOCATE (ri_data%img_to_RI_cell)
322 IF (
ALLOCATED(ri_data%RI_cell_to_img))
DEALLOCATE (ri_data%RI_cell_to_img)
324 IF (
ALLOCATED(ri_data%kp_mat_2c_pot))
THEN
325 DO j = 1,
SIZE(ri_data%kp_mat_2c_pot, 2)
326 DO i = 1,
SIZE(ri_data%kp_mat_2c_pot, 1)
330 DEALLOCATE (ri_data%kp_mat_2c_pot)
333 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
334 DO i = 1,
SIZE(ri_data%kp_t_3c_int)
335 CALL dbt_destroy(ri_data%kp_t_3c_int(i))
337 DEALLOCATE (ri_data%kp_t_3c_int)
340 IF (
ALLOCATED(ri_data%t_2c_inv))
THEN
341 DO j = 1,
SIZE(ri_data%t_2c_inv, 2)
342 DO i = 1,
SIZE(ri_data%t_2c_inv, 1)
343 CALL dbt_destroy(ri_data%t_2c_inv(i, j))
346 DEALLOCATE (ri_data%t_2c_inv)
349 IF (
ALLOCATED(ri_data%t_2c_int))
THEN
350 DO j = 1,
SIZE(ri_data%t_2c_int, 2)
351 DO i = 1,
SIZE(ri_data%t_2c_int, 1)
352 CALL dbt_destroy(ri_data%t_2c_int(i, j))
355 DEALLOCATE (ri_data%t_2c_int)
358 IF (
ALLOCATED(ri_data%t_2c_pot))
THEN
359 DO j = 1,
SIZE(ri_data%t_2c_pot, 2)
360 DO i = 1,
SIZE(ri_data%t_2c_pot, 1)
361 CALL dbt_destroy(ri_data%t_2c_pot(i, j))
364 DEALLOCATE (ri_data%t_2c_pot)
367 IF (
ALLOCATED(ri_data%t_3c_int_ctr_1))
THEN
368 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_1, 2)
369 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_1, 1)
370 CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
373 DEALLOCATE (ri_data%t_3c_int_ctr_1)
376 IF (
ALLOCATED(ri_data%t_3c_int_ctr_2))
THEN
377 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_2, 2)
378 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_2, 1)
379 CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
382 DEALLOCATE (ri_data%t_3c_int_ctr_2)
385 IF (
ALLOCATED(ri_data%rho_ao_t))
THEN
386 DO j = 1,
SIZE(ri_data%rho_ao_t, 2)
387 DO i = 1,
SIZE(ri_data%rho_ao_t, 1)
388 CALL dbt_destroy(ri_data%rho_ao_t(i, j))
391 DEALLOCATE (ri_data%rho_ao_t)
394 IF (
ALLOCATED(ri_data%ks_t))
THEN
395 DO j = 1,
SIZE(ri_data%ks_t, 2)
396 DO i = 1,
SIZE(ri_data%ks_t, 1)
397 CALL dbt_destroy(ri_data%ks_t(i, j))
400 DEALLOCATE (ri_data%ks_t)
403 END SUBROUTINE cleanup_kp
412 SUBROUTINE print_progress_bar(b_img, nimg, iprint, ri_data)
413 INTEGER,
INTENT(IN) :: b_img, nimg
414 INTEGER,
INTENT(INOUT) :: iprint
415 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
417 CHARACTER(LEN=default_string_length) :: bar
420 IF (ri_data%unit_nr > 0)
THEN
422 WRITE (ri_data%unit_nr,
'(/T6,A)', advance=
"no")
'[-'
425 IF (b_img > iprint*nimg/71)
THEN
426 rep = max(1, 71/nimg)
427 bar = repeat(
"-", rep)
428 WRITE (ri_data%unit_nr,
'(A)', advance=
"no") trim(bar)
432 IF (b_img == nimg)
THEN
433 rep = max(0, 1 + 71 - iprint*rep)
434 bar = repeat(
"-", rep)
435 WRITE (ri_data%unit_nr,
'(A,A)') trim(bar),
'-]'
440 END SUBROUTINE print_progress_bar
454 geometry_did_change, nspins, hf_fraction)
458 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
459 REAL(kind=
dp),
INTENT(OUT) :: ehfx
461 LOGICAL,
INTENT(IN) :: geometry_did_change
462 INTEGER,
INTENT(IN) :: nspins
463 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
465 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_ks_kp'
467 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_spin, iatom, &
468 iblk, igroup, iprint, jatom, mb_img, n_batch_nze, n_nze, natom, ngroups, nimg, nimg_nze
469 INTEGER(int_8) :: mem, nflop, nze
470 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges_at, batch_ranges_nze, &
472 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
473 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: sparsity_pattern
474 LOGICAL :: estimate_mem, print_progress, use_delta_p
475 REAL(
dp) :: etmp,
fac, occ, pfac, pref, t1, t2, t3, &
478 TYPE(
dbcsr_type) :: ks_desymm, rho_desymm, tmp
479 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
481 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ks_t_split, t_2c_ao_tmp, t_2c_work, &
482 t_3c_int, t_3c_work_2, t_3c_work_3
483 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: ks_t, ks_t_sub, t_3c_apc, t_3c_apc_sub
487 NULLIFY (para_env, para_env_sub, blacs_env_sub, hfx_section, dbcsr_template, print_section)
491 CALL timeset(routinen, handle)
493 CALL get_qs_env(qs_env, para_env=para_env, natom=natom)
495 IF (nspins == 1)
THEN
496 fac = 0.5_dp*hf_fraction
498 fac = 1.0_dp*hf_fraction
505 ri_data%kp_stack_size = batch_size
506 ri_data%kp_ngroups = ngroups
508 IF (geometry_did_change)
THEN
509 CALL hfx_ri_pre_scf_kp(ks_matrix(1, 1)%matrix, ri_data, qs_env)
512 nimg_nze = ri_data%nimg_nze
535 IF (.NOT. use_delta_p) pfac = 0.0_dp
536 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, pfac, ri_data, qs_env)
540 DO i_spin = 1, nspins
547 IF (n_nze == nspins)
THEN
548 cpwarn(
"It is highly recommended to restart from a converged GGA K-point calculations.")
551 ALLOCATE (ks_t(nspins, nimg))
553 DO i_spin = 1, nspins
554 CALL dbt_create(ri_data%ks_t(1, 1), ks_t(i_spin, i_img))
558 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
559 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
567 ALLOCATE (t_3c_apc(nspins, nimg))
569 DO i_spin = 1, nspins
570 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
573 CALL contract_pmat_3c(t_3c_apc, ri_data%rho_ao_t, ri_data, qs_env)
575 IF (mod(para_env%num_pe, ngroups) /= 0)
THEN
576 cpwarn(
"KP_NGROUPS must be an integer divisor of the total number of MPI ranks. It was set to 1.")
580 IF ((mod(ngroups, natom) /= 0) .AND. (mod(natom, ngroups) /= 0) .AND. geometry_did_change)
THEN
581 IF (ngroups > 1)
THEN
582 cpwarn(
"Better load balancing is reached if NGROUPS is a multiple/divisor of the number of atoms")
585 group_size = para_env%num_pe/ngroups
586 igroup = para_env%mepos/group_size
588 ALLOCATE (para_env_sub)
589 CALL para_env_sub%from_split(para_env, igroup)
593 ALLOCATE (sparsity_pattern(natom, natom, nimg))
594 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
595 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
598 ALLOCATE (mat_2c_pot(nimg), ks_t_sub(nspins, nimg), t_2c_ao_tmp(1), ks_t_split(2), t_2c_work(3))
599 CALL get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
600 group_size, ngroups, para_env, para_env_sub, ri_data)
602 ALLOCATE (t_3c_int(nimg), t_3c_apc_sub(nspins, nimg), t_3c_work_2(3), t_3c_work_3(3))
603 CALL get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
604 group_size, ngroups, para_env, para_env_sub, ri_data)
608 ALLOCATE (batch_ranges_at(natom + 1))
609 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
611 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
612 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
614 batch_ranges_at(iatom) = iblk
618 n_batch_nze = nimg_nze/batch_size
619 IF (
modulo(nimg_nze, batch_size) /= 0) n_batch_nze = n_batch_nze + 1
620 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
621 DO i_batch = 1, n_batch_nze
622 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
624 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
630 ALLOCATE (iapc_pairs(nimg, 2))
631 IF (estimate_mem .AND. geometry_did_change)
THEN
633 CALL get_iapc_pairs(iapc_pairs, 1, ri_data, qs_env)
634 CALL fill_3c_stack(t_3c_work_3(1), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
635 filter_at=1, filter_dim=2, idx_to_at=idx_to_at_ao, &
636 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
637 CALL fill_3c_stack(t_3c_work_3(2), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
638 filter_at=1, filter_dim=2, idx_to_at=idx_to_at_ao, &
639 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
640 CALL fill_3c_stack(t_3c_work_2(1), t_3c_apc_sub(1, :), iapc_pairs(:, 2), 3, &
641 ri_data, filter_at=1, filter_dim=1, idx_to_at=idx_to_at_ao, &
642 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
643 CALL fill_3c_stack(t_3c_work_2(2), t_3c_apc_sub(1, :), iapc_pairs(:, 2), 3, &
644 ri_data, filter_at=1, filter_dim=1, idx_to_at=idx_to_at_ao, &
645 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
646 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, 1, 1, 1, ri_data, qs_env, &
647 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
648 dbcsr_template=dbcsr_template)
650 CALL para_env%max(mem)
651 CALL dbt_clear(t_3c_work_2(1))
652 CALL dbt_clear(t_3c_work_2(2))
653 CALL dbt_clear(t_3c_work_3(1))
654 CALL dbt_clear(t_3c_work_3(2))
655 CALL dbt_clear(t_2c_work(1))
657 IF (ri_data%unit_nr > 0)
THEN
658 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I14)") &
659 "KP-HFX_RI_INFO| Estimated peak memory usage per MPI rank (MiB):", mem/(1024*1024)
664 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
665 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
666 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
667 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
671 ri_data%kp_cost(:, :, :) = 0.0_dp
673 IF (print_progress)
CALL print_progress_bar(b_img, nimg, iprint, ri_data)
674 CALL dbt_batched_contract_init(ks_t_split(1))
675 CALL dbt_batched_contract_init(ks_t_split(2))
678 IF (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup) cycle
680 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp
686 CALL timeset(routinen//
"_2c", handle2)
687 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
688 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
689 dbcsr_template=dbcsr_template)
690 CALL dbt_copy(t_2c_work(1), t_2c_work(2), move_data=.true.)
691 CALL dbt_filter(t_2c_work(2), ri_data%filter_eps)
692 CALL timestop(handle2)
694 CALL dbt_batched_contract_init(t_2c_work(2))
695 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env)
696 CALL timeset(routinen//
"_3c", handle2)
699 DO i_batch = 1, n_batch_nze
700 CALL fill_3c_stack(t_3c_work_3(3), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
701 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
702 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
703 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1), move_data=.true.)
705 CALL dbt_contract(1.0_dp, t_2c_work(2), t_3c_work_3(1), &
706 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
707 contract_1=[2], notcontract_1=[1], &
708 contract_2=[1], notcontract_2=[2, 3], &
709 filter_eps=ri_data%filter_eps, flop=nflop)
710 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
711 CALL dbt_copy(t_3c_work_3(2), t_3c_work_2(2), order=[2, 1, 3], move_data=.true.)
712 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1))
716 DO i_spin = 1, nspins
717 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
718 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
719 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
723 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
724 CALL dbt_contract(-pref*
fac, t_3c_work_2(1), t_3c_work_2(2), &
725 1.0_dp, ks_t_split(i_spin), map_1=[1], map_2=[2], &
726 contract_1=[2, 3], notcontract_1=[1], &
727 contract_2=[2, 3], notcontract_2=[1], &
728 filter_eps=ri_data%filter_eps, &
729 move_data=i_spin == nspins, flop=nflop)
730 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
733 CALL timestop(handle2)
734 CALL dbt_batched_contract_finalize(t_2c_work(2))
737 ri_data%kp_cost(iatom, jatom, b_img) = t4 - t3
740 CALL dbt_batched_contract_finalize(ks_t_split(1))
741 CALL dbt_batched_contract_finalize(ks_t_split(2))
743 DO i_spin = 1, nspins
744 CALL dbt_copy(ks_t_split(i_spin), t_2c_ao_tmp(1), move_data=.true.)
745 CALL dbt_copy(t_2c_ao_tmp(1), ks_t_sub(i_spin, b_img), summation=.true.)
748 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
749 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
750 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
751 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
753 CALL para_env%sum(ri_data%dbcsr_nflop)
754 CALL para_env%sum(ri_data%kp_cost)
756 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
759 CALL gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
763 CALL dbt_copy(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img), move_data=.true.)
767 CALL dbt_destroy(t_2c_ao_tmp(1))
768 CALL dbt_destroy(ks_t_split(1))
769 CALL dbt_destroy(ks_t_split(2))
770 CALL dbt_destroy(t_2c_work(1))
771 CALL dbt_destroy(t_2c_work(2))
772 CALL dbt_destroy(t_3c_work_2(1))
773 CALL dbt_destroy(t_3c_work_2(2))
774 CALL dbt_destroy(t_3c_work_2(3))
775 CALL dbt_destroy(t_3c_work_3(1))
776 CALL dbt_destroy(t_3c_work_3(2))
777 CALL dbt_destroy(t_3c_work_3(3))
779 CALL dbt_destroy(t_3c_int(i_img))
781 DO i_spin = 1, nspins
782 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
783 CALL dbt_destroy(ks_t_sub(i_spin, i_img))
786 IF (
ASSOCIATED(dbcsr_template))
THEN
788 DEALLOCATE (dbcsr_template)
793 CALL para_env_sub%free()
794 DEALLOCATE (para_env_sub)
799 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
800 DO i_spin = 1, nspins
802 CALL dbt_copy(ks_t(i_spin, b_img), ri_data%ks_t(i_spin, b_img), summation=.true.)
805 mb_img = get_opp_index(b_img, qs_env)
806 IF (mb_img > 0 .AND. mb_img <= nimg)
THEN
807 CALL dbt_copy(ks_t(i_spin, mb_img), ri_data%ks_t(i_spin, b_img), order=[2, 1], summation=.true.)
812 DO i_spin = 1, nspins
813 CALL dbt_destroy(ks_t(i_spin, b_img))
818 CALL dbt_create(ri_data%ks_t(1, 1), t_2c_ao_tmp(1))
819 CALL dbcsr_create(tmp, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
820 CALL dbcsr_create(ks_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
821 CALL dbcsr_create(rho_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
824 DO i_spin = 1, nspins
825 CALL dbt_filter(ri_data%ks_t(i_spin, i_img), ri_data%filter_eps)
826 CALL dbt_copy(ri_data%ks_t(i_spin, i_img), t_2c_ao_tmp(1))
827 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), ks_desymm)
828 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), tmp)
829 CALL dbcsr_add(ks_matrix(i_spin, i_img)%matrix, tmp, 1.0_dp, 1.0_dp)
831 CALL dbt_copy(ri_data%rho_ao_t(i_spin, i_img), t_2c_ao_tmp(1))
832 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), rho_desymm)
834 CALL dbcsr_dot(ks_desymm, rho_desymm, etmp)
835 ehfx = ehfx + 0.5_dp*etmp
837 IF (.NOT. use_delta_p)
CALL dbt_clear(ri_data%ks_t(i_spin, i_img))
843 CALL dbt_destroy(t_2c_ao_tmp(1))
845 CALL timestop(handle)
864 INTEGER,
INTENT(IN) :: nspins
865 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
867 LOGICAL,
INTENT(IN),
OPTIONAL :: use_virial
869 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_forces_kp'
871 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_loop, i_spin, &
872 i_xyz, iatom, iblk, igroup, j_xyz, jatom, k_xyz, n_batch, natom, ngroups, nimg, nimg_nze
873 INTEGER(int_8) :: nflop, nze
874 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, batch_ranges_at, &
875 batch_ranges_nze, dist1, dist2, &
876 i_images, idx_to_at_ao, idx_to_at_ri, &
878 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
879 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: force_pattern, sparsity_pattern
880 INTEGER,
DIMENSION(2, 1) :: bounds_iat, bounds_jat
881 LOGICAL :: use_virial_prv
882 REAL(
dp) ::
fac, occ, pref, t1, t2
883 REAL(
dp),
DIMENSION(3, 3) :: work_virial
887 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
888 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_pot, mat_der_pot_sub
890 TYPE(dbt_type) :: t_2c_r, t_2c_r_split
891 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_bint, t_2c_binv, t_2c_der_pot, &
892 t_2c_inv, t_2c_metric, t_2c_work, &
893 t_3c_der_stack, t_3c_work_2, &
895 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
896 t_2c_der_metric_sub, t_3c_apc, t_3c_apc_sub, t_3c_der_ao, t_3c_der_ao_sub, t_3c_der_ri, &
904 NULLIFY (para_env, para_env_sub, hfx_section, blacs_env_sub, dbcsr_template, force, atomic_kind_set, &
905 virial, particle_set, cell)
907 CALL timeset(routinen, handle)
909 use_virial_prv = .false.
910 IF (
PRESENT(use_virial)) use_virial_prv = use_virial
912 IF (nspins == 1)
THEN
913 fac = 0.5_dp*hf_fraction
915 fac = 1.0_dp*hf_fraction
918 CALL get_qs_env(qs_env, natom=natom, para_env=para_env, force=force, cell=cell, virial=virial, &
919 atomic_kind_set=atomic_kind_set, particle_set=particle_set)
922 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
923 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
925 ALLOCATE (idx_to_at_ri(
SIZE(ri_data%bsizes_RI_split)))
926 CALL get_idx_to_atom(idx_to_at_ri, ri_data%bsizes_RI_split, ri_data%bsizes_RI)
929 ALLOCATE (t_3c_der_ri(nimg, 3), t_3c_der_ao(nimg, 3), mat_der_pot(nimg, 3), t_2c_der_metric(natom, 3))
933 CALL precalc_derivatives(t_3c_der_ri, t_3c_der_ao, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
936 ALLOCATE (rho_ao_t(nspins, nimg))
938 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
940 DEALLOCATE (dist1, dist2)
941 IF (nspins == 2)
CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(2, 1))
943 DO i_spin = 1, nspins
944 CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(i_spin, i_img))
947 CALL get_pmat_images(rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
950 ALLOCATE (t_3c_apc(nspins, nimg))
952 DO i_spin = 1, nspins
953 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
956 CALL contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
961 group_size = para_env%num_pe/ngroups
962 igroup = para_env%mepos/group_size
964 ALLOCATE (para_env_sub)
965 CALL para_env_sub%from_split(para_env, igroup)
969 ALLOCATE (sparsity_pattern(natom, natom, nimg))
970 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
971 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
974 ALLOCATE (t_2c_inv(natom), mat_2c_pot(nimg), rho_ao_t_sub(nspins, nimg), t_2c_work(5), &
975 t_2c_der_metric_sub(natom, 3), mat_der_pot_sub(nimg, 3), t_2c_bint(natom), &
976 t_2c_metric(natom), t_2c_binv(natom))
977 CALL get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
978 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
979 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
980 CALL dbt_create(t_2c_work(1), t_2c_r)
981 CALL dbt_create(t_2c_work(5), t_2c_r_split)
983 ALLOCATE (t_2c_der_pot(3))
985 CALL dbt_create(t_2c_r, t_2c_der_pot(i_xyz))
989 ALLOCATE (t_3c_work_2(3), t_3c_work_3(4), t_3c_der_stack(6), t_3c_der_ao_sub(nimg, 3), &
990 t_3c_der_ri_sub(nimg, 3), t_3c_apc_sub(nspins, nimg))
991 CALL get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_ao, t_3c_der_ao_sub, &
992 t_3c_der_ri, t_3c_der_ri_sub, t_3c_apc, t_3c_apc_sub, t_3c_der_stack, &
993 group_size, ngroups, para_env, para_env_sub, ri_data)
996 ALLOCATE (batch_ranges_at(natom + 1))
997 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
999 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
1000 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
1002 batch_ranges_at(iatom) = iblk
1006 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
1007 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
1008 CALL dbt_batched_contract_init(t_3c_work_3(3), batch_range_2=batch_ranges_at)
1009 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
1010 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
1013 nimg_nze = ri_data%nimg_nze
1014 batch_size = ri_data%kp_stack_size
1015 n_batch = nimg_nze/batch_size
1016 IF (
modulo(nimg_nze, batch_size) /= 0) n_batch = n_batch + 1
1017 ALLOCATE (batch_ranges_nze(n_batch + 1))
1018 DO i_batch = 1, n_batch
1019 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
1021 batch_ranges_nze(n_batch + 1) = nimg_nze + 1
1026 CALL dbt_create(t_2c_inv(iatom), t_2c_binv(iatom))
1027 CALL dbt_copy(t_2c_inv(iatom), t_2c_binv(iatom))
1028 CALL apply_bump(t_2c_binv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
1029 CALL apply_bump(t_2c_inv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1033 work_virial = 0.0_dp
1034 ALLOCATE (iapc_pairs(nimg, 2), i_images(nimg))
1035 ALLOCATE (force_pattern(natom, natom, nimg))
1036 force_pattern(:, :, :) = -1
1045 IF (i_loop == 1 .AND. (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup)) cycle
1046 IF (i_loop == 2 .AND. (.NOT. force_pattern(iatom, jatom, b_img) == igroup)) cycle
1049 CALL timeset(routinen//
"_2c_1", handle2)
1050 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
1051 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
1052 dbcsr_template=dbcsr_template)
1053 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_inv(jatom), &
1054 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1055 contract_1=[2], notcontract_1=[1], &
1056 contract_2=[1], notcontract_2=[2], &
1057 filter_eps=ri_data%filter_eps, flop=nflop)
1058 CALL dbt_copy(t_2c_work(2), t_2c_work(5), move_data=.true.)
1059 CALL dbt_filter(t_2c_work(5), ri_data%filter_eps)
1060 CALL timestop(handle2)
1062 CALL timeset(routinen//
"_3c", handle2)
1063 bounds_iat(:, 1) = [sum(ri_data%bsizes_AO(1:iatom - 1)) + 1, sum(ri_data%bsizes_AO(1:iatom))]
1064 bounds_jat(:, 1) = [sum(ri_data%bsizes_AO(1:jatom - 1)) + 1, sum(ri_data%bsizes_AO(1:jatom))]
1065 CALL dbt_clear(t_2c_r_split)
1067 DO i_spin = 1, nspins
1068 CALL dbt_batched_contract_init(rho_ao_t_sub(i_spin, b_img))
1071 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env, i_images)
1072 DO i_batch = 1, n_batch
1076 CALL dbt_clear(t_3c_der_stack(i_xyz))
1077 CALL fill_3c_stack(t_3c_der_stack(i_xyz), t_3c_der_ri_sub(:, i_xyz), &
1078 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
1079 filter_dim=2, idx_to_at=idx_to_at_ao, &
1080 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1082 CALL dbt_clear(t_3c_der_stack(3 + i_xyz))
1083 CALL fill_3c_stack(t_3c_der_stack(3 + i_xyz), t_3c_der_ao_sub(:, i_xyz), &
1084 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
1085 filter_dim=2, idx_to_at=idx_to_at_ao, &
1086 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1089 DO i_spin = 1, nspins
1091 CALL dbt_clear(t_3c_work_2(3))
1092 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
1093 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
1094 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1097 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
1101 CALL dbt_contract(1.0_dp, rho_ao_t_sub(i_spin, b_img), t_3c_work_2(1), &
1102 0.0_dp, t_3c_work_2(2), map_1=[1], map_2=[2, 3], &
1103 contract_1=[1], notcontract_1=[2], &
1104 contract_2=[1], notcontract_2=[2, 3], &
1105 bounds_1=bounds_iat, bounds_2=bounds_jat, &
1106 filter_eps=ri_data%filter_eps, flop=nflop)
1107 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1113 CALL dbt_copy(t_3c_work_2(2), t_3c_work_3(1), order=[2, 1, 3], move_data=.true.)
1114 CALL dbt_batched_contract_init(t_2c_work(5))
1115 CALL dbt_contract(1.0_dp, t_2c_work(5), t_3c_work_3(1), &
1116 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
1117 contract_1=[1], notcontract_1=[2], &
1118 contract_2=[1], notcontract_2=[2, 3], &
1119 filter_eps=ri_data%filter_eps, flop=nflop)
1120 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1121 CALL dbt_batched_contract_finalize(t_2c_work(5))
1124 CALL dbt_copy(t_3c_work_3(2), t_3c_work_3(4), move_data=.true.)
1125 IF (use_virial_prv)
THEN
1127 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1128 idx_to_at_ri, idx_to_at_ao, i_images, &
1129 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1130 ri_data, qs_env, work_virial, cell, particle_set)
1133 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1134 idx_to_at_ri, idx_to_at_ao, i_images, &
1135 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1138 CALL dbt_clear(t_3c_work_3(4))
1142 IF (i_loop == 2) cycle
1145 CALL fill_3c_stack(t_3c_work_3(4), ri_data%kp_t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
1146 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
1147 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1148 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(3), move_data=.true.)
1150 CALL dbt_batched_contract_init(t_2c_r_split)
1151 CALL dbt_contract(1.0_dp, t_3c_work_3(1), t_3c_work_3(3), &
1152 1.0_dp, t_2c_r_split, map_1=[1], map_2=[2], &
1153 contract_1=[2, 3], notcontract_1=[1], &
1154 contract_2=[2, 3], notcontract_2=[1], &
1155 filter_eps=ri_data%filter_eps, flop=nflop)
1156 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1157 CALL dbt_batched_contract_finalize(t_2c_r_split)
1158 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(1))
1161 DO i_spin = 1, nspins
1162 CALL dbt_batched_contract_finalize(rho_ao_t_sub(i_spin, b_img))
1164 CALL timestop(handle2)
1166 IF (i_loop == 2) cycle
1168 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp*pref
1170 CALL timeset(routinen//
"_2c_2", handle2)
1172 CALL dbt_copy(t_2c_r_split, t_2c_r, move_data=.true.)
1174 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
1175 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
1176 dbcsr_template=dbcsr_template)
1189 CALL get_ext_2c_int(t_2c_der_pot(i_xyz), mat_der_pot_sub(:, i_xyz), iatom, jatom, &
1190 b_img, ri_data, qs_env, blacs_env_ext=blacs_env_sub, &
1191 para_env_ext=para_env_sub, dbcsr_template=dbcsr_template)
1194 IF (use_virial_prv)
THEN
1195 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1196 b_img, pref, ri_data, qs_env, work_virial, cell, particle_set)
1198 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1199 b_img, pref, ri_data, qs_env)
1203 CALL dbt_clear(t_2c_der_pot(i_xyz))
1207 CALL dbt_contract(1.0_dp, t_2c_metric(iatom), t_2c_r, &
1208 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1209 contract_1=[2], notcontract_1=[1], &
1210 contract_2=[1], notcontract_2=[2], &
1211 filter_eps=ri_data%filter_eps, flop=nflop)
1212 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1213 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_work(1), &
1214 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1215 contract_1=[2], notcontract_1=[1], &
1216 contract_2=[2], notcontract_2=[1], &
1217 filter_eps=ri_data%filter_eps, flop=nflop)
1218 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1222 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(iatom), &
1223 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1224 contract_1=[2], notcontract_1=[1], &
1225 contract_2=[1], notcontract_2=[2], &
1226 filter_eps=ri_data%filter_eps, flop=nflop)
1227 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1229 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(3), &
1230 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1231 contract_1=[1], notcontract_1=[2], &
1232 contract_2=[1], notcontract_2=[2], &
1233 filter_eps=ri_data%filter_eps, flop=nflop)
1234 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1236 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1237 CALL get_2c_bump_forces(force, t_2c_work(4), iatom, atom_of_kind, kind_of, pref, &
1238 ri_data, qs_env, work_virial)
1241 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(2), &
1242 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1243 contract_1=[1], notcontract_1=[2], &
1244 contract_2=[1], notcontract_2=[2], &
1245 filter_eps=ri_data%filter_eps, flop=nflop)
1246 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1248 IF (use_virial_prv)
THEN
1249 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1250 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1251 diag=.true., offdiag=.false.)
1253 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1254 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1258 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1259 CALL apply_bump(t_2c_work(2), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1261 IF (use_virial_prv)
THEN
1262 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1263 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1264 diag=.false., offdiag=.true.)
1266 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1267 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1272 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(iatom), &
1273 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1274 contract_1=[2], notcontract_1=[1], &
1275 contract_2=[1], notcontract_2=[2], &
1276 filter_eps=ri_data%filter_eps, flop=nflop)
1277 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1279 CALL dbt_contract(1.0_dp, t_2c_bint(iatom), t_2c_work(4), &
1280 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1281 contract_1=[1], notcontract_1=[2], &
1282 contract_2=[1], notcontract_2=[2], &
1283 filter_eps=ri_data%filter_eps, flop=nflop)
1284 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1286 CALL get_2c_bump_forces(force, t_2c_work(2), iatom, atom_of_kind, kind_of, -pref, &
1287 ri_data, qs_env, work_virial)
1290 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_r, &
1291 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1292 contract_1=[1], notcontract_1=[2], &
1293 contract_2=[1], notcontract_2=[2], &
1294 filter_eps=ri_data%filter_eps, flop=nflop)
1295 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1297 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_metric(jatom), &
1298 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1299 contract_1=[2], notcontract_1=[1], &
1300 contract_2=[1], notcontract_2=[2], &
1301 filter_eps=ri_data%filter_eps, flop=nflop)
1302 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1306 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(jatom), &
1307 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1308 contract_1=[2], notcontract_1=[1], &
1309 contract_2=[1], notcontract_2=[2], &
1310 filter_eps=ri_data%filter_eps, flop=nflop)
1311 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1313 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(3), &
1314 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1315 contract_1=[1], notcontract_1=[2], &
1316 contract_2=[1], notcontract_2=[2], &
1317 filter_eps=ri_data%filter_eps, flop=nflop)
1318 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1320 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1321 CALL get_2c_bump_forces(force, t_2c_work(4), jatom, atom_of_kind, kind_of, pref, &
1322 ri_data, qs_env, work_virial)
1325 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(2), &
1326 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1327 contract_1=[1], notcontract_1=[2], &
1328 contract_2=[1], notcontract_2=[2], &
1329 filter_eps=ri_data%filter_eps, flop=nflop)
1330 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1332 IF (use_virial_prv)
THEN
1333 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1334 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1335 diag=.true., offdiag=.false.)
1337 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1338 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1342 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1343 CALL apply_bump(t_2c_work(2), jatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1345 IF (use_virial_prv)
THEN
1346 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1347 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1348 diag=.false., offdiag=.true.)
1350 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1351 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1356 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(jatom), &
1357 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1358 contract_1=[2], notcontract_1=[1], &
1359 contract_2=[1], notcontract_2=[2], &
1360 filter_eps=ri_data%filter_eps, flop=nflop)
1361 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1363 CALL dbt_contract(1.0_dp, t_2c_bint(jatom), t_2c_work(4), &
1364 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1365 contract_1=[1], notcontract_1=[2], &
1366 contract_2=[1], notcontract_2=[2], &
1367 filter_eps=ri_data%filter_eps, flop=nflop)
1368 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1370 CALL get_2c_bump_forces(force, t_2c_work(2), jatom, atom_of_kind, kind_of, -pref, &
1371 ri_data, qs_env, work_virial)
1373 CALL timestop(handle2)
1378 IF (i_loop == 1)
THEN
1379 CALL update_pattern_to_forces(force_pattern, sparsity_pattern, ngroups, ri_data, qs_env)
1383 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
1384 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
1385 CALL dbt_batched_contract_finalize(t_3c_work_3(3))
1386 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
1387 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
1389 IF (use_virial_prv)
THEN
1393 virial%pv_fock_4c(i_xyz, j_xyz) = virial%pv_fock_4c(i_xyz, j_xyz) &
1394 + work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1402 CALL para_env_sub%free()
1403 DEALLOCATE (para_env_sub)
1405 CALL para_env%sync()
1407 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
1410 IF (
ASSOCIATED(dbcsr_template))
THEN
1412 DEALLOCATE (dbcsr_template)
1414 CALL dbt_destroy(t_2c_r)
1415 CALL dbt_destroy(t_2c_r_split)
1416 CALL dbt_destroy(t_2c_work(1))
1417 CALL dbt_destroy(t_2c_work(2))
1418 CALL dbt_destroy(t_2c_work(3))
1419 CALL dbt_destroy(t_2c_work(4))
1420 CALL dbt_destroy(t_2c_work(5))
1421 CALL dbt_destroy(t_3c_work_2(1))
1422 CALL dbt_destroy(t_3c_work_2(2))
1423 CALL dbt_destroy(t_3c_work_2(3))
1424 CALL dbt_destroy(t_3c_work_3(1))
1425 CALL dbt_destroy(t_3c_work_3(2))
1426 CALL dbt_destroy(t_3c_work_3(3))
1427 CALL dbt_destroy(t_3c_work_3(4))
1428 CALL dbt_destroy(t_3c_der_stack(1))
1429 CALL dbt_destroy(t_3c_der_stack(2))
1430 CALL dbt_destroy(t_3c_der_stack(3))
1431 CALL dbt_destroy(t_3c_der_stack(4))
1432 CALL dbt_destroy(t_3c_der_stack(5))
1433 CALL dbt_destroy(t_3c_der_stack(6))
1435 CALL dbt_destroy(t_2c_der_pot(i_xyz))
1438 CALL dbt_destroy(t_2c_inv(iatom))
1439 CALL dbt_destroy(t_2c_binv(iatom))
1440 CALL dbt_destroy(t_2c_bint(iatom))
1441 CALL dbt_destroy(t_2c_metric(iatom))
1443 CALL dbt_destroy(t_2c_der_metric_sub(iatom, i_xyz))
1448 DO i_spin = 1, nspins
1449 CALL dbt_destroy(rho_ao_t_sub(i_spin, i_img))
1450 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
1455 CALL dbt_destroy(t_3c_der_ri_sub(i_img, i_xyz))
1456 CALL dbt_destroy(t_3c_der_ao_sub(i_img, i_xyz))
1461 CALL timestop(handle)
1476 SUBROUTINE apply_bump(t_2c_inout, atom_i, ri_data, qs_env, from_left, from_right, debump)
1477 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_inout
1478 INTEGER,
INTENT(IN) :: atom_i
1481 LOGICAL,
INTENT(IN),
OPTIONAL :: from_left, from_right, debump
1483 INTEGER :: i_img, i_ri, iatom, ind(2), j_img, j_ri, &
1484 jatom, natom, nblks(2), nimg, nkind
1485 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1486 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1487 LOGICAL :: found, my_debump, my_left, my_right
1488 REAL(
dp) :: bval, r0, r1, ri(3), rj(3), rref(3), &
1490 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1492 TYPE(dbt_iterator_type) :: iter
1495 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1497 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1499 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1500 kpoints=kpoints, particle_set=particle_set)
1501 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1504 IF (
PRESENT(debump)) my_debump = debump
1507 IF (
PRESENT(from_left)) my_left = from_left
1510 IF (
PRESENT(from_right)) my_right = from_right
1511 cpassert(my_left .OR. my_right)
1513 CALL dbt_get_info(t_2c_inout, nblks_total=nblks)
1514 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1515 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1520 r1 = ri_data%kp_RI_range
1521 r0 = ri_data%kp_bump_rad
1522 rref =
pbc(particle_set(atom_i)%r, cell)
1527 CALL dbt_iterator_start(iter, t_2c_inout)
1528 DO WHILE (dbt_iterator_blocks_left(iter))
1529 CALL dbt_iterator_next_block(iter, ind)
1530 CALL dbt_get_block(t_2c_inout, ind, blk, found)
1531 IF (.NOT. found) cycle
1533 i_ri = (ind(1) - 1)/natom + 1
1534 i_img = ri_data%RI_cell_to_img(i_ri)
1535 iatom = ind(1) - (i_ri - 1)*natom
1538 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
1540 j_ri = (ind(2) - 1)/natom + 1
1541 j_img = ri_data%RI_cell_to_img(j_ri)
1542 jatom = ind(2) - (j_ri - 1)*natom
1545 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1547 IF (.NOT. my_debump)
THEN
1548 IF (my_left) blk(:, :) = blk(:, :)*bump(norm2(ri - rref), r0, r1)
1549 IF (my_right) blk(:, :) = blk(:, :)*bump(norm2(rj - rref), r0, r1)
1553 bval = bump(norm2(ri - rref), r0, r1)
1554 IF (my_left .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1555 bval = bump(norm2(rj - rref), r0, r1)
1556 IF (my_right .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1559 CALL dbt_put_block(t_2c_inout, ind, shape(blk), blk)
1563 CALL dbt_iterator_stop(iter)
1565 CALL dbt_filter(t_2c_inout, ri_data%filter_eps)
1567 END SUBROUTINE apply_bump
1581 SUBROUTINE get_2c_bump_forces(force, t_2c_in, atom_i, atom_of_kind, kind_of, pref, ri_data, &
1582 qs_env, work_virial)
1584 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_in
1585 INTEGER,
INTENT(IN) :: atom_i
1586 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
1587 REAL(
dp),
INTENT(IN) :: pref
1590 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT) :: work_virial
1592 INTEGER :: i, i_img, i_ri, i_xyz, iat_of_kind, iatom, ikind, ind(2), j_img, j_ri, j_xyz, &
1593 jat_of_kind, jatom, jkind, natom, nblks(2), nimg, nkind
1594 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1595 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1597 REAL(
dp) :: new_force, r0, r1, ri(3), rj(3), &
1598 rref(3), scoord(3), x
1599 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1601 TYPE(dbt_iterator_type) :: iter
1604 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1606 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1608 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1609 kpoints=kpoints, particle_set=particle_set)
1610 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1612 CALL dbt_get_info(t_2c_in, nblks_total=nblks)
1613 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1614 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1619 r1 = ri_data%kp_RI_range
1620 r0 = ri_data%kp_bump_rad
1621 rref =
pbc(particle_set(atom_i)%r, cell)
1623 iat_of_kind = atom_of_kind(atom_i)
1624 ikind = kind_of(atom_i)
1630 CALL dbt_iterator_start(iter, t_2c_in)
1631 DO WHILE (dbt_iterator_blocks_left(iter))
1632 CALL dbt_iterator_next_block(iter, ind)
1633 IF (ind(1) /= ind(2)) cycle
1635 CALL dbt_get_block(t_2c_in, ind, blk, found)
1636 IF (.NOT. found) cycle
1639 j_ri = (ind(2) - 1)/natom + 1
1640 j_img = ri_data%RI_cell_to_img(j_ri)
1641 jatom = ind(2) - (j_ri - 1)*natom
1642 jat_of_kind = atom_of_kind(jatom)
1643 jkind = kind_of(jatom)
1646 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1647 x = norm2(rj - rref)
1648 IF (x < r0 .OR. x > r1) cycle
1651 DO i = 1,
SIZE(blk, 1)
1652 new_force = new_force + blk(i, i)
1654 new_force = pref*new_force*dbump(x, r0, r1)
1660 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) + &
1661 new_force*(rj(i_xyz) - rref(i_xyz))/x
1667 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1668 + new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1673 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) - &
1674 new_force*(rj(i_xyz) - rref(i_xyz))/x
1680 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1681 - new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1687 CALL dbt_iterator_stop(iter)
1690 END SUBROUTINE get_2c_bump_forces
1699 FUNCTION bump(x, r0, r1)
RESULT(b)
1700 REAL(
dp),
INTENT(IN) :: x, r0, r1
1708 r = (x - r0)/(r1 - r0)
1709 b = -6.0_dp*r**5 + 15.0_dp*r**4 - 10.0_dp*r**3 + 1.0_dp
1710 IF (x >= r1) b = 0.0_dp
1711 IF (x <= r0) b = 1.0_dp
1722 FUNCTION dbump(x, r0, r1)
RESULT(b)
1723 REAL(
dp),
INTENT(IN) :: x, r0, r1
1728 r = (x - r0)/(r1 - r0)
1729 b = (-30.0_dp*r**4 + 60.0_dp*r**3 - 30.0_dp*r**2)/(r1 - r0)
1730 IF (x >= r1) b = 0.0_dp
1731 IF (x <= r0) b = 0.0_dp
1742 FUNCTION get_apc_index_from_ib(i_index, b_index, qs_env)
RESULT(apc_index)
1743 INTEGER,
INTENT(IN) :: i_index, b_index
1745 INTEGER :: apc_index
1747 INTEGER,
DIMENSION(3) :: cell_apc
1748 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1749 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1753 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1756 cell_apc(:) = index_to_cell(:, i_index) + index_to_cell(:, b_index)
1758 IF (any([cell_apc(1), cell_apc(2), cell_apc(3)] < lbound(cell_to_index)) .OR. &
1759 any([cell_apc(1), cell_apc(2), cell_apc(3)] > ubound(cell_to_index)))
THEN
1763 apc_index = cell_to_index(cell_apc(1), cell_apc(2), cell_apc(3))
1766 END FUNCTION get_apc_index_from_ib
1775 FUNCTION get_apc_index(a_index, c_index, qs_env)
RESULT(i_index)
1776 INTEGER,
INTENT(IN) :: a_index, c_index
1780 INTEGER,
DIMENSION(3) :: cell_i
1781 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1782 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1786 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1788 cell_i(:) = index_to_cell(:, a_index) + index_to_cell(:, c_index)
1790 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1791 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1795 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1798 END FUNCTION get_apc_index
1807 FUNCTION get_i_index(apc_index, b_index, qs_env)
RESULT(i_index)
1808 INTEGER,
INTENT(IN) :: apc_index, b_index
1812 INTEGER,
DIMENSION(3) :: cell_i
1813 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1814 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1818 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1820 cell_i(:) = index_to_cell(:, apc_index) - index_to_cell(:, b_index)
1822 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1823 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1827 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1830 END FUNCTION get_i_index
1841 SUBROUTINE get_ac_pairs(ac_pairs, apc_index, ri_data, qs_env)
1842 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: ac_pairs
1843 INTEGER,
INTENT(IN) :: apc_index
1847 INTEGER :: a_index, actual_img, c_index, nimg
1849 nimg =
SIZE(ac_pairs, 1)
1854 DO a_index = 1, nimg
1855 actual_img = ri_data%idx_to_img(a_index)
1857 c_index = get_i_index(apc_index, actual_img, qs_env)
1858 ac_pairs(a_index, 1) = a_index
1859 ac_pairs(a_index, 2) = c_index
1863 END SUBROUTINE get_ac_pairs
1875 SUBROUTINE get_iapc_pairs(iapc_pairs, b_index, ri_data, qs_env, actual_i_img)
1876 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: iapc_pairs
1877 INTEGER,
INTENT(IN) :: b_index
1880 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: actual_i_img
1882 INTEGER :: actual_img, apc_index, i_index, nimg
1884 nimg =
SIZE(iapc_pairs, 1)
1885 IF (
PRESENT(actual_i_img)) actual_i_img(:) = 0
1887 iapc_pairs(:, :) = 0
1890 DO i_index = 1, nimg
1891 actual_img = ri_data%idx_to_img(i_index)
1892 apc_index = get_apc_index_from_ib(actual_img, b_index, qs_env)
1893 IF (apc_index == 0) cycle
1894 iapc_pairs(i_index, 1) = i_index
1895 iapc_pairs(i_index, 2) = apc_index
1896 IF (
PRESENT(actual_i_img)) actual_i_img(i_index) = actual_img
1899 END SUBROUTINE get_iapc_pairs
1908 FUNCTION get_opp_index(a_index, qs_env)
RESULT(opp_index)
1909 INTEGER,
INTENT(IN) :: a_index
1911 INTEGER :: opp_index
1913 INTEGER,
DIMENSION(3) :: opp_cell
1914 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1915 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1918 NULLIFY (kpoints, cell_to_index, index_to_cell)
1921 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1923 opp_cell(:) = -index_to_cell(:, a_index)
1925 IF (any([opp_cell(1), opp_cell(2), opp_cell(3)] < lbound(cell_to_index)) .OR. &
1926 any([opp_cell(1), opp_cell(2), opp_cell(3)] > ubound(cell_to_index)))
THEN
1930 opp_index = cell_to_index(opp_cell(1), opp_cell(2), opp_cell(3))
1933 END FUNCTION get_opp_index
1944 SUBROUTINE get_pmat_images(rho_ao_t, rho_ao, scale_prev_p, ri_data, qs_env)
1945 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t
1946 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao
1947 REAL(
dp),
INTENT(IN) :: scale_prev_p
1951 INTEGER :: cell_j(3), i_img, i_spin, iatom, icol, &
1952 irow, j_img, jatom, mi_img, mj_img, &
1954 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1957 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock, pblock_desymm
1958 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, rho_desymm
1959 TYPE(dbt_type) :: tmp
1963 DIMENSION(:),
POINTER :: nl_iterator
1965 POINTER :: sab_nl, sab_nl_nosym
1968 NULLIFY (rho_desymm, kpoints, sab_nl_nosym, scf_env, matrix_ks, dft_control, &
1969 sab_nl, nl_iterator, cell_to_index, pblock, pblock_desymm)
1971 CALL get_qs_env(qs_env, kpoints=kpoints, scf_env=scf_env, matrix_ks_kp=matrix_ks, dft_control=dft_control)
1972 CALL get_kpoint_info(kpoints, sab_nl_nosym=sab_nl_nosym, cell_to_index=cell_to_index, sab_nl=sab_nl)
1974 IF (dft_control%do_admm)
THEN
1975 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks)
1978 nspins =
SIZE(matrix_ks, 1)
1981 ALLOCATE (rho_desymm(nspins, nimg))
1983 DO i_spin = 1, nspins
1984 ALLOCATE (rho_desymm(i_spin, i_img)%matrix)
1985 CALL dbcsr_create(rho_desymm(i_spin, i_img)%matrix, template=matrix_ks(i_spin, i_img)%matrix, &
1986 matrix_type=dbcsr_type_no_symmetry)
1990 CALL dbt_create(rho_desymm(1, 1)%matrix, tmp)
1997 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1998 IF (j_img > nimg .OR. j_img < 1) cycle
2001 IF (iatom == jatom)
fac = 0.5_dp
2002 mj_img = get_opp_index(j_img, qs_env)
2004 IF (mj_img == 0)
fac = 1.0_dp
2008 IF (iatom > jatom)
THEN
2014 DO i_spin = 1, nspins
2016 IF (.NOT. found) cycle
2019 CALL dbcsr_get_block_p(rho_desymm(i_spin, j_img)%matrix, iatom, jatom, pblock_desymm, found)
2020 IF (.NOT. found) cycle
2022 IF (iatom > jatom)
THEN
2023 pblock_desymm(:, :) =
fac*transpose(pblock(:, :))
2025 pblock_desymm(:, :) =
fac*pblock(:, :)
2032 DO i_spin = 1, nspins
2033 CALL dbt_scale(rho_ao_t(i_spin, i_img), scale_prev_p)
2035 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, i_img)%matrix, tmp)
2036 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), summation=.true., move_data=.true.)
2039 mi_img = get_opp_index(i_img, qs_env)
2040 IF (mi_img > 0 .AND. mi_img <= nimg)
THEN
2041 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, mi_img)%matrix, tmp)
2042 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), order=[2, 1], summation=.true., move_data=.true.)
2044 CALL dbt_filter(rho_ao_t(i_spin, i_img), ri_data%filter_eps)
2049 DO i_spin = 1, nspins
2051 DEALLOCATE (rho_desymm(i_spin, i_img)%matrix)
2055 CALL dbt_destroy(tmp)
2056 DEALLOCATE (rho_desymm)
2058 END SUBROUTINE get_pmat_images
2077 SUBROUTINE get_ext_2c_int(t_2c_pot, mat_orig, atom_i, atom_j, img_b, ri_data, qs_env, do_inverse, &
2078 para_env_ext, blacs_env_ext, dbcsr_template, off_diagonal, skip_inverse)
2079 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_pot
2080 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_orig
2081 INTEGER,
INTENT(IN) :: atom_i, atom_j, img_b
2084 LOGICAL,
INTENT(IN),
OPTIONAL :: do_inverse
2087 TYPE(
dbcsr_type),
OPTIONAL,
POINTER :: dbcsr_template
2088 LOGICAL,
INTENT(IN),
OPTIONAL :: off_diagonal, skip_inverse
2090 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_ext_2c_int'
2092 INTEGER :: group, handle, handle2, i_img, i_ri, iatom, iblk, ikind, img_tot, j_img, j_ri, &
2093 jatom, jblk, jkind, n_dependent, natom, nblks_ri, nimg, nkind
2094 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2
2095 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: present_atoms_i, present_atoms_j
2096 INTEGER,
DIMENSION(3) :: cell_b, cell_i, cell_j, cell_tot
2097 INTEGER,
DIMENSION(:),
POINTER :: col_dist, col_dist_ext, ri_blk_size_ext, &
2098 row_dist, row_dist_ext
2099 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell, pgrid
2100 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
2101 LOGICAL :: do_inverse_prv, found, my_offd, &
2102 skip_inverse_prv, use_template
2103 REAL(
dp) :: bfac, dij, r0, r1, threshold
2104 REAL(
dp),
DIMENSION(3) :: ri, rij, rj, rref, scoord
2105 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock
2110 TYPE(
dbcsr_type) :: work, work_tight, work_tight_inv
2111 TYPE(dbt_type) :: t_2c_tmp
2114 DIMENSION(:),
TARGET :: basis_set_ri
2118 DIMENSION(:),
POINTER :: nl_iterator
2122 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2124 NULLIFY (qs_kind_set, nl_2c, nl_iterator, cell, kpoints, cell_to_index, index_to_cell, dist_2d, &
2125 para_env, pblock, blacs_env, particle_set, col_dist, row_dist, pgrid, &
2126 col_dist_ext, row_dist_ext)
2128 CALL timeset(routinen, handle)
2133 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
2134 kpoints=kpoints, para_env=para_env, blacs_env=blacs_env, particle_set=particle_set)
2135 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
2137 do_inverse_prv = .false.
2138 IF (
PRESENT(do_inverse)) do_inverse_prv = do_inverse
2139 IF (do_inverse_prv)
THEN
2140 cpassert(atom_i == atom_j)
2143 skip_inverse_prv = .false.
2144 IF (
PRESENT(skip_inverse)) skip_inverse_prv = skip_inverse
2147 IF (
PRESENT(off_diagonal)) my_offd = off_diagonal
2149 IF (
PRESENT(para_env_ext)) para_env => para_env_ext
2150 IF (
PRESENT(blacs_env_ext)) blacs_env => blacs_env_ext
2152 nimg =
SIZE(mat_orig)
2154 CALL timeset(routinen//
"_nl_iter", handle2)
2157 ALLOCATE (dist1(natom), dist2(natom))
2159 dist1(iatom) = mod(iatom, blacs_env%num_pe(1))
2160 dist2(iatom) = mod(iatom, blacs_env%num_pe(2))
2164 ALLOCATE (basis_set_ri(nkind))
2168 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
2170 ALLOCATE (present_atoms_i(natom, nimg), present_atoms_j(natom, nimg))
2176 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij, cell=cell_j, &
2177 ikind=ikind, jkind=jkind)
2181 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2182 IF (j_img > nimg .OR. j_img < 1) cycle
2184 IF (iatom == atom_i .AND. dij <= ri_data%kp_RI_range) present_atoms_i(jatom, j_img) = 1
2185 IF (iatom == atom_j .AND. dij <= ri_data%kp_RI_range) present_atoms_j(jatom, j_img) = 1
2190 CALL timestop(handle2)
2192 CALL para_env%sum(present_atoms_i)
2193 CALL para_env%sum(present_atoms_j)
2197 use_template = .false.
2198 IF (
PRESENT(dbcsr_template))
THEN
2199 IF (
ASSOCIATED(dbcsr_template)) use_template = .true.
2202 IF (use_template)
THEN
2207 ALLOCATE (row_dist_ext(ri_data%ncell_RI*natom), col_dist_ext(ri_data%ncell_RI*natom))
2208 ALLOCATE (ri_blk_size_ext(ri_data%ncell_RI*natom))
2209 DO i_ri = 1, ri_data%ncell_RI
2210 row_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = row_dist(:)
2211 col_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = col_dist(:)
2212 ri_blk_size_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2216 row_dist=row_dist_ext, col_dist=col_dist_ext)
2217 CALL dbcsr_create(work, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2218 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2220 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2222 IF (
PRESENT(dbcsr_template))
THEN
2223 ALLOCATE (dbcsr_template)
2228 cell_b(:) = index_to_cell(:, img_b)
2230 i_ri = ri_data%img_to_RI_cell(i_img)
2231 IF (i_ri == 0) cycle
2232 cell_i(:) = index_to_cell(:, i_img)
2234 j_ri = ri_data%img_to_RI_cell(j_img)
2235 IF (j_ri == 0) cycle
2236 cell_j(:) = index_to_cell(:, j_img)
2237 cell_tot = cell_j - cell_i + cell_b
2239 IF (any([cell_tot(1), cell_tot(2), cell_tot(3)] < lbound(cell_to_index)) .OR. &
2240 any([cell_tot(1), cell_tot(2), cell_tot(3)] > ubound(cell_to_index))) cycle
2241 img_tot = cell_to_index(cell_tot(1), cell_tot(2), cell_tot(3))
2242 IF (img_tot > nimg .OR. img_tot < 1) cycle
2247 IF (present_atoms_i(iatom, i_img) == 0) cycle
2248 IF (present_atoms_j(jatom, j_img) == 0) cycle
2249 IF (my_offd .AND. (i_ri - 1)*natom + iatom == (j_ri - 1)*natom + jatom) cycle
2252 IF (.NOT. found) cycle
2254 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2263 IF (do_inverse_prv)
THEN
2265 r1 = ri_data%kp_RI_range
2266 r0 = ri_data%kp_bump_rad
2269 nblks_ri = sum(present_atoms_i)
2270 ALLOCATE (col_dist_ext(nblks_ri), row_dist_ext(nblks_ri), ri_blk_size_ext(nblks_ri))
2273 i_ri = ri_data%img_to_RI_cell(i_img)
2274 IF (i_ri == 0) cycle
2276 IF (present_atoms_i(iatom, i_img) == 0) cycle
2278 col_dist_ext(iblk) = col_dist(iatom)
2279 row_dist_ext(iblk) = row_dist(iatom)
2280 ri_blk_size_ext(iblk) = ri_data%bsizes_RI(iatom)
2285 row_dist=row_dist_ext, col_dist=col_dist_ext)
2286 CALL dbcsr_create(work_tight, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2287 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2288 CALL dbcsr_create(work_tight_inv, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2289 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2291 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2295 rref =
pbc(particle_set(atom_i)%r, cell)
2299 i_ri = ri_data%img_to_RI_cell(i_img)
2300 IF (i_ri == 0) cycle
2302 IF (present_atoms_i(iatom, i_img) == 0) cycle
2306 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
2310 j_ri = ri_data%img_to_RI_cell(j_img)
2311 IF (j_ri == 0) cycle
2313 IF (present_atoms_j(jatom, j_img) == 0) cycle
2317 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
2319 CALL dbcsr_get_block_p(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock, found)
2320 IF (.NOT. found) cycle
2323 IF (iblk /= jblk) bfac = bump(norm2(ri - rref), r0, r1)*bump(norm2(rj - rref), r0, r1)
2332 IF (.NOT. skip_inverse_prv)
THEN
2333 SELECT CASE (ri_data%t2c_method)
2335 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
2336 CALL invert_hotelling(work_tight_inv, work_tight, threshold=threshold, silent=.false.)
2341 uplo_to_full=.true.)
2344 CALL cp_dbcsr_power(work_tight_inv, -1.0_dp, ri_data%eps_eigval, n_dependent, &
2345 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
2356 i_ri = ri_data%img_to_RI_cell(i_img)
2357 IF (i_ri == 0) cycle
2359 IF (present_atoms_i(iatom, i_img) == 0) cycle
2364 j_ri = ri_data%img_to_RI_cell(j_img)
2365 IF (j_ri == 0) cycle
2367 IF (present_atoms_j(jatom, j_img) == 0) cycle
2371 IF (.NOT. found) cycle
2373 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2384 CALL dbt_create(work, t_2c_tmp)
2385 CALL dbt_copy_matrix_to_tensor(work, t_2c_tmp)
2386 CALL dbt_copy(t_2c_tmp, t_2c_pot, move_data=.true.)
2387 CALL dbt_filter(t_2c_pot, ri_data%filter_eps)
2389 CALL dbt_destroy(t_2c_tmp)
2392 CALL timestop(handle)
2394 END SUBROUTINE get_ext_2c_int
2404 SUBROUTINE contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
2405 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, rho_ao_t
2409 CHARACTER(len=*),
PARAMETER :: routinen =
'contract_pmat_3c'
2411 INTEGER :: apc_img, b_img, batch_size, handle, &
2412 i_batch, i_img, i_spin,
idx, j_batch, &
2413 n_batch_img, n_batch_nze, nimg, &
2415 INTEGER(int_8) :: nflop, nze
2416 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: apc_filter, batch_ranges_img, &
2417 batch_ranges_nze, int_indices
2418 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: ac_pairs, iapc_pairs
2419 REAL(
dp) :: occ, t1, t2
2420 TYPE(dbt_type) :: t_3c_tmp
2421 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ints_stack, res_stack, rho_stack
2424 CALL timeset(routinen, handle)
2426 CALL get_qs_env(qs_env, dft_control=dft_control)
2429 nimg_nze = ri_data%nimg_nze
2430 nspins = dft_control%nspins
2432 CALL dbt_create(t_3c_apc(1, 1), t_3c_tmp)
2434 batch_size = ri_data%kp_stack_size
2436 ALLOCATE (apc_filter(nimg), iapc_pairs(nimg, 2))
2439 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env)
2440 DO i_img = 1, nimg_nze
2441 idx = iapc_pairs(i_img, 2)
2442 IF (
idx < 1 .OR.
idx > nimg) cycle
2448 n_batch_img = nimg/batch_size
2449 IF (
modulo(nimg, batch_size) /= 0) n_batch_img = n_batch_img + 1
2450 ALLOCATE (batch_ranges_img(n_batch_img + 1))
2451 DO i_batch = 1, n_batch_img
2452 batch_ranges_img(i_batch) = (i_batch - 1)*batch_size + 1
2454 batch_ranges_img(n_batch_img + 1) = nimg + 1
2457 n_batch_nze = nimg_nze/batch_size
2458 IF (
modulo(nimg_nze, batch_size) /= 0) n_batch_nze = n_batch_nze + 1
2459 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
2460 DO i_batch = 1, n_batch_nze
2461 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
2463 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
2466 ALLOCATE (rho_stack(2), ints_stack(2), res_stack(2))
2467 CALL get_stack_tensors(res_stack, rho_stack, ints_stack, rho_ao_t(1, 1), &
2468 ri_data%t_3c_int_ctr_1(1, 1), batch_size, ri_data, qs_env)
2470 ALLOCATE (ac_pairs(nimg, 2), int_indices(nimg_nze))
2471 DO i_img = 1, nimg_nze
2472 int_indices(i_img) = i_img
2476 DO j_batch = 1, n_batch_nze
2478 CALL fill_3c_stack(ints_stack(1), ri_data%t_3c_int_ctr_1(1, :), int_indices, 3, ri_data, &
2479 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)])
2480 CALL dbt_copy(ints_stack(1), ints_stack(2), move_data=.true.)
2482 DO i_spin = 1, nspins
2483 DO i_batch = 1, n_batch_img
2485 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2486 IF (apc_filter(apc_img) == 0) cycle
2487 CALL get_ac_pairs(ac_pairs, apc_img, ri_data, qs_env)
2488 CALL fill_2c_stack(rho_stack(1), rho_ao_t(i_spin, :), ac_pairs(:, 2), 1, ri_data, &
2489 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)], &
2490 shift=apc_img - batch_ranges_img(i_batch) + 1)
2495 CALL dbt_copy(rho_stack(1), rho_stack(2), move_data=.true.)
2498 CALL dbt_batched_contract_init(rho_stack(2))
2499 CALL dbt_contract(1.0_dp, ints_stack(2), rho_stack(2), &
2500 0.0_dp, res_stack(2), map_1=[1, 2], map_2=[3], &
2501 contract_1=[3], notcontract_1=[1, 2], &
2502 contract_2=[1], notcontract_2=[2], &
2503 filter_eps=ri_data%filter_eps, flop=nflop)
2504 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2505 CALL dbt_batched_contract_finalize(rho_stack(2))
2506 CALL dbt_copy(res_stack(2), res_stack(1), move_data=.true.)
2508 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2510 IF (apc_filter(apc_img) == 0) cycle
2511 CALL unstack_t_3c_apc(t_3c_tmp, res_stack(1), apc_img - batch_ranges_img(i_batch) + 1)
2512 CALL dbt_copy(t_3c_tmp, t_3c_apc(i_spin, apc_img), summation=.true., move_data=.true.)
2518 DEALLOCATE (batch_ranges_img)
2519 DEALLOCATE (batch_ranges_nze)
2521 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
2523 CALL dbt_destroy(rho_stack(1))
2524 CALL dbt_destroy(rho_stack(2))
2525 CALL dbt_destroy(ints_stack(1))
2526 CALL dbt_destroy(ints_stack(2))
2527 CALL dbt_destroy(res_stack(1))
2528 CALL dbt_destroy(res_stack(2))
2529 CALL dbt_destroy(t_3c_tmp)
2531 CALL timestop(handle)
2533 END SUBROUTINE contract_pmat_3c
2541 SUBROUTINE precontract_3c_ints(t_3c_int, ri_data, qs_env)
2542 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_int
2546 CHARACTER(len=*),
PARAMETER :: routinen =
'precontract_3c_ints'
2548 INTEGER :: batch_size, handle, i_batch, i_img, &
2549 i_ri, iatom, is, n_batch, natom, &
2550 nblks, nblks_3c(3), nimg
2551 INTEGER(int_8) :: nflop
2552 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges, bsizes_ri_ext, bsizes_ri_ext_split, &
2553 bsizes_stack, dist1, dist2, dist3, dist_stack3, idx_to_at_ao, int_indices
2554 TYPE(dbt_distribution_type) :: t_dist
2555 TYPE(dbt_type) :: t_2c_ri_tmp(2), t_3c_tmp(3)
2557 CALL timeset(routinen, handle)
2562 ALLOCATE (int_indices(nimg))
2564 int_indices(i_img) = i_img
2567 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
2568 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
2570 nblks =
SIZE(ri_data%bsizes_RI_split)
2571 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
2572 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
2573 DO i_ri = 1, ri_data%ncell_RI
2574 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2575 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
2578 bsizes_ri_ext, bsizes_ri_ext, &
2580 DEALLOCATE (dist1, dist2)
2582 bsizes_ri_ext_split, bsizes_ri_ext_split, &
2584 DEALLOCATE (dist1, dist2)
2587 batch_size = ri_data%kp_stack_size
2588 n_batch = nimg/batch_size
2589 IF (
modulo(nimg, batch_size) /= 0) n_batch = n_batch + 1
2590 ALLOCATE (batch_ranges(n_batch + 1))
2591 DO i_batch = 1, n_batch
2592 batch_ranges(i_batch) = (i_batch - 1)*batch_size + 1
2594 batch_ranges(n_batch + 1) = nimg + 1
2596 nblks =
SIZE(ri_data%bsizes_AO_split)
2597 ALLOCATE (bsizes_stack(batch_size*nblks))
2598 DO is = 1, batch_size
2599 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
2602 CALL dbt_get_info(t_3c_int(1, 1), nblks_total=nblks_3c)
2603 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)), dist_stack3(batch_size*nblks_3c(3)))
2604 CALL dbt_get_info(t_3c_int(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, proc_dist_3=dist3)
2605 DO is = 1, batch_size
2606 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
2609 CALL dbt_distribution_new(t_dist, ri_data%pgrid, dist1, dist2, dist_stack3)
2610 CALL dbt_create(t_3c_tmp(1),
"ints_stack", t_dist, [1], [2, 3], bsizes_ri_ext_split, &
2611 ri_data%bsizes_AO_split, bsizes_stack)
2612 CALL dbt_distribution_destroy(t_dist)
2613 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
2615 CALL dbt_create(t_3c_tmp(1), t_3c_tmp(2))
2616 CALL dbt_create(t_3c_int(1, 1), t_3c_tmp(3))
2619 CALL dbt_copy(ri_data%t_2c_inv(1, iatom), t_2c_ri_tmp(1))
2620 CALL apply_bump(t_2c_ri_tmp(1), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
2621 CALL dbt_copy(t_2c_ri_tmp(1), t_2c_ri_tmp(2), move_data=.true.)
2623 CALL dbt_batched_contract_init(t_2c_ri_tmp(2))
2624 DO i_batch = 1, n_batch
2626 CALL fill_3c_stack(t_3c_tmp(1), t_3c_int(1, :), int_indices, 3, ri_data, &
2627 img_bounds=[batch_ranges(i_batch), batch_ranges(i_batch + 1)], &
2628 filter_at=iatom, filter_dim=2, idx_to_at=idx_to_at_ao)
2630 CALL dbt_contract(1.0_dp, t_2c_ri_tmp(2), t_3c_tmp(1), &
2631 0.0_dp, t_3c_tmp(2), map_1=[1], map_2=[2, 3], &
2632 contract_1=[2], notcontract_1=[1], &
2633 contract_2=[1], notcontract_2=[2, 3], &
2634 filter_eps=ri_data%filter_eps, flop=nflop)
2635 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2637 DO i_img = batch_ranges(i_batch), batch_ranges(i_batch + 1) - 1
2638 CALL unstack_t_3c_apc(t_3c_tmp(3), t_3c_tmp(2), i_img - batch_ranges(i_batch) + 1)
2639 CALL dbt_copy(t_3c_tmp(3), ri_data%t_3c_int_ctr_1(1, i_img), summation=.true., &
2640 order=[2, 1, 3], move_data=.true.)
2642 CALL dbt_clear(t_3c_tmp(1))
2644 CALL dbt_batched_contract_finalize(t_2c_ri_tmp(2))
2647 CALL dbt_destroy(t_2c_ri_tmp(1))
2648 CALL dbt_destroy(t_2c_ri_tmp(2))
2649 CALL dbt_destroy(t_3c_tmp(1))
2650 CALL dbt_destroy(t_3c_tmp(2))
2651 CALL dbt_destroy(t_3c_tmp(3))
2654 CALL dbt_destroy(t_3c_int(1, i_img))
2657 CALL timestop(handle)
2659 END SUBROUTINE precontract_3c_ints
2670 SUBROUTINE copy_2c_to_subgroup(t2c_sub, t2c_main, group_size, ngroups, para_env)
2671 TYPE(dbt_type),
INTENT(INOUT) :: t2c_sub, t2c_main
2672 INTEGER,
INTENT(IN) :: group_size, ngroups
2675 INTEGER :: batch_size, i, i_batch, i_msg, iblk, &
2676 igroup, iproc, ir, is, jblk, n_batch, &
2678 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2
2679 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: block_dest, block_source
2680 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: current_dest
2681 INTEGER,
DIMENSION(2) :: ind, nblks
2683 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
2684 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2685 TYPE(dbt_iterator_type) :: iter
2686 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2692 CALL dbt_get_info(t2c_main, nblks_total=nblks)
2695 ALLOCATE (block_source(nblks(1), nblks(2)))
2699 CALL dbt_iterator_start(iter, t2c_main)
2700 DO WHILE (dbt_iterator_blocks_left(iter))
2701 CALL dbt_iterator_next_block(iter, ind)
2702 CALL dbt_get_block(t2c_main, ind, blk, found)
2703 IF (.NOT. found) cycle
2705 block_source(ind(1), ind(2)) = para_env%mepos
2710 CALL dbt_iterator_stop(iter)
2713 CALL para_env%sum(nocc)
2714 CALL para_env%sum(block_source)
2715 block_source = block_source + para_env%num_pe - 1
2716 IF (nocc == 0)
RETURN
2719 igroup = para_env%mepos/group_size
2720 ALLOCATE (block_dest(nblks(1), nblks(2)))
2722 DO jblk = 1, nblks(2)
2723 DO iblk = 1, nblks(1)
2724 IF (block_source(iblk, jblk) == -1) cycle
2726 CALL dbt_get_stored_coordinates(t2c_sub, [iblk, jblk], iproc)
2727 block_dest(iblk, jblk) = igroup*group_size + iproc
2731 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)))
2732 CALL dbt_get_info(t2c_main, blk_size_1=bsizes1, blk_size_2=bsizes2)
2734 ALLOCATE (current_dest(nblks(1), nblks(2), 0:ngroups - 1))
2735 DO igroup = 0, ngroups - 1
2737 current_dest(:, :, igroup) = block_dest(:, :)
2738 CALL para_env%bcast(current_dest(:, :, igroup), source=igroup*group_size)
2742 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2743 n_batch = (nocc*ngroups)/batch_size
2744 IF (
modulo(nocc*ngroups, batch_size) /= 0) n_batch = n_batch + 1
2746 DO i_batch = 1, n_batch
2748 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2749 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2753 DO jblk = 1, nblks(2)
2754 DO iblk = 1, nblks(1)
2755 DO igroup = 0, ngroups - 1
2756 IF (block_source(iblk, jblk) == -1) cycle
2759 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2762 tag = i_msg - (i_batch - 1)*batch_size
2765 IF (para_env%mepos == block_source(iblk, jblk))
THEN
2766 CALL dbt_get_block(t2c_main, [iblk, jblk], blk, found)
2770 IF (block_source(iblk, jblk) == current_dest(iblk, jblk, igroup))
THEN
2771 IF (found)
CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2773 IF (para_env%mepos == block_source(iblk, jblk) .AND. found)
THEN
2774 ALLOCATE (send_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2775 send_buff(tag)%array(:, :) = blk(:, :)
2777 CALL para_env%isend(msgin=send_buff(tag)%array, dest=current_dest(iblk, jblk, igroup), &
2778 request=send_req(is), tag=tag)
2781 IF (para_env%mepos == current_dest(iblk, jblk, igroup))
THEN
2782 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2784 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk), &
2785 request=recv_req(ir), tag=tag)
2789 IF (found)
DEALLOCATE (blk)
2797 DO i = 1, batch_size
2798 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
2803 DO jblk = 1, nblks(2)
2804 DO iblk = 1, nblks(1)
2805 DO igroup = 0, ngroups - 1
2806 IF (block_source(iblk, jblk) == -1) cycle
2809 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2812 tag = i_msg - (i_batch - 1)*batch_size
2814 IF (para_env%mepos == current_dest(iblk, jblk, igroup) .AND. &
2815 block_source(iblk, jblk) /= current_dest(iblk, jblk, igroup))
THEN
2817 ALLOCATE (blk(bsizes1(iblk), bsizes2(jblk)))
2818 blk(:, :) = recv_buff(tag)%array(:, :)
2819 CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2827 DO i = 1, batch_size
2828 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
2830 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
2832 CALL dbt_finalize(t2c_sub)
2834 END SUBROUTINE copy_2c_to_subgroup
2845 SUBROUTINE get_3c_subgroup_dest(subgroup_dest, t3c_sub, t3c_main, group_size, ngroups, para_env)
2846 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :), &
2847 INTENT(INOUT) :: subgroup_dest
2848 TYPE(dbt_type),
INTENT(INOUT) :: t3c_sub, t3c_main
2849 INTEGER,
INTENT(IN) :: group_size, ngroups
2852 INTEGER :: iblk, igroup, iproc, jblk, kblk
2853 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: block_dest
2854 INTEGER,
DIMENSION(3) :: nblks
2856 CALL dbt_get_info(t3c_main, nblks_total=nblks)
2859 igroup = para_env%mepos/group_size
2860 ALLOCATE (block_dest(nblks(1), nblks(2), nblks(3)))
2861 DO kblk = 1, nblks(3)
2862 DO jblk = 1, nblks(2)
2863 DO iblk = 1, nblks(1)
2864 CALL dbt_get_stored_coordinates(t3c_sub, [iblk, jblk, kblk], iproc)
2865 block_dest(iblk, jblk, kblk) = igroup*group_size + iproc
2870 ALLOCATE (subgroup_dest(nblks(1), nblks(2), nblks(3), ngroups))
2871 DO igroup = 0, ngroups - 1
2873 subgroup_dest(:, :, :, igroup + 1) = block_dest(:, :, :)
2874 CALL para_env%bcast(subgroup_dest(:, :, :, igroup + 1), source=igroup*group_size)
2877 END SUBROUTINE get_3c_subgroup_dest
2891 SUBROUTINE copy_3c_to_subgroup(t3c_sub, t3c_main, ngroups, para_env, subgroup_dest, &
2892 iatom_to_subgroup, dim_at, idx_to_at)
2893 TYPE(dbt_type),
INTENT(INOUT) :: t3c_sub, t3c_main
2894 INTEGER,
INTENT(IN) :: ngroups
2896 INTEGER,
DIMENSION(:, :, :, :),
INTENT(IN) :: subgroup_dest
2898 INTENT(INOUT),
OPTIONAL :: iatom_to_subgroup
2899 INTEGER,
INTENT(IN),
OPTIONAL :: dim_at
2900 INTEGER,
DIMENSION(:),
OPTIONAL :: idx_to_at
2902 INTEGER :: batch_size, i, i_batch, i_msg, iatom, &
2903 iblk, igroup, ir, is, isbuff, jblk, &
2904 kblk, n_batch, nocc, tag
2905 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2, bsizes3
2906 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: block_source
2907 INTEGER,
DIMENSION(3) :: ind, nblks
2908 LOGICAL :: filter_at, found
2909 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
2910 TYPE(
cp_3d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2911 TYPE(dbt_iterator_type) :: iter
2912 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2918 CALL dbt_get_info(t3c_main, nblks_total=nblks)
2922 IF (
PRESENT(iatom_to_subgroup) .AND.
PRESENT(dim_at) .AND.
PRESENT(idx_to_at))
THEN
2924 cpassert(nblks(dim_at) ==
SIZE(idx_to_at))
2928 ALLOCATE (block_source(nblks(1), nblks(2), nblks(3)))
2932 CALL dbt_iterator_start(iter, t3c_main)
2933 DO WHILE (dbt_iterator_blocks_left(iter))
2934 CALL dbt_iterator_next_block(iter, ind)
2935 CALL dbt_get_block(t3c_main, ind, blk, found)
2936 IF (.NOT. found) cycle
2938 block_source(ind(1), ind(2), ind(3)) = para_env%mepos
2943 CALL dbt_iterator_stop(iter)
2946 CALL para_env%sum(nocc)
2947 CALL para_env%sum(block_source)
2948 block_source = block_source + para_env%num_pe - 1
2949 IF (nocc == 0)
RETURN
2951 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)), bsizes3(nblks(3)))
2952 CALL dbt_get_info(t3c_main, blk_size_1=bsizes1, blk_size_2=bsizes2, blk_size_3=bsizes3)
2955 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2956 n_batch = (nocc*ngroups)/batch_size
2957 IF (
modulo(nocc*ngroups, batch_size) /= 0) n_batch = n_batch + 1
2959 DO i_batch = 1, n_batch
2961 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2962 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2967 DO kblk = 1, nblks(3)
2968 DO jblk = 1, nblks(2)
2969 DO iblk = 1, nblks(1)
2970 IF (block_source(iblk, jblk, kblk) == -1) cycle
2973 IF (para_env%mepos == block_source(iblk, jblk, kblk))
THEN
2974 CALL dbt_get_block(t3c_main, [iblk, jblk, kblk], blk, found)
2977 ALLOCATE (send_buff(isbuff)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2981 DO igroup = 0, ngroups - 1
2984 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2987 tag = i_msg - (i_batch - 1)*batch_size
2990 ind(:) = [iblk, jblk, kblk]
2991 iatom = idx_to_at(ind(dim_at))
2992 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
2996 IF (block_source(iblk, jblk, kblk) == subgroup_dest(iblk, jblk, kblk, igroup + 1))
THEN
2997 IF (found)
CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(blk), blk)
2999 IF (para_env%mepos == block_source(iblk, jblk, kblk) .AND. found)
THEN
3000 send_buff(isbuff)%array(:, :, :) = blk(:, :, :)
3002 CALL para_env%isend(msgin=send_buff(isbuff)%array, &
3003 dest=subgroup_dest(iblk, jblk, kblk, igroup + 1), &
3004 request=send_req(is), tag=tag)
3007 IF (para_env%mepos == subgroup_dest(iblk, jblk, kblk, igroup + 1))
THEN
3008 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
3010 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk, kblk), &
3011 request=recv_req(ir), tag=tag)
3016 IF (found)
DEALLOCATE (blk)
3024 DO kblk = 1, nblks(3)
3025 DO jblk = 1, nblks(2)
3026 DO iblk = 1, nblks(1)
3027 DO igroup = 0, ngroups - 1
3028 IF (block_source(iblk, jblk, kblk) == -1) cycle
3031 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
3034 tag = i_msg - (i_batch - 1)*batch_size
3037 ind(:) = [iblk, jblk, kblk]
3038 iatom = idx_to_at(ind(dim_at))
3039 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
3042 IF (para_env%mepos == subgroup_dest(iblk, jblk, kblk, igroup + 1) .AND. &
3043 block_source(iblk, jblk, kblk) /= subgroup_dest(iblk, jblk, kblk, igroup + 1))
THEN
3047 CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(recv_buff(tag)%array), recv_buff(tag)%array)
3056 DO i = 1, batch_size
3057 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
3058 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
3060 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
3062 CALL dbt_finalize(t3c_sub)
3064 END SUBROUTINE copy_3c_to_subgroup
3076 SUBROUTINE gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
3077 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t, ks_t_sub
3078 INTEGER,
INTENT(IN) :: group_size
3079 INTEGER,
DIMENSION(:, :, :),
INTENT(IN) :: sparsity_pattern
3083 CHARACTER(len=*),
PARAMETER :: routinen =
'gather_ks_matrix'
3085 INTEGER :: b_img, dest, handle, i, i_spin, iatom, &
3086 igroup, ir, is, jatom, n_mess, natom, &
3087 nimg, nspins, source, tag
3089 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
3090 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
3091 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
3093 CALL timeset(routinen, handle)
3095 nimg =
SIZE(sparsity_pattern, 3)
3096 natom =
SIZE(sparsity_pattern, 2)
3097 nspins =
SIZE(ks_t, 1)
3101 DO i_spin = 1, nspins
3104 IF (sparsity_pattern(iatom, jatom, b_img) > -1) n_mess = n_mess + 1
3109 ALLOCATE (send_buff(n_mess), recv_buff(n_mess))
3110 ALLOCATE (send_req(n_mess), recv_req(n_mess))
3116 DO i_spin = 1, nspins
3119 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3124 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3125 CALL dbt_get_stored_coordinates(ks_t_sub(i_spin, b_img), [iatom, jatom], source)
3126 igroup = sparsity_pattern(iatom, jatom, b_img)
3127 source = source + igroup*group_size
3128 IF (para_env%mepos == source)
THEN
3129 CALL dbt_get_block(ks_t_sub(i_spin, b_img), [iatom, jatom], blk, found)
3130 IF (source == dest)
THEN
3131 IF (found)
CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3133 ALLOCATE (send_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3134 send_buff(n_mess)%array(:, :) = 0.0_dp
3136 send_buff(n_mess)%array(:, :) = blk(:, :)
3139 CALL para_env%isend(msgin=send_buff(n_mess)%array, dest=dest, &
3140 request=send_req(is), tag=tag)
3146 IF (para_env%mepos == dest .AND. source /= dest)
THEN
3147 ALLOCATE (recv_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3149 CALL para_env%irecv(msgout=recv_buff(n_mess)%array, source=source, &
3150 request=recv_req(ir), tag=tag)
3161 DO i_spin = 1, nspins
3164 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3167 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3168 IF (para_env%mepos == dest)
THEN
3169 IF (.NOT.
ASSOCIATED(recv_buff(n_mess)%array)) cycle
3170 ALLOCATE (blk(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3171 blk(:, :) = recv_buff(n_mess)%array(:, :)
3172 CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3181 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
3182 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
3184 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
3187 CALL timestop(handle)
3189 END SUBROUTINE gather_ks_matrix
3204 SUBROUTINE get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
3205 group_size, ngroups, para_env, para_env_sub, ri_data)
3206 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3207 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work, t_2c_ao_tmp, ks_t_split
3208 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t_sub
3209 INTEGER,
INTENT(IN) :: group_size, ngroups
3213 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_tensors'
3215 INTEGER :: handle, i, i_img, i_ri, i_spin, iproc, &
3216 j, natom, nblks, nimg, nspins
3217 INTEGER(int_8) :: nze
3218 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3220 INTEGER,
DIMENSION(2) :: pdims_2d
3221 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3222 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3225 TYPE(dbt_pgrid_type) :: pgrid_2d
3226 TYPE(dbt_type) :: work, work_sub
3228 CALL timeset(routinen, handle)
3232 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3234 natom =
SIZE(ri_data%bsizes_RI)
3235 nblks =
SIZE(ri_data%bsizes_RI_split)
3236 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3237 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3238 DO i_ri = 1, ri_data%ncell_RI
3239 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3240 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3245 bsizes_ri_ext, bsizes_ri_ext, &
3247 DEALLOCATE (dist1, dist2)
3250 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3252 DEALLOCATE (dist1, dist2)
3256 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3258 DEALLOCATE (dist1, dist2)
3259 CALL dbt_create(ks_t_split(1), ks_t_split(2))
3262 ri_data%bsizes_AO, ri_data%bsizes_AO, &
3264 DEALLOCATE (dist1, dist2)
3266 nspins =
SIZE(ks_t_sub, 1)
3267 nimg =
SIZE(ks_t_sub, 2)
3269 DO i_spin = 1, nspins
3270 CALL dbt_create(t_2c_ao_tmp(1), ks_t_sub(i_spin, i_img))
3277 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3279 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3281 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3283 DO i = 0, pdims_2d(1) - 1
3284 DO j = 0, pdims_2d(2) - 1
3285 dbcsr_pgrid(i, j) = iproc
3291 ALLOCATE (col_dist(natom), row_dist(natom))
3292 row_dist(:) = dist1(:)
3293 col_dist(:) = dist2(:)
3295 ALLOCATE (ri_blk_size(natom))
3296 ri_blk_size(:) = ri_data%bsizes_RI(:)
3299 row_dist=row_dist, col_dist=col_dist)
3300 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3301 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3304 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3305 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3309 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3310 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3311 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3312 CALL dbt_clear(work_sub)
3315 CALL dbt_destroy(work)
3316 CALL dbt_destroy(work_sub)
3317 CALL dbt_pgrid_destroy(pgrid_2d)
3319 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3320 CALL timestop(handle)
3322 END SUBROUTINE get_subgroup_2c_tensors
3337 SUBROUTINE get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
3338 group_size, ngroups, para_env, para_env_sub, ri_data)
3339 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_int, t_3c_work_2, t_3c_work_3
3340 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, t_3c_apc_sub
3341 INTEGER,
INTENT(IN) :: group_size, ngroups
3345 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_tensors'
3347 INTEGER :: batch_size, bo(2), handle, handle2, &
3348 i_blk, i_img, i_ri, i_spin, ib, natom, &
3349 nblks_ao, nblks_ri, nimg, nspins
3350 INTEGER(int_8) :: nze
3351 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3352 bsizes_stack, bsizes_tmp, dist1, &
3353 dist2, dist3, dist_stack, idx_to_at
3354 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :) :: subgroup_dest
3355 INTEGER,
DIMENSION(3) :: pdims
3357 TYPE(dbt_distribution_type) :: t_dist
3358 TYPE(dbt_pgrid_type) :: pgrid
3359 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3361 CALL timeset(routinen, handle)
3363 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3364 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3365 DO i_ri = 1, ri_data%ncell_RI
3366 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3370 natom =
SIZE(ri_data%bsizes_RI)
3372 ALLOCATE (bsizes_tmp(nblks_ri))
3373 DO i_blk = 1, nblks_ri
3374 bo =
get_limit(natom, nblks_ri, i_blk - 1)
3375 bsizes_tmp(i_blk) = sum(ri_data%bsizes_RI(bo(1):bo(2)))
3377 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3378 DO i_ri = 1, ri_data%ncell_RI
3379 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = bsizes_tmp(:)
3382 batch_size = ri_data%kp_stack_size
3383 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3384 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3385 DO ib = 1, batch_size
3386 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3390 natom =
SIZE(ri_data%bsizes_RI)
3392 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3393 tensor_dims=[
SIZE(bsizes_ri_ext_split), 1, batch_size*
SIZE(ri_data%bsizes_AO_split)])
3397 pgrid, bsizes_ri_ext_split, ri_data%bsizes_AO_split, &
3398 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(RI | AO AO)")
3399 nimg =
SIZE(t_3c_int)
3401 CALL dbt_create(t_3c_int(1), t_3c_int(i_img))
3405 ALLOCATE (dist_stack(batch_size*nblks_ao))
3406 DO ib = 1, batch_size
3407 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3410 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3411 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3412 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3413 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3414 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3415 CALL dbt_distribution_destroy(t_dist)
3416 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3420 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3421 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3422 DEALLOCATE (dist1, dist2, dist3)
3425 ri_data%pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3426 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3427 DEALLOCATE (dist1, dist2, dist3)
3429 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3430 group_size, ngroups, para_env)
3433 CALL timeset(routinen//
"_ints", handle2)
3434 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
3436 CALL dbt_copy(ri_data%kp_t_3c_int(i_img), t_3c_int(i_img), move_data=.true.)
3439 ALLOCATE (ri_data%kp_t_3c_int(nimg))
3441 CALL dbt_create(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img))
3444 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, i_img), work_atom_block, order=[2, 1, 3])
3445 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3446 ngroups, para_env, subgroup_dest)
3447 CALL dbt_copy(work_atom_block_sub, t_3c_int(i_img), move_data=.true.)
3448 CALL dbt_filter(t_3c_int(i_img), ri_data%filter_eps)
3451 CALL timestop(handle2)
3452 CALL dbt_pgrid_destroy(pgrid)
3453 CALL dbt_destroy(work_atom_block)
3454 CALL dbt_destroy(work_atom_block_sub)
3455 DEALLOCATE (subgroup_dest)
3459 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3460 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3464 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3465 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3466 DEALLOCATE (dist1, dist2, dist3)
3469 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3470 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3471 DEALLOCATE (dist1, dist2, dist3)
3473 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3474 group_size, ngroups, para_env)
3478 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3479 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3482 ALLOCATE (dist_stack(batch_size*nblks_ao))
3483 DO ib = 1, batch_size
3484 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3487 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3488 CALL dbt_create(t_3c_work_2(1),
"work_2_stack", t_dist, [1], [2, 3], &
3489 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3490 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3491 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3492 CALL dbt_distribution_destroy(t_dist)
3493 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3496 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3498 nspins =
SIZE(t_3c_apc, 1)
3499 CALL timeset(routinen//
"_apc", handle2)
3501 DO i_spin = 1, nspins
3502 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3505 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3506 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, ngroups, para_env, &
3507 subgroup_dest, ri_data%iatom_to_subgroup, 1, idx_to_at)
3508 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3509 CALL dbt_filter(t_3c_apc_sub(i_spin, i_img), ri_data%filter_eps)
3511 DO i_spin = 1, nspins
3512 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3515 CALL timestop(handle2)
3516 CALL dbt_pgrid_destroy(pgrid)
3517 CALL dbt_destroy(tmp)
3518 CALL dbt_destroy(work_atom_block)
3519 CALL dbt_destroy(work_atom_block_sub)
3521 CALL timestop(handle)
3523 END SUBROUTINE get_subgroup_3c_tensors
3545 SUBROUTINE get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
3546 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
3547 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
3548 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_inv, t_2c_bint, t_2c_metric
3549 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3550 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work
3551 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
3553 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot, mat_der_pot_sub
3554 INTEGER,
INTENT(IN) :: group_size, ngroups
3558 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_derivs'
3560 INTEGER :: handle, i, i_img, i_ri, i_spin, i_xyz, &
3561 iatom, iproc, j, natom, nblks, nimg, &
3563 INTEGER(int_8) :: nze
3564 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3566 INTEGER,
DIMENSION(2) :: pdims_2d
3567 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3568 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3571 TYPE(dbt_pgrid_type) :: pgrid_2d
3572 TYPE(dbt_type) :: work, work_sub
3574 CALL timeset(routinen, handle)
3579 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3581 natom =
SIZE(ri_data%bsizes_RI)
3582 nblks =
SIZE(ri_data%bsizes_RI_split)
3583 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3584 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3585 DO i_ri = 1, ri_data%ncell_RI
3586 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3587 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3592 bsizes_ri_ext, bsizes_ri_ext, &
3594 DEALLOCATE (dist1, dist2)
3596 CALL dbt_create(t_2c_inv(1), t_2c_bint(1))
3597 CALL dbt_create(t_2c_inv(1), t_2c_metric(1))
3599 CALL dbt_create(t_2c_inv(1), t_2c_inv(iatom))
3600 CALL dbt_create(t_2c_inv(1), t_2c_bint(iatom))
3601 CALL dbt_create(t_2c_inv(1), t_2c_metric(iatom))
3603 CALL dbt_create(t_2c_inv(1), t_2c_work(1))
3604 CALL dbt_create(t_2c_inv(1), t_2c_work(2))
3605 CALL dbt_create(t_2c_inv(1), t_2c_work(3))
3606 CALL dbt_create(t_2c_inv(1), t_2c_work(4))
3609 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3611 DEALLOCATE (dist1, dist2)
3615 CALL copy_2c_to_subgroup(t_2c_inv(iatom), ri_data%t_2c_inv(1, iatom), group_size, ngroups, para_env)
3616 CALL copy_2c_to_subgroup(t_2c_bint(iatom), ri_data%t_2c_int(1, iatom), group_size, ngroups, para_env)
3617 CALL copy_2c_to_subgroup(t_2c_metric(iatom), ri_data%t_2c_pot(1, iatom), group_size, ngroups, para_env)
3623 CALL dbt_create(t_2c_inv(1), t_2c_der_metric_sub(iatom, i_xyz))
3624 CALL copy_2c_to_subgroup(t_2c_der_metric_sub(iatom, i_xyz), t_2c_der_metric(iatom, i_xyz), &
3625 group_size, ngroups, para_env)
3626 CALL dbt_destroy(t_2c_der_metric(iatom, i_xyz))
3632 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3634 DEALLOCATE (dist1, dist2)
3635 nspins =
SIZE(rho_ao_t, 1)
3636 nimg =
SIZE(rho_ao_t, 2)
3639 DO i_spin = 1, nspins
3640 IF (.NOT. (i_img == 1 .AND. i_spin == 1)) &
3641 CALL dbt_create(rho_ao_t_sub(1, 1), rho_ao_t_sub(i_spin, i_img))
3642 CALL copy_2c_to_subgroup(rho_ao_t_sub(i_spin, i_img), rho_ao_t(i_spin, i_img), &
3643 group_size, ngroups, para_env)
3644 CALL dbt_destroy(rho_ao_t(i_spin, i_img))
3650 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3652 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3654 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3656 DO i = 0, pdims_2d(1) - 1
3657 DO j = 0, pdims_2d(2) - 1
3658 dbcsr_pgrid(i, j) = iproc
3664 ALLOCATE (col_dist(natom), row_dist(natom))
3665 row_dist(:) = dist1(:)
3666 col_dist(:) = dist2(:)
3668 ALLOCATE (ri_blk_size(natom))
3669 ri_blk_size(:) = ri_data%bsizes_RI(:)
3672 row_dist=row_dist, col_dist=col_dist)
3673 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3674 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3678 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3679 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3683 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3684 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3685 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3686 CALL dbt_clear(work_sub)
3692 CALL dbcsr_create(mat_der_pot_sub(i_img, i_xyz), template=mat_2c_pot(1))
3693 CALL dbt_copy_matrix_to_tensor(mat_der_pot(i_img, i_xyz), work)
3698 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3699 CALL dbt_copy_tensor_to_matrix(work_sub, mat_der_pot_sub(i_img, i_xyz))
3700 CALL dbcsr_filter(mat_der_pot_sub(i_img, i_xyz), ri_data%filter_eps)
3701 CALL dbt_clear(work_sub)
3705 CALL dbt_destroy(work)
3706 CALL dbt_destroy(work_sub)
3707 CALL dbt_pgrid_destroy(pgrid_2d)
3709 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3711 CALL timestop(handle)
3713 END SUBROUTINE get_subgroup_2c_derivs
3733 SUBROUTINE get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_AO, t_3c_der_AO_sub, &
3734 t_3c_der_RI, t_3c_der_RI_sub, t_3c_apc, t_3c_apc_sub, &
3735 t_3c_der_stack, group_size, ngroups, para_env, para_env_sub, &
3737 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_work_2, t_3c_work_3
3738 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ao, t_3c_der_ao_sub, &
3739 t_3c_der_ri, t_3c_der_ri_sub, &
3740 t_3c_apc, t_3c_apc_sub
3741 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_der_stack
3742 INTEGER,
INTENT(IN) :: group_size, ngroups
3746 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_derivs'
3748 INTEGER :: batch_size, handle, i_img, i_ri, i_spin, &
3749 i_xyz, ib, nblks_ao, nblks_ri, nimg, &
3751 INTEGER(int_8) :: nze
3752 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3753 bsizes_stack, dist1, dist2, dist3, &
3754 dist_stack, idx_to_at
3755 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :) :: subgroup_dest
3757 TYPE(dbt_distribution_type) :: t_dist
3758 TYPE(dbt_pgrid_type) :: pgrid
3759 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3761 CALL timeset(routinen, handle)
3764 nblks_ri =
SIZE(ri_data%bsizes_RI)
3765 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3766 DO i_ri = 1, ri_data%ncell_RI
3767 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
3770 CALL dbt_get_info(ri_data%kp_t_3c_int(1), pdims=pdims)
3771 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3774 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3775 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3776 DEALLOCATE (dist1, dist2, dist3)
3779 ri_data%pgrid_2, bsizes_ri_ext, ri_data%bsizes_AO, &
3780 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3781 DEALLOCATE (dist1, dist2, dist3)
3782 CALL dbt_pgrid_destroy(pgrid)
3784 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3785 group_size, ngroups, para_env)
3791 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ao_sub(i_img, i_xyz))
3795 CALL dbt_copy(t_3c_der_ao(i_img, i_xyz), work_atom_block, move_data=.true.)
3796 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3797 ngroups, para_env, subgroup_dest)
3798 CALL dbt_copy(work_atom_block_sub, t_3c_der_ao_sub(i_img, i_xyz), move_data=.true.)
3799 CALL dbt_filter(t_3c_der_ao_sub(i_img, i_xyz), ri_data%filter_eps)
3803 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ri_sub(i_img, i_xyz))
3807 CALL dbt_copy(t_3c_der_ri(i_img, i_xyz), work_atom_block, move_data=.true.)
3808 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3809 ngroups, para_env, subgroup_dest)
3810 CALL dbt_copy(work_atom_block_sub, t_3c_der_ri_sub(i_img, i_xyz), move_data=.true.)
3811 CALL dbt_filter(t_3c_der_ri_sub(i_img, i_xyz), ri_data%filter_eps)
3815 CALL dbt_destroy(t_3c_der_ri(i_img, i_xyz))
3816 CALL dbt_destroy(t_3c_der_ao(i_img, i_xyz))
3819 CALL dbt_destroy(work_atom_block_sub)
3820 CALL dbt_destroy(work_atom_block)
3821 DEALLOCATE (subgroup_dest)
3824 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3825 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3826 DO i_ri = 1, ri_data%ncell_RI
3827 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3831 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3832 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3835 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3836 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3837 DEALLOCATE (dist1, dist2, dist3)
3840 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3841 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3842 DEALLOCATE (dist1, dist2, dist3)
3845 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3846 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3847 DEALLOCATE (dist1, dist2, dist3)
3849 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3850 group_size, ngroups, para_env)
3852 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3854 nspins =
SIZE(t_3c_apc, 1)
3856 DO i_spin = 1, nspins
3857 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3860 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3861 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, ngroups, para_env, &
3862 subgroup_dest, ri_data%iatom_to_subgroup, 1, idx_to_at)
3863 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3864 CALL dbt_filter(t_3c_apc_sub(i_spin, i_img), ri_data%filter_eps)
3866 DO i_spin = 1, nspins
3867 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3870 CALL dbt_destroy(tmp)
3871 CALL dbt_destroy(work_atom_block)
3872 CALL dbt_destroy(work_atom_block_sub)
3873 CALL dbt_pgrid_destroy(pgrid)
3876 batch_size = ri_data%kp_stack_size
3877 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3878 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3879 DO ib = 1, batch_size
3880 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3883 ALLOCATE (dist1(ri_data%ncell_RI*nblks_ri), dist2(nblks_ao), dist3(nblks_ao))
3884 CALL dbt_get_info(ri_data%kp_t_3c_int(1), proc_dist_1=dist1, proc_dist_2=dist2, &
3885 proc_dist_3=dist3, pdims=pdims)
3887 ALLOCATE (dist_stack(batch_size*nblks_ao))
3888 DO ib = 1, batch_size
3889 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3892 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3893 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3894 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3895 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3896 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3897 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3898 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(4))
3899 CALL dbt_distribution_destroy(t_dist)
3900 CALL dbt_pgrid_destroy(pgrid)
3901 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3904 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(1))
3905 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(2))
3906 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(3))
3907 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(4))
3908 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(5))
3909 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(6))
3912 ALLOCATE (dist1(nblks_ao), dist2(ri_data%ncell_RI*nblks_ri), dist3(nblks_ao))
3913 CALL dbt_get_info(t_3c_apc_sub(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, &
3914 proc_dist_3=dist3, pdims=pdims)
3916 ALLOCATE (dist_stack(batch_size*nblks_ao))
3917 DO ib = 1, batch_size
3918 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3921 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3922 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3923 CALL dbt_create(t_3c_work_2(1),
"work_3_stack", t_dist, [1], [2, 3], &
3924 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3925 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3926 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3927 CALL dbt_distribution_destroy(t_dist)
3928 CALL dbt_pgrid_destroy(pgrid)
3929 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3931 CALL timestop(handle)
3933 END SUBROUTINE get_subgroup_3c_derivs
3941 SUBROUTINE reorder_3c_ints(t_3c_ints, ri_data)
3942 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_ints
3945 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_ints'
3947 INTEGER :: handle, i_img,
idx, idx_empty, idx_full, &
3949 INTEGER(int_8) :: nze
3951 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
3953 CALL timeset(routinen, handle)
3956 ALLOCATE (t_3c_tmp(nimg))
3958 CALL dbt_create(t_3c_ints(i_img), t_3c_tmp(i_img))
3959 CALL dbt_copy(t_3c_ints(i_img), t_3c_tmp(i_img), move_data=.true.)
3964 ALLOCATE (ri_data%idx_to_img(nimg))
3966 idx_empty = nimg + 1
3971 idx_empty = idx_empty - 1
3972 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_empty), move_data=.true.)
3973 ri_data%idx_to_img(idx_empty) = i_img
3975 idx_full = idx_full + 1
3976 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_full), move_data=.true.)
3977 ri_data%idx_to_img(idx_full) = i_img
3979 CALL dbt_destroy(t_3c_tmp(i_img))
3983 ri_data%nimg_nze = idx_full
3985 ALLOCATE (ri_data%img_to_idx(nimg))
3987 ri_data%img_to_idx(ri_data%idx_to_img(
idx)) =
idx
3990 CALL timestop(handle)
3992 END SUBROUTINE reorder_3c_ints
4000 SUBROUTINE reorder_3c_derivs(t_3c_derivs, ri_data)
4001 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_derivs
4004 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_derivs'
4006 INTEGER :: handle, i_img, i_xyz,
idx, nimg
4007 INTEGER(int_8) :: nze
4009 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
4011 CALL timeset(routinen, handle)
4014 ALLOCATE (t_3c_tmp(nimg))
4016 CALL dbt_create(t_3c_derivs(1, 1), t_3c_tmp(i_img))
4021 CALL dbt_copy(t_3c_derivs(i_img, i_xyz), t_3c_tmp(i_img), move_data=.true.)
4024 idx = ri_data%img_to_idx(i_img)
4025 CALL dbt_copy(t_3c_tmp(i_img), t_3c_derivs(
idx, i_xyz), move_data=.true.)
4027 IF (nze > 0) ri_data%nimg_nze = max(
idx, ri_data%nimg_nze)
4032 CALL dbt_destroy(t_3c_tmp(i_img))
4035 CALL timestop(handle)
4037 END SUBROUTINE reorder_3c_derivs
4045 SUBROUTINE get_sparsity_pattern(pattern, ri_data, qs_env)
4046 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: pattern
4050 INTEGER :: iatom, j_img, jatom, mj_img, natom, nimg
4051 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bins
4052 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: tmp_pattern
4053 INTEGER,
DIMENSION(3) :: cell_j
4054 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
4055 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
4060 DIMENSION(:),
POINTER :: nl_iterator
4064 NULLIFY (nl_2c, nl_iterator, kpoints, cell_to_index, dft_control, index_to_cell, para_env)
4066 CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, para_env=para_env, natom=natom)
4067 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell, sab_nl=nl_2c)
4070 pattern(:, :, :) = 0
4077 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4078 IF (j_img > nimg .OR. j_img < 1) cycle
4080 mj_img = get_opp_index(j_img, qs_env)
4081 IF (mj_img > nimg .OR. mj_img < 1) cycle
4083 IF (ri_data%present_images(j_img) == 0) cycle
4085 pattern(iatom, jatom, j_img) = 1
4096 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4097 IF (j_img > nimg .OR. j_img < 1) cycle
4099 mj_img = get_opp_index(j_img, qs_env)
4100 IF (mj_img <= nimg .AND. mj_img > 0) cycle
4102 IF (ri_data%present_images(j_img) == 0) cycle
4104 pattern(iatom, jatom, j_img) = 1
4108 CALL para_env%sum(pattern)
4113 IF (pattern(iatom, iatom, j_img) /= 0)
THEN
4114 mj_img = get_opp_index(j_img, qs_env)
4115 IF (mj_img > nimg .OR. mj_img < 1) cycle
4116 pattern(iatom, iatom, mj_img) = 0
4123 ALLOCATE (bins(natom))
4126 ALLOCATE (tmp_pattern(natom, natom, nimg))
4127 tmp_pattern(:, :, :) = 0
4131 IF (pattern(iatom, jatom, j_img) == 0) cycle
4132 mj_img = get_opp_index(j_img, qs_env)
4135 IF (mj_img > nimg .OR. mj_img < 1)
THEN
4137 bins(iatom) = bins(iatom) + 1
4138 tmp_pattern(iatom, jatom, j_img) = 1
4141 IF (bins(iatom) > bins(jatom))
THEN
4142 bins(jatom) = bins(jatom) + 1
4143 tmp_pattern(jatom, iatom, mj_img) = 1
4145 bins(iatom) = bins(iatom) + 1
4146 tmp_pattern(iatom, jatom, j_img) = 1
4154 pattern(:, :, :) = tmp_pattern(:, :, :) - 1
4156 END SUBROUTINE get_sparsity_pattern
4166 SUBROUTINE get_sub_dist(sparsity_pattern, ngroups, ri_data)
4167 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: sparsity_pattern
4168 INTEGER,
INTENT(IN) :: ngroups
4171 INTEGER :: b_img, ctr, iat, iatom, igroup, jatom, &
4173 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: max_at_per_group
4175 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4177 natom =
SIZE(sparsity_pattern, 2)
4178 nimg =
SIZE(sparsity_pattern, 3)
4183 IF (.NOT.
ALLOCATED(ri_data%iatom_to_subgroup))
THEN
4184 ALLOCATE (ri_data%iatom_to_subgroup(natom), max_at_per_group(ngroups))
4186 NULLIFY (ri_data%iatom_to_subgroup(iatom)%array)
4187 ALLOCATE (ri_data%iatom_to_subgroup(iatom)%array(ngroups))
4188 ri_data%iatom_to_subgroup(iatom)%array(:) = .false.
4192 IF (ub*ngroups < natom) ub = ub + 1
4193 max_at_per_group(:) = max(1, ub)
4198 DO WHILE (
modulo(sum(max_at_per_group), natom) /= 0)
4199 igroup =
modulo(ctr, ngroups) + 1
4200 max_at_per_group(igroup) = max_at_per_group(igroup) + 1
4205 DO igroup = 1, ngroups
4206 DO iat = 1, max_at_per_group(igroup)
4207 iatom =
modulo(ctr, natom) + 1
4208 ri_data%iatom_to_subgroup(iatom)%array(igroup) = .true.
4214 ALLOCATE (bins(ngroups))
4219 IF (sparsity_pattern(iatom, jatom, b_img) == -1) cycle
4220 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4223 IF (any(ri_data%kp_cost > epsilon(0.0_dp)))
THEN
4224 cost = ri_data%kp_cost(iatom, jatom, b_img)
4226 cost = real(ri_data%bsizes_AO(iatom)*ri_data%bsizes_AO(jatom),
dp)
4228 bins(igroup + 1) = bins(igroup + 1) + cost
4229 sparsity_pattern(iatom, jatom, b_img) = igroup
4234 END SUBROUTINE get_sub_dist
4245 SUBROUTINE update_pattern_to_forces(force_pattern, scf_pattern, ngroups, ri_data, qs_env)
4246 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: force_pattern, scf_pattern
4247 INTEGER,
INTENT(IN) :: ngroups
4251 INTEGER :: b_img, iatom, igroup, jatom, mb_img, &
4253 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4255 natom =
SIZE(scf_pattern, 2)
4256 nimg =
SIZE(scf_pattern, 3)
4258 ALLOCATE (bins(ngroups))
4262 mb_img = get_opp_index(b_img, qs_env)
4266 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4269 IF (scf_pattern(iatom, jatom, b_img) > -1) cycle
4272 IF (mb_img > 0 .AND. mb_img <= nimg)
THEN
4273 IF (scf_pattern(jatom, iatom, mb_img) == -1) cycle
4274 bins(igroup + 1) = bins(igroup + 1) + ri_data%kp_cost(jatom, iatom, mb_img)
4275 force_pattern(iatom, jatom, b_img) = igroup
4281 END SUBROUTINE update_pattern_to_forces
4289 SUBROUTINE get_kp_and_ri_images(ri_data, qs_env)
4293 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_kp_and_ri_images'
4295 CHARACTER(LEN=512) :: warning_msg
4296 INTEGER :: cell_j(3), cell_k(3), handle, i_img, iatom, ikind, j_img, jatom, jcell, katom, &
4297 kcell, kp_index_lbounds(3), kp_index_ubounds(3), natom, ngroups, nimg, nkind, pcoord(3), &
4299 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_ri, &
4300 nri_per_atom, present_img, ri_cells
4301 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
4302 REAL(
dp) :: bump_fact, dij, dik, image_range, &
4303 ri_range, rij(3), rik(3)
4304 TYPE(dbt_type) :: t_dummy
4309 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4316 DIMENSION(:),
POINTER :: nl_iterator
4320 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4323 NULLIFY (qs_kind_set, dist_2d, nl_2c, nl_iterator, dft_control, &
4324 particle_set, kpoints, para_env, cell_to_index, hfx_section)
4326 CALL timeset(routinen, handle)
4328 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, &
4329 dft_control=dft_control, particle_set=particle_set, kpoints=kpoints, &
4330 para_env=para_env, natom=natom)
4331 nimg = dft_control%nimages
4333 kp_index_lbounds = lbound(cell_to_index)
4334 kp_index_ubounds = ubound(cell_to_index)
4339 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4346 CALL erfc_cutoff(ri_data%eps_schwarz, ri_data%hfx_pot%omega, ri_data%hfx_pot%cutoff_radius)
4347 WRITE (warning_msg,
'(A)') &
4348 "The SHORTANGE HFX potential typically extends over many periodic images, "// &
4349 "possibly slowing down the calculation. Consider using the TRUNCATED "// &
4350 "potential for better computational performance."
4355 ri_data%kp_RI_range = 0.0_dp
4356 ri_data%kp_image_range = 0.0_dp
4361 ri_data%kp_RI_range = max(ri_range, ri_data%kp_RI_range)
4365 CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, kind_radius=image_range)
4368 ri_data%kp_image_range = max(image_range, ri_data%kp_image_range)
4372 ri_data%kp_bump_rad = bump_fact*ri_data%kp_RI_range
4378 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
4380 ALLOCATE (present_img(nimg))
4389 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4390 IF (j_img > nimg .OR. j_img < 1) cycle
4392 IF (dij > ri_data%kp_image_range) cycle
4394 ri_data%nimg = max(j_img, ri_data%nimg)
4395 present_img(j_img) = 1
4400 CALL para_env%max(ri_data%nimg)
4401 IF (ri_data%nimg > nimg) &
4402 cpabort(
"Make sure the smallest exponent of the RI-HFX basis is larger than that of the ORB basis.")
4405 CALL para_env%sum(present_img)
4406 ALLOCATE (ri_data%present_images(ri_data%nimg))
4407 ri_data%present_images = 0
4408 DO i_img = 1, ri_data%nimg
4409 IF (present_img(i_img) > 0) ri_data%present_images(i_img) = 1
4413 ri_data%pgrid, ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4414 map1=[1, 2], map2=[3], name=
"(AO AO | RI)")
4416 CALL dbt_mp_environ_pgrid(ri_data%pgrid, pdims, pcoord)
4417 CALL mp_comm_t3c%create(ri_data%pgrid%mp_comm_2d, 3, pdims)
4419 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4420 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4421 CALL dbt_destroy(t_dummy)
4426 ri_data%ri_metric,
"HFX_3c_nl", qs_env, op_pos=2, sym_ij=.false., &
4429 ALLOCATE (ri_cells(nimg))
4432 ALLOCATE (nri_per_atom(natom))
4438 iatom=iatom, jatom=jatom, katom=katom)
4441 IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
4442 any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
4444 jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4445 IF (jcell > nimg .OR. jcell < 1) cycle
4447 IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
4448 any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
4450 kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
4451 IF (kcell > nimg .OR. kcell < 1) cycle
4453 IF (dik > ri_data%kp_RI_range) cycle
4456 IF (jcell == 1 .AND. iatom == jatom) nri_per_atom(iatom) = nri_per_atom(iatom) + ri_data%bsizes_RI(katom)
4460 CALL para_env%sum(ri_cells)
4461 CALL para_env%sum(nri_per_atom)
4463 ALLOCATE (ri_data%img_to_RI_cell(nimg))
4464 ri_data%ncell_RI = 0
4465 ri_data%img_to_RI_cell = 0
4467 IF (ri_cells(i_img) > 0)
THEN
4468 ri_data%ncell_RI = ri_data%ncell_RI + 1
4469 ri_data%img_to_RI_cell(i_img) = ri_data%ncell_RI
4473 ALLOCATE (ri_data%RI_cell_to_img(ri_data%ncell_RI))
4475 IF (ri_data%img_to_RI_cell(i_img) > 0) ri_data%RI_cell_to_img(ri_data%img_to_RI_cell(i_img)) = i_img
4479 IF (ri_data%unit_nr > 0)
THEN
4480 WRITE (ri_data%unit_nr, fmt=
"(/T3,A,I29)") &
4481 "KP-HFX_RI_INFO| Number of RI-KP parallel groups:", ngroups
4482 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I29)") &
4483 "KP-HFX_RI_INFO| Tensor stack size: ", ri_data%kp_stack_size
4484 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F31.3,A)") &
4485 "KP-HFX_RI_INFO| RI basis extension radius:", ri_data%kp_RI_range*
angstrom,
" Ang"
4486 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F12.3,A, F6.3, A)") &
4487 "KP-HFX_RI_INFO| RI basis bump factor and bump radius:", bump_fact,
" /", &
4488 ri_data%kp_bump_rad*
angstrom,
" Ang"
4489 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I16,A)") &
4490 "KP-HFX_RI_INFO| The extended RI bases cover up to ", ri_data%ncell_RI,
" unit cells"
4491 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I18)") &
4492 "KP-HFX_RI_INFO| Average number of sgf in extended RI bases:", sum(nri_per_atom)/natom
4493 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F13.3,A)") &
4494 "KP-HFX_RI_INFO| Consider all image cells within a radius of ", ri_data%kp_image_range*
angstrom,
" Ang"
4495 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I27/)") &
4496 "KP-HFX_RI_INFO| Number of image cells considered: ", ri_data%nimg
4500 CALL timestop(handle)
4502 END SUBROUTINE get_kp_and_ri_images
4517 SUBROUTINE get_stack_tensors(res_stack, rho_stack, ints_stack, rho_template, ints_template, &
4518 stack_size, ri_data, qs_env)
4519 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: res_stack, rho_stack, ints_stack
4520 TYPE(dbt_type),
INTENT(INOUT) :: rho_template, ints_template
4521 INTEGER,
INTENT(IN) :: stack_size
4525 INTEGER :: is, nblks, nblks_3c(3), pdims_3d(3)
4526 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_stack, dist1, &
4527 dist2, dist3, dist_stack1, &
4528 dist_stack2, dist_stack3
4529 TYPE(dbt_distribution_type) :: t_dist
4530 TYPE(dbt_pgrid_type) :: pgrid
4537 nblks =
SIZE(ri_data%bsizes_AO_split)
4538 ALLOCATE (bsizes_stack(stack_size*nblks))
4539 DO is = 1, stack_size
4540 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
4543 ALLOCATE (dist1(nblks), dist2(nblks), dist_stack1(stack_size*nblks), dist_stack2(stack_size*nblks))
4544 CALL dbt_get_info(rho_template, proc_dist_1=dist1, proc_dist_2=dist2)
4545 DO is = 1, stack_size
4546 dist_stack1((is - 1)*nblks + 1:is*nblks) = dist1(:)
4547 dist_stack2((is - 1)*nblks + 1:is*nblks) = dist2(:)
4552 CALL dbt_distribution_new(t_dist, ri_data%pgrid_2d, dist_stack1, dist_stack2)
4553 CALL dbt_create(rho_stack(1),
"RHO_stack", t_dist, [1], [2], bsizes_stack, bsizes_stack)
4554 CALL dbt_distribution_destroy(t_dist)
4555 DEALLOCATE (dist1, dist2, dist_stack1, dist_stack2)
4558 CALL create_2c_tensor(rho_stack(2), dist1, dist2, ri_data%pgrid_2d, bsizes_stack, bsizes_stack, name=
"RHO_stack")
4559 DEALLOCATE (dist1, dist2)
4561 CALL dbt_get_info(ints_template, nblks_total=nblks_3c)
4562 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)))
4563 ALLOCATE (dist_stack3(stack_size*nblks_3c(3)), bsizes_ri_ext(nblks_3c(2)))
4564 CALL dbt_get_info(ints_template, proc_dist_1=dist1, proc_dist_2=dist2, &
4565 proc_dist_3=dist3, blk_size_2=bsizes_ri_ext)
4566 DO is = 1, stack_size
4567 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
4571 CALL dbt_distribution_new(t_dist, ri_data%pgrid_1, dist1, dist2, dist_stack3)
4572 CALL dbt_create(ints_stack(1),
"ints_stack", t_dist, [1, 2], [3], ri_data%bsizes_AO_split, &
4573 bsizes_ri_ext, bsizes_stack)
4574 CALL dbt_distribution_destroy(t_dist)
4575 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
4579 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid, tensor_dims=[nblks_3c(1), nblks_3c(2), stack_size*nblks_3c(3)])
4580 CALL create_3c_tensor(ints_stack(2), dist1, dist2, dist3, pgrid, ri_data%bsizes_AO_split, &
4581 bsizes_ri_ext, bsizes_stack, [1, 2], [3], name=
"ints_stack")
4582 DEALLOCATE (dist1, dist2, dist3)
4583 CALL dbt_pgrid_destroy(pgrid)
4586 CALL dbt_create(ints_stack(1), res_stack(1))
4587 CALL dbt_create(ints_stack(2), res_stack(2))
4589 END SUBROUTINE get_stack_tensors
4603 SUBROUTINE fill_3c_stack(t_3c_stack, t_3c_in, images, stack_dim, ri_data, filter_at, filter_dim, &
4604 idx_to_at, img_bounds)
4605 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_stack
4606 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_in
4607 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4608 INTEGER,
INTENT(IN) :: stack_dim
4610 INTEGER,
INTENT(IN),
OPTIONAL :: filter_at, filter_dim
4611 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: idx_to_at
4612 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2)
4614 INTEGER :: dest(3), i_img,
idx, ind(3), lb, nblks, &
4616 LOGICAL :: do_filter, found
4617 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4618 TYPE(dbt_iterator_type) :: iter
4623 nblks =
SIZE(ri_data%bsizes_AO_split)
4626 IF (
PRESENT(filter_at) .AND.
PRESENT(filter_dim) .AND.
PRESENT(idx_to_at)) do_filter = .true.
4631 IF (
PRESENT(img_bounds))
THEN
4633 ub = img_bounds(2) - 1
4639 IF (i_img == 0 .OR. i_img > nimg) cycle
4644 CALL dbt_iterator_start(iter, t_3c_in(i_img))
4645 DO WHILE (dbt_iterator_blocks_left(iter))
4646 CALL dbt_iterator_next_block(iter, ind)
4647 CALL dbt_get_block(t_3c_in(i_img), ind, blk, found)
4648 IF (.NOT. found) cycle
4651 IF (.NOT. idx_to_at(ind(filter_dim)) == filter_at) cycle
4654 IF (stack_dim == 1)
THEN
4655 dest = [(
idx - offset - 1)*nblks + ind(1), ind(2), ind(3)]
4656 ELSE IF (stack_dim == 2)
THEN
4657 dest = [ind(1), (
idx - offset - 1)*nblks + ind(2), ind(3)]
4659 dest = [ind(1), ind(2), (
idx - offset - 1)*nblks + ind(3)]
4662 CALL dbt_put_block(t_3c_stack, dest, shape(blk), blk)
4665 CALL dbt_iterator_stop(iter)
4668 CALL dbt_finalize(t_3c_stack)
4670 END SUBROUTINE fill_3c_stack
4682 SUBROUTINE fill_2c_stack(t_2c_stack, t_2c_in, images, stack_dim, ri_data, img_bounds, shift)
4683 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_stack
4684 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_in
4685 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4686 INTEGER,
INTENT(IN) :: stack_dim
4688 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2), shift
4690 INTEGER :: dest(2), i_img,
idx, ind(2), lb, &
4691 my_shift, nblks, nimg, offset, ub
4693 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
4694 TYPE(dbt_iterator_type) :: iter
4699 nblks =
SIZE(ri_data%bsizes_AO_split)
4704 IF (
PRESENT(img_bounds))
THEN
4706 ub = img_bounds(2) - 1
4711 IF (
PRESENT(shift)) my_shift = shift
4715 IF (i_img == 0 .OR. i_img > nimg) cycle
4719 CALL dbt_iterator_start(iter, t_2c_in(i_img))
4720 DO WHILE (dbt_iterator_blocks_left(iter))
4721 CALL dbt_iterator_next_block(iter, ind)
4722 CALL dbt_get_block(t_2c_in(i_img), ind, blk, found)
4723 IF (.NOT. found) cycle
4725 IF (stack_dim == 1)
THEN
4726 dest = [(
idx - offset - 1)*nblks + ind(1), (my_shift - 1)*nblks + ind(2)]
4728 dest = [(my_shift - 1)*nblks + ind(1), (
idx - offset - 1)*nblks + ind(2)]
4731 CALL dbt_put_block(t_2c_stack, dest, shape(blk), blk)
4734 CALL dbt_iterator_stop(iter)
4737 CALL dbt_finalize(t_2c_stack)
4739 END SUBROUTINE fill_2c_stack
4747 SUBROUTINE unstack_t_3c_apc(t_3c_apc, t_stacked, idx)
4748 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_apc, t_stacked
4749 INTEGER,
INTENT(IN) ::
idx
4751 INTEGER :: current_idx
4752 INTEGER,
DIMENSION(3) :: ind, nblks_3c
4754 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4755 TYPE(dbt_iterator_type) :: iter
4758 CALL dbt_get_info(t_3c_apc, nblks_total=nblks_3c)
4761 CALL dbt_iterator_start(iter, t_stacked)
4762 DO WHILE (dbt_iterator_blocks_left(iter))
4763 CALL dbt_iterator_next_block(iter, ind)
4766 current_idx = (ind(3) - 1)/nblks_3c(3) + 1
4767 IF (.NOT.
idx == current_idx) cycle
4769 CALL dbt_get_block(t_stacked, ind, blk, found)
4770 IF (.NOT. found) cycle
4772 CALL dbt_put_block(t_3c_apc, [ind(1), ind(2), ind(3) - (
idx - 1)*nblks_3c(3)], shape(blk), blk)
4775 CALL dbt_iterator_stop(iter)
4778 END SUBROUTINE unstack_t_3c_apc
4788 SUBROUTINE get_atom_3c_ints(t_3c_at, t_3c_ints, iatom, dim_at, idx_to_at)
4789 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_at, t_3c_ints
4790 INTEGER,
INTENT(IN) :: iatom, dim_at
4791 INTEGER,
DIMENSION(:),
INTENT(IN) :: idx_to_at
4793 INTEGER,
DIMENSION(3) :: ind
4795 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4796 TYPE(dbt_iterator_type) :: iter
4799 CALL dbt_iterator_start(iter, t_3c_ints)
4800 DO WHILE (dbt_iterator_blocks_left(iter))
4801 CALL dbt_iterator_next_block(iter, ind)
4802 IF (.NOT. idx_to_at(ind(dim_at)) == iatom) cycle
4804 CALL dbt_get_block(t_3c_ints, ind, blk, found)
4805 IF (.NOT. found) cycle
4807 CALL dbt_put_block(t_3c_at, ind, shape(blk), blk)
4810 CALL dbt_iterator_stop(iter)
4812 CALL dbt_finalize(t_3c_at)
4814 END SUBROUTINE get_atom_3c_ints
4825 SUBROUTINE precalc_derivatives(t_3c_der_RI, t_3c_der_AO, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
4826 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ri, t_3c_der_ao
4827 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot
4828 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_2c_der_metric
4832 CHARACTER(LEN=*),
PARAMETER :: routinen =
'precalc_derivatives'
4834 INTEGER :: handle, handle2, i_img, i_mem, i_ri, &
4835 i_xyz, iatom, n_mem, natom, nblks_ri, &
4836 ncell_ri, nimg, nkind, nthreads
4837 INTEGER(int_8) :: nze
4838 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, dist_ao_1, &
4839 dist_ao_2, dist_ri, dist_ri_ext, dummy_end, dummy_start, end_blocks, start_blocks
4840 INTEGER,
DIMENSION(3) :: pcoord, pdims
4841 INTEGER,
DIMENSION(:),
POINTER :: col_bsize, row_bsize
4845 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_metric
4846 TYPE(dbt_distribution_type) :: t_dist
4847 TYPE(dbt_pgrid_type) :: pgrid
4848 TYPE(dbt_type) :: t_3c_template
4849 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
4854 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4861 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4863 NULLIFY (qs_kind_set, dist_2d, nl_2c, particle_set, dft_control, para_env, row_bsize, col_bsize)
4865 CALL timeset(routinen, handle)
4867 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, natom=natom, &
4868 particle_set=particle_set, dft_control=dft_control, para_env=para_env)
4871 ncell_ri = ri_data%ncell_RI
4873 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4883 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[max(1, natom/(ri_data%n_mem*nthreads)), natom, natom])
4885 CALL create_3c_tensor(t_3c_template, dist_ao_1, dist_ao_2, dist_ri, pgrid, &
4886 ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4887 map1=[1, 2], map2=[3], name=
"tmp")
4888 CALL dbt_destroy(t_3c_template)
4891 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
4892 ALLOCATE (dist_ri_ext(natom*ncell_ri))
4893 ALLOCATE (bsizes_ri_ext(natom*ncell_ri))
4894 ALLOCATE (bsizes_ri_ext_split(nblks_ri*ncell_ri))
4895 DO i_ri = 1, ncell_ri
4896 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
4897 dist_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = dist_ri(:)
4898 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
4901 CALL dbt_distribution_new(t_dist, pgrid, dist_ao_1, dist_ao_2, dist_ri_ext)
4902 CALL dbt_create(t_3c_template,
"KP_3c_der", t_dist, [1, 2], [3], &
4903 ri_data%bsizes_AO, ri_data%bsizes_AO, bsizes_ri_ext)
4904 CALL dbt_distribution_destroy(t_dist)
4906 ALLOCATE (t_3c_der_ri_prv(nimg, 1, 3), t_3c_der_ao_prv(nimg, 1, 3))
4909 CALL dbt_create(t_3c_template, t_3c_der_ri_prv(i_img, 1, i_xyz))
4910 CALL dbt_create(t_3c_template, t_3c_der_ao_prv(i_img, 1, i_xyz))
4913 CALL dbt_destroy(t_3c_template)
4915 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
4916 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
4918 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4919 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4920 CALL dbt_pgrid_destroy(pgrid)
4923 "HFX_3c_nl", qs_env, op_pos=2, sym_jk=.false., own_dist=.true.)
4925 n_mem = ri_data%n_mem
4927 start_blocks, end_blocks)
4928 DEALLOCATE (dummy_start, dummy_end)
4930 CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, ri_data%pgrid_2, &
4931 bsizes_ri_ext_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
4932 map1=[1], map2=[2, 3], name=
"der (RI | AO AO)")
4935 CALL dbt_create(t_3c_template, t_3c_der_ri(i_img, i_xyz))
4936 CALL dbt_create(t_3c_template, t_3c_der_ao(i_img, i_xyz))
4942 nl_3c, basis_set_ao, basis_set_ao, basis_set_ri, &
4943 ri_data%ri_metric, der_eps=ri_data%eps_schwarz_forces, op_pos=2, &
4944 do_kpoints=.true., do_hfx_kpoints=.true., &
4945 bounds_k=[start_blocks(i_mem), end_blocks(i_mem)], &
4946 ri_range=ri_data%kp_RI_range, img_to_ri_cell=ri_data%img_to_RI_cell)
4948 CALL timeset(routinen//
"_cpy", handle2)
4955 CALL dbt_copy(t_3c_der_ao_prv(i_img, 1, i_xyz), t_3c_template, &
4956 order=[3, 2, 1], move_data=.true.)
4957 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4958 CALL dbt_copy(t_3c_template, t_3c_der_ao(i_img, i_xyz), &
4959 order=[1, 3, 2], move_data=.true., summation=.true.)
4965 CALL dbt_copy(t_3c_der_ri_prv(i_img, 1, i_xyz), t_3c_template, &
4966 order=[3, 2, 1], move_data=.true.)
4967 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4968 CALL dbt_copy(t_3c_template, t_3c_der_ri(i_img, i_xyz), &
4969 order=[1, 3, 2], move_data=.true., summation=.true.)
4973 CALL timestop(handle2)
4975 CALL dbt_destroy(t_3c_template)
4980 CALL dbt_destroy(t_3c_der_ri_prv(i_img, 1, i_xyz))
4981 CALL dbt_destroy(t_3c_der_ao_prv(i_img, 1, i_xyz))
4984 DEALLOCATE (t_3c_der_ri_prv, t_3c_der_ao_prv)
4987 CALL reorder_3c_derivs(t_3c_der_ri, ri_data)
4988 CALL reorder_3c_derivs(t_3c_der_ao, ri_data)
4990 CALL timeset(routinen//
"_2c", handle2)
4993 ALLOCATE (row_bsize(
SIZE(ri_data%bsizes_RI)))
4994 ALLOCATE (col_bsize(
SIZE(ri_data%bsizes_RI)))
4995 row_bsize(:) = ri_data%bsizes_RI
4996 col_bsize(:) = ri_data%bsizes_RI
4998 CALL dbcsr_create(dbcsr_template,
"2c_der", dbcsr_dist, dbcsr_type_no_symmetry, &
4999 row_bsize, col_bsize)
5001 DEALLOCATE (col_bsize, row_bsize)
5003 ALLOCATE (mat_der_metric(nimg, 3))
5006 CALL dbcsr_create(mat_der_pot(i_img, i_xyz), template=dbcsr_template)
5007 CALL dbcsr_create(mat_der_metric(i_img, i_xyz), template=dbcsr_template)
5014 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
5016 basis_set_ri, basis_set_ri, ri_data%hfx_pot, do_kpoints=.true.)
5021 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
5023 basis_set_ri, basis_set_ri, ri_data%ri_metric, do_kpoints=.true.)
5029 CALL dbt_create(ri_data%t_2c_inv(1, 1), t_2c_der_metric(iatom, i_xyz))
5030 CALL get_ext_2c_int(t_2c_der_metric(iatom, i_xyz), mat_der_metric(:, i_xyz), &
5031 iatom, iatom, 1, ri_data, qs_env)
5037 CALL timestop(handle2)
5039 CALL timestop(handle)
5041 END SUBROUTINE precalc_derivatives
5062 SUBROUTINE get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, img, pref, &
5063 ri_data, qs_env, work_virial, cell, particle_set, diag, offdiag)
5066 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_contr
5067 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_der
5068 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
5069 INTEGER,
INTENT(IN) :: img
5070 REAL(
dp),
INTENT(IN) :: pref
5073 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
5074 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
5076 POINTER :: particle_set
5077 LOGICAL,
INTENT(IN),
OPTIONAL ::
diag, offdiag
5079 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_2c_der_force'
5081 INTEGER :: handle, i_img, i_ri, i_xyz, iat, &
5082 iat_of_kind, ikind, j_img, j_ri, &
5083 j_xyz, jat, jat_of_kind, jkind, natom
5084 INTEGER,
DIMENSION(2) :: ind
5085 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
5086 LOGICAL :: found, my_diag, my_offdiag, use_virial
5087 REAL(
dp) :: new_force
5088 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :),
TARGET :: contr_blk, der_blk
5089 REAL(
dp),
DIMENSION(3) :: scoord
5090 TYPE(dbt_iterator_type) :: iter
5093 NULLIFY (kpoints, index_to_cell)
5098 CALL timeset(routinen, handle)
5100 use_virial = .false.
5101 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
5106 my_offdiag = .false.
5107 IF (
PRESENT(
diag)) my_offdiag = offdiag
5109 CALL get_qs_env(qs_env, kpoints=kpoints, natom=natom)
5118 CALL dbt_iterator_start(iter, t_2c_der(i_xyz))
5119 DO WHILE (dbt_iterator_blocks_left(iter))
5120 CALL dbt_iterator_next_block(iter, ind)
5123 IF ((my_diag .AND. .NOT. my_offdiag) .OR. (.NOT. my_diag .AND. my_offdiag))
THEN
5124 IF (my_diag .AND. (ind(1) /= ind(2))) cycle
5125 IF (my_offdiag .AND. (ind(1) == ind(2))) cycle
5128 CALL dbt_get_block(t_2c_der(i_xyz), ind, der_blk, found)
5130 CALL dbt_get_block(t_2c_contr, ind, contr_blk, found)
5136 new_force = pref*sum(der_blk(:, :)*contr_blk(:, :))
5138 i_ri = (ind(1) - 1)/natom + 1
5139 i_img = ri_data%RI_cell_to_img(i_ri)
5140 iat = ind(1) - (i_ri - 1)*natom
5141 iat_of_kind = atom_of_kind(iat)
5142 ikind = kind_of(iat)
5144 j_ri = (ind(2) - 1)/natom + 1
5145 j_img = ri_data%RI_cell_to_img(j_ri)
5146 jat = ind(2) - (j_ri - 1)*natom
5147 jat_of_kind = atom_of_kind(jat)
5148 jkind = kind_of(jat)
5152 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5155 IF (use_virial)
THEN
5158 scoord(:) = scoord(:) + real(index_to_cell(:, i_img),
dp)
5162 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5168 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5171 IF (use_virial)
THEN
5174 scoord(:) = scoord(:) + real(index_to_cell(:, j_img) + index_to_cell(:, img),
dp)
5178 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - new_force*scoord(j_xyz)
5182 DEALLOCATE (contr_blk)
5185 DEALLOCATE (der_blk)
5187 CALL dbt_iterator_stop(iter)
5191 CALL timestop(handle)
5217 idx_to_at_RI, idx_to_at_AO, i_images, lb_img, pref, &
5218 ri_data, qs_env, work_virial, cell, particle_set)
5221 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_contr
5222 TYPE(dbt_type),
DIMENSION(3),
INTENT(INOUT) :: t_3c_der_1, t_3c_der_2
5223 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of, idx_to_at_ri, &
5224 idx_to_at_ao, i_images
5225 INTEGER,
INTENT(IN) :: lb_img
5226 REAL(
dp),
INTENT(IN) :: pref
5229 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
5230 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
5232 POINTER :: particle_set
5234 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_force_from_3c_trace'
5236 INTEGER :: handle, i_ri, i_xyz, iat, iat_of_kind,
idx, ikind, j_xyz, jat, jat_of_kind, &
5237 jkind, kat, kat_of_kind, kkind, nblks_ao, nblks_ri, ri_img
5238 INTEGER,
DIMENSION(3) :: ind
5239 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
5240 LOGICAL :: found, found_1, found_2, use_virial
5241 REAL(
dp) :: new_force
5242 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :),
TARGET :: contr_blk, der_blk_1, der_blk_2, &
5244 REAL(
dp),
DIMENSION(3) :: scoord
5245 TYPE(dbt_iterator_type) :: iter
5248 NULLIFY (kpoints, index_to_cell)
5250 CALL timeset(routinen, handle)
5255 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
5256 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
5258 use_virial = .false.
5259 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
5266 CALL dbt_iterator_start(iter, t_3c_contr)
5267 DO WHILE (dbt_iterator_blocks_left(iter))
5268 CALL dbt_iterator_next_block(iter, ind)
5270 CALL dbt_get_block(t_3c_contr, ind, contr_blk, found)
5274 CALL dbt_get_block(t_3c_der_1(i_xyz), ind, der_blk_1, found_1)
5275 IF (.NOT. found_1)
THEN
5276 DEALLOCATE (der_blk_1)
5277 ALLOCATE (der_blk_1(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5278 der_blk_1(:, :, :) = 0.0_dp
5280 CALL dbt_get_block(t_3c_der_2(i_xyz), ind, der_blk_2, found_2)
5281 IF (.NOT. found_2)
THEN
5282 DEALLOCATE (der_blk_2)
5283 ALLOCATE (der_blk_2(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5284 der_blk_2(:, :, :) = 0.0_dp
5287 ALLOCATE (der_blk_3(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5288 der_blk_3(:, :, :) = -(der_blk_1(:, :, :) + der_blk_2(:, :, :))
5294 new_force = pref*sum(der_blk_1(:, :, :)*contr_blk(:, :, :))
5296 i_ri = (ind(1) - 1)/nblks_ri + 1
5297 ri_img = ri_data%RI_cell_to_img(i_ri)
5298 iat = idx_to_at_ri(ind(1) - (i_ri - 1)*nblks_ri)
5299 iat_of_kind = atom_of_kind(iat)
5300 ikind = kind_of(iat)
5303 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5306 IF (use_virial)
THEN
5309 scoord(:) = scoord(:) + real(index_to_cell(:, ri_img),
dp)
5313 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5318 new_force = pref*sum(der_blk_2(:, :, :)*contr_blk(:, :, :))
5319 jat = idx_to_at_ao(ind(2))
5320 jat_of_kind = atom_of_kind(jat)
5321 jkind = kind_of(jat)
5324 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5327 IF (use_virial)
THEN
5333 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5339 new_force = pref*sum(der_blk_3(:, :, :)*contr_blk(:, :, :))
5340 idx = (ind(3) - 1)/nblks_ao + 1
5341 kat = idx_to_at_ao(ind(3) - (
idx - 1)*nblks_ao)
5342 kat_of_kind = atom_of_kind(kat)
5343 kkind = kind_of(kat)
5346 force(kkind)%fock_4c(i_xyz, kat_of_kind) = force(kkind)%fock_4c(i_xyz, kat_of_kind) &
5349 IF (use_virial)
THEN
5351 scoord(:) = scoord(:) + real(index_to_cell(:, i_images(lb_img - 1 +
idx)),
dp)
5355 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5359 DEALLOCATE (der_blk_1, der_blk_2, der_blk_3)
5361 DEALLOCATE (contr_blk)
5364 CALL dbt_iterator_stop(iter)
5366 CALL timestop(handle)
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Types and set/get functions for auxiliary density matrix methods.
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind 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)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public bussy2024
Handles all functions related to the CELL.
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_distribution_release(dist)
...
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_clear(matrix)
...
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_power(matrix, exponent, threshold, n_dependent, para_env, blacs_env, verbose, eigenvectors, eigenvalues)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
This is the start of a dbt_api, all publically needed functions are exported here....
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
subroutine, public distribution_2d_create(distribution_2d, blacs_env, local_rows_ptr, n_local_rows, local_cols_ptr, row_distribution_ptr, col_distribution_ptr, n_local_cols, n_row_distribution, n_col_distribution)
initializes the distribution_2d
subroutine, public distribution_2d_release(distribution_2d)
...
RI-methods for HFX and K-points. \auhtor Augustin Bussy (01.2023)
subroutine, public hfx_ri_update_forces_kp(qs_env, ri_data, nspins, hf_fraction, rho_ao, use_virial)
Update the K-points RI-HFX forces.
subroutine, public hfx_ri_update_ks_kp(qs_env, ri_data, ks_matrix, ehfx, rho_ao, geometry_did_change, nspins, hf_fraction)
Update the KS matrices for each real-space image.
subroutine, public get_force_from_3c_trace(force, t_3c_contr, t_3c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, deriv_dim)
This routines calculates the force contribution from a trace over 3D tensors, i.e....
subroutine, public get_idx_to_atom(idx_to_at, bsizes_split, bsizes_orig)
a small utility function that returns the atom corresponding to a block of a split tensor
subroutine, public hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_ri, t_2c_int_pot, t_3c_int, do_kpoints)
Calculate 2-center and 3-center integrals.
subroutine, public get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, do_ovlp)
Update the forces due to the derivative of the a 2-center product d/dR (Q|R)
Types and set/get functions for HFX.
Routines useful for iterative matrix calculations.
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
integer, parameter, public default_string_length
Types and basic routines needed for a kpoint calculation.
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.
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
real(kind=dp), parameter, public cutoff_screen_factor
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Collection of simple mathematical functions and subroutines.
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public angstrom
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Define the quickstep kind type and their sub types.
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
module that contains the definitions of the scf types
Utility methods to build 3-center integral tensors of various types.
subroutine, public distribution_3d_create(dist_3d, dist1, dist2, dist3, nkind, particle_set, mp_comm_3d, own_comm)
Create a 3d distribution.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
subroutine, public create_tensor_batches(sizes, nbatches, starts_array, ends_array, starts_array_block, ends_array_block)
...
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Utility methods to build 3-center integral tensors of various types.
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.
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...
recursive integer function, public neighbor_list_3c_iterate(iterator)
Iterate 3c-nl iterator.
subroutine, public neighbor_list_3c_iterator_destroy(iterator)
Destroy 3c-nl iterator.
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
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.
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
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.
subroutine, public neighbor_list_3c_iterator_create(iterator, ijk_nl)
Create a 3-center neighborlist iterator.
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.
All kind of helpful little routines.
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
represent a pointer to a 1d array
represent a pointer to a 2d array
represent a pointer to a 3d array
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
distributes pairs on a 2d grid of processors
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.