46 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
47 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
48 dbt_distribution_destroy, dbt_distribution_new, dbt_distribution_type, dbt_filter, &
49 dbt_finalize, dbt_get_block, dbt_get_info, dbt_get_stored_coordinates, &
50 dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, &
51 dbt_iterator_type, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
52 dbt_pgrid_type, dbt_put_block, dbt_scale, dbt_type
114#include "./base/base_uses.f90"
123 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'hfx_ri_kp'
133 SUBROUTINE adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
134 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
135 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
136 TYPE(qs_environment_type),
POINTER :: qs_env
138 INTEGER :: i_img, i_RI, i_spin, iatom, natom, &
139 nblks_RI, nimg, nkind, nspins
140 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_RI_ext, dist1, dist2, dist3
141 TYPE(dft_control_type),
POINTER :: dft_control
142 TYPE(mp_para_env_type),
POINTER :: para_env
144 NULLIFY (dft_control, para_env)
151 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, para_env=para_env, nkind=nkind)
155 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
156 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
157 DO i_ri = 1, ri_data%ncell_RI
158 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
161 ALLOCATE (ri_data%t_3c_int_ctr_1(1, nimg))
163 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
164 ri_data%bsizes_AO_split, [1, 2], [3], name=
"(AO RI | AO)")
167 CALL dbt_create(ri_data%t_3c_int_ctr_1(1, 1), ri_data%t_3c_int_ctr_1(1, i_img))
169 DEALLOCATE (dist1, dist2, dist3)
171 ALLOCATE (ri_data%t_3c_int_ctr_2(1, 1))
173 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
174 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
175 DEALLOCATE (dist1, dist2, dist3)
178 DEALLOCATE (bsizes_ri_ext)
179 nblks_ri =
SIZE(ri_data%bsizes_RI)
180 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
181 DO i_ri = 1, ri_data%ncell_RI
182 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
185 ALLOCATE (ri_data%t_2c_inv(1, natom), ri_data%t_2c_int(1, natom), ri_data%t_2c_pot(1, natom))
186 CALL create_2c_tensor(ri_data%t_2c_inv(1, 1), dist1, dist2, ri_data%pgrid_2d, &
187 bsizes_ri_ext, bsizes_ri_ext, &
189 DEALLOCATE (dist1, dist2)
190 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, 1))
191 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, 1))
193 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_inv(1, iatom))
194 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, iatom))
195 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, iatom))
198 ALLOCATE (ri_data%kp_cost(natom, natom, nimg))
199 ri_data%kp_cost = 0.0_dp
202 nspins = dft_control%nspins
203 ALLOCATE (ri_data%rho_ao_t(nspins, nimg), ri_data%ks_t(nspins, nimg))
204 CALL create_2c_tensor(ri_data%rho_ao_t(1, 1), dist1, dist2, ri_data%pgrid_2d, &
205 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
207 DEALLOCATE (dist1, dist2)
209 CALL dbt_create(dbcsr_template, ri_data%ks_t(1, 1))
211 IF (nspins == 2)
THEN
212 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
213 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
216 DO i_spin = 1, nspins
217 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(i_spin, i_img))
218 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(i_spin, i_img))
222 END SUBROUTINE adapt_ri_data_to_kp
230 SUBROUTINE hfx_ri_pre_scf_kp(dbcsr_template, ri_data, qs_env)
231 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
232 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
233 TYPE(qs_environment_type),
POINTER :: qs_env
235 CHARACTER(LEN=*),
PARAMETER :: routineN =
'hfx_ri_pre_scf_kp'
237 INTEGER :: handle, i_img, iatom, natom, nimg, nkind
238 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_op_pot, t_2c_op_RI
239 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_int
240 TYPE(dft_control_type),
POINTER :: dft_control
242 NULLIFY (dft_control)
244 CALL timeset(routinen, handle)
246 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, nkind=nkind)
248 CALL cleanup_kp(ri_data)
251 IF (ri_data%flavor /=
ri_pmat) cpabort(
"K-points RI-HFX only with RHO flavor")
252 IF (ri_data%same_op) ri_data%same_op = .false.
253 IF (abs(ri_data%eps_pgf_orb - dft_control%qs_control%eps_pgf_orb) > 1.0e-16_dp)
THEN
254 cpabort(
"RI%EPS_PGF_ORB and QS%EPS_PGF_ORB must be identical for RI-HFX k-points")
257 CALL get_kp_and_ri_images(ri_data, qs_env)
261 ALLOCATE (t_2c_op_pot(nimg), t_2c_op_ri(nimg))
262 ALLOCATE (t_3c_int(1, nimg))
266 CALL adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
271 CALL get_ext_2c_int(ri_data%t_2c_inv(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
275 CALL get_ext_2c_int(ri_data%t_2c_int(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, &
276 qs_env, off_diagonal=.true.)
277 CALL apply_bump(ri_data%t_2c_int(1, iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
280 CALL get_ext_2c_int(ri_data%t_2c_pot(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
281 do_inverse=.true., skip_inverse=.true.)
282 CALL apply_bump(ri_data%t_2c_pot(1, iatom), iatom, ri_data, qs_env, from_left=.true., &
283 from_right=.true., debump=.true.)
291 ALLOCATE (ri_data%kp_mat_2c_pot(1, nimg))
293 CALL dbcsr_create(ri_data%kp_mat_2c_pot(1, i_img), template=t_2c_op_pot(i_img))
294 CALL dbcsr_copy(ri_data%kp_mat_2c_pot(1, i_img), t_2c_op_pot(i_img))
299 CALL reorder_3c_ints(t_3c_int(1, :), ri_data)
303 CALL precontract_3c_ints(t_3c_int, ri_data, qs_env)
305 CALL timestop(handle)
307 END SUBROUTINE hfx_ri_pre_scf_kp
313 SUBROUTINE cleanup_kp(ri_data)
314 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
318 IF (
ALLOCATED(ri_data%kp_cost))
DEALLOCATE (ri_data%kp_cost)
319 IF (
ALLOCATED(ri_data%idx_to_img))
DEALLOCATE (ri_data%idx_to_img)
320 IF (
ALLOCATED(ri_data%img_to_idx))
DEALLOCATE (ri_data%img_to_idx)
321 IF (
ALLOCATED(ri_data%present_images))
DEALLOCATE (ri_data%present_images)
322 IF (
ALLOCATED(ri_data%img_to_RI_cell))
DEALLOCATE (ri_data%img_to_RI_cell)
323 IF (
ALLOCATED(ri_data%RI_cell_to_img))
DEALLOCATE (ri_data%RI_cell_to_img)
325 IF (
ALLOCATED(ri_data%kp_mat_2c_pot))
THEN
326 DO j = 1,
SIZE(ri_data%kp_mat_2c_pot, 2)
327 DO i = 1,
SIZE(ri_data%kp_mat_2c_pot, 1)
331 DEALLOCATE (ri_data%kp_mat_2c_pot)
334 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
335 DO i = 1,
SIZE(ri_data%kp_t_3c_int)
336 CALL dbt_destroy(ri_data%kp_t_3c_int(i))
338 DEALLOCATE (ri_data%kp_t_3c_int)
341 IF (
ALLOCATED(ri_data%t_2c_inv))
THEN
342 DO j = 1,
SIZE(ri_data%t_2c_inv, 2)
343 DO i = 1,
SIZE(ri_data%t_2c_inv, 1)
344 CALL dbt_destroy(ri_data%t_2c_inv(i, j))
347 DEALLOCATE (ri_data%t_2c_inv)
350 IF (
ALLOCATED(ri_data%t_2c_int))
THEN
351 DO j = 1,
SIZE(ri_data%t_2c_int, 2)
352 DO i = 1,
SIZE(ri_data%t_2c_int, 1)
353 CALL dbt_destroy(ri_data%t_2c_int(i, j))
356 DEALLOCATE (ri_data%t_2c_int)
359 IF (
ALLOCATED(ri_data%t_2c_pot))
THEN
360 DO j = 1,
SIZE(ri_data%t_2c_pot, 2)
361 DO i = 1,
SIZE(ri_data%t_2c_pot, 1)
362 CALL dbt_destroy(ri_data%t_2c_pot(i, j))
365 DEALLOCATE (ri_data%t_2c_pot)
368 IF (
ALLOCATED(ri_data%t_3c_int_ctr_1))
THEN
369 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_1, 2)
370 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_1, 1)
371 CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
374 DEALLOCATE (ri_data%t_3c_int_ctr_1)
377 IF (
ALLOCATED(ri_data%t_3c_int_ctr_2))
THEN
378 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_2, 2)
379 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_2, 1)
380 CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
383 DEALLOCATE (ri_data%t_3c_int_ctr_2)
386 IF (
ALLOCATED(ri_data%rho_ao_t))
THEN
387 DO j = 1,
SIZE(ri_data%rho_ao_t, 2)
388 DO i = 1,
SIZE(ri_data%rho_ao_t, 1)
389 CALL dbt_destroy(ri_data%rho_ao_t(i, j))
392 DEALLOCATE (ri_data%rho_ao_t)
395 IF (
ALLOCATED(ri_data%ks_t))
THEN
396 DO j = 1,
SIZE(ri_data%ks_t, 2)
397 DO i = 1,
SIZE(ri_data%ks_t, 1)
398 CALL dbt_destroy(ri_data%ks_t(i, j))
401 DEALLOCATE (ri_data%ks_t)
404 END SUBROUTINE cleanup_kp
413 SUBROUTINE print_progress_bar(b_img, nimg, iprint, ri_data)
414 INTEGER,
INTENT(IN) :: b_img, nimg
415 INTEGER,
INTENT(INOUT) :: iprint
416 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
418 CHARACTER(LEN=default_string_length) :: bar
421 IF (ri_data%unit_nr > 0)
THEN
423 WRITE (ri_data%unit_nr,
'(/T6,A)', advance=
"no")
'[-'
426 IF (b_img > iprint*nimg/71)
THEN
427 rep = max(1, 71/nimg)
428 bar = repeat(
"-", rep)
429 WRITE (ri_data%unit_nr,
'(A)', advance=
"no") trim(bar)
433 IF (b_img == nimg)
THEN
434 rep = max(0, 1 + 71 - iprint*rep)
435 bar = repeat(
"-", rep)
436 WRITE (ri_data%unit_nr,
'(A,A)') trim(bar),
'-]'
441 END SUBROUTINE print_progress_bar
455 geometry_did_change, nspins, hf_fraction)
459 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
460 REAL(kind=
dp),
INTENT(OUT) :: ehfx
462 LOGICAL,
INTENT(IN) :: geometry_did_change
463 INTEGER,
INTENT(IN) :: nspins
464 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
466 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_ks_kp'
468 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_spin, iatom, &
469 iblk, igroup, iprint, jatom, mb_img, n_batch_nze, n_nze, natom, ngroups, nimg, nimg_nze
470 INTEGER(int_8) :: mem, nflop, nze
471 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges_at, batch_ranges_nze, &
473 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
474 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: sparsity_pattern
475 LOGICAL :: estimate_mem, print_progress, use_delta_p
476 REAL(
dp) :: etmp,
fac, occ, pfac, pref, t1, t2, t3, &
479 TYPE(
dbcsr_type) :: ks_desymm, rho_desymm, tmp
480 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
482 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ks_t_split, t_2c_ao_tmp, t_2c_work, &
483 t_3c_int, t_3c_work_2, t_3c_work_3
484 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: ks_t, ks_t_sub, t_3c_apc, t_3c_apc_sub
488 NULLIFY (para_env, para_env_sub, blacs_env_sub, hfx_section, dbcsr_template, print_section)
492 CALL timeset(routinen, handle)
494 CALL get_qs_env(qs_env, para_env=para_env, natom=natom)
496 IF (nspins == 1)
THEN
497 fac = 0.5_dp*hf_fraction
499 fac = 1.0_dp*hf_fraction
506 ri_data%kp_stack_size = batch_size
507 ri_data%kp_ngroups = ngroups
509 IF (geometry_did_change)
THEN
510 CALL hfx_ri_pre_scf_kp(ks_matrix(1, 1)%matrix, ri_data, qs_env)
513 nimg_nze = ri_data%nimg_nze
536 IF (.NOT. use_delta_p) pfac = 0.0_dp
537 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, pfac, ri_data, qs_env)
541 DO i_spin = 1, nspins
548 IF (n_nze == nspins)
THEN
549 cpwarn(
"It is highly recommended to restart from a converged GGA K-point calculations.")
552 ALLOCATE (ks_t(nspins, nimg))
554 DO i_spin = 1, nspins
555 CALL dbt_create(ri_data%ks_t(1, 1), ks_t(i_spin, i_img))
559 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
560 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
568 ALLOCATE (t_3c_apc(nspins, nimg))
570 DO i_spin = 1, nspins
571 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
574 CALL contract_pmat_3c(t_3c_apc, ri_data%rho_ao_t, ri_data, qs_env)
576 IF (mod(para_env%num_pe, ngroups) /= 0)
THEN
577 cpwarn(
"KP_NGROUPS must be an integer divisor of the total number of MPI ranks. It was set to 1.")
581 IF ((mod(ngroups, natom) /= 0) .AND. (mod(natom, ngroups) /= 0) .AND. geometry_did_change)
THEN
582 IF (ngroups > 1)
THEN
583 cpwarn(
"Better load balancing is reached if NGROUPS is a multiple/divisor of the number of atoms")
586 group_size = para_env%num_pe/ngroups
587 igroup = para_env%mepos/group_size
589 ALLOCATE (para_env_sub)
590 CALL para_env_sub%from_split(para_env, igroup)
594 ALLOCATE (sparsity_pattern(natom, natom, nimg))
595 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
596 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
599 ALLOCATE (mat_2c_pot(nimg), ks_t_sub(nspins, nimg), t_2c_ao_tmp(1), ks_t_split(2), t_2c_work(3))
600 CALL get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
601 group_size, ngroups, para_env, para_env_sub, ri_data)
603 ALLOCATE (t_3c_int(nimg), t_3c_apc_sub(nspins, nimg), t_3c_work_2(3), t_3c_work_3(3))
604 CALL get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
605 group_size, ngroups, para_env, para_env_sub, ri_data)
609 ALLOCATE (batch_ranges_at(natom + 1))
610 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
612 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
613 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
615 batch_ranges_at(iatom) = iblk
619 n_batch_nze = nimg_nze/batch_size
620 IF (
modulo(nimg_nze, batch_size) /= 0) n_batch_nze = n_batch_nze + 1
621 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
622 DO i_batch = 1, n_batch_nze
623 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
625 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
631 ALLOCATE (iapc_pairs(nimg, 2))
632 IF (estimate_mem .AND. geometry_did_change)
THEN
634 CALL get_iapc_pairs(iapc_pairs, 1, ri_data, qs_env)
635 CALL fill_3c_stack(t_3c_work_3(1), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
636 filter_at=1, filter_dim=2, idx_to_at=idx_to_at_ao, &
637 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
638 CALL fill_3c_stack(t_3c_work_3(2), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
639 filter_at=1, filter_dim=2, idx_to_at=idx_to_at_ao, &
640 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
641 CALL fill_3c_stack(t_3c_work_2(1), t_3c_apc_sub(1, :), iapc_pairs(:, 2), 3, &
642 ri_data, filter_at=1, filter_dim=1, idx_to_at=idx_to_at_ao, &
643 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
644 CALL fill_3c_stack(t_3c_work_2(2), t_3c_apc_sub(1, :), iapc_pairs(:, 2), 3, &
645 ri_data, filter_at=1, filter_dim=1, idx_to_at=idx_to_at_ao, &
646 img_bounds=[batch_ranges_nze(1), batch_ranges_nze(2)])
647 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, 1, 1, 1, ri_data, qs_env, &
648 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
649 dbcsr_template=dbcsr_template)
651 CALL para_env%max(mem)
652 CALL dbt_clear(t_3c_work_2(1))
653 CALL dbt_clear(t_3c_work_2(2))
654 CALL dbt_clear(t_3c_work_3(1))
655 CALL dbt_clear(t_3c_work_3(2))
656 CALL dbt_clear(t_2c_work(1))
658 IF (ri_data%unit_nr > 0)
THEN
659 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I14)") &
660 "KP-HFX_RI_INFO| Estimated peak memory usage per MPI rank (MiB):", mem/(1024*1024)
665 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
666 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
667 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
668 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
672 ri_data%kp_cost(:, :, :) = 0.0_dp
674 IF (print_progress)
CALL print_progress_bar(b_img, nimg, iprint, ri_data)
675 CALL dbt_batched_contract_init(ks_t_split(1))
676 CALL dbt_batched_contract_init(ks_t_split(2))
679 IF (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup) cycle
681 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp
687 CALL timeset(routinen//
"_2c", handle2)
688 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
689 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
690 dbcsr_template=dbcsr_template)
691 CALL dbt_copy(t_2c_work(1), t_2c_work(2), move_data=.true.)
692 CALL dbt_filter(t_2c_work(2), ri_data%filter_eps)
693 CALL timestop(handle2)
695 CALL dbt_batched_contract_init(t_2c_work(2))
696 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env)
697 CALL timeset(routinen//
"_3c", handle2)
700 DO i_batch = 1, n_batch_nze
701 CALL fill_3c_stack(t_3c_work_3(3), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
702 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
703 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
704 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1), move_data=.true.)
706 CALL dbt_contract(1.0_dp, t_2c_work(2), t_3c_work_3(1), &
707 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
708 contract_1=[2], notcontract_1=[1], &
709 contract_2=[1], notcontract_2=[2, 3], &
710 filter_eps=ri_data%filter_eps, flop=nflop)
711 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
712 CALL dbt_copy(t_3c_work_3(2), t_3c_work_2(2), order=[2, 1, 3], move_data=.true.)
713 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1))
717 DO i_spin = 1, nspins
718 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
719 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
720 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
724 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
725 CALL dbt_contract(-pref*
fac, t_3c_work_2(1), t_3c_work_2(2), &
726 1.0_dp, ks_t_split(i_spin), map_1=[1], map_2=[2], &
727 contract_1=[2, 3], notcontract_1=[1], &
728 contract_2=[2, 3], notcontract_2=[1], &
729 filter_eps=ri_data%filter_eps, &
730 move_data=i_spin == nspins, flop=nflop)
731 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
734 CALL timestop(handle2)
735 CALL dbt_batched_contract_finalize(t_2c_work(2))
738 ri_data%kp_cost(iatom, jatom, b_img) = t4 - t3
741 CALL dbt_batched_contract_finalize(ks_t_split(1))
742 CALL dbt_batched_contract_finalize(ks_t_split(2))
744 DO i_spin = 1, nspins
745 CALL dbt_copy(ks_t_split(i_spin), t_2c_ao_tmp(1), move_data=.true.)
746 CALL dbt_copy(t_2c_ao_tmp(1), ks_t_sub(i_spin, b_img), summation=.true.)
749 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
750 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
751 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
752 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
754 CALL para_env%sum(ri_data%dbcsr_nflop)
755 CALL para_env%sum(ri_data%kp_cost)
757 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
760 CALL gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
764 CALL dbt_copy(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img), move_data=.true.)
768 CALL dbt_destroy(t_2c_ao_tmp(1))
769 CALL dbt_destroy(ks_t_split(1))
770 CALL dbt_destroy(ks_t_split(2))
771 CALL dbt_destroy(t_2c_work(1))
772 CALL dbt_destroy(t_2c_work(2))
773 CALL dbt_destroy(t_3c_work_2(1))
774 CALL dbt_destroy(t_3c_work_2(2))
775 CALL dbt_destroy(t_3c_work_2(3))
776 CALL dbt_destroy(t_3c_work_3(1))
777 CALL dbt_destroy(t_3c_work_3(2))
778 CALL dbt_destroy(t_3c_work_3(3))
780 CALL dbt_destroy(t_3c_int(i_img))
782 DO i_spin = 1, nspins
783 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
784 CALL dbt_destroy(ks_t_sub(i_spin, i_img))
787 IF (
ASSOCIATED(dbcsr_template))
THEN
789 DEALLOCATE (dbcsr_template)
794 CALL para_env_sub%free()
795 DEALLOCATE (para_env_sub)
800 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
801 DO i_spin = 1, nspins
803 CALL dbt_copy(ks_t(i_spin, b_img), ri_data%ks_t(i_spin, b_img), summation=.true.)
806 mb_img = get_opp_index(b_img, qs_env)
807 IF (mb_img > 0 .AND. mb_img <= nimg)
THEN
808 CALL dbt_copy(ks_t(i_spin, mb_img), ri_data%ks_t(i_spin, b_img), order=[2, 1], summation=.true.)
813 DO i_spin = 1, nspins
814 CALL dbt_destroy(ks_t(i_spin, b_img))
819 CALL dbt_create(ri_data%ks_t(1, 1), t_2c_ao_tmp(1))
820 CALL dbcsr_create(tmp, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
821 CALL dbcsr_create(ks_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
822 CALL dbcsr_create(rho_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
825 DO i_spin = 1, nspins
826 CALL dbt_filter(ri_data%ks_t(i_spin, i_img), ri_data%filter_eps)
827 CALL dbt_copy(ri_data%ks_t(i_spin, i_img), t_2c_ao_tmp(1))
828 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), ks_desymm)
829 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), tmp)
830 CALL dbcsr_add(ks_matrix(i_spin, i_img)%matrix, tmp, 1.0_dp, 1.0_dp)
832 CALL dbt_copy(ri_data%rho_ao_t(i_spin, i_img), t_2c_ao_tmp(1))
833 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), rho_desymm)
835 CALL dbcsr_dot(ks_desymm, rho_desymm, etmp)
836 ehfx = ehfx + 0.5_dp*etmp
838 IF (.NOT. use_delta_p)
CALL dbt_clear(ri_data%ks_t(i_spin, i_img))
844 CALL dbt_destroy(t_2c_ao_tmp(1))
846 CALL timestop(handle)
865 INTEGER,
INTENT(IN) :: nspins
866 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
868 LOGICAL,
INTENT(IN),
OPTIONAL :: use_virial
870 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_forces_kp'
872 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_loop, i_spin, &
873 i_xyz, iatom, iblk, igroup, j_xyz, jatom, k_xyz, n_batch, natom, ngroups, nimg, nimg_nze
874 INTEGER(int_8) :: nflop, nze
875 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, batch_ranges_at, &
876 batch_ranges_nze, dist1, dist2, &
877 i_images, idx_to_at_ao, idx_to_at_ri, &
879 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
880 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: force_pattern, sparsity_pattern
881 INTEGER,
DIMENSION(2, 1) :: bounds_iat, bounds_jat
882 LOGICAL :: use_virial_prv
883 REAL(
dp) ::
fac, occ, pref, t1, t2
884 REAL(
dp),
DIMENSION(3, 3) :: work_virial
888 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
889 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_pot, mat_der_pot_sub
891 TYPE(dbt_type) :: t_2c_r, t_2c_r_split
892 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_bint, t_2c_binv, t_2c_der_pot, &
893 t_2c_inv, t_2c_metric, t_2c_work, &
894 t_3c_der_stack, t_3c_work_2, &
896 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
897 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, &
905 NULLIFY (para_env, para_env_sub, hfx_section, blacs_env_sub, dbcsr_template, force, atomic_kind_set, &
906 virial, particle_set, cell)
908 CALL timeset(routinen, handle)
910 use_virial_prv = .false.
911 IF (
PRESENT(use_virial)) use_virial_prv = use_virial
913 IF (nspins == 1)
THEN
914 fac = 0.5_dp*hf_fraction
916 fac = 1.0_dp*hf_fraction
919 CALL get_qs_env(qs_env, natom=natom, para_env=para_env, force=force, cell=cell, virial=virial, &
920 atomic_kind_set=atomic_kind_set, particle_set=particle_set)
923 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
924 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
926 ALLOCATE (idx_to_at_ri(
SIZE(ri_data%bsizes_RI_split)))
927 CALL get_idx_to_atom(idx_to_at_ri, ri_data%bsizes_RI_split, ri_data%bsizes_RI)
930 ALLOCATE (t_3c_der_ri(nimg, 3), t_3c_der_ao(nimg, 3), mat_der_pot(nimg, 3), t_2c_der_metric(natom, 3))
934 CALL precalc_derivatives(t_3c_der_ri, t_3c_der_ao, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
937 ALLOCATE (rho_ao_t(nspins, nimg))
939 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
941 DEALLOCATE (dist1, dist2)
942 IF (nspins == 2)
CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(2, 1))
944 DO i_spin = 1, nspins
945 CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(i_spin, i_img))
948 CALL get_pmat_images(rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
951 ALLOCATE (t_3c_apc(nspins, nimg))
953 DO i_spin = 1, nspins
954 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
957 CALL contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
962 group_size = para_env%num_pe/ngroups
963 igroup = para_env%mepos/group_size
965 ALLOCATE (para_env_sub)
966 CALL para_env_sub%from_split(para_env, igroup)
970 ALLOCATE (sparsity_pattern(natom, natom, nimg))
971 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
972 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
975 ALLOCATE (t_2c_inv(natom), mat_2c_pot(nimg), rho_ao_t_sub(nspins, nimg), t_2c_work(5), &
976 t_2c_der_metric_sub(natom, 3), mat_der_pot_sub(nimg, 3), t_2c_bint(natom), &
977 t_2c_metric(natom), t_2c_binv(natom))
978 CALL get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
979 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
980 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
981 CALL dbt_create(t_2c_work(1), t_2c_r)
982 CALL dbt_create(t_2c_work(5), t_2c_r_split)
984 ALLOCATE (t_2c_der_pot(3))
986 CALL dbt_create(t_2c_r, t_2c_der_pot(i_xyz))
990 ALLOCATE (t_3c_work_2(3), t_3c_work_3(4), t_3c_der_stack(6), t_3c_der_ao_sub(nimg, 3), &
991 t_3c_der_ri_sub(nimg, 3), t_3c_apc_sub(nspins, nimg))
992 CALL get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_ao, t_3c_der_ao_sub, &
993 t_3c_der_ri, t_3c_der_ri_sub, t_3c_apc, t_3c_apc_sub, t_3c_der_stack, &
994 group_size, ngroups, para_env, para_env_sub, ri_data)
997 ALLOCATE (batch_ranges_at(natom + 1))
998 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
1000 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
1001 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
1003 batch_ranges_at(iatom) = iblk
1007 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
1008 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
1009 CALL dbt_batched_contract_init(t_3c_work_3(3), batch_range_2=batch_ranges_at)
1010 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
1011 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
1014 nimg_nze = ri_data%nimg_nze
1015 batch_size = ri_data%kp_stack_size
1016 n_batch = nimg_nze/batch_size
1017 IF (
modulo(nimg_nze, batch_size) /= 0) n_batch = n_batch + 1
1018 ALLOCATE (batch_ranges_nze(n_batch + 1))
1019 DO i_batch = 1, n_batch
1020 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
1022 batch_ranges_nze(n_batch + 1) = nimg_nze + 1
1027 CALL dbt_create(t_2c_inv(iatom), t_2c_binv(iatom))
1028 CALL dbt_copy(t_2c_inv(iatom), t_2c_binv(iatom))
1029 CALL apply_bump(t_2c_binv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
1030 CALL apply_bump(t_2c_inv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1034 work_virial = 0.0_dp
1035 ALLOCATE (iapc_pairs(nimg, 2), i_images(nimg))
1036 ALLOCATE (force_pattern(natom, natom, nimg))
1037 force_pattern(:, :, :) = -1
1046 IF (i_loop == 1 .AND. (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup)) cycle
1047 IF (i_loop == 2 .AND. (.NOT. force_pattern(iatom, jatom, b_img) == igroup)) cycle
1050 CALL timeset(routinen//
"_2c_1", handle2)
1051 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
1052 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
1053 dbcsr_template=dbcsr_template)
1054 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_inv(jatom), &
1055 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1056 contract_1=[2], notcontract_1=[1], &
1057 contract_2=[1], notcontract_2=[2], &
1058 filter_eps=ri_data%filter_eps, flop=nflop)
1059 CALL dbt_copy(t_2c_work(2), t_2c_work(5), move_data=.true.)
1060 CALL dbt_filter(t_2c_work(5), ri_data%filter_eps)
1061 CALL timestop(handle2)
1063 CALL timeset(routinen//
"_3c", handle2)
1064 bounds_iat(:, 1) = [sum(ri_data%bsizes_AO(1:iatom - 1)) + 1, sum(ri_data%bsizes_AO(1:iatom))]
1065 bounds_jat(:, 1) = [sum(ri_data%bsizes_AO(1:jatom - 1)) + 1, sum(ri_data%bsizes_AO(1:jatom))]
1066 CALL dbt_clear(t_2c_r_split)
1068 DO i_spin = 1, nspins
1069 CALL dbt_batched_contract_init(rho_ao_t_sub(i_spin, b_img))
1072 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env, i_images)
1073 DO i_batch = 1, n_batch
1077 CALL dbt_clear(t_3c_der_stack(i_xyz))
1078 CALL fill_3c_stack(t_3c_der_stack(i_xyz), t_3c_der_ri_sub(:, i_xyz), &
1079 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
1080 filter_dim=2, idx_to_at=idx_to_at_ao, &
1081 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1083 CALL dbt_clear(t_3c_der_stack(3 + i_xyz))
1084 CALL fill_3c_stack(t_3c_der_stack(3 + i_xyz), t_3c_der_ao_sub(:, i_xyz), &
1085 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
1086 filter_dim=2, idx_to_at=idx_to_at_ao, &
1087 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1090 DO i_spin = 1, nspins
1092 CALL dbt_clear(t_3c_work_2(3))
1093 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
1094 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
1095 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1098 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
1102 CALL dbt_contract(1.0_dp, rho_ao_t_sub(i_spin, b_img), t_3c_work_2(1), &
1103 0.0_dp, t_3c_work_2(2), map_1=[1], map_2=[2, 3], &
1104 contract_1=[1], notcontract_1=[2], &
1105 contract_2=[1], notcontract_2=[2, 3], &
1106 bounds_1=bounds_iat, bounds_2=bounds_jat, &
1107 filter_eps=ri_data%filter_eps, flop=nflop)
1108 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1114 CALL dbt_copy(t_3c_work_2(2), t_3c_work_3(1), order=[2, 1, 3], move_data=.true.)
1115 CALL dbt_batched_contract_init(t_2c_work(5))
1116 CALL dbt_contract(1.0_dp, t_2c_work(5), t_3c_work_3(1), &
1117 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
1118 contract_1=[1], notcontract_1=[2], &
1119 contract_2=[1], notcontract_2=[2, 3], &
1120 filter_eps=ri_data%filter_eps, flop=nflop)
1121 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1122 CALL dbt_batched_contract_finalize(t_2c_work(5))
1125 CALL dbt_copy(t_3c_work_3(2), t_3c_work_3(4), move_data=.true.)
1126 IF (use_virial_prv)
THEN
1128 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1129 idx_to_at_ri, idx_to_at_ao, i_images, &
1130 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1131 ri_data, qs_env, work_virial, cell, particle_set)
1134 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1135 idx_to_at_ri, idx_to_at_ao, i_images, &
1136 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1139 CALL dbt_clear(t_3c_work_3(4))
1143 IF (i_loop == 2) cycle
1146 CALL fill_3c_stack(t_3c_work_3(4), ri_data%kp_t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
1147 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
1148 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1149 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(3), move_data=.true.)
1151 CALL dbt_batched_contract_init(t_2c_r_split)
1152 CALL dbt_contract(1.0_dp, t_3c_work_3(1), t_3c_work_3(3), &
1153 1.0_dp, t_2c_r_split, map_1=[1], map_2=[2], &
1154 contract_1=[2, 3], notcontract_1=[1], &
1155 contract_2=[2, 3], notcontract_2=[1], &
1156 filter_eps=ri_data%filter_eps, flop=nflop)
1157 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1158 CALL dbt_batched_contract_finalize(t_2c_r_split)
1159 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(1))
1162 DO i_spin = 1, nspins
1163 CALL dbt_batched_contract_finalize(rho_ao_t_sub(i_spin, b_img))
1165 CALL timestop(handle2)
1167 IF (i_loop == 2) cycle
1169 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp*pref
1171 CALL timeset(routinen//
"_2c_2", handle2)
1173 CALL dbt_copy(t_2c_r_split, t_2c_r, move_data=.true.)
1175 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
1176 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
1177 dbcsr_template=dbcsr_template)
1190 CALL get_ext_2c_int(t_2c_der_pot(i_xyz), mat_der_pot_sub(:, i_xyz), iatom, jatom, &
1191 b_img, ri_data, qs_env, blacs_env_ext=blacs_env_sub, &
1192 para_env_ext=para_env_sub, dbcsr_template=dbcsr_template)
1195 IF (use_virial_prv)
THEN
1196 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1197 b_img, pref, ri_data, qs_env, work_virial, cell, particle_set)
1199 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1200 b_img, pref, ri_data, qs_env)
1204 CALL dbt_clear(t_2c_der_pot(i_xyz))
1208 CALL dbt_contract(1.0_dp, t_2c_metric(iatom), t_2c_r, &
1209 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1210 contract_1=[2], notcontract_1=[1], &
1211 contract_2=[1], notcontract_2=[2], &
1212 filter_eps=ri_data%filter_eps, flop=nflop)
1213 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1214 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_work(1), &
1215 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1216 contract_1=[2], notcontract_1=[1], &
1217 contract_2=[2], notcontract_2=[1], &
1218 filter_eps=ri_data%filter_eps, flop=nflop)
1219 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1223 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(iatom), &
1224 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1225 contract_1=[2], notcontract_1=[1], &
1226 contract_2=[1], notcontract_2=[2], &
1227 filter_eps=ri_data%filter_eps, flop=nflop)
1228 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1230 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(3), &
1231 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1232 contract_1=[1], notcontract_1=[2], &
1233 contract_2=[1], notcontract_2=[2], &
1234 filter_eps=ri_data%filter_eps, flop=nflop)
1235 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1237 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1238 CALL get_2c_bump_forces(force, t_2c_work(4), iatom, atom_of_kind, kind_of, pref, &
1239 ri_data, qs_env, work_virial)
1242 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(2), &
1243 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1244 contract_1=[1], notcontract_1=[2], &
1245 contract_2=[1], notcontract_2=[2], &
1246 filter_eps=ri_data%filter_eps, flop=nflop)
1247 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1249 IF (use_virial_prv)
THEN
1250 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1251 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1252 diag=.true., offdiag=.false.)
1254 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1255 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1259 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1260 CALL apply_bump(t_2c_work(2), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1262 IF (use_virial_prv)
THEN
1263 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1264 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1265 diag=.false., offdiag=.true.)
1267 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1268 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1273 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(iatom), &
1274 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1275 contract_1=[2], notcontract_1=[1], &
1276 contract_2=[1], notcontract_2=[2], &
1277 filter_eps=ri_data%filter_eps, flop=nflop)
1278 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1280 CALL dbt_contract(1.0_dp, t_2c_bint(iatom), t_2c_work(4), &
1281 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1282 contract_1=[1], notcontract_1=[2], &
1283 contract_2=[1], notcontract_2=[2], &
1284 filter_eps=ri_data%filter_eps, flop=nflop)
1285 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1287 CALL get_2c_bump_forces(force, t_2c_work(2), iatom, atom_of_kind, kind_of, -pref, &
1288 ri_data, qs_env, work_virial)
1291 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_r, &
1292 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1293 contract_1=[1], notcontract_1=[2], &
1294 contract_2=[1], notcontract_2=[2], &
1295 filter_eps=ri_data%filter_eps, flop=nflop)
1296 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1298 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_metric(jatom), &
1299 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1300 contract_1=[2], notcontract_1=[1], &
1301 contract_2=[1], notcontract_2=[2], &
1302 filter_eps=ri_data%filter_eps, flop=nflop)
1303 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1307 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(jatom), &
1308 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1309 contract_1=[2], notcontract_1=[1], &
1310 contract_2=[1], notcontract_2=[2], &
1311 filter_eps=ri_data%filter_eps, flop=nflop)
1312 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1314 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(3), &
1315 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1316 contract_1=[1], notcontract_1=[2], &
1317 contract_2=[1], notcontract_2=[2], &
1318 filter_eps=ri_data%filter_eps, flop=nflop)
1319 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1321 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1322 CALL get_2c_bump_forces(force, t_2c_work(4), jatom, atom_of_kind, kind_of, pref, &
1323 ri_data, qs_env, work_virial)
1326 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(2), &
1327 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1328 contract_1=[1], notcontract_1=[2], &
1329 contract_2=[1], notcontract_2=[2], &
1330 filter_eps=ri_data%filter_eps, flop=nflop)
1331 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1333 IF (use_virial_prv)
THEN
1334 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1335 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1336 diag=.true., offdiag=.false.)
1338 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1339 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1343 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1344 CALL apply_bump(t_2c_work(2), jatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1346 IF (use_virial_prv)
THEN
1347 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1348 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1349 diag=.false., offdiag=.true.)
1351 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1352 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1357 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(jatom), &
1358 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1359 contract_1=[2], notcontract_1=[1], &
1360 contract_2=[1], notcontract_2=[2], &
1361 filter_eps=ri_data%filter_eps, flop=nflop)
1362 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1364 CALL dbt_contract(1.0_dp, t_2c_bint(jatom), t_2c_work(4), &
1365 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1366 contract_1=[1], notcontract_1=[2], &
1367 contract_2=[1], notcontract_2=[2], &
1368 filter_eps=ri_data%filter_eps, flop=nflop)
1369 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1371 CALL get_2c_bump_forces(force, t_2c_work(2), jatom, atom_of_kind, kind_of, -pref, &
1372 ri_data, qs_env, work_virial)
1374 CALL timestop(handle2)
1379 IF (i_loop == 1)
THEN
1380 CALL update_pattern_to_forces(force_pattern, sparsity_pattern, ngroups, ri_data, qs_env)
1384 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
1385 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
1386 CALL dbt_batched_contract_finalize(t_3c_work_3(3))
1387 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
1388 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
1390 IF (use_virial_prv)
THEN
1394 virial%pv_fock_4c(i_xyz, j_xyz) = virial%pv_fock_4c(i_xyz, j_xyz) &
1395 + work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1403 CALL para_env_sub%free()
1404 DEALLOCATE (para_env_sub)
1406 CALL para_env%sync()
1408 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
1411 IF (
ASSOCIATED(dbcsr_template))
THEN
1413 DEALLOCATE (dbcsr_template)
1415 CALL dbt_destroy(t_2c_r)
1416 CALL dbt_destroy(t_2c_r_split)
1417 CALL dbt_destroy(t_2c_work(1))
1418 CALL dbt_destroy(t_2c_work(2))
1419 CALL dbt_destroy(t_2c_work(3))
1420 CALL dbt_destroy(t_2c_work(4))
1421 CALL dbt_destroy(t_2c_work(5))
1422 CALL dbt_destroy(t_3c_work_2(1))
1423 CALL dbt_destroy(t_3c_work_2(2))
1424 CALL dbt_destroy(t_3c_work_2(3))
1425 CALL dbt_destroy(t_3c_work_3(1))
1426 CALL dbt_destroy(t_3c_work_3(2))
1427 CALL dbt_destroy(t_3c_work_3(3))
1428 CALL dbt_destroy(t_3c_work_3(4))
1429 CALL dbt_destroy(t_3c_der_stack(1))
1430 CALL dbt_destroy(t_3c_der_stack(2))
1431 CALL dbt_destroy(t_3c_der_stack(3))
1432 CALL dbt_destroy(t_3c_der_stack(4))
1433 CALL dbt_destroy(t_3c_der_stack(5))
1434 CALL dbt_destroy(t_3c_der_stack(6))
1436 CALL dbt_destroy(t_2c_der_pot(i_xyz))
1439 CALL dbt_destroy(t_2c_inv(iatom))
1440 CALL dbt_destroy(t_2c_binv(iatom))
1441 CALL dbt_destroy(t_2c_bint(iatom))
1442 CALL dbt_destroy(t_2c_metric(iatom))
1444 CALL dbt_destroy(t_2c_der_metric_sub(iatom, i_xyz))
1449 DO i_spin = 1, nspins
1450 CALL dbt_destroy(rho_ao_t_sub(i_spin, i_img))
1451 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
1456 CALL dbt_destroy(t_3c_der_ri_sub(i_img, i_xyz))
1457 CALL dbt_destroy(t_3c_der_ao_sub(i_img, i_xyz))
1462 CALL timestop(handle)
1477 SUBROUTINE apply_bump(t_2c_inout, atom_i, ri_data, qs_env, from_left, from_right, debump)
1478 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_inout
1479 INTEGER,
INTENT(IN) :: atom_i
1482 LOGICAL,
INTENT(IN),
OPTIONAL :: from_left, from_right, debump
1484 INTEGER :: i_img, i_ri, iatom, ind(2), j_img, j_ri, &
1485 jatom, natom, nblks(2), nimg, nkind
1486 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1487 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1488 LOGICAL :: found, my_debump, my_left, my_right
1489 REAL(
dp) :: bval, r0, r1, ri(3), rj(3), rref(3), &
1491 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1493 TYPE(dbt_iterator_type) :: iter
1496 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1498 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1500 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1501 kpoints=kpoints, particle_set=particle_set)
1502 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1505 IF (
PRESENT(debump)) my_debump = debump
1508 IF (
PRESENT(from_left)) my_left = from_left
1511 IF (
PRESENT(from_right)) my_right = from_right
1512 cpassert(my_left .OR. my_right)
1514 CALL dbt_get_info(t_2c_inout, nblks_total=nblks)
1515 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1516 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1521 r1 = ri_data%kp_RI_range
1522 r0 = ri_data%kp_bump_rad
1523 rref =
pbc(particle_set(atom_i)%r, cell)
1528 CALL dbt_iterator_start(iter, t_2c_inout)
1529 DO WHILE (dbt_iterator_blocks_left(iter))
1530 CALL dbt_iterator_next_block(iter, ind)
1531 CALL dbt_get_block(t_2c_inout, ind, blk, found)
1532 IF (.NOT. found) cycle
1534 i_ri = (ind(1) - 1)/natom + 1
1535 i_img = ri_data%RI_cell_to_img(i_ri)
1536 iatom = ind(1) - (i_ri - 1)*natom
1539 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
1541 j_ri = (ind(2) - 1)/natom + 1
1542 j_img = ri_data%RI_cell_to_img(j_ri)
1543 jatom = ind(2) - (j_ri - 1)*natom
1546 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1548 IF (.NOT. my_debump)
THEN
1549 IF (my_left) blk(:, :) = blk(:, :)*bump(norm2(ri - rref), r0, r1)
1550 IF (my_right) blk(:, :) = blk(:, :)*bump(norm2(rj - rref), r0, r1)
1554 bval = bump(norm2(ri - rref), r0, r1)
1555 IF (my_left .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1556 bval = bump(norm2(rj - rref), r0, r1)
1557 IF (my_right .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1560 CALL dbt_put_block(t_2c_inout, ind, shape(blk), blk)
1564 CALL dbt_iterator_stop(iter)
1566 CALL dbt_filter(t_2c_inout, ri_data%filter_eps)
1568 END SUBROUTINE apply_bump
1582 SUBROUTINE get_2c_bump_forces(force, t_2c_in, atom_i, atom_of_kind, kind_of, pref, ri_data, &
1583 qs_env, work_virial)
1585 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_in
1586 INTEGER,
INTENT(IN) :: atom_i
1587 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
1588 REAL(
dp),
INTENT(IN) :: pref
1591 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT) :: work_virial
1593 INTEGER :: i, i_img, i_ri, i_xyz, iat_of_kind, iatom, ikind, ind(2), j_img, j_ri, j_xyz, &
1594 jat_of_kind, jatom, jkind, natom, nblks(2), nimg, nkind
1595 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1596 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1598 REAL(
dp) :: new_force, r0, r1, ri(3), rj(3), &
1599 rref(3), scoord(3), x
1600 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1602 TYPE(dbt_iterator_type) :: iter
1605 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1607 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1609 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1610 kpoints=kpoints, particle_set=particle_set)
1611 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1613 CALL dbt_get_info(t_2c_in, nblks_total=nblks)
1614 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1615 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1620 r1 = ri_data%kp_RI_range
1621 r0 = ri_data%kp_bump_rad
1622 rref =
pbc(particle_set(atom_i)%r, cell)
1624 iat_of_kind = atom_of_kind(atom_i)
1625 ikind = kind_of(atom_i)
1631 CALL dbt_iterator_start(iter, t_2c_in)
1632 DO WHILE (dbt_iterator_blocks_left(iter))
1633 CALL dbt_iterator_next_block(iter, ind)
1634 IF (ind(1) /= ind(2)) cycle
1636 CALL dbt_get_block(t_2c_in, ind, blk, found)
1637 IF (.NOT. found) cycle
1640 j_ri = (ind(2) - 1)/natom + 1
1641 j_img = ri_data%RI_cell_to_img(j_ri)
1642 jatom = ind(2) - (j_ri - 1)*natom
1643 jat_of_kind = atom_of_kind(jatom)
1644 jkind = kind_of(jatom)
1647 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1648 x = norm2(rj - rref)
1649 IF (x < r0 .OR. x > r1) cycle
1652 DO i = 1,
SIZE(blk, 1)
1653 new_force = new_force + blk(i, i)
1655 new_force = pref*new_force*dbump(x, r0, r1)
1661 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) + &
1662 new_force*(rj(i_xyz) - rref(i_xyz))/x
1668 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1669 + new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1674 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) - &
1675 new_force*(rj(i_xyz) - rref(i_xyz))/x
1681 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1682 - new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1688 CALL dbt_iterator_stop(iter)
1691 END SUBROUTINE get_2c_bump_forces
1700 FUNCTION bump(x, r0, r1)
RESULT(b)
1701 REAL(
dp),
INTENT(IN) :: x, r0, r1
1709 r = (x - r0)/(r1 - r0)
1710 b = -6.0_dp*r**5 + 15.0_dp*r**4 - 10.0_dp*r**3 + 1.0_dp
1711 IF (x >= r1) b = 0.0_dp
1712 IF (x <= r0) b = 1.0_dp
1723 FUNCTION dbump(x, r0, r1)
RESULT(b)
1724 REAL(
dp),
INTENT(IN) :: x, r0, r1
1729 r = (x - r0)/(r1 - r0)
1730 b = (-30.0_dp*r**4 + 60.0_dp*r**3 - 30.0_dp*r**2)/(r1 - r0)
1731 IF (x >= r1) b = 0.0_dp
1732 IF (x <= r0) b = 0.0_dp
1743 FUNCTION get_apc_index_from_ib(i_index, b_index, qs_env)
RESULT(apc_index)
1744 INTEGER,
INTENT(IN) :: i_index, b_index
1746 INTEGER :: apc_index
1748 INTEGER,
DIMENSION(3) :: cell_apc
1749 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1750 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1754 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1757 cell_apc(:) = index_to_cell(:, i_index) + index_to_cell(:, b_index)
1759 IF (any([cell_apc(1), cell_apc(2), cell_apc(3)] < lbound(cell_to_index)) .OR. &
1760 any([cell_apc(1), cell_apc(2), cell_apc(3)] > ubound(cell_to_index)))
THEN
1764 apc_index = cell_to_index(cell_apc(1), cell_apc(2), cell_apc(3))
1767 END FUNCTION get_apc_index_from_ib
1776 FUNCTION get_apc_index(a_index, c_index, qs_env)
RESULT(i_index)
1777 INTEGER,
INTENT(IN) :: a_index, c_index
1781 INTEGER,
DIMENSION(3) :: cell_i
1782 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1783 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1787 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1789 cell_i(:) = index_to_cell(:, a_index) + index_to_cell(:, c_index)
1791 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1792 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1796 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1799 END FUNCTION get_apc_index
1808 FUNCTION get_i_index(apc_index, b_index, qs_env)
RESULT(i_index)
1809 INTEGER,
INTENT(IN) :: apc_index, b_index
1813 INTEGER,
DIMENSION(3) :: cell_i
1814 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1815 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1819 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1821 cell_i(:) = index_to_cell(:, apc_index) - index_to_cell(:, b_index)
1823 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1824 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1828 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1831 END FUNCTION get_i_index
1842 SUBROUTINE get_ac_pairs(ac_pairs, apc_index, ri_data, qs_env)
1843 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: ac_pairs
1844 INTEGER,
INTENT(IN) :: apc_index
1848 INTEGER :: a_index, actual_img, c_index, nimg
1850 nimg =
SIZE(ac_pairs, 1)
1855 DO a_index = 1, nimg
1856 actual_img = ri_data%idx_to_img(a_index)
1858 c_index = get_i_index(apc_index, actual_img, qs_env)
1859 ac_pairs(a_index, 1) = a_index
1860 ac_pairs(a_index, 2) = c_index
1864 END SUBROUTINE get_ac_pairs
1876 SUBROUTINE get_iapc_pairs(iapc_pairs, b_index, ri_data, qs_env, actual_i_img)
1877 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: iapc_pairs
1878 INTEGER,
INTENT(IN) :: b_index
1881 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: actual_i_img
1883 INTEGER :: actual_img, apc_index, i_index, nimg
1885 nimg =
SIZE(iapc_pairs, 1)
1886 IF (
PRESENT(actual_i_img)) actual_i_img(:) = 0
1888 iapc_pairs(:, :) = 0
1891 DO i_index = 1, nimg
1892 actual_img = ri_data%idx_to_img(i_index)
1893 apc_index = get_apc_index_from_ib(actual_img, b_index, qs_env)
1894 IF (apc_index == 0) cycle
1895 iapc_pairs(i_index, 1) = i_index
1896 iapc_pairs(i_index, 2) = apc_index
1897 IF (
PRESENT(actual_i_img)) actual_i_img(i_index) = actual_img
1900 END SUBROUTINE get_iapc_pairs
1909 FUNCTION get_opp_index(a_index, qs_env)
RESULT(opp_index)
1910 INTEGER,
INTENT(IN) :: a_index
1912 INTEGER :: opp_index
1914 INTEGER,
DIMENSION(3) :: opp_cell
1915 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1916 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1919 NULLIFY (kpoints, cell_to_index, index_to_cell)
1922 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1924 opp_cell(:) = -index_to_cell(:, a_index)
1926 IF (any([opp_cell(1), opp_cell(2), opp_cell(3)] < lbound(cell_to_index)) .OR. &
1927 any([opp_cell(1), opp_cell(2), opp_cell(3)] > ubound(cell_to_index)))
THEN
1931 opp_index = cell_to_index(opp_cell(1), opp_cell(2), opp_cell(3))
1934 END FUNCTION get_opp_index
1945 SUBROUTINE get_pmat_images(rho_ao_t, rho_ao, scale_prev_p, ri_data, qs_env)
1946 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t
1947 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao
1948 REAL(
dp),
INTENT(IN) :: scale_prev_p
1952 INTEGER :: cell_j(3), i_img, i_spin, iatom, icol, &
1953 irow, j_img, jatom, mi_img, mj_img, &
1955 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1958 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock, pblock_desymm
1959 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, rho_desymm
1960 TYPE(dbt_type) :: tmp
1964 DIMENSION(:),
POINTER :: nl_iterator
1966 POINTER :: sab_nl, sab_nl_nosym
1969 NULLIFY (rho_desymm, kpoints, sab_nl_nosym, scf_env, matrix_ks, dft_control, &
1970 sab_nl, nl_iterator, cell_to_index, pblock, pblock_desymm)
1972 CALL get_qs_env(qs_env, kpoints=kpoints, scf_env=scf_env, matrix_ks_kp=matrix_ks, dft_control=dft_control)
1973 CALL get_kpoint_info(kpoints, sab_nl_nosym=sab_nl_nosym, cell_to_index=cell_to_index, sab_nl=sab_nl)
1975 IF (dft_control%do_admm)
THEN
1976 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks)
1979 nspins =
SIZE(matrix_ks, 1)
1982 ALLOCATE (rho_desymm(nspins, nimg))
1984 DO i_spin = 1, nspins
1985 ALLOCATE (rho_desymm(i_spin, i_img)%matrix)
1986 CALL dbcsr_create(rho_desymm(i_spin, i_img)%matrix, template=matrix_ks(i_spin, i_img)%matrix, &
1987 matrix_type=dbcsr_type_no_symmetry)
1991 CALL dbt_create(rho_desymm(1, 1)%matrix, tmp)
1998 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1999 IF (j_img > nimg .OR. j_img < 1) cycle
2002 IF (iatom == jatom)
fac = 0.5_dp
2003 mj_img = get_opp_index(j_img, qs_env)
2005 IF (mj_img == 0)
fac = 1.0_dp
2009 IF (iatom > jatom)
THEN
2015 DO i_spin = 1, nspins
2017 IF (.NOT. found) cycle
2020 CALL dbcsr_get_block_p(rho_desymm(i_spin, j_img)%matrix, iatom, jatom, pblock_desymm, found)
2021 IF (.NOT. found) cycle
2023 IF (iatom > jatom)
THEN
2024 pblock_desymm(:, :) =
fac*transpose(pblock(:, :))
2026 pblock_desymm(:, :) =
fac*pblock(:, :)
2033 DO i_spin = 1, nspins
2034 CALL dbt_scale(rho_ao_t(i_spin, i_img), scale_prev_p)
2036 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, i_img)%matrix, tmp)
2037 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), summation=.true., move_data=.true.)
2040 mi_img = get_opp_index(i_img, qs_env)
2041 IF (mi_img > 0 .AND. mi_img <= nimg)
THEN
2042 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, mi_img)%matrix, tmp)
2043 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), order=[2, 1], summation=.true., move_data=.true.)
2045 CALL dbt_filter(rho_ao_t(i_spin, i_img), ri_data%filter_eps)
2050 DO i_spin = 1, nspins
2052 DEALLOCATE (rho_desymm(i_spin, i_img)%matrix)
2056 CALL dbt_destroy(tmp)
2057 DEALLOCATE (rho_desymm)
2059 END SUBROUTINE get_pmat_images
2078 SUBROUTINE get_ext_2c_int(t_2c_pot, mat_orig, atom_i, atom_j, img_b, ri_data, qs_env, do_inverse, &
2079 para_env_ext, blacs_env_ext, dbcsr_template, off_diagonal, skip_inverse)
2080 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_pot
2081 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_orig
2082 INTEGER,
INTENT(IN) :: atom_i, atom_j, img_b
2085 LOGICAL,
INTENT(IN),
OPTIONAL :: do_inverse
2088 TYPE(
dbcsr_type),
OPTIONAL,
POINTER :: dbcsr_template
2089 LOGICAL,
INTENT(IN),
OPTIONAL :: off_diagonal, skip_inverse
2091 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_ext_2c_int'
2093 INTEGER :: group, handle, handle2, i_img, i_ri, iatom, iblk, ikind, img_tot, j_img, j_ri, &
2094 jatom, jblk, jkind, n_dependent, natom, nblks_ri, nimg, nkind
2095 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2
2096 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: present_atoms_i, present_atoms_j
2097 INTEGER,
DIMENSION(3) :: cell_b, cell_i, cell_j, cell_tot
2098 INTEGER,
DIMENSION(:),
POINTER :: col_dist, col_dist_ext, ri_blk_size_ext, &
2099 row_dist, row_dist_ext
2100 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell, pgrid
2101 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
2102 LOGICAL :: do_inverse_prv, found, my_offd, &
2103 skip_inverse_prv, use_template
2104 REAL(
dp) :: bfac, dij, r0, r1, threshold
2105 REAL(
dp),
DIMENSION(3) :: ri, rij, rj, rref, scoord
2106 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock
2111 TYPE(
dbcsr_type) :: work, work_tight, work_tight_inv
2112 TYPE(dbt_type) :: t_2c_tmp
2115 DIMENSION(:),
TARGET :: basis_set_ri
2119 DIMENSION(:),
POINTER :: nl_iterator
2123 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2125 NULLIFY (qs_kind_set, nl_2c, nl_iterator, cell, kpoints, cell_to_index, index_to_cell, dist_2d, &
2126 para_env, pblock, blacs_env, particle_set, col_dist, row_dist, pgrid, &
2127 col_dist_ext, row_dist_ext)
2129 CALL timeset(routinen, handle)
2134 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
2135 kpoints=kpoints, para_env=para_env, blacs_env=blacs_env, particle_set=particle_set)
2136 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
2138 do_inverse_prv = .false.
2139 IF (
PRESENT(do_inverse)) do_inverse_prv = do_inverse
2140 IF (do_inverse_prv)
THEN
2141 cpassert(atom_i == atom_j)
2144 skip_inverse_prv = .false.
2145 IF (
PRESENT(skip_inverse)) skip_inverse_prv = skip_inverse
2148 IF (
PRESENT(off_diagonal)) my_offd = off_diagonal
2150 IF (
PRESENT(para_env_ext)) para_env => para_env_ext
2151 IF (
PRESENT(blacs_env_ext)) blacs_env => blacs_env_ext
2153 nimg =
SIZE(mat_orig)
2155 CALL timeset(routinen//
"_nl_iter", handle2)
2158 ALLOCATE (dist1(natom), dist2(natom))
2160 dist1(iatom) = mod(iatom, blacs_env%num_pe(1))
2161 dist2(iatom) = mod(iatom, blacs_env%num_pe(2))
2165 ALLOCATE (basis_set_ri(nkind))
2169 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
2171 ALLOCATE (present_atoms_i(natom, nimg), present_atoms_j(natom, nimg))
2177 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij, cell=cell_j, &
2178 ikind=ikind, jkind=jkind)
2182 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2183 IF (j_img > nimg .OR. j_img < 1) cycle
2185 IF (iatom == atom_i .AND. dij <= ri_data%kp_RI_range) present_atoms_i(jatom, j_img) = 1
2186 IF (iatom == atom_j .AND. dij <= ri_data%kp_RI_range) present_atoms_j(jatom, j_img) = 1
2191 CALL timestop(handle2)
2193 CALL para_env%sum(present_atoms_i)
2194 CALL para_env%sum(present_atoms_j)
2198 use_template = .false.
2199 IF (
PRESENT(dbcsr_template))
THEN
2200 IF (
ASSOCIATED(dbcsr_template)) use_template = .true.
2203 IF (use_template)
THEN
2208 ALLOCATE (row_dist_ext(ri_data%ncell_RI*natom), col_dist_ext(ri_data%ncell_RI*natom))
2209 ALLOCATE (ri_blk_size_ext(ri_data%ncell_RI*natom))
2210 DO i_ri = 1, ri_data%ncell_RI
2211 row_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = row_dist(:)
2212 col_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = col_dist(:)
2213 ri_blk_size_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2217 row_dist=row_dist_ext, col_dist=col_dist_ext)
2218 CALL dbcsr_create(work, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2219 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2221 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2223 IF (
PRESENT(dbcsr_template))
THEN
2224 ALLOCATE (dbcsr_template)
2229 cell_b(:) = index_to_cell(:, img_b)
2231 i_ri = ri_data%img_to_RI_cell(i_img)
2232 IF (i_ri == 0) cycle
2233 cell_i(:) = index_to_cell(:, i_img)
2235 j_ri = ri_data%img_to_RI_cell(j_img)
2236 IF (j_ri == 0) cycle
2237 cell_j(:) = index_to_cell(:, j_img)
2238 cell_tot = cell_j - cell_i + cell_b
2240 IF (any([cell_tot(1), cell_tot(2), cell_tot(3)] < lbound(cell_to_index)) .OR. &
2241 any([cell_tot(1), cell_tot(2), cell_tot(3)] > ubound(cell_to_index))) cycle
2242 img_tot = cell_to_index(cell_tot(1), cell_tot(2), cell_tot(3))
2243 IF (img_tot > nimg .OR. img_tot < 1) cycle
2248 IF (present_atoms_i(iatom, i_img) == 0) cycle
2249 IF (present_atoms_j(jatom, j_img) == 0) cycle
2250 IF (my_offd .AND. (i_ri - 1)*natom + iatom == (j_ri - 1)*natom + jatom) cycle
2253 IF (.NOT. found) cycle
2255 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2264 IF (do_inverse_prv)
THEN
2266 r1 = ri_data%kp_RI_range
2267 r0 = ri_data%kp_bump_rad
2270 nblks_ri = sum(present_atoms_i)
2271 ALLOCATE (col_dist_ext(nblks_ri), row_dist_ext(nblks_ri), ri_blk_size_ext(nblks_ri))
2274 i_ri = ri_data%img_to_RI_cell(i_img)
2275 IF (i_ri == 0) cycle
2277 IF (present_atoms_i(iatom, i_img) == 0) cycle
2279 col_dist_ext(iblk) = col_dist(iatom)
2280 row_dist_ext(iblk) = row_dist(iatom)
2281 ri_blk_size_ext(iblk) = ri_data%bsizes_RI(iatom)
2286 row_dist=row_dist_ext, col_dist=col_dist_ext)
2287 CALL dbcsr_create(work_tight, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2288 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2289 CALL dbcsr_create(work_tight_inv, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2290 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2292 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2296 rref =
pbc(particle_set(atom_i)%r, cell)
2300 i_ri = ri_data%img_to_RI_cell(i_img)
2301 IF (i_ri == 0) cycle
2303 IF (present_atoms_i(iatom, i_img) == 0) cycle
2307 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
2311 j_ri = ri_data%img_to_RI_cell(j_img)
2312 IF (j_ri == 0) cycle
2314 IF (present_atoms_j(jatom, j_img) == 0) cycle
2318 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
2320 CALL dbcsr_get_block_p(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock, found)
2321 IF (.NOT. found) cycle
2324 IF (iblk /= jblk) bfac = bump(norm2(ri - rref), r0, r1)*bump(norm2(rj - rref), r0, r1)
2333 IF (.NOT. skip_inverse_prv)
THEN
2334 SELECT CASE (ri_data%t2c_method)
2336 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
2337 CALL invert_hotelling(work_tight_inv, work_tight, threshold=threshold, silent=.false.)
2342 uplo_to_full=.true.)
2345 CALL cp_dbcsr_power(work_tight_inv, -1.0_dp, ri_data%eps_eigval, n_dependent, &
2346 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
2357 i_ri = ri_data%img_to_RI_cell(i_img)
2358 IF (i_ri == 0) cycle
2360 IF (present_atoms_i(iatom, i_img) == 0) cycle
2365 j_ri = ri_data%img_to_RI_cell(j_img)
2366 IF (j_ri == 0) cycle
2368 IF (present_atoms_j(jatom, j_img) == 0) cycle
2372 IF (.NOT. found) cycle
2374 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2385 CALL dbt_create(work, t_2c_tmp)
2386 CALL dbt_copy_matrix_to_tensor(work, t_2c_tmp)
2387 CALL dbt_copy(t_2c_tmp, t_2c_pot, move_data=.true.)
2388 CALL dbt_filter(t_2c_pot, ri_data%filter_eps)
2390 CALL dbt_destroy(t_2c_tmp)
2393 CALL timestop(handle)
2395 END SUBROUTINE get_ext_2c_int
2405 SUBROUTINE contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
2406 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, rho_ao_t
2410 CHARACTER(len=*),
PARAMETER :: routinen =
'contract_pmat_3c'
2412 INTEGER :: apc_img, b_img, batch_size, handle, &
2413 i_batch, i_img, i_spin,
idx, j_batch, &
2414 n_batch_img, n_batch_nze, nimg, &
2416 INTEGER(int_8) :: nflop, nze
2417 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: apc_filter, batch_ranges_img, &
2418 batch_ranges_nze, int_indices
2419 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: ac_pairs, iapc_pairs
2420 REAL(
dp) :: occ, t1, t2
2421 TYPE(dbt_type) :: t_3c_tmp
2422 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ints_stack, res_stack, rho_stack
2425 CALL timeset(routinen, handle)
2427 CALL get_qs_env(qs_env, dft_control=dft_control)
2430 nimg_nze = ri_data%nimg_nze
2431 nspins = dft_control%nspins
2433 CALL dbt_create(t_3c_apc(1, 1), t_3c_tmp)
2435 batch_size = ri_data%kp_stack_size
2437 ALLOCATE (apc_filter(nimg), iapc_pairs(nimg, 2))
2440 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env)
2441 DO i_img = 1, nimg_nze
2442 idx = iapc_pairs(i_img, 2)
2443 IF (
idx < 1 .OR.
idx > nimg) cycle
2449 n_batch_img = nimg/batch_size
2450 IF (
modulo(nimg, batch_size) /= 0) n_batch_img = n_batch_img + 1
2451 ALLOCATE (batch_ranges_img(n_batch_img + 1))
2452 DO i_batch = 1, n_batch_img
2453 batch_ranges_img(i_batch) = (i_batch - 1)*batch_size + 1
2455 batch_ranges_img(n_batch_img + 1) = nimg + 1
2458 n_batch_nze = nimg_nze/batch_size
2459 IF (
modulo(nimg_nze, batch_size) /= 0) n_batch_nze = n_batch_nze + 1
2460 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
2461 DO i_batch = 1, n_batch_nze
2462 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
2464 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
2467 ALLOCATE (rho_stack(2), ints_stack(2), res_stack(2))
2468 CALL get_stack_tensors(res_stack, rho_stack, ints_stack, rho_ao_t(1, 1), &
2469 ri_data%t_3c_int_ctr_1(1, 1), batch_size, ri_data, qs_env)
2471 ALLOCATE (ac_pairs(nimg, 2), int_indices(nimg_nze))
2472 DO i_img = 1, nimg_nze
2473 int_indices(i_img) = i_img
2477 DO j_batch = 1, n_batch_nze
2479 CALL fill_3c_stack(ints_stack(1), ri_data%t_3c_int_ctr_1(1, :), int_indices, 3, ri_data, &
2480 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)])
2481 CALL dbt_copy(ints_stack(1), ints_stack(2), move_data=.true.)
2483 DO i_spin = 1, nspins
2484 DO i_batch = 1, n_batch_img
2486 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2487 IF (apc_filter(apc_img) == 0) cycle
2488 CALL get_ac_pairs(ac_pairs, apc_img, ri_data, qs_env)
2489 CALL fill_2c_stack(rho_stack(1), rho_ao_t(i_spin, :), ac_pairs(:, 2), 1, ri_data, &
2490 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)], &
2491 shift=apc_img - batch_ranges_img(i_batch) + 1)
2496 CALL dbt_copy(rho_stack(1), rho_stack(2), move_data=.true.)
2499 CALL dbt_batched_contract_init(rho_stack(2))
2500 CALL dbt_contract(1.0_dp, ints_stack(2), rho_stack(2), &
2501 0.0_dp, res_stack(2), map_1=[1, 2], map_2=[3], &
2502 contract_1=[3], notcontract_1=[1, 2], &
2503 contract_2=[1], notcontract_2=[2], &
2504 filter_eps=ri_data%filter_eps, flop=nflop)
2505 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2506 CALL dbt_batched_contract_finalize(rho_stack(2))
2507 CALL dbt_copy(res_stack(2), res_stack(1), move_data=.true.)
2509 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2511 IF (apc_filter(apc_img) == 0) cycle
2512 CALL unstack_t_3c_apc(t_3c_tmp, res_stack(1), apc_img - batch_ranges_img(i_batch) + 1)
2513 CALL dbt_copy(t_3c_tmp, t_3c_apc(i_spin, apc_img), summation=.true., move_data=.true.)
2519 DEALLOCATE (batch_ranges_img)
2520 DEALLOCATE (batch_ranges_nze)
2522 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
2524 CALL dbt_destroy(rho_stack(1))
2525 CALL dbt_destroy(rho_stack(2))
2526 CALL dbt_destroy(ints_stack(1))
2527 CALL dbt_destroy(ints_stack(2))
2528 CALL dbt_destroy(res_stack(1))
2529 CALL dbt_destroy(res_stack(2))
2530 CALL dbt_destroy(t_3c_tmp)
2532 CALL timestop(handle)
2534 END SUBROUTINE contract_pmat_3c
2542 SUBROUTINE precontract_3c_ints(t_3c_int, ri_data, qs_env)
2543 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_int
2547 CHARACTER(len=*),
PARAMETER :: routinen =
'precontract_3c_ints'
2549 INTEGER :: batch_size, handle, i_batch, i_img, &
2550 i_ri, iatom, is, n_batch, natom, &
2551 nblks, nblks_3c(3), nimg
2552 INTEGER(int_8) :: nflop
2553 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges, bsizes_ri_ext, bsizes_ri_ext_split, &
2554 bsizes_stack, dist1, dist2, dist3, dist_stack3, idx_to_at_ao, int_indices
2555 TYPE(dbt_distribution_type) :: t_dist
2556 TYPE(dbt_type) :: t_2c_ri_tmp(2), t_3c_tmp(3)
2558 CALL timeset(routinen, handle)
2563 ALLOCATE (int_indices(nimg))
2565 int_indices(i_img) = i_img
2568 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
2569 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
2571 nblks =
SIZE(ri_data%bsizes_RI_split)
2572 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
2573 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
2574 DO i_ri = 1, ri_data%ncell_RI
2575 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2576 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
2579 bsizes_ri_ext, bsizes_ri_ext, &
2581 DEALLOCATE (dist1, dist2)
2583 bsizes_ri_ext_split, bsizes_ri_ext_split, &
2585 DEALLOCATE (dist1, dist2)
2588 batch_size = ri_data%kp_stack_size
2589 n_batch = nimg/batch_size
2590 IF (
modulo(nimg, batch_size) /= 0) n_batch = n_batch + 1
2591 ALLOCATE (batch_ranges(n_batch + 1))
2592 DO i_batch = 1, n_batch
2593 batch_ranges(i_batch) = (i_batch - 1)*batch_size + 1
2595 batch_ranges(n_batch + 1) = nimg + 1
2597 nblks =
SIZE(ri_data%bsizes_AO_split)
2598 ALLOCATE (bsizes_stack(batch_size*nblks))
2599 DO is = 1, batch_size
2600 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
2603 CALL dbt_get_info(t_3c_int(1, 1), nblks_total=nblks_3c)
2604 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)), dist_stack3(batch_size*nblks_3c(3)))
2605 CALL dbt_get_info(t_3c_int(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, proc_dist_3=dist3)
2606 DO is = 1, batch_size
2607 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
2610 CALL dbt_distribution_new(t_dist, ri_data%pgrid, dist1, dist2, dist_stack3)
2611 CALL dbt_create(t_3c_tmp(1),
"ints_stack", t_dist, [1], [2, 3], bsizes_ri_ext_split, &
2612 ri_data%bsizes_AO_split, bsizes_stack)
2613 CALL dbt_distribution_destroy(t_dist)
2614 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
2616 CALL dbt_create(t_3c_tmp(1), t_3c_tmp(2))
2617 CALL dbt_create(t_3c_int(1, 1), t_3c_tmp(3))
2620 CALL dbt_copy(ri_data%t_2c_inv(1, iatom), t_2c_ri_tmp(1))
2621 CALL apply_bump(t_2c_ri_tmp(1), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
2622 CALL dbt_copy(t_2c_ri_tmp(1), t_2c_ri_tmp(2), move_data=.true.)
2624 CALL dbt_batched_contract_init(t_2c_ri_tmp(2))
2625 DO i_batch = 1, n_batch
2627 CALL fill_3c_stack(t_3c_tmp(1), t_3c_int(1, :), int_indices, 3, ri_data, &
2628 img_bounds=[batch_ranges(i_batch), batch_ranges(i_batch + 1)], &
2629 filter_at=iatom, filter_dim=2, idx_to_at=idx_to_at_ao)
2631 CALL dbt_contract(1.0_dp, t_2c_ri_tmp(2), t_3c_tmp(1), &
2632 0.0_dp, t_3c_tmp(2), map_1=[1], map_2=[2, 3], &
2633 contract_1=[2], notcontract_1=[1], &
2634 contract_2=[1], notcontract_2=[2, 3], &
2635 filter_eps=ri_data%filter_eps, flop=nflop)
2636 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2638 DO i_img = batch_ranges(i_batch), batch_ranges(i_batch + 1) - 1
2639 CALL unstack_t_3c_apc(t_3c_tmp(3), t_3c_tmp(2), i_img - batch_ranges(i_batch) + 1)
2640 CALL dbt_copy(t_3c_tmp(3), ri_data%t_3c_int_ctr_1(1, i_img), summation=.true., &
2641 order=[2, 1, 3], move_data=.true.)
2643 CALL dbt_clear(t_3c_tmp(1))
2645 CALL dbt_batched_contract_finalize(t_2c_ri_tmp(2))
2648 CALL dbt_destroy(t_2c_ri_tmp(1))
2649 CALL dbt_destroy(t_2c_ri_tmp(2))
2650 CALL dbt_destroy(t_3c_tmp(1))
2651 CALL dbt_destroy(t_3c_tmp(2))
2652 CALL dbt_destroy(t_3c_tmp(3))
2655 CALL dbt_destroy(t_3c_int(1, i_img))
2658 CALL timestop(handle)
2660 END SUBROUTINE precontract_3c_ints
2671 SUBROUTINE copy_2c_to_subgroup(t2c_sub, t2c_main, group_size, ngroups, para_env)
2672 TYPE(dbt_type),
INTENT(INOUT) :: t2c_sub, t2c_main
2673 INTEGER,
INTENT(IN) :: group_size, ngroups
2676 INTEGER :: batch_size, i, i_batch, i_msg, iblk, &
2677 igroup, iproc, ir, is, jblk, n_batch, &
2679 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2
2680 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: block_dest, block_source
2681 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: current_dest
2682 INTEGER,
DIMENSION(2) :: ind, nblks
2684 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
2685 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2686 TYPE(dbt_iterator_type) :: iter
2687 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2693 CALL dbt_get_info(t2c_main, nblks_total=nblks)
2696 ALLOCATE (block_source(nblks(1), nblks(2)))
2700 CALL dbt_iterator_start(iter, t2c_main)
2701 DO WHILE (dbt_iterator_blocks_left(iter))
2702 CALL dbt_iterator_next_block(iter, ind)
2703 CALL dbt_get_block(t2c_main, ind, blk, found)
2704 IF (.NOT. found) cycle
2706 block_source(ind(1), ind(2)) = para_env%mepos
2711 CALL dbt_iterator_stop(iter)
2714 CALL para_env%sum(nocc)
2715 CALL para_env%sum(block_source)
2716 block_source = block_source + para_env%num_pe - 1
2717 IF (nocc == 0)
RETURN
2720 igroup = para_env%mepos/group_size
2721 ALLOCATE (block_dest(nblks(1), nblks(2)))
2723 DO jblk = 1, nblks(2)
2724 DO iblk = 1, nblks(1)
2725 IF (block_source(iblk, jblk) == -1) cycle
2727 CALL dbt_get_stored_coordinates(t2c_sub, [iblk, jblk], iproc)
2728 block_dest(iblk, jblk) = igroup*group_size + iproc
2732 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)))
2733 CALL dbt_get_info(t2c_main, blk_size_1=bsizes1, blk_size_2=bsizes2)
2735 ALLOCATE (current_dest(nblks(1), nblks(2), 0:ngroups - 1))
2736 DO igroup = 0, ngroups - 1
2738 current_dest(:, :, igroup) = block_dest(:, :)
2739 CALL para_env%bcast(current_dest(:, :, igroup), source=igroup*group_size)
2743 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2744 n_batch = (nocc*ngroups)/batch_size
2745 IF (
modulo(nocc*ngroups, batch_size) /= 0) n_batch = n_batch + 1
2747 DO i_batch = 1, n_batch
2749 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2750 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2754 DO jblk = 1, nblks(2)
2755 DO iblk = 1, nblks(1)
2756 DO igroup = 0, ngroups - 1
2757 IF (block_source(iblk, jblk) == -1) cycle
2760 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2763 tag = i_msg - (i_batch - 1)*batch_size
2766 IF (para_env%mepos == block_source(iblk, jblk))
THEN
2767 CALL dbt_get_block(t2c_main, [iblk, jblk], blk, found)
2771 IF (block_source(iblk, jblk) == current_dest(iblk, jblk, igroup))
THEN
2772 IF (found)
CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2774 IF (para_env%mepos == block_source(iblk, jblk) .AND. found)
THEN
2775 ALLOCATE (send_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2776 send_buff(tag)%array(:, :) = blk(:, :)
2778 CALL para_env%isend(msgin=send_buff(tag)%array, dest=current_dest(iblk, jblk, igroup), &
2779 request=send_req(is), tag=tag)
2782 IF (para_env%mepos == current_dest(iblk, jblk, igroup))
THEN
2783 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2785 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk), &
2786 request=recv_req(ir), tag=tag)
2790 IF (found)
DEALLOCATE (blk)
2798 DO i = 1, batch_size
2799 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
2804 DO jblk = 1, nblks(2)
2805 DO iblk = 1, nblks(1)
2806 DO igroup = 0, ngroups - 1
2807 IF (block_source(iblk, jblk) == -1) cycle
2810 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2813 tag = i_msg - (i_batch - 1)*batch_size
2815 IF (para_env%mepos == current_dest(iblk, jblk, igroup) .AND. &
2816 block_source(iblk, jblk) /= current_dest(iblk, jblk, igroup))
THEN
2818 ALLOCATE (blk(bsizes1(iblk), bsizes2(jblk)))
2819 blk(:, :) = recv_buff(tag)%array(:, :)
2820 CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2828 DO i = 1, batch_size
2829 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
2831 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
2833 CALL dbt_finalize(t2c_sub)
2835 END SUBROUTINE copy_2c_to_subgroup
2846 SUBROUTINE get_3c_subgroup_dest(subgroup_dest, t3c_sub, t3c_main, group_size, ngroups, para_env)
2847 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :), &
2848 INTENT(INOUT) :: subgroup_dest
2849 TYPE(dbt_type),
INTENT(INOUT) :: t3c_sub, t3c_main
2850 INTEGER,
INTENT(IN) :: group_size, ngroups
2853 INTEGER :: iblk, igroup, iproc, jblk, kblk
2854 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: block_dest
2855 INTEGER,
DIMENSION(3) :: nblks
2857 CALL dbt_get_info(t3c_main, nblks_total=nblks)
2860 igroup = para_env%mepos/group_size
2861 ALLOCATE (block_dest(nblks(1), nblks(2), nblks(3)))
2862 DO kblk = 1, nblks(3)
2863 DO jblk = 1, nblks(2)
2864 DO iblk = 1, nblks(1)
2865 CALL dbt_get_stored_coordinates(t3c_sub, [iblk, jblk, kblk], iproc)
2866 block_dest(iblk, jblk, kblk) = igroup*group_size + iproc
2871 ALLOCATE (subgroup_dest(nblks(1), nblks(2), nblks(3), ngroups))
2872 DO igroup = 0, ngroups - 1
2874 subgroup_dest(:, :, :, igroup + 1) = block_dest(:, :, :)
2875 CALL para_env%bcast(subgroup_dest(:, :, :, igroup + 1), source=igroup*group_size)
2878 END SUBROUTINE get_3c_subgroup_dest
2892 SUBROUTINE copy_3c_to_subgroup(t3c_sub, t3c_main, ngroups, para_env, subgroup_dest, &
2893 iatom_to_subgroup, dim_at, idx_to_at)
2894 TYPE(dbt_type),
INTENT(INOUT) :: t3c_sub, t3c_main
2895 INTEGER,
INTENT(IN) :: ngroups
2897 INTEGER,
DIMENSION(:, :, :, :),
INTENT(IN) :: subgroup_dest
2899 INTENT(INOUT),
OPTIONAL :: iatom_to_subgroup
2900 INTEGER,
INTENT(IN),
OPTIONAL :: dim_at
2901 INTEGER,
DIMENSION(:),
OPTIONAL :: idx_to_at
2903 INTEGER :: batch_size, i, i_batch, i_msg, iatom, &
2904 iblk, igroup, ir, is, isbuff, jblk, &
2905 kblk, n_batch, nocc, tag
2906 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2, bsizes3
2907 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: block_source
2908 INTEGER,
DIMENSION(3) :: ind, nblks
2909 LOGICAL :: filter_at, found
2910 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
2911 TYPE(
cp_3d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2912 TYPE(dbt_iterator_type) :: iter
2913 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2919 CALL dbt_get_info(t3c_main, nblks_total=nblks)
2923 IF (
PRESENT(iatom_to_subgroup) .AND.
PRESENT(dim_at) .AND.
PRESENT(idx_to_at))
THEN
2925 cpassert(nblks(dim_at) ==
SIZE(idx_to_at))
2929 ALLOCATE (block_source(nblks(1), nblks(2), nblks(3)))
2933 CALL dbt_iterator_start(iter, t3c_main)
2934 DO WHILE (dbt_iterator_blocks_left(iter))
2935 CALL dbt_iterator_next_block(iter, ind)
2936 CALL dbt_get_block(t3c_main, ind, blk, found)
2937 IF (.NOT. found) cycle
2939 block_source(ind(1), ind(2), ind(3)) = para_env%mepos
2944 CALL dbt_iterator_stop(iter)
2947 CALL para_env%sum(nocc)
2948 CALL para_env%sum(block_source)
2949 block_source = block_source + para_env%num_pe - 1
2950 IF (nocc == 0)
RETURN
2952 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)), bsizes3(nblks(3)))
2953 CALL dbt_get_info(t3c_main, blk_size_1=bsizes1, blk_size_2=bsizes2, blk_size_3=bsizes3)
2956 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2957 n_batch = (nocc*ngroups)/batch_size
2958 IF (
modulo(nocc*ngroups, batch_size) /= 0) n_batch = n_batch + 1
2960 DO i_batch = 1, n_batch
2962 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2963 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2968 DO kblk = 1, nblks(3)
2969 DO jblk = 1, nblks(2)
2970 DO iblk = 1, nblks(1)
2971 IF (block_source(iblk, jblk, kblk) == -1) cycle
2974 IF (para_env%mepos == block_source(iblk, jblk, kblk))
THEN
2975 CALL dbt_get_block(t3c_main, [iblk, jblk, kblk], blk, found)
2978 ALLOCATE (send_buff(isbuff)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2982 DO igroup = 0, ngroups - 1
2985 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2988 tag = i_msg - (i_batch - 1)*batch_size
2991 ind(:) = [iblk, jblk, kblk]
2992 iatom = idx_to_at(ind(dim_at))
2993 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
2997 IF (block_source(iblk, jblk, kblk) == subgroup_dest(iblk, jblk, kblk, igroup + 1))
THEN
2998 IF (found)
CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(blk), blk)
3000 IF (para_env%mepos == block_source(iblk, jblk, kblk) .AND. found)
THEN
3001 send_buff(isbuff)%array(:, :, :) = blk(:, :, :)
3003 CALL para_env%isend(msgin=send_buff(isbuff)%array, &
3004 dest=subgroup_dest(iblk, jblk, kblk, igroup + 1), &
3005 request=send_req(is), tag=tag)
3008 IF (para_env%mepos == subgroup_dest(iblk, jblk, kblk, igroup + 1))
THEN
3009 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
3011 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk, kblk), &
3012 request=recv_req(ir), tag=tag)
3017 IF (found)
DEALLOCATE (blk)
3025 DO kblk = 1, nblks(3)
3026 DO jblk = 1, nblks(2)
3027 DO iblk = 1, nblks(1)
3028 DO igroup = 0, ngroups - 1
3029 IF (block_source(iblk, jblk, kblk) == -1) cycle
3032 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
3035 tag = i_msg - (i_batch - 1)*batch_size
3038 ind(:) = [iblk, jblk, kblk]
3039 iatom = idx_to_at(ind(dim_at))
3040 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
3043 IF (para_env%mepos == subgroup_dest(iblk, jblk, kblk, igroup + 1) .AND. &
3044 block_source(iblk, jblk, kblk) /= subgroup_dest(iblk, jblk, kblk, igroup + 1))
THEN
3048 CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(recv_buff(tag)%array), recv_buff(tag)%array)
3057 DO i = 1, batch_size
3058 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
3059 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
3061 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
3063 CALL dbt_finalize(t3c_sub)
3065 END SUBROUTINE copy_3c_to_subgroup
3077 SUBROUTINE gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
3078 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t, ks_t_sub
3079 INTEGER,
INTENT(IN) :: group_size
3080 INTEGER,
DIMENSION(:, :, :),
INTENT(IN) :: sparsity_pattern
3084 CHARACTER(len=*),
PARAMETER :: routinen =
'gather_ks_matrix'
3086 INTEGER :: b_img, dest, handle, i, i_spin, iatom, &
3087 igroup, ir, is, jatom, n_mess, natom, &
3088 nimg, nspins, source, tag
3090 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
3091 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
3092 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
3094 CALL timeset(routinen, handle)
3096 nimg =
SIZE(sparsity_pattern, 3)
3097 natom =
SIZE(sparsity_pattern, 2)
3098 nspins =
SIZE(ks_t, 1)
3102 DO i_spin = 1, nspins
3105 IF (sparsity_pattern(iatom, jatom, b_img) > -1) n_mess = n_mess + 1
3110 ALLOCATE (send_buff(n_mess), recv_buff(n_mess))
3111 ALLOCATE (send_req(n_mess), recv_req(n_mess))
3117 DO i_spin = 1, nspins
3120 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3125 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3126 CALL dbt_get_stored_coordinates(ks_t_sub(i_spin, b_img), [iatom, jatom], source)
3127 igroup = sparsity_pattern(iatom, jatom, b_img)
3128 source = source + igroup*group_size
3129 IF (para_env%mepos == source)
THEN
3130 CALL dbt_get_block(ks_t_sub(i_spin, b_img), [iatom, jatom], blk, found)
3131 IF (source == dest)
THEN
3132 IF (found)
CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3134 ALLOCATE (send_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3135 send_buff(n_mess)%array(:, :) = 0.0_dp
3137 send_buff(n_mess)%array(:, :) = blk(:, :)
3140 CALL para_env%isend(msgin=send_buff(n_mess)%array, dest=dest, &
3141 request=send_req(is), tag=tag)
3147 IF (para_env%mepos == dest .AND. source /= dest)
THEN
3148 ALLOCATE (recv_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3150 CALL para_env%irecv(msgout=recv_buff(n_mess)%array, source=source, &
3151 request=recv_req(ir), tag=tag)
3162 DO i_spin = 1, nspins
3165 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3168 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3169 IF (para_env%mepos == dest)
THEN
3170 IF (.NOT.
ASSOCIATED(recv_buff(n_mess)%array)) cycle
3171 ALLOCATE (blk(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3172 blk(:, :) = recv_buff(n_mess)%array(:, :)
3173 CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3182 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
3183 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
3185 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
3188 CALL timestop(handle)
3190 END SUBROUTINE gather_ks_matrix
3205 SUBROUTINE get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
3206 group_size, ngroups, para_env, para_env_sub, ri_data)
3207 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3208 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work, t_2c_ao_tmp, ks_t_split
3209 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t_sub
3210 INTEGER,
INTENT(IN) :: group_size, ngroups
3214 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_tensors'
3216 INTEGER :: handle, i, i_img, i_ri, i_spin, iproc, &
3217 j, natom, nblks, nimg, nspins
3218 INTEGER(int_8) :: nze
3219 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3221 INTEGER,
DIMENSION(2) :: pdims_2d
3222 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3223 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3226 TYPE(dbt_pgrid_type) :: pgrid_2d
3227 TYPE(dbt_type) :: work, work_sub
3229 CALL timeset(routinen, handle)
3233 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3235 natom =
SIZE(ri_data%bsizes_RI)
3236 nblks =
SIZE(ri_data%bsizes_RI_split)
3237 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3238 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3239 DO i_ri = 1, ri_data%ncell_RI
3240 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3241 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3246 bsizes_ri_ext, bsizes_ri_ext, &
3248 DEALLOCATE (dist1, dist2)
3251 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3253 DEALLOCATE (dist1, dist2)
3257 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3259 DEALLOCATE (dist1, dist2)
3260 CALL dbt_create(ks_t_split(1), ks_t_split(2))
3263 ri_data%bsizes_AO, ri_data%bsizes_AO, &
3265 DEALLOCATE (dist1, dist2)
3267 nspins =
SIZE(ks_t_sub, 1)
3268 nimg =
SIZE(ks_t_sub, 2)
3270 DO i_spin = 1, nspins
3271 CALL dbt_create(t_2c_ao_tmp(1), ks_t_sub(i_spin, i_img))
3278 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3280 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3282 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3284 DO i = 0, pdims_2d(1) - 1
3285 DO j = 0, pdims_2d(2) - 1
3286 dbcsr_pgrid(i, j) = iproc
3292 ALLOCATE (col_dist(natom), row_dist(natom))
3293 row_dist(:) = dist1(:)
3294 col_dist(:) = dist2(:)
3296 ALLOCATE (ri_blk_size(natom))
3297 ri_blk_size(:) = ri_data%bsizes_RI(:)
3300 row_dist=row_dist, col_dist=col_dist)
3301 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3302 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3305 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3306 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3310 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3311 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3312 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3313 CALL dbt_clear(work_sub)
3316 CALL dbt_destroy(work)
3317 CALL dbt_destroy(work_sub)
3318 CALL dbt_pgrid_destroy(pgrid_2d)
3320 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3321 CALL timestop(handle)
3323 END SUBROUTINE get_subgroup_2c_tensors
3338 SUBROUTINE get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
3339 group_size, ngroups, para_env, para_env_sub, ri_data)
3340 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_int, t_3c_work_2, t_3c_work_3
3341 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, t_3c_apc_sub
3342 INTEGER,
INTENT(IN) :: group_size, ngroups
3346 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_tensors'
3348 INTEGER :: batch_size, bo(2), handle, handle2, &
3349 i_blk, i_img, i_ri, i_spin, ib, natom, &
3350 nblks_ao, nblks_ri, nimg, nspins
3351 INTEGER(int_8) :: nze
3352 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3353 bsizes_stack, bsizes_tmp, dist1, &
3354 dist2, dist3, dist_stack, idx_to_at
3355 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :) :: subgroup_dest
3356 INTEGER,
DIMENSION(3) :: pdims
3358 TYPE(dbt_distribution_type) :: t_dist
3359 TYPE(dbt_pgrid_type) :: pgrid
3360 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3362 CALL timeset(routinen, handle)
3364 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3365 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3366 DO i_ri = 1, ri_data%ncell_RI
3367 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3371 natom =
SIZE(ri_data%bsizes_RI)
3373 ALLOCATE (bsizes_tmp(nblks_ri))
3374 DO i_blk = 1, nblks_ri
3375 bo =
get_limit(natom, nblks_ri, i_blk - 1)
3376 bsizes_tmp(i_blk) = sum(ri_data%bsizes_RI(bo(1):bo(2)))
3378 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3379 DO i_ri = 1, ri_data%ncell_RI
3380 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = bsizes_tmp(:)
3383 batch_size = ri_data%kp_stack_size
3384 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3385 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3386 DO ib = 1, batch_size
3387 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3391 natom =
SIZE(ri_data%bsizes_RI)
3393 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3394 tensor_dims=[
SIZE(bsizes_ri_ext_split), 1, batch_size*
SIZE(ri_data%bsizes_AO_split)])
3398 pgrid, bsizes_ri_ext_split, ri_data%bsizes_AO_split, &
3399 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(RI | AO AO)")
3400 nimg =
SIZE(t_3c_int)
3402 CALL dbt_create(t_3c_int(1), t_3c_int(i_img))
3406 ALLOCATE (dist_stack(batch_size*nblks_ao))
3407 DO ib = 1, batch_size
3408 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3411 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3412 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3413 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3414 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3415 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3416 CALL dbt_distribution_destroy(t_dist)
3417 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3421 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3422 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3423 DEALLOCATE (dist1, dist2, dist3)
3426 ri_data%pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3427 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3428 DEALLOCATE (dist1, dist2, dist3)
3430 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3431 group_size, ngroups, para_env)
3434 CALL timeset(routinen//
"_ints", handle2)
3435 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
3437 CALL dbt_copy(ri_data%kp_t_3c_int(i_img), t_3c_int(i_img), move_data=.true.)
3440 ALLOCATE (ri_data%kp_t_3c_int(nimg))
3442 CALL dbt_create(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img))
3445 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, i_img), work_atom_block, order=[2, 1, 3])
3446 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3447 ngroups, para_env, subgroup_dest)
3448 CALL dbt_copy(work_atom_block_sub, t_3c_int(i_img), move_data=.true.)
3449 CALL dbt_filter(t_3c_int(i_img), ri_data%filter_eps)
3452 CALL timestop(handle2)
3453 CALL dbt_pgrid_destroy(pgrid)
3454 CALL dbt_destroy(work_atom_block)
3455 CALL dbt_destroy(work_atom_block_sub)
3456 DEALLOCATE (subgroup_dest)
3460 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3461 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3465 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3466 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3467 DEALLOCATE (dist1, dist2, dist3)
3470 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3471 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3472 DEALLOCATE (dist1, dist2, dist3)
3474 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3475 group_size, ngroups, para_env)
3479 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3480 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3483 ALLOCATE (dist_stack(batch_size*nblks_ao))
3484 DO ib = 1, batch_size
3485 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3488 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3489 CALL dbt_create(t_3c_work_2(1),
"work_2_stack", t_dist, [1], [2, 3], &
3490 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3491 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3492 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3493 CALL dbt_distribution_destroy(t_dist)
3494 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3497 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3499 nspins =
SIZE(t_3c_apc, 1)
3500 CALL timeset(routinen//
"_apc", handle2)
3502 DO i_spin = 1, nspins
3503 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3506 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3507 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, ngroups, para_env, &
3508 subgroup_dest, ri_data%iatom_to_subgroup, 1, idx_to_at)
3509 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3510 CALL dbt_filter(t_3c_apc_sub(i_spin, i_img), ri_data%filter_eps)
3512 DO i_spin = 1, nspins
3513 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3516 CALL timestop(handle2)
3517 CALL dbt_pgrid_destroy(pgrid)
3518 CALL dbt_destroy(tmp)
3519 CALL dbt_destroy(work_atom_block)
3520 CALL dbt_destroy(work_atom_block_sub)
3522 CALL timestop(handle)
3524 END SUBROUTINE get_subgroup_3c_tensors
3546 SUBROUTINE get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
3547 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
3548 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
3549 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_inv, t_2c_bint, t_2c_metric
3550 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3551 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work
3552 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
3554 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot, mat_der_pot_sub
3555 INTEGER,
INTENT(IN) :: group_size, ngroups
3559 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_derivs'
3561 INTEGER :: handle, i, i_img, i_ri, i_spin, i_xyz, &
3562 iatom, iproc, j, natom, nblks, nimg, &
3564 INTEGER(int_8) :: nze
3565 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3567 INTEGER,
DIMENSION(2) :: pdims_2d
3568 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3569 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3572 TYPE(dbt_pgrid_type) :: pgrid_2d
3573 TYPE(dbt_type) :: work, work_sub
3575 CALL timeset(routinen, handle)
3580 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3582 natom =
SIZE(ri_data%bsizes_RI)
3583 nblks =
SIZE(ri_data%bsizes_RI_split)
3584 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3585 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3586 DO i_ri = 1, ri_data%ncell_RI
3587 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3588 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3593 bsizes_ri_ext, bsizes_ri_ext, &
3595 DEALLOCATE (dist1, dist2)
3597 CALL dbt_create(t_2c_inv(1), t_2c_bint(1))
3598 CALL dbt_create(t_2c_inv(1), t_2c_metric(1))
3600 CALL dbt_create(t_2c_inv(1), t_2c_inv(iatom))
3601 CALL dbt_create(t_2c_inv(1), t_2c_bint(iatom))
3602 CALL dbt_create(t_2c_inv(1), t_2c_metric(iatom))
3604 CALL dbt_create(t_2c_inv(1), t_2c_work(1))
3605 CALL dbt_create(t_2c_inv(1), t_2c_work(2))
3606 CALL dbt_create(t_2c_inv(1), t_2c_work(3))
3607 CALL dbt_create(t_2c_inv(1), t_2c_work(4))
3610 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3612 DEALLOCATE (dist1, dist2)
3616 CALL copy_2c_to_subgroup(t_2c_inv(iatom), ri_data%t_2c_inv(1, iatom), group_size, ngroups, para_env)
3617 CALL copy_2c_to_subgroup(t_2c_bint(iatom), ri_data%t_2c_int(1, iatom), group_size, ngroups, para_env)
3618 CALL copy_2c_to_subgroup(t_2c_metric(iatom), ri_data%t_2c_pot(1, iatom), group_size, ngroups, para_env)
3624 CALL dbt_create(t_2c_inv(1), t_2c_der_metric_sub(iatom, i_xyz))
3625 CALL copy_2c_to_subgroup(t_2c_der_metric_sub(iatom, i_xyz), t_2c_der_metric(iatom, i_xyz), &
3626 group_size, ngroups, para_env)
3627 CALL dbt_destroy(t_2c_der_metric(iatom, i_xyz))
3633 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3635 DEALLOCATE (dist1, dist2)
3636 nspins =
SIZE(rho_ao_t, 1)
3637 nimg =
SIZE(rho_ao_t, 2)
3640 DO i_spin = 1, nspins
3641 IF (.NOT. (i_img == 1 .AND. i_spin == 1))
THEN
3642 CALL dbt_create(rho_ao_t_sub(1, 1), rho_ao_t_sub(i_spin, i_img))
3644 CALL copy_2c_to_subgroup(rho_ao_t_sub(i_spin, i_img), rho_ao_t(i_spin, i_img), &
3645 group_size, ngroups, para_env)
3646 CALL dbt_destroy(rho_ao_t(i_spin, i_img))
3652 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3654 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3656 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3658 DO i = 0, pdims_2d(1) - 1
3659 DO j = 0, pdims_2d(2) - 1
3660 dbcsr_pgrid(i, j) = iproc
3666 ALLOCATE (col_dist(natom), row_dist(natom))
3667 row_dist(:) = dist1(:)
3668 col_dist(:) = dist2(:)
3670 ALLOCATE (ri_blk_size(natom))
3671 ri_blk_size(:) = ri_data%bsizes_RI(:)
3674 row_dist=row_dist, col_dist=col_dist)
3675 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3676 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3680 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3681 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3685 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3686 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3687 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3688 CALL dbt_clear(work_sub)
3694 CALL dbcsr_create(mat_der_pot_sub(i_img, i_xyz), template=mat_2c_pot(1))
3695 CALL dbt_copy_matrix_to_tensor(mat_der_pot(i_img, i_xyz), work)
3700 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3701 CALL dbt_copy_tensor_to_matrix(work_sub, mat_der_pot_sub(i_img, i_xyz))
3702 CALL dbcsr_filter(mat_der_pot_sub(i_img, i_xyz), ri_data%filter_eps)
3703 CALL dbt_clear(work_sub)
3707 CALL dbt_destroy(work)
3708 CALL dbt_destroy(work_sub)
3709 CALL dbt_pgrid_destroy(pgrid_2d)
3711 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3713 CALL timestop(handle)
3715 END SUBROUTINE get_subgroup_2c_derivs
3735 SUBROUTINE get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_AO, t_3c_der_AO_sub, &
3736 t_3c_der_RI, t_3c_der_RI_sub, t_3c_apc, t_3c_apc_sub, &
3737 t_3c_der_stack, group_size, ngroups, para_env, para_env_sub, &
3739 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_work_2, t_3c_work_3
3740 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ao, t_3c_der_ao_sub, &
3741 t_3c_der_ri, t_3c_der_ri_sub, &
3742 t_3c_apc, t_3c_apc_sub
3743 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_der_stack
3744 INTEGER,
INTENT(IN) :: group_size, ngroups
3748 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_derivs'
3750 INTEGER :: batch_size, handle, i_img, i_ri, i_spin, &
3751 i_xyz, ib, nblks_ao, nblks_ri, nimg, &
3753 INTEGER(int_8) :: nze
3754 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3755 bsizes_stack, dist1, dist2, dist3, &
3756 dist_stack, idx_to_at
3757 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :) :: subgroup_dest
3759 TYPE(dbt_distribution_type) :: t_dist
3760 TYPE(dbt_pgrid_type) :: pgrid
3761 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3763 CALL timeset(routinen, handle)
3766 nblks_ri =
SIZE(ri_data%bsizes_RI)
3767 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3768 DO i_ri = 1, ri_data%ncell_RI
3769 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
3772 CALL dbt_get_info(ri_data%kp_t_3c_int(1), pdims=pdims)
3773 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3776 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3777 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3778 DEALLOCATE (dist1, dist2, dist3)
3781 ri_data%pgrid_2, bsizes_ri_ext, ri_data%bsizes_AO, &
3782 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3783 DEALLOCATE (dist1, dist2, dist3)
3784 CALL dbt_pgrid_destroy(pgrid)
3786 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3787 group_size, ngroups, para_env)
3793 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ao_sub(i_img, i_xyz))
3797 CALL dbt_copy(t_3c_der_ao(i_img, i_xyz), work_atom_block, move_data=.true.)
3798 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3799 ngroups, para_env, subgroup_dest)
3800 CALL dbt_copy(work_atom_block_sub, t_3c_der_ao_sub(i_img, i_xyz), move_data=.true.)
3801 CALL dbt_filter(t_3c_der_ao_sub(i_img, i_xyz), ri_data%filter_eps)
3805 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ri_sub(i_img, i_xyz))
3809 CALL dbt_copy(t_3c_der_ri(i_img, i_xyz), work_atom_block, move_data=.true.)
3810 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3811 ngroups, para_env, subgroup_dest)
3812 CALL dbt_copy(work_atom_block_sub, t_3c_der_ri_sub(i_img, i_xyz), move_data=.true.)
3813 CALL dbt_filter(t_3c_der_ri_sub(i_img, i_xyz), ri_data%filter_eps)
3817 CALL dbt_destroy(t_3c_der_ri(i_img, i_xyz))
3818 CALL dbt_destroy(t_3c_der_ao(i_img, i_xyz))
3821 CALL dbt_destroy(work_atom_block_sub)
3822 CALL dbt_destroy(work_atom_block)
3823 DEALLOCATE (subgroup_dest)
3826 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3827 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3828 DO i_ri = 1, ri_data%ncell_RI
3829 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3833 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3834 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3837 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3838 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3839 DEALLOCATE (dist1, dist2, dist3)
3842 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3843 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3844 DEALLOCATE (dist1, dist2, dist3)
3847 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3848 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3849 DEALLOCATE (dist1, dist2, dist3)
3851 CALL get_3c_subgroup_dest(subgroup_dest, work_atom_block_sub, work_atom_block, &
3852 group_size, ngroups, para_env)
3854 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3856 nspins =
SIZE(t_3c_apc, 1)
3858 DO i_spin = 1, nspins
3859 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3862 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3863 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, ngroups, para_env, &
3864 subgroup_dest, ri_data%iatom_to_subgroup, 1, idx_to_at)
3865 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3866 CALL dbt_filter(t_3c_apc_sub(i_spin, i_img), ri_data%filter_eps)
3868 DO i_spin = 1, nspins
3869 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3872 CALL dbt_destroy(tmp)
3873 CALL dbt_destroy(work_atom_block)
3874 CALL dbt_destroy(work_atom_block_sub)
3875 CALL dbt_pgrid_destroy(pgrid)
3878 batch_size = ri_data%kp_stack_size
3879 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3880 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3881 DO ib = 1, batch_size
3882 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3885 ALLOCATE (dist1(ri_data%ncell_RI*nblks_ri), dist2(nblks_ao), dist3(nblks_ao))
3886 CALL dbt_get_info(ri_data%kp_t_3c_int(1), proc_dist_1=dist1, proc_dist_2=dist2, &
3887 proc_dist_3=dist3, pdims=pdims)
3889 ALLOCATE (dist_stack(batch_size*nblks_ao))
3890 DO ib = 1, batch_size
3891 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3894 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3895 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3896 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3897 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3898 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3899 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3900 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(4))
3901 CALL dbt_distribution_destroy(t_dist)
3902 CALL dbt_pgrid_destroy(pgrid)
3903 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3906 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(1))
3907 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(2))
3908 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(3))
3909 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(4))
3910 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(5))
3911 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(6))
3914 ALLOCATE (dist1(nblks_ao), dist2(ri_data%ncell_RI*nblks_ri), dist3(nblks_ao))
3915 CALL dbt_get_info(t_3c_apc_sub(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, &
3916 proc_dist_3=dist3, pdims=pdims)
3918 ALLOCATE (dist_stack(batch_size*nblks_ao))
3919 DO ib = 1, batch_size
3920 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3923 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3924 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3925 CALL dbt_create(t_3c_work_2(1),
"work_3_stack", t_dist, [1], [2, 3], &
3926 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3927 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3928 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3929 CALL dbt_distribution_destroy(t_dist)
3930 CALL dbt_pgrid_destroy(pgrid)
3931 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3933 CALL timestop(handle)
3935 END SUBROUTINE get_subgroup_3c_derivs
3943 SUBROUTINE reorder_3c_ints(t_3c_ints, ri_data)
3944 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_ints
3947 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_ints'
3949 INTEGER :: handle, i_img,
idx, idx_empty, idx_full, &
3951 INTEGER(int_8) :: nze
3953 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
3955 CALL timeset(routinen, handle)
3958 ALLOCATE (t_3c_tmp(nimg))
3960 CALL dbt_create(t_3c_ints(i_img), t_3c_tmp(i_img))
3961 CALL dbt_copy(t_3c_ints(i_img), t_3c_tmp(i_img), move_data=.true.)
3966 ALLOCATE (ri_data%idx_to_img(nimg))
3968 idx_empty = nimg + 1
3973 idx_empty = idx_empty - 1
3974 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_empty), move_data=.true.)
3975 ri_data%idx_to_img(idx_empty) = i_img
3977 idx_full = idx_full + 1
3978 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_full), move_data=.true.)
3979 ri_data%idx_to_img(idx_full) = i_img
3981 CALL dbt_destroy(t_3c_tmp(i_img))
3985 ri_data%nimg_nze = idx_full
3987 ALLOCATE (ri_data%img_to_idx(nimg))
3989 ri_data%img_to_idx(ri_data%idx_to_img(
idx)) =
idx
3992 CALL timestop(handle)
3994 END SUBROUTINE reorder_3c_ints
4002 SUBROUTINE reorder_3c_derivs(t_3c_derivs, ri_data)
4003 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_derivs
4006 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_derivs'
4008 INTEGER :: handle, i_img, i_xyz,
idx, nimg
4009 INTEGER(int_8) :: nze
4011 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
4013 CALL timeset(routinen, handle)
4016 ALLOCATE (t_3c_tmp(nimg))
4018 CALL dbt_create(t_3c_derivs(1, 1), t_3c_tmp(i_img))
4023 CALL dbt_copy(t_3c_derivs(i_img, i_xyz), t_3c_tmp(i_img), move_data=.true.)
4026 idx = ri_data%img_to_idx(i_img)
4027 CALL dbt_copy(t_3c_tmp(i_img), t_3c_derivs(
idx, i_xyz), move_data=.true.)
4029 IF (nze > 0) ri_data%nimg_nze = max(
idx, ri_data%nimg_nze)
4034 CALL dbt_destroy(t_3c_tmp(i_img))
4037 CALL timestop(handle)
4039 END SUBROUTINE reorder_3c_derivs
4047 SUBROUTINE get_sparsity_pattern(pattern, ri_data, qs_env)
4048 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: pattern
4052 INTEGER :: iatom, j_img, jatom, mj_img, natom, nimg
4053 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bins
4054 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: tmp_pattern
4055 INTEGER,
DIMENSION(3) :: cell_j
4056 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
4057 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
4062 DIMENSION(:),
POINTER :: nl_iterator
4066 NULLIFY (nl_2c, nl_iterator, kpoints, cell_to_index, dft_control, index_to_cell, para_env)
4068 CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, para_env=para_env, natom=natom)
4069 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell, sab_nl=nl_2c)
4072 pattern(:, :, :) = 0
4079 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4080 IF (j_img > nimg .OR. j_img < 1) cycle
4082 mj_img = get_opp_index(j_img, qs_env)
4083 IF (mj_img > nimg .OR. mj_img < 1) cycle
4085 IF (ri_data%present_images(j_img) == 0) cycle
4087 pattern(iatom, jatom, j_img) = 1
4098 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4099 IF (j_img > nimg .OR. j_img < 1) cycle
4101 mj_img = get_opp_index(j_img, qs_env)
4102 IF (mj_img <= nimg .AND. mj_img > 0) cycle
4104 IF (ri_data%present_images(j_img) == 0) cycle
4106 pattern(iatom, jatom, j_img) = 1
4110 CALL para_env%sum(pattern)
4115 IF (pattern(iatom, iatom, j_img) /= 0)
THEN
4116 mj_img = get_opp_index(j_img, qs_env)
4117 IF (mj_img > nimg .OR. mj_img < 1) cycle
4118 pattern(iatom, iatom, mj_img) = 0
4125 ALLOCATE (bins(natom))
4128 ALLOCATE (tmp_pattern(natom, natom, nimg))
4129 tmp_pattern(:, :, :) = 0
4133 IF (pattern(iatom, jatom, j_img) == 0) cycle
4134 mj_img = get_opp_index(j_img, qs_env)
4137 IF (mj_img > nimg .OR. mj_img < 1)
THEN
4139 bins(iatom) = bins(iatom) + 1
4140 tmp_pattern(iatom, jatom, j_img) = 1
4143 IF (bins(iatom) > bins(jatom))
THEN
4144 bins(jatom) = bins(jatom) + 1
4145 tmp_pattern(jatom, iatom, mj_img) = 1
4147 bins(iatom) = bins(iatom) + 1
4148 tmp_pattern(iatom, jatom, j_img) = 1
4156 pattern(:, :, :) = tmp_pattern(:, :, :) - 1
4158 END SUBROUTINE get_sparsity_pattern
4168 SUBROUTINE get_sub_dist(sparsity_pattern, ngroups, ri_data)
4169 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: sparsity_pattern
4170 INTEGER,
INTENT(IN) :: ngroups
4173 INTEGER :: b_img, ctr, iat, iatom, igroup, jatom, &
4175 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: max_at_per_group
4177 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4179 natom =
SIZE(sparsity_pattern, 2)
4180 nimg =
SIZE(sparsity_pattern, 3)
4185 IF (.NOT.
ALLOCATED(ri_data%iatom_to_subgroup))
THEN
4186 ALLOCATE (ri_data%iatom_to_subgroup(natom), max_at_per_group(ngroups))
4188 NULLIFY (ri_data%iatom_to_subgroup(iatom)%array)
4189 ALLOCATE (ri_data%iatom_to_subgroup(iatom)%array(ngroups))
4190 ri_data%iatom_to_subgroup(iatom)%array(:) = .false.
4194 IF (ub*ngroups < natom) ub = ub + 1
4195 max_at_per_group(:) = max(1, ub)
4200 DO WHILE (
modulo(sum(max_at_per_group), natom) /= 0)
4201 igroup =
modulo(ctr, ngroups) + 1
4202 max_at_per_group(igroup) = max_at_per_group(igroup) + 1
4207 DO igroup = 1, ngroups
4208 DO iat = 1, max_at_per_group(igroup)
4209 iatom =
modulo(ctr, natom) + 1
4210 ri_data%iatom_to_subgroup(iatom)%array(igroup) = .true.
4216 ALLOCATE (bins(ngroups))
4221 IF (sparsity_pattern(iatom, jatom, b_img) == -1) cycle
4222 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4225 IF (any(ri_data%kp_cost > epsilon(0.0_dp)))
THEN
4226 cost = ri_data%kp_cost(iatom, jatom, b_img)
4228 cost = real(ri_data%bsizes_AO(iatom)*ri_data%bsizes_AO(jatom),
dp)
4230 bins(igroup + 1) = bins(igroup + 1) + cost
4231 sparsity_pattern(iatom, jatom, b_img) = igroup
4236 END SUBROUTINE get_sub_dist
4247 SUBROUTINE update_pattern_to_forces(force_pattern, scf_pattern, ngroups, ri_data, qs_env)
4248 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: force_pattern, scf_pattern
4249 INTEGER,
INTENT(IN) :: ngroups
4253 INTEGER :: b_img, iatom, igroup, jatom, mb_img, &
4255 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4257 natom =
SIZE(scf_pattern, 2)
4258 nimg =
SIZE(scf_pattern, 3)
4260 ALLOCATE (bins(ngroups))
4264 mb_img = get_opp_index(b_img, qs_env)
4268 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4271 IF (scf_pattern(iatom, jatom, b_img) > -1) cycle
4274 IF (mb_img > 0 .AND. mb_img <= nimg)
THEN
4275 IF (scf_pattern(jatom, iatom, mb_img) == -1) cycle
4276 bins(igroup + 1) = bins(igroup + 1) + ri_data%kp_cost(jatom, iatom, mb_img)
4277 force_pattern(iatom, jatom, b_img) = igroup
4283 END SUBROUTINE update_pattern_to_forces
4291 SUBROUTINE get_kp_and_ri_images(ri_data, qs_env)
4295 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_kp_and_ri_images'
4297 CHARACTER(LEN=512) :: warning_msg
4298 INTEGER :: cell_j(3), cell_k(3), handle, i_img, iatom, ikind, j_img, jatom, jcell, katom, &
4299 kcell, kp_index_lbounds(3), kp_index_ubounds(3), natom, ngroups, nimg, nkind, pcoord(3), &
4301 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_ri, &
4302 nri_per_atom, present_img, ri_cells
4303 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
4304 REAL(
dp) :: bump_fact, dij, dik, image_range, &
4305 ri_range, rij(3), rik(3)
4306 TYPE(dbt_type) :: t_dummy
4311 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4318 DIMENSION(:),
POINTER :: nl_iterator
4322 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4325 NULLIFY (qs_kind_set, dist_2d, nl_2c, nl_iterator, dft_control, &
4326 particle_set, kpoints, para_env, cell_to_index, hfx_section)
4328 CALL timeset(routinen, handle)
4330 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, &
4331 dft_control=dft_control, particle_set=particle_set, kpoints=kpoints, &
4332 para_env=para_env, natom=natom)
4333 nimg = dft_control%nimages
4335 kp_index_lbounds = lbound(cell_to_index)
4336 kp_index_ubounds = ubound(cell_to_index)
4341 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4348 CALL erfc_cutoff(ri_data%eps_schwarz, ri_data%hfx_pot%omega, ri_data%hfx_pot%cutoff_radius)
4349 WRITE (warning_msg,
'(A)') &
4350 "The SHORTANGE HFX potential typically extends over many periodic images, "// &
4351 "possibly slowing down the calculation. Consider using the TRUNCATED "// &
4352 "potential for better computational performance."
4357 ri_data%kp_RI_range = 0.0_dp
4358 ri_data%kp_image_range = 0.0_dp
4363 ri_data%kp_RI_range = max(ri_range, ri_data%kp_RI_range)
4367 CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, kind_radius=image_range)
4370 ri_data%kp_image_range = max(image_range, ri_data%kp_image_range)
4374 ri_data%kp_bump_rad = bump_fact*ri_data%kp_RI_range
4380 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
4382 ALLOCATE (present_img(nimg))
4391 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4392 IF (j_img > nimg .OR. j_img < 1) cycle
4394 IF (dij > ri_data%kp_image_range) cycle
4396 ri_data%nimg = max(j_img, ri_data%nimg)
4397 present_img(j_img) = 1
4402 CALL para_env%max(ri_data%nimg)
4403 IF (ri_data%nimg > nimg)
THEN
4404 cpabort(
"Make sure the smallest exponent of the RI-HFX basis is larger than that of the ORB basis.")
4408 CALL para_env%sum(present_img)
4409 ALLOCATE (ri_data%present_images(ri_data%nimg))
4410 ri_data%present_images = 0
4411 DO i_img = 1, ri_data%nimg
4412 IF (present_img(i_img) > 0) ri_data%present_images(i_img) = 1
4416 ri_data%pgrid, ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4417 map1=[1, 2], map2=[3], name=
"(AO AO | RI)")
4419 CALL dbt_mp_environ_pgrid(ri_data%pgrid, pdims, pcoord)
4420 CALL mp_comm_t3c%create(ri_data%pgrid%mp_comm_2d, 3, pdims)
4422 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4423 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4424 CALL dbt_destroy(t_dummy)
4429 ri_data%ri_metric,
"HFX_3c_nl", qs_env, op_pos=2, sym_ij=.false., &
4432 ALLOCATE (ri_cells(nimg))
4435 ALLOCATE (nri_per_atom(natom))
4441 iatom=iatom, jatom=jatom, katom=katom)
4444 IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
4445 any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
4447 jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4448 IF (jcell > nimg .OR. jcell < 1) cycle
4450 IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
4451 any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
4453 kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
4454 IF (kcell > nimg .OR. kcell < 1) cycle
4456 IF (dik > ri_data%kp_RI_range) cycle
4459 IF (jcell == 1 .AND. iatom == jatom) nri_per_atom(iatom) = nri_per_atom(iatom) + ri_data%bsizes_RI(katom)
4463 CALL para_env%sum(ri_cells)
4464 CALL para_env%sum(nri_per_atom)
4466 ALLOCATE (ri_data%img_to_RI_cell(nimg))
4467 ri_data%ncell_RI = 0
4468 ri_data%img_to_RI_cell = 0
4470 IF (ri_cells(i_img) > 0)
THEN
4471 ri_data%ncell_RI = ri_data%ncell_RI + 1
4472 ri_data%img_to_RI_cell(i_img) = ri_data%ncell_RI
4476 ALLOCATE (ri_data%RI_cell_to_img(ri_data%ncell_RI))
4478 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
4482 IF (ri_data%unit_nr > 0)
THEN
4483 WRITE (ri_data%unit_nr, fmt=
"(/T3,A,I29)") &
4484 "KP-HFX_RI_INFO| Number of RI-KP parallel groups:", ngroups
4485 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I29)") &
4486 "KP-HFX_RI_INFO| Tensor stack size: ", ri_data%kp_stack_size
4487 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F31.3,A)") &
4488 "KP-HFX_RI_INFO| RI basis extension radius:", ri_data%kp_RI_range*
angstrom,
" Ang"
4489 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F12.3,A, F6.3, A)") &
4490 "KP-HFX_RI_INFO| RI basis bump factor and bump radius:", bump_fact,
" /", &
4491 ri_data%kp_bump_rad*
angstrom,
" Ang"
4492 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I16,A)") &
4493 "KP-HFX_RI_INFO| The extended RI bases cover up to ", ri_data%ncell_RI,
" unit cells"
4494 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I18)") &
4495 "KP-HFX_RI_INFO| Average number of sgf in extended RI bases:", sum(nri_per_atom)/natom
4496 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F13.3,A)") &
4497 "KP-HFX_RI_INFO| Consider all image cells within a radius of ", ri_data%kp_image_range*
angstrom,
" Ang"
4498 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I27/)") &
4499 "KP-HFX_RI_INFO| Number of image cells considered: ", ri_data%nimg
4503 CALL timestop(handle)
4505 END SUBROUTINE get_kp_and_ri_images
4520 SUBROUTINE get_stack_tensors(res_stack, rho_stack, ints_stack, rho_template, ints_template, &
4521 stack_size, ri_data, qs_env)
4522 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: res_stack, rho_stack, ints_stack
4523 TYPE(dbt_type),
INTENT(INOUT) :: rho_template, ints_template
4524 INTEGER,
INTENT(IN) :: stack_size
4528 INTEGER :: is, nblks, nblks_3c(3), pdims_3d(3)
4529 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_stack, dist1, &
4530 dist2, dist3, dist_stack1, &
4531 dist_stack2, dist_stack3
4532 TYPE(dbt_distribution_type) :: t_dist
4533 TYPE(dbt_pgrid_type) :: pgrid
4540 nblks =
SIZE(ri_data%bsizes_AO_split)
4541 ALLOCATE (bsizes_stack(stack_size*nblks))
4542 DO is = 1, stack_size
4543 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
4546 ALLOCATE (dist1(nblks), dist2(nblks), dist_stack1(stack_size*nblks), dist_stack2(stack_size*nblks))
4547 CALL dbt_get_info(rho_template, proc_dist_1=dist1, proc_dist_2=dist2)
4548 DO is = 1, stack_size
4549 dist_stack1((is - 1)*nblks + 1:is*nblks) = dist1(:)
4550 dist_stack2((is - 1)*nblks + 1:is*nblks) = dist2(:)
4555 CALL dbt_distribution_new(t_dist, ri_data%pgrid_2d, dist_stack1, dist_stack2)
4556 CALL dbt_create(rho_stack(1),
"RHO_stack", t_dist, [1], [2], bsizes_stack, bsizes_stack)
4557 CALL dbt_distribution_destroy(t_dist)
4558 DEALLOCATE (dist1, dist2, dist_stack1, dist_stack2)
4561 CALL create_2c_tensor(rho_stack(2), dist1, dist2, ri_data%pgrid_2d, bsizes_stack, bsizes_stack, name=
"RHO_stack")
4562 DEALLOCATE (dist1, dist2)
4564 CALL dbt_get_info(ints_template, nblks_total=nblks_3c)
4565 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)))
4566 ALLOCATE (dist_stack3(stack_size*nblks_3c(3)), bsizes_ri_ext(nblks_3c(2)))
4567 CALL dbt_get_info(ints_template, proc_dist_1=dist1, proc_dist_2=dist2, &
4568 proc_dist_3=dist3, blk_size_2=bsizes_ri_ext)
4569 DO is = 1, stack_size
4570 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
4574 CALL dbt_distribution_new(t_dist, ri_data%pgrid_1, dist1, dist2, dist_stack3)
4575 CALL dbt_create(ints_stack(1),
"ints_stack", t_dist, [1, 2], [3], ri_data%bsizes_AO_split, &
4576 bsizes_ri_ext, bsizes_stack)
4577 CALL dbt_distribution_destroy(t_dist)
4578 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
4582 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid, tensor_dims=[nblks_3c(1), nblks_3c(2), stack_size*nblks_3c(3)])
4583 CALL create_3c_tensor(ints_stack(2), dist1, dist2, dist3, pgrid, ri_data%bsizes_AO_split, &
4584 bsizes_ri_ext, bsizes_stack, [1, 2], [3], name=
"ints_stack")
4585 DEALLOCATE (dist1, dist2, dist3)
4586 CALL dbt_pgrid_destroy(pgrid)
4589 CALL dbt_create(ints_stack(1), res_stack(1))
4590 CALL dbt_create(ints_stack(2), res_stack(2))
4592 END SUBROUTINE get_stack_tensors
4606 SUBROUTINE fill_3c_stack(t_3c_stack, t_3c_in, images, stack_dim, ri_data, filter_at, filter_dim, &
4607 idx_to_at, img_bounds)
4608 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_stack
4609 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_in
4610 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4611 INTEGER,
INTENT(IN) :: stack_dim
4613 INTEGER,
INTENT(IN),
OPTIONAL :: filter_at, filter_dim
4614 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: idx_to_at
4615 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2)
4617 INTEGER :: dest(3), i_img,
idx, ind(3), lb, nblks, &
4619 LOGICAL :: do_filter, found
4620 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4621 TYPE(dbt_iterator_type) :: iter
4626 nblks =
SIZE(ri_data%bsizes_AO_split)
4629 IF (
PRESENT(filter_at) .AND.
PRESENT(filter_dim) .AND.
PRESENT(idx_to_at)) do_filter = .true.
4634 IF (
PRESENT(img_bounds))
THEN
4636 ub = img_bounds(2) - 1
4642 IF (i_img == 0 .OR. i_img > nimg) cycle
4647 CALL dbt_iterator_start(iter, t_3c_in(i_img))
4648 DO WHILE (dbt_iterator_blocks_left(iter))
4649 CALL dbt_iterator_next_block(iter, ind)
4650 CALL dbt_get_block(t_3c_in(i_img), ind, blk, found)
4651 IF (.NOT. found) cycle
4654 IF (.NOT. idx_to_at(ind(filter_dim)) == filter_at) cycle
4657 IF (stack_dim == 1)
THEN
4658 dest = [(
idx - offset - 1)*nblks + ind(1), ind(2), ind(3)]
4659 ELSE IF (stack_dim == 2)
THEN
4660 dest = [ind(1), (
idx - offset - 1)*nblks + ind(2), ind(3)]
4662 dest = [ind(1), ind(2), (
idx - offset - 1)*nblks + ind(3)]
4665 CALL dbt_put_block(t_3c_stack, dest, shape(blk), blk)
4668 CALL dbt_iterator_stop(iter)
4671 CALL dbt_finalize(t_3c_stack)
4673 END SUBROUTINE fill_3c_stack
4685 SUBROUTINE fill_2c_stack(t_2c_stack, t_2c_in, images, stack_dim, ri_data, img_bounds, shift)
4686 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_stack
4687 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_in
4688 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4689 INTEGER,
INTENT(IN) :: stack_dim
4691 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2), shift
4693 INTEGER :: dest(2), i_img,
idx, ind(2), lb, &
4694 my_shift, nblks, nimg, offset, ub
4696 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
4697 TYPE(dbt_iterator_type) :: iter
4702 nblks =
SIZE(ri_data%bsizes_AO_split)
4707 IF (
PRESENT(img_bounds))
THEN
4709 ub = img_bounds(2) - 1
4714 IF (
PRESENT(shift)) my_shift = shift
4718 IF (i_img == 0 .OR. i_img > nimg) cycle
4722 CALL dbt_iterator_start(iter, t_2c_in(i_img))
4723 DO WHILE (dbt_iterator_blocks_left(iter))
4724 CALL dbt_iterator_next_block(iter, ind)
4725 CALL dbt_get_block(t_2c_in(i_img), ind, blk, found)
4726 IF (.NOT. found) cycle
4728 IF (stack_dim == 1)
THEN
4729 dest = [(
idx - offset - 1)*nblks + ind(1), (my_shift - 1)*nblks + ind(2)]
4731 dest = [(my_shift - 1)*nblks + ind(1), (
idx - offset - 1)*nblks + ind(2)]
4734 CALL dbt_put_block(t_2c_stack, dest, shape(blk), blk)
4737 CALL dbt_iterator_stop(iter)
4740 CALL dbt_finalize(t_2c_stack)
4742 END SUBROUTINE fill_2c_stack
4750 SUBROUTINE unstack_t_3c_apc(t_3c_apc, t_stacked, idx)
4751 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_apc, t_stacked
4752 INTEGER,
INTENT(IN) ::
idx
4754 INTEGER :: current_idx
4755 INTEGER,
DIMENSION(3) :: ind, nblks_3c
4757 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4758 TYPE(dbt_iterator_type) :: iter
4761 CALL dbt_get_info(t_3c_apc, nblks_total=nblks_3c)
4764 CALL dbt_iterator_start(iter, t_stacked)
4765 DO WHILE (dbt_iterator_blocks_left(iter))
4766 CALL dbt_iterator_next_block(iter, ind)
4769 current_idx = (ind(3) - 1)/nblks_3c(3) + 1
4770 IF (.NOT.
idx == current_idx) cycle
4772 CALL dbt_get_block(t_stacked, ind, blk, found)
4773 IF (.NOT. found) cycle
4775 CALL dbt_put_block(t_3c_apc, [ind(1), ind(2), ind(3) - (
idx - 1)*nblks_3c(3)], shape(blk), blk)
4778 CALL dbt_iterator_stop(iter)
4781 END SUBROUTINE unstack_t_3c_apc
4791 SUBROUTINE get_atom_3c_ints(t_3c_at, t_3c_ints, iatom, dim_at, idx_to_at)
4792 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_at, t_3c_ints
4793 INTEGER,
INTENT(IN) :: iatom, dim_at
4794 INTEGER,
DIMENSION(:),
INTENT(IN) :: idx_to_at
4796 INTEGER,
DIMENSION(3) :: ind
4798 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4799 TYPE(dbt_iterator_type) :: iter
4802 CALL dbt_iterator_start(iter, t_3c_ints)
4803 DO WHILE (dbt_iterator_blocks_left(iter))
4804 CALL dbt_iterator_next_block(iter, ind)
4805 IF (.NOT. idx_to_at(ind(dim_at)) == iatom) cycle
4807 CALL dbt_get_block(t_3c_ints, ind, blk, found)
4808 IF (.NOT. found) cycle
4810 CALL dbt_put_block(t_3c_at, ind, shape(blk), blk)
4813 CALL dbt_iterator_stop(iter)
4815 CALL dbt_finalize(t_3c_at)
4817 END SUBROUTINE get_atom_3c_ints
4828 SUBROUTINE precalc_derivatives(t_3c_der_RI, t_3c_der_AO, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
4829 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ri, t_3c_der_ao
4830 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot
4831 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_2c_der_metric
4835 CHARACTER(LEN=*),
PARAMETER :: routinen =
'precalc_derivatives'
4837 INTEGER :: handle, handle2, i_img, i_mem, i_ri, &
4838 i_xyz, iatom, n_mem, natom, nblks_ri, &
4839 ncell_ri, nimg, nkind, nthreads
4840 INTEGER(int_8) :: nze
4841 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, dist_ao_1, &
4842 dist_ao_2, dist_ri, dist_ri_ext, dummy_end, dummy_start, end_blocks, start_blocks
4843 INTEGER,
DIMENSION(3) :: pcoord, pdims
4844 INTEGER,
DIMENSION(:),
POINTER :: col_bsize, row_bsize
4848 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_metric
4849 TYPE(dbt_distribution_type) :: t_dist
4850 TYPE(dbt_pgrid_type) :: pgrid
4851 TYPE(dbt_type) :: t_3c_template
4852 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
4857 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4864 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4866 NULLIFY (qs_kind_set, dist_2d, nl_2c, particle_set, dft_control, para_env, row_bsize, col_bsize)
4868 CALL timeset(routinen, handle)
4870 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, natom=natom, &
4871 particle_set=particle_set, dft_control=dft_control, para_env=para_env)
4874 ncell_ri = ri_data%ncell_RI
4876 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4886 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[max(1, natom/(ri_data%n_mem*nthreads)), natom, natom])
4888 CALL create_3c_tensor(t_3c_template, dist_ao_1, dist_ao_2, dist_ri, pgrid, &
4889 ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4890 map1=[1, 2], map2=[3], name=
"tmp")
4891 CALL dbt_destroy(t_3c_template)
4894 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
4895 ALLOCATE (dist_ri_ext(natom*ncell_ri))
4896 ALLOCATE (bsizes_ri_ext(natom*ncell_ri))
4897 ALLOCATE (bsizes_ri_ext_split(nblks_ri*ncell_ri))
4898 DO i_ri = 1, ncell_ri
4899 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
4900 dist_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = dist_ri(:)
4901 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
4904 CALL dbt_distribution_new(t_dist, pgrid, dist_ao_1, dist_ao_2, dist_ri_ext)
4905 CALL dbt_create(t_3c_template,
"KP_3c_der", t_dist, [1, 2], [3], &
4906 ri_data%bsizes_AO, ri_data%bsizes_AO, bsizes_ri_ext)
4907 CALL dbt_distribution_destroy(t_dist)
4909 ALLOCATE (t_3c_der_ri_prv(nimg, 1, 3), t_3c_der_ao_prv(nimg, 1, 3))
4912 CALL dbt_create(t_3c_template, t_3c_der_ri_prv(i_img, 1, i_xyz))
4913 CALL dbt_create(t_3c_template, t_3c_der_ao_prv(i_img, 1, i_xyz))
4916 CALL dbt_destroy(t_3c_template)
4918 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
4919 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
4921 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4922 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4923 CALL dbt_pgrid_destroy(pgrid)
4926 "HFX_3c_nl", qs_env, op_pos=2, sym_jk=.false., own_dist=.true.)
4928 n_mem = ri_data%n_mem
4930 start_blocks, end_blocks)
4931 DEALLOCATE (dummy_start, dummy_end)
4933 CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, ri_data%pgrid_2, &
4934 bsizes_ri_ext_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
4935 map1=[1], map2=[2, 3], name=
"der (RI | AO AO)")
4938 CALL dbt_create(t_3c_template, t_3c_der_ri(i_img, i_xyz))
4939 CALL dbt_create(t_3c_template, t_3c_der_ao(i_img, i_xyz))
4945 nl_3c, basis_set_ao, basis_set_ao, basis_set_ri, &
4946 ri_data%ri_metric, der_eps=ri_data%eps_schwarz_forces, op_pos=2, &
4947 do_kpoints=.true., do_hfx_kpoints=.true., &
4948 bounds_k=[start_blocks(i_mem), end_blocks(i_mem)], &
4949 ri_range=ri_data%kp_RI_range, img_to_ri_cell=ri_data%img_to_RI_cell)
4951 CALL timeset(routinen//
"_cpy", handle2)
4958 CALL dbt_copy(t_3c_der_ao_prv(i_img, 1, i_xyz), t_3c_template, &
4959 order=[3, 2, 1], move_data=.true.)
4960 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4961 CALL dbt_copy(t_3c_template, t_3c_der_ao(i_img, i_xyz), &
4962 order=[1, 3, 2], move_data=.true., summation=.true.)
4968 CALL dbt_copy(t_3c_der_ri_prv(i_img, 1, i_xyz), t_3c_template, &
4969 order=[3, 2, 1], move_data=.true.)
4970 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4971 CALL dbt_copy(t_3c_template, t_3c_der_ri(i_img, i_xyz), &
4972 order=[1, 3, 2], move_data=.true., summation=.true.)
4976 CALL timestop(handle2)
4978 CALL dbt_destroy(t_3c_template)
4983 CALL dbt_destroy(t_3c_der_ri_prv(i_img, 1, i_xyz))
4984 CALL dbt_destroy(t_3c_der_ao_prv(i_img, 1, i_xyz))
4987 DEALLOCATE (t_3c_der_ri_prv, t_3c_der_ao_prv)
4990 CALL reorder_3c_derivs(t_3c_der_ri, ri_data)
4991 CALL reorder_3c_derivs(t_3c_der_ao, ri_data)
4993 CALL timeset(routinen//
"_2c", handle2)
4996 ALLOCATE (row_bsize(
SIZE(ri_data%bsizes_RI)))
4997 ALLOCATE (col_bsize(
SIZE(ri_data%bsizes_RI)))
4998 row_bsize(:) = ri_data%bsizes_RI
4999 col_bsize(:) = ri_data%bsizes_RI
5001 CALL dbcsr_create(dbcsr_template,
"2c_der", dbcsr_dist, dbcsr_type_no_symmetry, &
5002 row_bsize, col_bsize)
5004 DEALLOCATE (col_bsize, row_bsize)
5006 ALLOCATE (mat_der_metric(nimg, 3))
5009 CALL dbcsr_create(mat_der_pot(i_img, i_xyz), template=dbcsr_template)
5010 CALL dbcsr_create(mat_der_metric(i_img, i_xyz), template=dbcsr_template)
5017 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
5019 basis_set_ri, basis_set_ri, ri_data%hfx_pot, do_kpoints=.true.)
5024 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
5026 basis_set_ri, basis_set_ri, ri_data%ri_metric, do_kpoints=.true.)
5032 CALL dbt_create(ri_data%t_2c_inv(1, 1), t_2c_der_metric(iatom, i_xyz))
5033 CALL get_ext_2c_int(t_2c_der_metric(iatom, i_xyz), mat_der_metric(:, i_xyz), &
5034 iatom, iatom, 1, ri_data, qs_env)
5040 CALL timestop(handle2)
5042 CALL timestop(handle)
5044 END SUBROUTINE precalc_derivatives
5065 SUBROUTINE get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, img, pref, &
5066 ri_data, qs_env, work_virial, cell, particle_set, diag, offdiag)
5069 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_contr
5070 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_der
5071 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
5072 INTEGER,
INTENT(IN) :: img
5073 REAL(
dp),
INTENT(IN) :: pref
5076 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
5077 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
5079 POINTER :: particle_set
5080 LOGICAL,
INTENT(IN),
OPTIONAL ::
diag, offdiag
5082 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_2c_der_force'
5084 INTEGER :: handle, i_img, i_ri, i_xyz, iat, &
5085 iat_of_kind, ikind, j_img, j_ri, &
5086 j_xyz, jat, jat_of_kind, jkind, natom
5087 INTEGER,
DIMENSION(2) :: ind
5088 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
5089 LOGICAL :: found, my_diag, my_offdiag, use_virial
5090 REAL(
dp) :: new_force
5091 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :),
TARGET :: contr_blk, der_blk
5092 REAL(
dp),
DIMENSION(3) :: scoord
5093 TYPE(dbt_iterator_type) :: iter
5096 NULLIFY (kpoints, index_to_cell)
5101 CALL timeset(routinen, handle)
5103 use_virial = .false.
5104 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
5109 my_offdiag = .false.
5110 IF (
PRESENT(
diag)) my_offdiag = offdiag
5112 CALL get_qs_env(qs_env, kpoints=kpoints, natom=natom)
5121 CALL dbt_iterator_start(iter, t_2c_der(i_xyz))
5122 DO WHILE (dbt_iterator_blocks_left(iter))
5123 CALL dbt_iterator_next_block(iter, ind)
5126 IF ((my_diag .AND. .NOT. my_offdiag) .OR. (.NOT. my_diag .AND. my_offdiag))
THEN
5127 IF (my_diag .AND. (ind(1) /= ind(2))) cycle
5128 IF (my_offdiag .AND. (ind(1) == ind(2))) cycle
5131 CALL dbt_get_block(t_2c_der(i_xyz), ind, der_blk, found)
5133 CALL dbt_get_block(t_2c_contr, ind, contr_blk, found)
5139 new_force = pref*sum(der_blk(:, :)*contr_blk(:, :))
5141 i_ri = (ind(1) - 1)/natom + 1
5142 i_img = ri_data%RI_cell_to_img(i_ri)
5143 iat = ind(1) - (i_ri - 1)*natom
5144 iat_of_kind = atom_of_kind(iat)
5145 ikind = kind_of(iat)
5147 j_ri = (ind(2) - 1)/natom + 1
5148 j_img = ri_data%RI_cell_to_img(j_ri)
5149 jat = ind(2) - (j_ri - 1)*natom
5150 jat_of_kind = atom_of_kind(jat)
5151 jkind = kind_of(jat)
5155 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5158 IF (use_virial)
THEN
5161 scoord(:) = scoord(:) + real(index_to_cell(:, i_img),
dp)
5165 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5171 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5174 IF (use_virial)
THEN
5177 scoord(:) = scoord(:) + real(index_to_cell(:, j_img) + index_to_cell(:, img),
dp)
5181 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - new_force*scoord(j_xyz)
5185 DEALLOCATE (contr_blk)
5188 DEALLOCATE (der_blk)
5190 CALL dbt_iterator_stop(iter)
5194 CALL timestop(handle)
5220 idx_to_at_RI, idx_to_at_AO, i_images, lb_img, pref, &
5221 ri_data, qs_env, work_virial, cell, particle_set)
5224 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_contr
5225 TYPE(dbt_type),
DIMENSION(3),
INTENT(INOUT) :: t_3c_der_1, t_3c_der_2
5226 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of, idx_to_at_ri, &
5227 idx_to_at_ao, i_images
5228 INTEGER,
INTENT(IN) :: lb_img
5229 REAL(
dp),
INTENT(IN) :: pref
5232 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
5233 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
5235 POINTER :: particle_set
5237 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_force_from_3c_trace'
5239 INTEGER :: handle, i_ri, i_xyz, iat, iat_of_kind,
idx, ikind, j_xyz, jat, jat_of_kind, &
5240 jkind, kat, kat_of_kind, kkind, nblks_ao, nblks_ri, ri_img
5241 INTEGER,
DIMENSION(3) :: ind
5242 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
5243 LOGICAL :: found, found_1, found_2, use_virial
5244 REAL(
dp) :: new_force
5245 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :),
TARGET :: contr_blk, der_blk_1, der_blk_2, &
5247 REAL(
dp),
DIMENSION(3) :: scoord
5248 TYPE(dbt_iterator_type) :: iter
5251 NULLIFY (kpoints, index_to_cell)
5253 CALL timeset(routinen, handle)
5258 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
5259 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
5261 use_virial = .false.
5262 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
5269 CALL dbt_iterator_start(iter, t_3c_contr)
5270 DO WHILE (dbt_iterator_blocks_left(iter))
5271 CALL dbt_iterator_next_block(iter, ind)
5273 CALL dbt_get_block(t_3c_contr, ind, contr_blk, found)
5277 CALL dbt_get_block(t_3c_der_1(i_xyz), ind, der_blk_1, found_1)
5278 IF (.NOT. found_1)
THEN
5279 DEALLOCATE (der_blk_1)
5280 ALLOCATE (der_blk_1(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5281 der_blk_1(:, :, :) = 0.0_dp
5283 CALL dbt_get_block(t_3c_der_2(i_xyz), ind, der_blk_2, found_2)
5284 IF (.NOT. found_2)
THEN
5285 DEALLOCATE (der_blk_2)
5286 ALLOCATE (der_blk_2(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5287 der_blk_2(:, :, :) = 0.0_dp
5290 ALLOCATE (der_blk_3(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5291 der_blk_3(:, :, :) = -(der_blk_1(:, :, :) + der_blk_2(:, :, :))
5297 new_force = pref*sum(der_blk_1(:, :, :)*contr_blk(:, :, :))
5299 i_ri = (ind(1) - 1)/nblks_ri + 1
5300 ri_img = ri_data%RI_cell_to_img(i_ri)
5301 iat = idx_to_at_ri(ind(1) - (i_ri - 1)*nblks_ri)
5302 iat_of_kind = atom_of_kind(iat)
5303 ikind = kind_of(iat)
5306 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5309 IF (use_virial)
THEN
5312 scoord(:) = scoord(:) + real(index_to_cell(:, ri_img),
dp)
5316 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5321 new_force = pref*sum(der_blk_2(:, :, :)*contr_blk(:, :, :))
5322 jat = idx_to_at_ao(ind(2))
5323 jat_of_kind = atom_of_kind(jat)
5324 jkind = kind_of(jat)
5327 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5330 IF (use_virial)
THEN
5336 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5342 new_force = pref*sum(der_blk_3(:, :, :)*contr_blk(:, :, :))
5343 idx = (ind(3) - 1)/nblks_ao + 1
5344 kat = idx_to_at_ao(ind(3) - (
idx - 1)*nblks_ao)
5345 kat_of_kind = atom_of_kind(kat)
5346 kkind = kind_of(kat)
5349 force(kkind)%fock_4c(i_xyz, kat_of_kind) = force(kkind)%fock_4c(i_xyz, kat_of_kind) &
5352 IF (use_virial)
THEN
5354 scoord(:) = scoord(:) + real(index_to_cell(:, i_images(lb_img - 1 +
idx)),
dp)
5358 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5362 DEALLOCATE (der_blk_1, der_blk_2, der_blk_3)
5364 DEALLOCATE (contr_blk)
5367 CALL dbt_iterator_stop(iter)
5369 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, ccon)
...
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)
...
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_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
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, ncgf)
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, mimic, 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, matrix_vhxc, 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, xcint_weights, 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.