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 .NE.
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) .NE. 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) .NE. 0) .AND. (mod(natom, ngroups) .NE. 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) .NE. 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 .LE. 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) .NE. 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) .NE. 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 .GE. r1) b = 0.0_dp
1711 IF (x .LE. 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 .GE. r1) b = 0.0_dp
1731 IF (x .LE. 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 .LE. 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 .LE. ri_data%kp_RI_range) present_atoms_i(jatom, j_img) = 1
2185 IF (iatom == atom_j .AND. dij .LE. 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 .NE. 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) .NE. 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) .NE. 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) .NE. 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) .NE. 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) .NE. 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) .NE. 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) .NE. 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 .NE. 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 .LE. 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) .NE. 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) .NE. 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 .LE. 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) .NE. 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.