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
112#include "./base/base_uses.f90"
121 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'hfx_ri_kp'
138 SUBROUTINE adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
139 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
140 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
141 TYPE(qs_environment_type),
POINTER :: qs_env
143 INTEGER :: i_img, i_RI, i_spin, iatom, natom, &
144 nblks_RI, nimg, nkind, nspins
145 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_RI_ext, dist1, dist2, dist3
146 TYPE(dft_control_type),
POINTER :: dft_control
147 TYPE(mp_para_env_type),
POINTER :: para_env
149 NULLIFY (dft_control, para_env)
156 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, para_env=para_env, nkind=nkind)
160 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
161 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
162 DO i_ri = 1, ri_data%ncell_RI
163 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
166 ALLOCATE (ri_data%t_3c_int_ctr_1(1, nimg))
168 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
169 ri_data%bsizes_AO_split, [1, 2], [3], name=
"(AO RI | AO)")
172 CALL dbt_create(ri_data%t_3c_int_ctr_1(1, 1), ri_data%t_3c_int_ctr_1(1, i_img))
174 DEALLOCATE (dist1, dist2, dist3)
176 ALLOCATE (ri_data%t_3c_int_ctr_2(1, 1))
178 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
179 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
180 DEALLOCATE (dist1, dist2, dist3)
183 DEALLOCATE (bsizes_ri_ext)
184 nblks_ri =
SIZE(ri_data%bsizes_RI)
185 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
186 DO i_ri = 1, ri_data%ncell_RI
187 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
190 ALLOCATE (ri_data%t_2c_inv(1, natom), ri_data%t_2c_int(1, natom), ri_data%t_2c_pot(1, natom))
191 CALL create_2c_tensor(ri_data%t_2c_inv(1, 1), dist1, dist2, ri_data%pgrid_2d, &
192 bsizes_ri_ext, bsizes_ri_ext, &
194 DEALLOCATE (dist1, dist2)
195 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, 1))
196 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, 1))
198 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_inv(1, iatom))
199 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, iatom))
200 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, iatom))
203 ALLOCATE (ri_data%kp_cost(natom, natom, nimg))
204 ri_data%kp_cost = 0.0_dp
207 nspins = dft_control%nspins
208 ALLOCATE (ri_data%rho_ao_t(nspins, nimg), ri_data%ks_t(nspins, nimg))
209 CALL create_2c_tensor(ri_data%rho_ao_t(1, 1), dist1, dist2, ri_data%pgrid_2d, &
210 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
212 DEALLOCATE (dist1, dist2)
214 CALL dbt_create(dbcsr_template, ri_data%ks_t(1, 1))
216 IF (nspins == 2)
THEN
217 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
218 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
221 DO i_spin = 1, nspins
222 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(i_spin, i_img))
223 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(i_spin, i_img))
227 END SUBROUTINE adapt_ri_data_to_kp
235 SUBROUTINE hfx_ri_pre_scf_kp(dbcsr_template, ri_data, qs_env)
236 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
237 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
238 TYPE(qs_environment_type),
POINTER :: qs_env
240 CHARACTER(LEN=*),
PARAMETER :: routineN =
'hfx_ri_pre_scf_kp'
242 INTEGER :: handle, i_img, iatom, natom, nimg, nkind
243 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_op_pot, t_2c_op_RI
244 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_int
245 TYPE(dft_control_type),
POINTER :: dft_control
247 NULLIFY (dft_control)
249 CALL timeset(routinen, handle)
251 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, nkind=nkind)
253 CALL cleanup_kp(ri_data)
256 IF (ri_data%flavor .NE.
ri_pmat) cpabort(
"K-points RI-HFX only with RHO flavor")
257 IF (ri_data%same_op) ri_data%same_op = .false.
258 IF (abs(ri_data%eps_pgf_orb - dft_control%qs_control%eps_pgf_orb) > 1.0e-16_dp) &
259 cpabort(
"RI%EPS_PGF_ORB and QS%EPS_PGF_ORB must be identical for RI-HFX k-points")
261 CALL get_kp_and_ri_images(ri_data, qs_env)
265 ALLOCATE (t_2c_op_pot(nimg), t_2c_op_ri(nimg))
266 ALLOCATE (t_3c_int(1, nimg))
270 CALL adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
275 CALL get_ext_2c_int(ri_data%t_2c_inv(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
279 CALL get_ext_2c_int(ri_data%t_2c_int(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, &
280 qs_env, off_diagonal=.true.)
281 CALL apply_bump(ri_data%t_2c_int(1, iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
284 CALL get_ext_2c_int(ri_data%t_2c_pot(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
285 do_inverse=.true., skip_inverse=.true.)
286 CALL apply_bump(ri_data%t_2c_pot(1, iatom), iatom, ri_data, qs_env, from_left=.true., &
287 from_right=.true., debump=.true.)
295 ALLOCATE (ri_data%kp_mat_2c_pot(1, nimg))
297 CALL dbcsr_create(ri_data%kp_mat_2c_pot(1, i_img), template=t_2c_op_pot(i_img))
298 CALL dbcsr_copy(ri_data%kp_mat_2c_pot(1, i_img), t_2c_op_pot(i_img))
304 CALL precontract_3c_ints(t_3c_int, ri_data, qs_env)
307 CALL reorder_3c_ints(ri_data%t_3c_int_ctr_1(1, :), ri_data)
309 CALL timestop(handle)
311 END SUBROUTINE hfx_ri_pre_scf_kp
317 SUBROUTINE cleanup_kp(ri_data)
318 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
322 IF (
ALLOCATED(ri_data%kp_cost))
DEALLOCATE (ri_data%kp_cost)
323 IF (
ALLOCATED(ri_data%idx_to_img))
DEALLOCATE (ri_data%idx_to_img)
324 IF (
ALLOCATED(ri_data%img_to_idx))
DEALLOCATE (ri_data%img_to_idx)
325 IF (
ALLOCATED(ri_data%present_images))
DEALLOCATE (ri_data%present_images)
326 IF (
ALLOCATED(ri_data%img_to_RI_cell))
DEALLOCATE (ri_data%img_to_RI_cell)
327 IF (
ALLOCATED(ri_data%RI_cell_to_img))
DEALLOCATE (ri_data%RI_cell_to_img)
329 IF (
ALLOCATED(ri_data%kp_mat_2c_pot))
THEN
330 DO j = 1,
SIZE(ri_data%kp_mat_2c_pot, 2)
331 DO i = 1,
SIZE(ri_data%kp_mat_2c_pot, 1)
335 DEALLOCATE (ri_data%kp_mat_2c_pot)
338 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
339 DO i = 1,
SIZE(ri_data%kp_t_3c_int)
340 CALL dbt_destroy(ri_data%kp_t_3c_int(i))
342 DEALLOCATE (ri_data%kp_t_3c_int)
345 IF (
ALLOCATED(ri_data%t_2c_inv))
THEN
346 DO j = 1,
SIZE(ri_data%t_2c_inv, 2)
347 DO i = 1,
SIZE(ri_data%t_2c_inv, 1)
348 CALL dbt_destroy(ri_data%t_2c_inv(i, j))
351 DEALLOCATE (ri_data%t_2c_inv)
354 IF (
ALLOCATED(ri_data%t_2c_int))
THEN
355 DO j = 1,
SIZE(ri_data%t_2c_int, 2)
356 DO i = 1,
SIZE(ri_data%t_2c_int, 1)
357 CALL dbt_destroy(ri_data%t_2c_int(i, j))
360 DEALLOCATE (ri_data%t_2c_int)
363 IF (
ALLOCATED(ri_data%t_2c_pot))
THEN
364 DO j = 1,
SIZE(ri_data%t_2c_pot, 2)
365 DO i = 1,
SIZE(ri_data%t_2c_pot, 1)
366 CALL dbt_destroy(ri_data%t_2c_pot(i, j))
369 DEALLOCATE (ri_data%t_2c_pot)
372 IF (
ALLOCATED(ri_data%t_3c_int_ctr_1))
THEN
373 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_1, 2)
374 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_1, 1)
375 CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
378 DEALLOCATE (ri_data%t_3c_int_ctr_1)
381 IF (
ALLOCATED(ri_data%t_3c_int_ctr_2))
THEN
382 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_2, 2)
383 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_2, 1)
384 CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
387 DEALLOCATE (ri_data%t_3c_int_ctr_2)
390 IF (
ALLOCATED(ri_data%rho_ao_t))
THEN
391 DO j = 1,
SIZE(ri_data%rho_ao_t, 2)
392 DO i = 1,
SIZE(ri_data%rho_ao_t, 1)
393 CALL dbt_destroy(ri_data%rho_ao_t(i, j))
396 DEALLOCATE (ri_data%rho_ao_t)
399 IF (
ALLOCATED(ri_data%ks_t))
THEN
400 DO j = 1,
SIZE(ri_data%ks_t, 2)
401 DO i = 1,
SIZE(ri_data%ks_t, 1)
402 CALL dbt_destroy(ri_data%ks_t(i, j))
405 DEALLOCATE (ri_data%ks_t)
408 END SUBROUTINE cleanup_kp
422 geometry_did_change, nspins, hf_fraction)
426 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
427 REAL(kind=
dp),
INTENT(OUT) :: ehfx
429 LOGICAL,
INTENT(IN) :: geometry_did_change
430 INTEGER,
INTENT(IN) :: nspins
431 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
433 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_ks_kp'
435 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_spin, iatom, &
436 iblk, igroup, jatom, mb_img, n_batch_nze, natom, ngroups, nimg, nimg_nze
437 INTEGER(int_8) :: nflop, nze
438 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges_at, batch_ranges_nze, &
440 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
441 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: sparsity_pattern
442 LOGICAL :: use_delta_p
443 REAL(
dp) :: etmp,
fac, occ, pfac, pref, t1, t2, t3, &
446 TYPE(
dbcsr_type) :: ks_desymm, rho_desymm, tmp
447 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
449 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ks_t_split, t_2c_ao_tmp, t_2c_work, &
450 t_3c_int, t_3c_work_2, t_3c_work_3
451 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: ks_t, ks_t_sub, t_3c_apc, t_3c_apc_sub
455 NULLIFY (para_env, para_env_sub, blacs_env_sub, hfx_section, dbcsr_template)
459 CALL timeset(routinen, handle)
461 CALL get_qs_env(qs_env, para_env=para_env, natom=natom)
463 IF (nspins == 1)
THEN
464 fac = 0.5_dp*hf_fraction
466 fac = 1.0_dp*hf_fraction
469 IF (geometry_did_change)
THEN
470 CALL hfx_ri_pre_scf_kp(ks_matrix(1, 1)%matrix, ri_data, qs_env)
473 nimg_nze = ri_data%nimg_nze
499 IF (.NOT. use_delta_p) pfac = 0.0_dp
500 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, pfac, ri_data, qs_env)
502 ALLOCATE (ks_t(nspins, nimg))
504 DO i_spin = 1, nspins
505 CALL dbt_create(ri_data%ks_t(1, 1), ks_t(i_spin, i_img))
509 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
510 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
518 ALLOCATE (t_3c_apc(nspins, nimg))
520 DO i_spin = 1, nspins
521 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
524 CALL contract_pmat_3c(t_3c_apc, ri_data%rho_ao_t, ri_data, qs_env)
529 ri_data%kp_stack_size = batch_size
531 IF (mod(para_env%num_pe, ngroups) .NE. 0)
THEN
532 cpwarn(
"KP_NGROUPS must be an integer divisor of the total number of MPI ranks. It was set to 1.")
536 IF ((mod(ngroups, natom) .NE. 0) .AND. (mod(natom, ngroups) .NE. 0) .AND. geometry_did_change)
THEN
537 IF (ngroups > 1)
THEN
538 cpwarn(
"Better load balancing is reached if NGROUPS is a multiple/divisor of the number of atoms")
541 group_size = para_env%num_pe/ngroups
542 igroup = para_env%mepos/group_size
544 ALLOCATE (para_env_sub)
545 CALL para_env_sub%from_split(para_env, igroup)
549 ALLOCATE (sparsity_pattern(natom, natom, nimg))
550 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
551 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
554 ALLOCATE (mat_2c_pot(nimg), ks_t_sub(nspins, nimg), t_2c_ao_tmp(1), ks_t_split(2), t_2c_work(3))
555 CALL get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
556 group_size, ngroups, para_env, para_env_sub, ri_data)
558 ALLOCATE (t_3c_int(nimg), t_3c_apc_sub(nspins, nimg), t_3c_work_2(3), t_3c_work_3(3))
559 CALL get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
560 group_size, ngroups, para_env, para_env_sub, ri_data)
564 ALLOCATE (batch_ranges_at(natom + 1))
565 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
567 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
568 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
570 batch_ranges_at(iatom) = iblk
574 n_batch_nze = nimg_nze/batch_size
575 IF (
modulo(nimg_nze, batch_size) .NE. 0) n_batch_nze = n_batch_nze + 1
576 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
577 DO i_batch = 1, n_batch_nze
578 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
580 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
582 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
583 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
584 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
585 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
588 ri_data%kp_cost(:, :, :) = 0.0_dp
589 ALLOCATE (iapc_pairs(nimg, 2))
591 CALL dbt_batched_contract_init(ks_t_split(1))
592 CALL dbt_batched_contract_init(ks_t_split(2))
595 IF (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup) cycle
597 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp
603 CALL timeset(routinen//
"_2c", handle2)
604 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
605 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
606 dbcsr_template=dbcsr_template)
607 CALL dbt_copy(t_2c_work(1), t_2c_work(2), move_data=.true.)
608 CALL dbt_filter(t_2c_work(2), ri_data%filter_eps)
609 CALL timestop(handle2)
611 CALL dbt_batched_contract_init(t_2c_work(2))
612 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env)
613 CALL timeset(routinen//
"_3c", handle2)
616 DO i_batch = 1, n_batch_nze
617 CALL fill_3c_stack(t_3c_work_3(3), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
618 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
619 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
620 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1), move_data=.true.)
622 CALL dbt_contract(1.0_dp, t_2c_work(2), t_3c_work_3(1), &
623 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
624 contract_1=[2], notcontract_1=[1], &
625 contract_2=[1], notcontract_2=[2, 3], &
626 filter_eps=ri_data%filter_eps, flop=nflop)
627 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
628 CALL dbt_copy(t_3c_work_3(2), t_3c_work_2(2), order=[2, 1, 3], move_data=.true.)
629 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1))
633 DO i_spin = 1, nspins
634 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
635 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
636 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
639 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
640 CALL dbt_contract(-pref*
fac, t_3c_work_2(1), t_3c_work_2(2), &
641 1.0_dp, ks_t_split(i_spin), map_1=[1], map_2=[2], &
642 contract_1=[2, 3], notcontract_1=[1], &
643 contract_2=[2, 3], notcontract_2=[1], &
644 filter_eps=ri_data%filter_eps, &
645 move_data=i_spin == nspins, flop=nflop)
646 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
649 CALL timestop(handle2)
650 CALL dbt_batched_contract_finalize(t_2c_work(2))
653 ri_data%kp_cost(iatom, jatom, b_img) = t4 - t3
656 CALL dbt_batched_contract_finalize(ks_t_split(1))
657 CALL dbt_batched_contract_finalize(ks_t_split(2))
659 DO i_spin = 1, nspins
660 CALL dbt_copy(ks_t_split(i_spin), t_2c_ao_tmp(1), move_data=.true.)
661 CALL dbt_copy(t_2c_ao_tmp(1), ks_t_sub(i_spin, b_img), summation=.true.)
664 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
665 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
666 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
667 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
669 CALL para_env%sum(ri_data%dbcsr_nflop)
670 CALL para_env%sum(ri_data%kp_cost)
672 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
675 CALL gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
679 CALL dbt_copy(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img), move_data=.true.)
683 CALL dbt_destroy(t_2c_ao_tmp(1))
684 CALL dbt_destroy(ks_t_split(1))
685 CALL dbt_destroy(ks_t_split(2))
686 CALL dbt_destroy(t_2c_work(1))
687 CALL dbt_destroy(t_2c_work(2))
688 CALL dbt_destroy(t_3c_work_2(1))
689 CALL dbt_destroy(t_3c_work_2(2))
690 CALL dbt_destroy(t_3c_work_2(3))
691 CALL dbt_destroy(t_3c_work_3(1))
692 CALL dbt_destroy(t_3c_work_3(2))
693 CALL dbt_destroy(t_3c_work_3(3))
695 CALL dbt_destroy(t_3c_int(i_img))
697 DO i_spin = 1, nspins
698 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
699 CALL dbt_destroy(ks_t_sub(i_spin, i_img))
702 IF (
ASSOCIATED(dbcsr_template))
THEN
704 DEALLOCATE (dbcsr_template)
709 CALL para_env_sub%free()
710 DEALLOCATE (para_env_sub)
715 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
716 DO i_spin = 1, nspins
718 CALL dbt_copy(ks_t(i_spin, b_img), ri_data%ks_t(i_spin, b_img), summation=.true.)
721 mb_img = get_opp_index(b_img, qs_env)
722 IF (mb_img > 0 .AND. mb_img .LE. nimg)
THEN
723 CALL dbt_copy(ks_t(i_spin, mb_img), ri_data%ks_t(i_spin, b_img), order=[2, 1], summation=.true.)
728 DO i_spin = 1, nspins
729 CALL dbt_destroy(ks_t(i_spin, b_img))
734 CALL dbt_create(ri_data%ks_t(1, 1), t_2c_ao_tmp(1))
735 CALL dbcsr_create(tmp, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
736 CALL dbcsr_create(ks_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
737 CALL dbcsr_create(rho_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
740 DO i_spin = 1, nspins
741 CALL dbt_filter(ri_data%ks_t(i_spin, i_img), ri_data%filter_eps)
742 CALL dbt_copy(ri_data%ks_t(i_spin, i_img), t_2c_ao_tmp(1))
743 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), ks_desymm)
744 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), tmp)
745 CALL dbcsr_add(ks_matrix(i_spin, i_img)%matrix, tmp, 1.0_dp, 1.0_dp)
747 CALL dbt_copy(ri_data%rho_ao_t(i_spin, i_img), t_2c_ao_tmp(1))
748 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), rho_desymm)
750 CALL dbcsr_dot(ks_desymm, rho_desymm, etmp)
751 ehfx = ehfx + 0.5_dp*etmp
753 IF (.NOT. use_delta_p)
CALL dbt_clear(ri_data%ks_t(i_spin, i_img))
759 CALL dbt_destroy(t_2c_ao_tmp(1))
761 CALL timestop(handle)
780 INTEGER,
INTENT(IN) :: nspins
781 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
783 LOGICAL,
INTENT(IN),
OPTIONAL :: use_virial
785 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_forces_kp'
787 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_loop, i_spin, &
788 i_xyz, iatom, iblk, igroup, j_xyz, jatom, k_xyz, n_batch, natom, ngroups, nimg, nimg_nze
789 INTEGER(int_8) :: nflop, nze
790 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, batch_ranges_at, &
791 batch_ranges_nze, dist1, dist2, &
792 i_images, idx_to_at_ao, idx_to_at_ri, &
794 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
795 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: force_pattern, sparsity_pattern
796 INTEGER,
DIMENSION(2, 1) :: bounds_iat, bounds_jat
797 LOGICAL :: use_virial_prv
798 REAL(
dp) ::
fac, occ, pref, t1, t2
799 REAL(
dp),
DIMENSION(3, 3) :: work_virial
803 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
804 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_pot, mat_der_pot_sub
806 TYPE(dbt_type) :: t_2c_r, t_2c_r_split
807 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_bint, t_2c_binv, t_2c_der_pot, &
808 t_2c_inv, t_2c_metric, t_2c_work, &
809 t_3c_der_stack, t_3c_work_2, &
811 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
812 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, &
820 NULLIFY (para_env, para_env_sub, hfx_section, blacs_env_sub, dbcsr_template, force, atomic_kind_set, &
821 virial, particle_set, cell)
823 CALL timeset(routinen, handle)
825 use_virial_prv = .false.
826 IF (
PRESENT(use_virial)) use_virial_prv = use_virial
828 IF (nspins == 1)
THEN
829 fac = 0.5_dp*hf_fraction
831 fac = 1.0_dp*hf_fraction
834 CALL get_qs_env(qs_env, natom=natom, para_env=para_env, force=force, cell=cell, virial=virial, &
835 atomic_kind_set=atomic_kind_set, particle_set=particle_set)
838 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
839 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
841 ALLOCATE (idx_to_at_ri(
SIZE(ri_data%bsizes_RI_split)))
842 CALL get_idx_to_atom(idx_to_at_ri, ri_data%bsizes_RI_split, ri_data%bsizes_RI)
845 ALLOCATE (t_3c_der_ri(nimg, 3), t_3c_der_ao(nimg, 3), mat_der_pot(nimg, 3), t_2c_der_metric(natom, 3))
849 CALL precalc_derivatives(t_3c_der_ri, t_3c_der_ao, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
852 ALLOCATE (rho_ao_t(nspins, nimg))
854 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
856 DEALLOCATE (dist1, dist2)
857 IF (nspins == 2)
CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(2, 1))
859 DO i_spin = 1, nspins
860 CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(i_spin, i_img))
863 CALL get_pmat_images(rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
866 ALLOCATE (t_3c_apc(nspins, nimg))
868 DO i_spin = 1, nspins
869 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
872 CALL contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
877 group_size = para_env%num_pe/ngroups
878 igroup = para_env%mepos/group_size
880 ALLOCATE (para_env_sub)
881 CALL para_env_sub%from_split(para_env, igroup)
885 ALLOCATE (sparsity_pattern(natom, natom, nimg))
886 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
887 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
890 ALLOCATE (t_2c_inv(natom), mat_2c_pot(nimg), rho_ao_t_sub(nspins, nimg), t_2c_work(5), &
891 t_2c_der_metric_sub(natom, 3), mat_der_pot_sub(nimg, 3), t_2c_bint(natom), &
892 t_2c_metric(natom), t_2c_binv(natom))
893 CALL get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
894 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
895 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
896 CALL dbt_create(t_2c_work(1), t_2c_r)
897 CALL dbt_create(t_2c_work(5), t_2c_r_split)
899 ALLOCATE (t_2c_der_pot(3))
901 CALL dbt_create(t_2c_r, t_2c_der_pot(i_xyz))
905 ALLOCATE (t_3c_work_2(3), t_3c_work_3(4), t_3c_der_stack(6), t_3c_der_ao_sub(nimg, 3), &
906 t_3c_der_ri_sub(nimg, 3), t_3c_apc_sub(nspins, nimg))
907 CALL get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_ao, t_3c_der_ao_sub, &
908 t_3c_der_ri, t_3c_der_ri_sub, t_3c_apc, t_3c_apc_sub, t_3c_der_stack, &
909 group_size, ngroups, para_env, para_env_sub, ri_data)
912 ALLOCATE (batch_ranges_at(natom + 1))
913 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
915 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
916 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
918 batch_ranges_at(iatom) = iblk
922 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
923 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
924 CALL dbt_batched_contract_init(t_3c_work_3(3), batch_range_2=batch_ranges_at)
925 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
926 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
929 nimg_nze = ri_data%nimg_nze
930 batch_size = ri_data%kp_stack_size
931 n_batch = nimg_nze/batch_size
932 IF (
modulo(nimg_nze, batch_size) .NE. 0) n_batch = n_batch + 1
933 ALLOCATE (batch_ranges_nze(n_batch + 1))
934 DO i_batch = 1, n_batch
935 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
937 batch_ranges_nze(n_batch + 1) = nimg_nze + 1
942 CALL dbt_create(t_2c_inv(iatom), t_2c_binv(iatom))
943 CALL dbt_copy(t_2c_inv(iatom), t_2c_binv(iatom))
944 CALL apply_bump(t_2c_binv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
945 CALL apply_bump(t_2c_inv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
950 ALLOCATE (iapc_pairs(nimg, 2), i_images(nimg))
951 ALLOCATE (force_pattern(natom, natom, nimg))
952 force_pattern(:, :, :) = -1
961 IF (i_loop == 1 .AND. (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup)) cycle
962 IF (i_loop == 2 .AND. (.NOT. force_pattern(iatom, jatom, b_img) == igroup)) cycle
965 CALL timeset(routinen//
"_2c_1", handle2)
966 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
967 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
968 dbcsr_template=dbcsr_template)
969 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_inv(jatom), &
970 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
971 contract_1=[2], notcontract_1=[1], &
972 contract_2=[1], notcontract_2=[2], &
973 filter_eps=ri_data%filter_eps, flop=nflop)
974 CALL dbt_copy(t_2c_work(2), t_2c_work(5), move_data=.true.)
975 CALL dbt_filter(t_2c_work(5), ri_data%filter_eps)
976 CALL timestop(handle2)
978 CALL timeset(routinen//
"_3c", handle2)
979 bounds_iat(:, 1) = [sum(ri_data%bsizes_AO(1:iatom - 1)) + 1, sum(ri_data%bsizes_AO(1:iatom))]
980 bounds_jat(:, 1) = [sum(ri_data%bsizes_AO(1:jatom - 1)) + 1, sum(ri_data%bsizes_AO(1:jatom))]
981 CALL dbt_clear(t_2c_r_split)
983 DO i_spin = 1, nspins
984 CALL dbt_batched_contract_init(rho_ao_t_sub(i_spin, b_img))
987 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env, i_images)
988 DO i_batch = 1, n_batch
992 CALL dbt_clear(t_3c_der_stack(i_xyz))
993 CALL fill_3c_stack(t_3c_der_stack(i_xyz), t_3c_der_ri_sub(:, i_xyz), &
994 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
995 filter_dim=2, idx_to_at=idx_to_at_ao, &
996 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
998 CALL dbt_clear(t_3c_der_stack(3 + i_xyz))
999 CALL fill_3c_stack(t_3c_der_stack(3 + i_xyz), t_3c_der_ao_sub(:, i_xyz), &
1000 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
1001 filter_dim=2, idx_to_at=idx_to_at_ao, &
1002 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1005 DO i_spin = 1, nspins
1007 CALL dbt_clear(t_3c_work_2(3))
1008 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
1009 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
1010 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1013 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
1017 CALL dbt_contract(1.0_dp, rho_ao_t_sub(i_spin, b_img), t_3c_work_2(1), &
1018 0.0_dp, t_3c_work_2(2), map_1=[1], map_2=[2, 3], &
1019 contract_1=[1], notcontract_1=[2], &
1020 contract_2=[1], notcontract_2=[2, 3], &
1021 bounds_1=bounds_iat, bounds_2=bounds_jat, &
1022 filter_eps=ri_data%filter_eps, flop=nflop)
1023 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1029 CALL dbt_copy(t_3c_work_2(2), t_3c_work_3(1), order=[2, 1, 3], move_data=.true.)
1030 CALL dbt_batched_contract_init(t_2c_work(5))
1031 CALL dbt_contract(1.0_dp, t_2c_work(5), t_3c_work_3(1), &
1032 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
1033 contract_1=[1], notcontract_1=[2], &
1034 contract_2=[1], notcontract_2=[2, 3], &
1035 filter_eps=ri_data%filter_eps, flop=nflop)
1036 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1037 CALL dbt_batched_contract_finalize(t_2c_work(5))
1040 CALL dbt_copy(t_3c_work_3(2), t_3c_work_3(4), move_data=.true.)
1041 IF (use_virial_prv)
THEN
1043 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1044 idx_to_at_ri, idx_to_at_ao, i_images, &
1045 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1046 ri_data, qs_env, work_virial, cell, particle_set)
1049 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1050 idx_to_at_ri, idx_to_at_ao, i_images, &
1051 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1054 CALL dbt_clear(t_3c_work_3(4))
1058 IF (i_loop == 2) cycle
1061 CALL fill_3c_stack(t_3c_work_3(4), ri_data%kp_t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
1062 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
1063 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1064 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(3), move_data=.true.)
1066 CALL dbt_batched_contract_init(t_2c_r_split)
1067 CALL dbt_contract(1.0_dp, t_3c_work_3(1), t_3c_work_3(3), &
1068 1.0_dp, t_2c_r_split, map_1=[1], map_2=[2], &
1069 contract_1=[2, 3], notcontract_1=[1], &
1070 contract_2=[2, 3], notcontract_2=[1], &
1071 filter_eps=ri_data%filter_eps, flop=nflop)
1072 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1073 CALL dbt_batched_contract_finalize(t_2c_r_split)
1074 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(1))
1077 DO i_spin = 1, nspins
1078 CALL dbt_batched_contract_finalize(rho_ao_t_sub(i_spin, b_img))
1080 CALL timestop(handle2)
1082 IF (i_loop == 2) cycle
1084 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp*pref
1086 CALL timeset(routinen//
"_2c_2", handle2)
1088 CALL dbt_copy(t_2c_r_split, t_2c_r, move_data=.true.)
1090 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
1091 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
1092 dbcsr_template=dbcsr_template)
1105 CALL get_ext_2c_int(t_2c_der_pot(i_xyz), mat_der_pot_sub(:, i_xyz), iatom, jatom, &
1106 b_img, ri_data, qs_env, blacs_env_ext=blacs_env_sub, &
1107 para_env_ext=para_env_sub, dbcsr_template=dbcsr_template)
1110 IF (use_virial_prv)
THEN
1111 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1112 b_img, pref, ri_data, qs_env, work_virial, cell, particle_set)
1114 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1115 b_img, pref, ri_data, qs_env)
1119 CALL dbt_clear(t_2c_der_pot(i_xyz))
1123 CALL dbt_contract(1.0_dp, t_2c_metric(iatom), t_2c_r, &
1124 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1125 contract_1=[2], notcontract_1=[1], &
1126 contract_2=[1], notcontract_2=[2], &
1127 filter_eps=ri_data%filter_eps, flop=nflop)
1128 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1129 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_work(1), &
1130 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1131 contract_1=[2], notcontract_1=[1], &
1132 contract_2=[2], notcontract_2=[1], &
1133 filter_eps=ri_data%filter_eps, flop=nflop)
1134 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1138 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(iatom), &
1139 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1140 contract_1=[2], notcontract_1=[1], &
1141 contract_2=[1], notcontract_2=[2], &
1142 filter_eps=ri_data%filter_eps, flop=nflop)
1143 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1145 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(3), &
1146 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1147 contract_1=[1], notcontract_1=[2], &
1148 contract_2=[1], notcontract_2=[2], &
1149 filter_eps=ri_data%filter_eps, flop=nflop)
1150 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1152 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1153 CALL get_2c_bump_forces(force, t_2c_work(4), iatom, atom_of_kind, kind_of, pref, &
1154 ri_data, qs_env, work_virial)
1157 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(2), &
1158 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1159 contract_1=[1], notcontract_1=[2], &
1160 contract_2=[1], notcontract_2=[2], &
1161 filter_eps=ri_data%filter_eps, flop=nflop)
1162 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1164 IF (use_virial_prv)
THEN
1165 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1166 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1167 diag=.true., offdiag=.false.)
1169 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1170 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1174 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1175 CALL apply_bump(t_2c_work(2), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1177 IF (use_virial_prv)
THEN
1178 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1179 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1180 diag=.false., offdiag=.true.)
1182 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1183 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1188 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(iatom), &
1189 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1190 contract_1=[2], notcontract_1=[1], &
1191 contract_2=[1], notcontract_2=[2], &
1192 filter_eps=ri_data%filter_eps, flop=nflop)
1193 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1195 CALL dbt_contract(1.0_dp, t_2c_bint(iatom), t_2c_work(4), &
1196 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1197 contract_1=[1], notcontract_1=[2], &
1198 contract_2=[1], notcontract_2=[2], &
1199 filter_eps=ri_data%filter_eps, flop=nflop)
1200 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1202 CALL get_2c_bump_forces(force, t_2c_work(2), iatom, atom_of_kind, kind_of, -pref, &
1203 ri_data, qs_env, work_virial)
1206 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_r, &
1207 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1208 contract_1=[1], notcontract_1=[2], &
1209 contract_2=[1], notcontract_2=[2], &
1210 filter_eps=ri_data%filter_eps, flop=nflop)
1211 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1213 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_metric(jatom), &
1214 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1215 contract_1=[2], notcontract_1=[1], &
1216 contract_2=[1], notcontract_2=[2], &
1217 filter_eps=ri_data%filter_eps, flop=nflop)
1218 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1222 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(jatom), &
1223 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1224 contract_1=[2], notcontract_1=[1], &
1225 contract_2=[1], notcontract_2=[2], &
1226 filter_eps=ri_data%filter_eps, flop=nflop)
1227 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1229 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(3), &
1230 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1231 contract_1=[1], notcontract_1=[2], &
1232 contract_2=[1], notcontract_2=[2], &
1233 filter_eps=ri_data%filter_eps, flop=nflop)
1234 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1236 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1237 CALL get_2c_bump_forces(force, t_2c_work(4), jatom, atom_of_kind, kind_of, pref, &
1238 ri_data, qs_env, work_virial)
1241 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(2), &
1242 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1243 contract_1=[1], notcontract_1=[2], &
1244 contract_2=[1], notcontract_2=[2], &
1245 filter_eps=ri_data%filter_eps, flop=nflop)
1246 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1248 IF (use_virial_prv)
THEN
1249 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1250 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1251 diag=.true., offdiag=.false.)
1253 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1254 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1258 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1259 CALL apply_bump(t_2c_work(2), jatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1261 IF (use_virial_prv)
THEN
1262 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1263 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1264 diag=.false., offdiag=.true.)
1266 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1267 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1272 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(jatom), &
1273 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1274 contract_1=[2], notcontract_1=[1], &
1275 contract_2=[1], notcontract_2=[2], &
1276 filter_eps=ri_data%filter_eps, flop=nflop)
1277 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1279 CALL dbt_contract(1.0_dp, t_2c_bint(jatom), t_2c_work(4), &
1280 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1281 contract_1=[1], notcontract_1=[2], &
1282 contract_2=[1], notcontract_2=[2], &
1283 filter_eps=ri_data%filter_eps, flop=nflop)
1284 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1286 CALL get_2c_bump_forces(force, t_2c_work(2), jatom, atom_of_kind, kind_of, -pref, &
1287 ri_data, qs_env, work_virial)
1289 CALL timestop(handle2)
1294 IF (i_loop == 1)
THEN
1295 CALL update_pattern_to_forces(force_pattern, sparsity_pattern, ngroups, ri_data, qs_env)
1299 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
1300 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
1301 CALL dbt_batched_contract_finalize(t_3c_work_3(3))
1302 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
1303 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
1305 IF (use_virial_prv)
THEN
1309 virial%pv_fock_4c(i_xyz, j_xyz) = virial%pv_fock_4c(i_xyz, j_xyz) &
1310 + work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1318 CALL para_env_sub%free()
1319 DEALLOCATE (para_env_sub)
1321 CALL para_env%sync()
1323 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
1326 IF (
ASSOCIATED(dbcsr_template))
THEN
1328 DEALLOCATE (dbcsr_template)
1330 CALL dbt_destroy(t_2c_r)
1331 CALL dbt_destroy(t_2c_r_split)
1332 CALL dbt_destroy(t_2c_work(1))
1333 CALL dbt_destroy(t_2c_work(2))
1334 CALL dbt_destroy(t_2c_work(3))
1335 CALL dbt_destroy(t_2c_work(4))
1336 CALL dbt_destroy(t_2c_work(5))
1337 CALL dbt_destroy(t_3c_work_2(1))
1338 CALL dbt_destroy(t_3c_work_2(2))
1339 CALL dbt_destroy(t_3c_work_2(3))
1340 CALL dbt_destroy(t_3c_work_3(1))
1341 CALL dbt_destroy(t_3c_work_3(2))
1342 CALL dbt_destroy(t_3c_work_3(3))
1343 CALL dbt_destroy(t_3c_work_3(4))
1344 CALL dbt_destroy(t_3c_der_stack(1))
1345 CALL dbt_destroy(t_3c_der_stack(2))
1346 CALL dbt_destroy(t_3c_der_stack(3))
1347 CALL dbt_destroy(t_3c_der_stack(4))
1348 CALL dbt_destroy(t_3c_der_stack(5))
1349 CALL dbt_destroy(t_3c_der_stack(6))
1351 CALL dbt_destroy(t_2c_der_pot(i_xyz))
1354 CALL dbt_destroy(t_2c_inv(iatom))
1355 CALL dbt_destroy(t_2c_binv(iatom))
1356 CALL dbt_destroy(t_2c_bint(iatom))
1357 CALL dbt_destroy(t_2c_metric(iatom))
1359 CALL dbt_destroy(t_2c_der_metric_sub(iatom, i_xyz))
1364 DO i_spin = 1, nspins
1365 CALL dbt_destroy(rho_ao_t_sub(i_spin, i_img))
1366 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
1371 CALL dbt_destroy(t_3c_der_ri_sub(i_img, i_xyz))
1372 CALL dbt_destroy(t_3c_der_ao_sub(i_img, i_xyz))
1377 CALL timestop(handle)
1392 SUBROUTINE apply_bump(t_2c_inout, atom_i, ri_data, qs_env, from_left, from_right, debump)
1393 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_inout
1394 INTEGER,
INTENT(IN) :: atom_i
1397 LOGICAL,
INTENT(IN),
OPTIONAL :: from_left, from_right, debump
1399 INTEGER :: i_img, i_ri, iatom, ind(2), j_img, j_ri, &
1400 jatom, natom, nblks(2), nimg, nkind
1401 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1402 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1403 LOGICAL :: found, my_debump, my_left, my_right
1404 REAL(
dp) :: bval, r0, r1, ri(3), rj(3), rref(3), &
1406 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1408 TYPE(dbt_iterator_type) :: iter
1411 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1413 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1415 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1416 kpoints=kpoints, particle_set=particle_set)
1417 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1420 IF (
PRESENT(debump)) my_debump = debump
1423 IF (
PRESENT(from_left)) my_left = from_left
1426 IF (
PRESENT(from_right)) my_right = from_right
1427 cpassert(my_left .OR. my_right)
1429 CALL dbt_get_info(t_2c_inout, nblks_total=nblks)
1430 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1431 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1436 r1 = ri_data%kp_RI_range
1437 r0 = ri_data%kp_bump_rad
1438 rref =
pbc(particle_set(atom_i)%r, cell)
1443 CALL dbt_iterator_start(iter, t_2c_inout)
1444 DO WHILE (dbt_iterator_blocks_left(iter))
1445 CALL dbt_iterator_next_block(iter, ind)
1446 CALL dbt_get_block(t_2c_inout, ind, blk, found)
1447 IF (.NOT. found) cycle
1449 i_ri = (ind(1) - 1)/natom + 1
1450 i_img = ri_data%RI_cell_to_img(i_ri)
1451 iatom = ind(1) - (i_ri - 1)*natom
1454 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
1456 j_ri = (ind(2) - 1)/natom + 1
1457 j_img = ri_data%RI_cell_to_img(j_ri)
1458 jatom = ind(2) - (j_ri - 1)*natom
1461 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1463 IF (.NOT. my_debump)
THEN
1464 IF (my_left) blk(:, :) = blk(:, :)*bump(norm2(ri - rref), r0, r1)
1465 IF (my_right) blk(:, :) = blk(:, :)*bump(norm2(rj - rref), r0, r1)
1469 bval = bump(norm2(ri - rref), r0, r1)
1470 IF (my_left .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1471 bval = bump(norm2(rj - rref), r0, r1)
1472 IF (my_right .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1475 CALL dbt_put_block(t_2c_inout, ind, shape(blk), blk)
1479 CALL dbt_iterator_stop(iter)
1481 CALL dbt_filter(t_2c_inout, ri_data%filter_eps)
1483 END SUBROUTINE apply_bump
1497 SUBROUTINE get_2c_bump_forces(force, t_2c_in, atom_i, atom_of_kind, kind_of, pref, ri_data, &
1498 qs_env, work_virial)
1500 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_in
1501 INTEGER,
INTENT(IN) :: atom_i
1502 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
1503 REAL(
dp),
INTENT(IN) :: pref
1506 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT) :: work_virial
1508 INTEGER :: i, i_img, i_ri, i_xyz, iat_of_kind, iatom, ikind, ind(2), j_img, j_ri, j_xyz, &
1509 jat_of_kind, jatom, jkind, natom, nblks(2), nimg, nkind
1510 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1511 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1513 REAL(
dp) :: new_force, r0, r1, ri(3), rj(3), &
1514 rref(3), scoord(3), x
1515 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1517 TYPE(dbt_iterator_type) :: iter
1520 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1522 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1524 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1525 kpoints=kpoints, particle_set=particle_set)
1526 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1528 CALL dbt_get_info(t_2c_in, nblks_total=nblks)
1529 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1530 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1535 r1 = ri_data%kp_RI_range
1536 r0 = ri_data%kp_bump_rad
1537 rref =
pbc(particle_set(atom_i)%r, cell)
1539 iat_of_kind = atom_of_kind(atom_i)
1540 ikind = kind_of(atom_i)
1546 CALL dbt_iterator_start(iter, t_2c_in)
1547 DO WHILE (dbt_iterator_blocks_left(iter))
1548 CALL dbt_iterator_next_block(iter, ind)
1549 IF (ind(1) .NE. ind(2)) cycle
1551 CALL dbt_get_block(t_2c_in, ind, blk, found)
1552 IF (.NOT. found) cycle
1555 j_ri = (ind(2) - 1)/natom + 1
1556 j_img = ri_data%RI_cell_to_img(j_ri)
1557 jatom = ind(2) - (j_ri - 1)*natom
1558 jat_of_kind = atom_of_kind(jatom)
1559 jkind = kind_of(jatom)
1562 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1563 x = norm2(rj - rref)
1564 IF (x < r0 .OR. x > r1) cycle
1567 DO i = 1,
SIZE(blk, 1)
1568 new_force = new_force + blk(i, i)
1570 new_force = pref*new_force*dbump(x, r0, r1)
1576 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) + &
1577 new_force*(rj(i_xyz) - rref(i_xyz))/x
1583 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1584 + new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1589 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) - &
1590 new_force*(rj(i_xyz) - rref(i_xyz))/x
1596 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1597 - new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1603 CALL dbt_iterator_stop(iter)
1606 END SUBROUTINE get_2c_bump_forces
1615 FUNCTION bump(x, r0, r1)
RESULT(b)
1616 REAL(
dp),
INTENT(IN) :: x, r0, r1
1624 r = (x - r0)/(r1 - r0)
1625 b = -6.0_dp*r**5 + 15.0_dp*r**4 - 10.0_dp*r**3 + 1.0_dp
1626 IF (x .GE. r1) b = 0.0_dp
1627 IF (x .LE. r0) b = 1.0_dp
1638 FUNCTION dbump(x, r0, r1)
RESULT(b)
1639 REAL(
dp),
INTENT(IN) :: x, r0, r1
1644 r = (x - r0)/(r1 - r0)
1645 b = (-30.0_dp*r**4 + 60.0_dp*r**3 - 30.0_dp*r**2)/(r1 - r0)
1646 IF (x .GE. r1) b = 0.0_dp
1647 IF (x .LE. r0) b = 0.0_dp
1658 FUNCTION get_apc_index_from_ib(i_index, b_index, qs_env)
RESULT(apc_index)
1659 INTEGER,
INTENT(IN) :: i_index, b_index
1661 INTEGER :: apc_index
1663 INTEGER,
DIMENSION(3) :: cell_apc
1664 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1665 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1669 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1672 cell_apc(:) = index_to_cell(:, i_index) + index_to_cell(:, b_index)
1674 IF (any([cell_apc(1), cell_apc(2), cell_apc(3)] < lbound(cell_to_index)) .OR. &
1675 any([cell_apc(1), cell_apc(2), cell_apc(3)] > ubound(cell_to_index)))
THEN
1679 apc_index = cell_to_index(cell_apc(1), cell_apc(2), cell_apc(3))
1682 END FUNCTION get_apc_index_from_ib
1691 FUNCTION get_apc_index(a_index, c_index, qs_env)
RESULT(i_index)
1692 INTEGER,
INTENT(IN) :: a_index, c_index
1696 INTEGER,
DIMENSION(3) :: cell_i
1697 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1698 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1702 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1704 cell_i(:) = index_to_cell(:, a_index) + index_to_cell(:, c_index)
1706 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1707 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1711 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1714 END FUNCTION get_apc_index
1723 FUNCTION get_i_index(apc_index, b_index, qs_env)
RESULT(i_index)
1724 INTEGER,
INTENT(IN) :: apc_index, b_index
1728 INTEGER,
DIMENSION(3) :: cell_i
1729 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1730 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1734 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1736 cell_i(:) = index_to_cell(:, apc_index) - index_to_cell(:, b_index)
1738 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1739 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1743 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1746 END FUNCTION get_i_index
1757 SUBROUTINE get_ac_pairs(ac_pairs, apc_index, ri_data, qs_env)
1758 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: ac_pairs
1759 INTEGER,
INTENT(IN) :: apc_index
1763 INTEGER :: a_index, actual_img, c_index, nimg
1765 nimg =
SIZE(ac_pairs, 1)
1770 DO a_index = 1, nimg
1771 actual_img = ri_data%idx_to_img(a_index)
1773 c_index = get_i_index(apc_index, actual_img, qs_env)
1774 ac_pairs(a_index, 1) = a_index
1775 ac_pairs(a_index, 2) = c_index
1779 END SUBROUTINE get_ac_pairs
1791 SUBROUTINE get_iapc_pairs(iapc_pairs, b_index, ri_data, qs_env, actual_i_img)
1792 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: iapc_pairs
1793 INTEGER,
INTENT(IN) :: b_index
1796 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: actual_i_img
1798 INTEGER :: actual_img, apc_index, i_index, nimg
1800 nimg =
SIZE(iapc_pairs, 1)
1801 IF (
PRESENT(actual_i_img)) actual_i_img(:) = 0
1803 iapc_pairs(:, :) = 0
1806 DO i_index = 1, nimg
1807 actual_img = ri_data%idx_to_img(i_index)
1808 apc_index = get_apc_index_from_ib(actual_img, b_index, qs_env)
1809 IF (apc_index == 0) cycle
1810 iapc_pairs(i_index, 1) = i_index
1811 iapc_pairs(i_index, 2) = apc_index
1812 IF (
PRESENT(actual_i_img)) actual_i_img(i_index) = actual_img
1815 END SUBROUTINE get_iapc_pairs
1824 FUNCTION get_opp_index(a_index, qs_env)
RESULT(opp_index)
1825 INTEGER,
INTENT(IN) :: a_index
1827 INTEGER :: opp_index
1829 INTEGER,
DIMENSION(3) :: opp_cell
1830 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1831 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1834 NULLIFY (kpoints, cell_to_index, index_to_cell)
1837 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1839 opp_cell(:) = -index_to_cell(:, a_index)
1841 IF (any([opp_cell(1), opp_cell(2), opp_cell(3)] < lbound(cell_to_index)) .OR. &
1842 any([opp_cell(1), opp_cell(2), opp_cell(3)] > ubound(cell_to_index)))
THEN
1846 opp_index = cell_to_index(opp_cell(1), opp_cell(2), opp_cell(3))
1849 END FUNCTION get_opp_index
1860 SUBROUTINE get_pmat_images(rho_ao_t, rho_ao, scale_prev_p, ri_data, qs_env)
1861 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t
1862 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao
1863 REAL(
dp),
INTENT(IN) :: scale_prev_p
1867 INTEGER :: cell_j(3), i_img, i_spin, iatom, icol, &
1868 irow, j_img, jatom, mi_img, mj_img, &
1870 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1873 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock, pblock_desymm
1874 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, rho_desymm
1875 TYPE(dbt_type) :: tmp
1879 DIMENSION(:),
POINTER :: nl_iterator
1881 POINTER :: sab_nl, sab_nl_nosym
1884 NULLIFY (rho_desymm, kpoints, sab_nl_nosym, scf_env, matrix_ks, dft_control, &
1885 sab_nl, nl_iterator, cell_to_index, pblock, pblock_desymm)
1887 CALL get_qs_env(qs_env, kpoints=kpoints, scf_env=scf_env, matrix_ks_kp=matrix_ks, dft_control=dft_control)
1888 CALL get_kpoint_info(kpoints, sab_nl_nosym=sab_nl_nosym, cell_to_index=cell_to_index, sab_nl=sab_nl)
1890 IF (dft_control%do_admm)
THEN
1891 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks)
1894 nspins =
SIZE(matrix_ks, 1)
1897 ALLOCATE (rho_desymm(nspins, nimg))
1899 DO i_spin = 1, nspins
1900 ALLOCATE (rho_desymm(i_spin, i_img)%matrix)
1901 CALL dbcsr_create(rho_desymm(i_spin, i_img)%matrix, template=matrix_ks(i_spin, i_img)%matrix, &
1902 matrix_type=dbcsr_type_no_symmetry)
1906 CALL dbt_create(rho_desymm(1, 1)%matrix, tmp)
1913 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1914 IF (j_img > nimg .OR. j_img < 1) cycle
1917 IF (iatom == jatom)
fac = 0.5_dp
1918 mj_img = get_opp_index(j_img, qs_env)
1920 IF (mj_img == 0)
fac = 1.0_dp
1924 IF (iatom > jatom)
THEN
1930 DO i_spin = 1, nspins
1932 IF (.NOT. found) cycle
1935 CALL dbcsr_get_block_p(rho_desymm(i_spin, j_img)%matrix, iatom, jatom, pblock_desymm, found)
1936 IF (.NOT. found) cycle
1938 IF (iatom > jatom)
THEN
1939 pblock_desymm(:, :) =
fac*transpose(pblock(:, :))
1941 pblock_desymm(:, :) =
fac*pblock(:, :)
1948 DO i_spin = 1, nspins
1949 CALL dbt_scale(rho_ao_t(i_spin, i_img), scale_prev_p)
1951 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, i_img)%matrix, tmp)
1952 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), summation=.true., move_data=.true.)
1955 mi_img = get_opp_index(i_img, qs_env)
1956 IF (mi_img > 0 .AND. mi_img .LE. nimg)
THEN
1957 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, mi_img)%matrix, tmp)
1958 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), order=[2, 1], summation=.true., move_data=.true.)
1960 CALL dbt_filter(rho_ao_t(i_spin, i_img), ri_data%filter_eps)
1965 DO i_spin = 1, nspins
1967 DEALLOCATE (rho_desymm(i_spin, i_img)%matrix)
1971 CALL dbt_destroy(tmp)
1972 DEALLOCATE (rho_desymm)
1974 END SUBROUTINE get_pmat_images
1993 SUBROUTINE get_ext_2c_int(t_2c_pot, mat_orig, atom_i, atom_j, img_b, ri_data, qs_env, do_inverse, &
1994 para_env_ext, blacs_env_ext, dbcsr_template, off_diagonal, skip_inverse)
1995 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_pot
1996 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_orig
1997 INTEGER,
INTENT(IN) :: atom_i, atom_j, img_b
2000 LOGICAL,
INTENT(IN),
OPTIONAL :: do_inverse
2003 TYPE(
dbcsr_type),
OPTIONAL,
POINTER :: dbcsr_template
2004 LOGICAL,
INTENT(IN),
OPTIONAL :: off_diagonal, skip_inverse
2006 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_ext_2c_int'
2008 INTEGER :: group, handle, handle2, i_img, i_ri, iatom, iblk, ikind, img_tot, j_img, j_ri, &
2009 jatom, jblk, jkind, n_dependent, natom, nblks_ri, nimg, nkind
2010 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2
2011 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: present_atoms_i, present_atoms_j
2012 INTEGER,
DIMENSION(3) :: cell_b, cell_i, cell_j, cell_tot
2013 INTEGER,
DIMENSION(:),
POINTER :: col_dist, col_dist_ext, ri_blk_size_ext, &
2014 row_dist, row_dist_ext
2015 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell, pgrid
2016 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
2017 LOGICAL :: do_inverse_prv, found, my_offd, &
2018 skip_inverse_prv, use_template
2019 REAL(
dp) :: bfac, dij, r0, r1, threshold
2020 REAL(
dp),
DIMENSION(3) :: ri, rij, rj, rref, scoord
2021 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock
2026 TYPE(
dbcsr_type) :: work, work_tight, work_tight_inv
2027 TYPE(dbt_type) :: t_2c_tmp
2030 DIMENSION(:),
TARGET :: basis_set_ri
2034 DIMENSION(:),
POINTER :: nl_iterator
2038 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2040 NULLIFY (qs_kind_set, nl_2c, nl_iterator, cell, kpoints, cell_to_index, index_to_cell, dist_2d, &
2041 para_env, pblock, blacs_env, particle_set, col_dist, row_dist, pgrid, &
2042 col_dist_ext, row_dist_ext)
2044 CALL timeset(routinen, handle)
2049 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
2050 kpoints=kpoints, para_env=para_env, blacs_env=blacs_env, particle_set=particle_set)
2051 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
2053 do_inverse_prv = .false.
2054 IF (
PRESENT(do_inverse)) do_inverse_prv = do_inverse
2055 IF (do_inverse_prv)
THEN
2056 cpassert(atom_i == atom_j)
2059 skip_inverse_prv = .false.
2060 IF (
PRESENT(skip_inverse)) skip_inverse_prv = skip_inverse
2063 IF (
PRESENT(off_diagonal)) my_offd = off_diagonal
2065 IF (
PRESENT(para_env_ext)) para_env => para_env_ext
2066 IF (
PRESENT(blacs_env_ext)) blacs_env => blacs_env_ext
2068 nimg =
SIZE(mat_orig)
2070 CALL timeset(routinen//
"_nl_iter", handle2)
2073 ALLOCATE (dist1(natom), dist2(natom))
2075 dist1(iatom) = mod(iatom, blacs_env%num_pe(1))
2076 dist2(iatom) = mod(iatom, blacs_env%num_pe(2))
2080 ALLOCATE (basis_set_ri(nkind))
2084 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
2086 ALLOCATE (present_atoms_i(natom, nimg), present_atoms_j(natom, nimg))
2092 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij, cell=cell_j, &
2093 ikind=ikind, jkind=jkind)
2097 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2098 IF (j_img > nimg .OR. j_img < 1) cycle
2100 IF (iatom == atom_i .AND. dij .LE. ri_data%kp_RI_range) present_atoms_i(jatom, j_img) = 1
2101 IF (iatom == atom_j .AND. dij .LE. ri_data%kp_RI_range) present_atoms_j(jatom, j_img) = 1
2106 CALL timestop(handle2)
2108 CALL para_env%sum(present_atoms_i)
2109 CALL para_env%sum(present_atoms_j)
2113 use_template = .false.
2114 IF (
PRESENT(dbcsr_template))
THEN
2115 IF (
ASSOCIATED(dbcsr_template)) use_template = .true.
2118 IF (use_template)
THEN
2123 ALLOCATE (row_dist_ext(ri_data%ncell_RI*natom), col_dist_ext(ri_data%ncell_RI*natom))
2124 ALLOCATE (ri_blk_size_ext(ri_data%ncell_RI*natom))
2125 DO i_ri = 1, ri_data%ncell_RI
2126 row_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = row_dist(:)
2127 col_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = col_dist(:)
2128 ri_blk_size_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2132 row_dist=row_dist_ext, col_dist=col_dist_ext)
2133 CALL dbcsr_create(work, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2134 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2136 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2138 IF (
PRESENT(dbcsr_template))
THEN
2139 ALLOCATE (dbcsr_template)
2144 cell_b(:) = index_to_cell(:, img_b)
2146 i_ri = ri_data%img_to_RI_cell(i_img)
2147 IF (i_ri == 0) cycle
2148 cell_i(:) = index_to_cell(:, i_img)
2150 j_ri = ri_data%img_to_RI_cell(j_img)
2151 IF (j_ri == 0) cycle
2152 cell_j(:) = index_to_cell(:, j_img)
2153 cell_tot = cell_j - cell_i + cell_b
2155 IF (any([cell_tot(1), cell_tot(2), cell_tot(3)] < lbound(cell_to_index)) .OR. &
2156 any([cell_tot(1), cell_tot(2), cell_tot(3)] > ubound(cell_to_index))) cycle
2157 img_tot = cell_to_index(cell_tot(1), cell_tot(2), cell_tot(3))
2158 IF (img_tot > nimg .OR. img_tot < 1) cycle
2163 IF (present_atoms_i(iatom, i_img) == 0) cycle
2164 IF (present_atoms_j(jatom, j_img) == 0) cycle
2165 IF (my_offd .AND. (i_ri - 1)*natom + iatom == (j_ri - 1)*natom + jatom) cycle
2168 IF (.NOT. found) cycle
2170 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2179 IF (do_inverse_prv)
THEN
2181 r1 = ri_data%kp_RI_range
2182 r0 = ri_data%kp_bump_rad
2185 nblks_ri = sum(present_atoms_i)
2186 ALLOCATE (col_dist_ext(nblks_ri), row_dist_ext(nblks_ri), ri_blk_size_ext(nblks_ri))
2189 i_ri = ri_data%img_to_RI_cell(i_img)
2190 IF (i_ri == 0) cycle
2192 IF (present_atoms_i(iatom, i_img) == 0) cycle
2194 col_dist_ext(iblk) = col_dist(iatom)
2195 row_dist_ext(iblk) = row_dist(iatom)
2196 ri_blk_size_ext(iblk) = ri_data%bsizes_RI(iatom)
2201 row_dist=row_dist_ext, col_dist=col_dist_ext)
2202 CALL dbcsr_create(work_tight, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2203 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2204 CALL dbcsr_create(work_tight_inv, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2205 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2207 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2211 rref =
pbc(particle_set(atom_i)%r, cell)
2215 i_ri = ri_data%img_to_RI_cell(i_img)
2216 IF (i_ri == 0) cycle
2218 IF (present_atoms_i(iatom, i_img) == 0) cycle
2222 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
2226 j_ri = ri_data%img_to_RI_cell(j_img)
2227 IF (j_ri == 0) cycle
2229 IF (present_atoms_j(jatom, j_img) == 0) cycle
2233 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
2235 CALL dbcsr_get_block_p(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock, found)
2236 IF (.NOT. found) cycle
2239 IF (iblk .NE. jblk) bfac = bump(norm2(ri - rref), r0, r1)*bump(norm2(rj - rref), r0, r1)
2248 IF (.NOT. skip_inverse_prv)
THEN
2249 SELECT CASE (ri_data%t2c_method)
2251 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
2252 CALL invert_hotelling(work_tight_inv, work_tight, threshold=threshold, silent=.false.)
2257 uplo_to_full=.true.)
2260 CALL cp_dbcsr_power(work_tight_inv, -1.0_dp, ri_data%eps_eigval, n_dependent, &
2261 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
2272 i_ri = ri_data%img_to_RI_cell(i_img)
2273 IF (i_ri == 0) cycle
2275 IF (present_atoms_i(iatom, i_img) == 0) cycle
2280 j_ri = ri_data%img_to_RI_cell(j_img)
2281 IF (j_ri == 0) cycle
2283 IF (present_atoms_j(jatom, j_img) == 0) cycle
2287 IF (.NOT. found) cycle
2289 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2300 CALL dbt_create(work, t_2c_tmp)
2301 CALL dbt_copy_matrix_to_tensor(work, t_2c_tmp)
2302 CALL dbt_copy(t_2c_tmp, t_2c_pot, move_data=.true.)
2303 CALL dbt_filter(t_2c_pot, ri_data%filter_eps)
2305 CALL dbt_destroy(t_2c_tmp)
2308 CALL timestop(handle)
2310 END SUBROUTINE get_ext_2c_int
2320 SUBROUTINE contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
2321 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, rho_ao_t
2325 CHARACTER(len=*),
PARAMETER :: routinen =
'contract_pmat_3c'
2327 INTEGER :: apc_img, batch_size, handle, i_batch, &
2328 i_img, i_spin, j_batch, n_batch_img, &
2329 n_batch_nze, nimg, nimg_nze, nspins
2330 INTEGER(int_8) :: nflop, nze
2331 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges_img, batch_ranges_nze, &
2333 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: ac_pairs
2334 REAL(
dp) :: occ, t1, t2
2335 TYPE(dbt_type) :: t_3c_tmp
2336 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ints_stack, res_stack, rho_stack
2339 CALL timeset(routinen, handle)
2341 CALL get_qs_env(qs_env, dft_control=dft_control)
2344 nimg_nze = ri_data%nimg_nze
2345 nspins = dft_control%nspins
2347 CALL dbt_create(t_3c_apc(1, 1), t_3c_tmp)
2349 batch_size = nimg/ri_data%n_mem
2352 n_batch_img = nimg/batch_size
2353 IF (
modulo(nimg, batch_size) .NE. 0) n_batch_img = n_batch_img + 1
2354 ALLOCATE (batch_ranges_img(n_batch_img + 1))
2355 DO i_batch = 1, n_batch_img
2356 batch_ranges_img(i_batch) = (i_batch - 1)*batch_size + 1
2358 batch_ranges_img(n_batch_img + 1) = nimg + 1
2361 n_batch_nze = nimg_nze/batch_size
2362 IF (
modulo(nimg_nze, batch_size) .NE. 0) n_batch_nze = n_batch_nze + 1
2363 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
2364 DO i_batch = 1, n_batch_nze
2365 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
2367 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
2370 ALLOCATE (rho_stack(2), ints_stack(2), res_stack(2))
2371 CALL get_stack_tensors(res_stack, rho_stack, ints_stack, rho_ao_t(1, 1), &
2372 ri_data%t_3c_int_ctr_1(1, 1), batch_size, ri_data, qs_env)
2374 ALLOCATE (ac_pairs(nimg, 2), int_indices(nimg_nze))
2375 DO i_img = 1, nimg_nze
2376 int_indices(i_img) = i_img
2380 DO j_batch = 1, n_batch_nze
2382 CALL fill_3c_stack(ints_stack(1), ri_data%t_3c_int_ctr_1(1, :), int_indices, 3, ri_data, &
2383 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)])
2384 CALL dbt_copy(ints_stack(1), ints_stack(2), move_data=.true.)
2386 DO i_spin = 1, nspins
2387 DO i_batch = 1, n_batch_img
2389 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2390 CALL get_ac_pairs(ac_pairs, apc_img, ri_data, qs_env)
2391 CALL fill_2c_stack(rho_stack(1), rho_ao_t(i_spin, :), ac_pairs(:, 2), 1, ri_data, &
2392 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)], &
2393 shift=apc_img - batch_ranges_img(i_batch) + 1)
2398 CALL dbt_copy(rho_stack(1), rho_stack(2), move_data=.true.)
2401 CALL dbt_batched_contract_init(rho_stack(2))
2402 CALL dbt_contract(1.0_dp, ints_stack(2), rho_stack(2), &
2403 0.0_dp, res_stack(2), map_1=[1, 2], map_2=[3], &
2404 contract_1=[3], notcontract_1=[1, 2], &
2405 contract_2=[1], notcontract_2=[2], &
2406 filter_eps=ri_data%filter_eps, flop=nflop)
2407 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2408 CALL dbt_batched_contract_finalize(rho_stack(2))
2409 CALL dbt_copy(res_stack(2), res_stack(1), move_data=.true.)
2411 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2413 CALL unstack_t_3c_apc(t_3c_tmp, res_stack(1), apc_img - batch_ranges_img(i_batch) + 1)
2414 CALL dbt_copy(t_3c_tmp, t_3c_apc(i_spin, apc_img), summation=.true., move_data=.true.)
2420 DEALLOCATE (batch_ranges_img)
2421 DEALLOCATE (batch_ranges_nze)
2423 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
2425 CALL dbt_destroy(rho_stack(1))
2426 CALL dbt_destroy(rho_stack(2))
2427 CALL dbt_destroy(ints_stack(1))
2428 CALL dbt_destroy(ints_stack(2))
2429 CALL dbt_destroy(res_stack(1))
2430 CALL dbt_destroy(res_stack(2))
2431 CALL dbt_destroy(t_3c_tmp)
2433 CALL timestop(handle)
2435 END SUBROUTINE contract_pmat_3c
2443 SUBROUTINE precontract_3c_ints(t_3c_int, ri_data, qs_env)
2444 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_int
2448 CHARACTER(len=*),
PARAMETER :: routinen =
'precontract_3c_ints'
2450 INTEGER :: batch_size, handle, i_batch, i_img, &
2451 i_ri, iatom, is, n_batch, natom, &
2452 nblks, nblks_3c(3), nimg
2453 INTEGER(int_8) :: nflop
2454 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges, bsizes_ri_ext, bsizes_ri_ext_split, &
2455 bsizes_stack, dist1, dist2, dist3, dist_stack3, idx_to_at_ao, int_indices
2456 TYPE(dbt_distribution_type) :: t_dist
2457 TYPE(dbt_type) :: t_2c_ri_tmp(2), t_3c_tmp(3)
2459 CALL timeset(routinen, handle)
2464 ALLOCATE (int_indices(nimg))
2466 int_indices(i_img) = i_img
2469 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
2470 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
2472 nblks =
SIZE(ri_data%bsizes_RI_split)
2473 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
2474 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
2475 DO i_ri = 1, ri_data%ncell_RI
2476 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2477 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
2480 bsizes_ri_ext, bsizes_ri_ext, &
2482 DEALLOCATE (dist1, dist2)
2484 bsizes_ri_ext_split, bsizes_ri_ext_split, &
2486 DEALLOCATE (dist1, dist2)
2489 batch_size = nimg/ri_data%n_mem
2490 n_batch = nimg/batch_size
2491 IF (
modulo(nimg, batch_size) .NE. 0) n_batch = n_batch + 1
2492 ALLOCATE (batch_ranges(n_batch + 1))
2493 DO i_batch = 1, n_batch
2494 batch_ranges(i_batch) = (i_batch - 1)*batch_size + 1
2496 batch_ranges(n_batch + 1) = nimg + 1
2498 nblks =
SIZE(ri_data%bsizes_AO_split)
2499 ALLOCATE (bsizes_stack(batch_size*nblks))
2500 DO is = 1, batch_size
2501 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
2504 CALL dbt_get_info(t_3c_int(1, 1), nblks_total=nblks_3c)
2505 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)), dist_stack3(batch_size*nblks_3c(3)))
2506 CALL dbt_get_info(t_3c_int(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, proc_dist_3=dist3)
2507 DO is = 1, batch_size
2508 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
2511 CALL dbt_distribution_new(t_dist, ri_data%pgrid, dist1, dist2, dist_stack3)
2512 CALL dbt_create(t_3c_tmp(1),
"ints_stack", t_dist, [1], [2, 3], bsizes_ri_ext_split, &
2513 ri_data%bsizes_AO_split, bsizes_stack)
2514 CALL dbt_distribution_destroy(t_dist)
2515 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
2517 CALL dbt_create(t_3c_tmp(1), t_3c_tmp(2))
2518 CALL dbt_create(t_3c_int(1, 1), t_3c_tmp(3))
2521 CALL dbt_copy(ri_data%t_2c_inv(1, iatom), t_2c_ri_tmp(1))
2522 CALL apply_bump(t_2c_ri_tmp(1), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
2523 CALL dbt_copy(t_2c_ri_tmp(1), t_2c_ri_tmp(2), move_data=.true.)
2525 CALL dbt_batched_contract_init(t_2c_ri_tmp(2))
2526 DO i_batch = 1, n_batch
2528 CALL fill_3c_stack(t_3c_tmp(1), t_3c_int(1, :), int_indices, 3, ri_data, &
2529 img_bounds=[batch_ranges(i_batch), batch_ranges(i_batch + 1)], &
2530 filter_at=iatom, filter_dim=2, idx_to_at=idx_to_at_ao)
2532 CALL dbt_contract(1.0_dp, t_2c_ri_tmp(2), t_3c_tmp(1), &
2533 0.0_dp, t_3c_tmp(2), map_1=[1], map_2=[2, 3], &
2534 contract_1=[2], notcontract_1=[1], &
2535 contract_2=[1], notcontract_2=[2, 3], &
2536 filter_eps=ri_data%filter_eps, flop=nflop)
2537 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2539 DO i_img = batch_ranges(i_batch), batch_ranges(i_batch + 1) - 1
2540 CALL unstack_t_3c_apc(t_3c_tmp(3), t_3c_tmp(2), i_img - batch_ranges(i_batch) + 1)
2541 CALL dbt_copy(t_3c_tmp(3), ri_data%t_3c_int_ctr_1(1, i_img), summation=.true., &
2542 order=[2, 1, 3], move_data=.true.)
2544 CALL dbt_clear(t_3c_tmp(1))
2546 CALL dbt_batched_contract_finalize(t_2c_ri_tmp(2))
2549 CALL dbt_destroy(t_2c_ri_tmp(1))
2550 CALL dbt_destroy(t_2c_ri_tmp(2))
2551 CALL dbt_destroy(t_3c_tmp(1))
2552 CALL dbt_destroy(t_3c_tmp(2))
2553 CALL dbt_destroy(t_3c_tmp(3))
2556 CALL dbt_destroy(t_3c_int(1, i_img))
2559 CALL timestop(handle)
2561 END SUBROUTINE precontract_3c_ints
2572 SUBROUTINE copy_2c_to_subgroup(t2c_sub, t2c_main, group_size, ngroups, para_env)
2573 TYPE(dbt_type),
INTENT(INOUT) :: t2c_sub, t2c_main
2574 INTEGER,
INTENT(IN) :: group_size, ngroups
2577 INTEGER :: batch_size, i, i_batch, i_msg, iblk, &
2578 igroup, iproc, ir, is, jblk, n_batch, &
2580 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2
2581 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: block_dest, block_source
2582 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: current_dest
2583 INTEGER,
DIMENSION(2) :: ind, nblks
2585 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
2586 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2587 TYPE(dbt_iterator_type) :: iter
2588 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2594 CALL dbt_get_info(t2c_main, nblks_total=nblks)
2597 ALLOCATE (block_source(nblks(1), nblks(2)))
2601 CALL dbt_iterator_start(iter, t2c_main)
2602 DO WHILE (dbt_iterator_blocks_left(iter))
2603 CALL dbt_iterator_next_block(iter, ind)
2604 CALL dbt_get_block(t2c_main, ind, blk, found)
2605 IF (.NOT. found) cycle
2607 block_source(ind(1), ind(2)) = para_env%mepos
2612 CALL dbt_iterator_stop(iter)
2615 CALL para_env%sum(nocc)
2616 CALL para_env%sum(block_source)
2617 block_source = block_source + para_env%num_pe - 1
2618 IF (nocc == 0)
RETURN
2621 igroup = para_env%mepos/group_size
2622 ALLOCATE (block_dest(nblks(1), nblks(2)))
2624 DO jblk = 1, nblks(2)
2625 DO iblk = 1, nblks(1)
2626 IF (block_source(iblk, jblk) == -1) cycle
2628 CALL dbt_get_stored_coordinates(t2c_sub, [iblk, jblk], iproc)
2629 block_dest(iblk, jblk) = igroup*group_size + iproc
2633 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)))
2634 CALL dbt_get_info(t2c_main, blk_size_1=bsizes1, blk_size_2=bsizes2)
2636 ALLOCATE (current_dest(nblks(1), nblks(2), 0:ngroups - 1))
2637 DO igroup = 0, ngroups - 1
2639 current_dest(:, :, igroup) = block_dest(:, :)
2640 CALL para_env%bcast(current_dest(:, :, igroup), source=igroup*group_size)
2644 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2645 n_batch = (nocc*ngroups)/batch_size
2646 IF (
modulo(nocc*ngroups, batch_size) .NE. 0) n_batch = n_batch + 1
2648 DO i_batch = 1, n_batch
2650 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2651 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2655 DO jblk = 1, nblks(2)
2656 DO iblk = 1, nblks(1)
2657 DO igroup = 0, ngroups - 1
2658 IF (block_source(iblk, jblk) == -1) cycle
2661 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2664 tag = i_msg - (i_batch - 1)*batch_size
2667 IF (para_env%mepos == block_source(iblk, jblk))
THEN
2668 CALL dbt_get_block(t2c_main, [iblk, jblk], blk, found)
2672 IF (block_source(iblk, jblk) == current_dest(iblk, jblk, igroup))
THEN
2673 IF (found)
CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2675 IF (para_env%mepos == block_source(iblk, jblk) .AND. found)
THEN
2676 ALLOCATE (send_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2677 send_buff(tag)%array(:, :) = blk(:, :)
2679 CALL para_env%isend(msgin=send_buff(tag)%array, dest=current_dest(iblk, jblk, igroup), &
2680 request=send_req(is), tag=tag)
2683 IF (para_env%mepos == current_dest(iblk, jblk, igroup))
THEN
2684 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2686 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk), &
2687 request=recv_req(ir), tag=tag)
2691 IF (found)
DEALLOCATE (blk)
2699 DO i = 1, batch_size
2700 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
2705 DO jblk = 1, nblks(2)
2706 DO iblk = 1, nblks(1)
2707 DO igroup = 0, ngroups - 1
2708 IF (block_source(iblk, jblk) == -1) cycle
2711 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2714 tag = i_msg - (i_batch - 1)*batch_size
2716 IF (para_env%mepos == current_dest(iblk, jblk, igroup) .AND. &
2717 block_source(iblk, jblk) .NE. current_dest(iblk, jblk, igroup))
THEN
2719 ALLOCATE (blk(bsizes1(iblk), bsizes2(jblk)))
2720 blk(:, :) = recv_buff(tag)%array(:, :)
2721 CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2729 DO i = 1, batch_size
2730 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
2732 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
2734 CALL dbt_finalize(t2c_sub)
2736 END SUBROUTINE copy_2c_to_subgroup
2750 SUBROUTINE copy_3c_to_subgroup(t3c_sub, t3c_main, group_size, ngroups, para_env, iatom_to_subgroup, &
2752 TYPE(dbt_type),
INTENT(INOUT) :: t3c_sub, t3c_main
2753 INTEGER,
INTENT(IN) :: group_size, ngroups
2756 INTENT(INOUT),
OPTIONAL :: iatom_to_subgroup
2757 INTEGER,
INTENT(IN),
OPTIONAL :: dim_at
2758 INTEGER,
DIMENSION(:),
OPTIONAL :: idx_to_at
2760 INTEGER :: batch_size, i, i_batch, i_msg, iatom, &
2761 iblk, igroup, iproc, ir, is, jblk, &
2762 kblk, n_batch, nocc, tag
2763 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2, bsizes3
2764 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: block_dest, block_source
2765 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :) :: current_dest
2766 INTEGER,
DIMENSION(3) :: ind, nblks
2767 LOGICAL :: filter_at, found
2768 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
2769 TYPE(
cp_3d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2770 TYPE(dbt_iterator_type) :: iter
2771 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2777 CALL dbt_get_info(t3c_main, nblks_total=nblks)
2781 IF (
PRESENT(iatom_to_subgroup) .AND.
PRESENT(dim_at) .AND.
PRESENT(idx_to_at))
THEN
2783 cpassert(nblks(dim_at) ==
SIZE(idx_to_at))
2787 ALLOCATE (block_source(nblks(1), nblks(2), nblks(3)))
2791 CALL dbt_iterator_start(iter, t3c_main)
2792 DO WHILE (dbt_iterator_blocks_left(iter))
2793 CALL dbt_iterator_next_block(iter, ind)
2794 CALL dbt_get_block(t3c_main, ind, blk, found)
2795 IF (.NOT. found) cycle
2797 block_source(ind(1), ind(2), ind(3)) = para_env%mepos
2802 CALL dbt_iterator_stop(iter)
2805 CALL para_env%sum(nocc)
2806 CALL para_env%sum(block_source)
2807 block_source = block_source + para_env%num_pe - 1
2808 IF (nocc == 0)
RETURN
2811 igroup = para_env%mepos/group_size
2812 ALLOCATE (block_dest(nblks(1), nblks(2), nblks(3)))
2814 DO kblk = 1, nblks(3)
2815 DO jblk = 1, nblks(2)
2816 DO iblk = 1, nblks(1)
2817 IF (block_source(iblk, jblk, kblk) == -1) cycle
2819 CALL dbt_get_stored_coordinates(t3c_sub, [iblk, jblk, kblk], iproc)
2820 block_dest(iblk, jblk, kblk) = igroup*group_size + iproc
2825 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)), bsizes3(nblks(3)))
2826 CALL dbt_get_info(t3c_main, blk_size_1=bsizes1, blk_size_2=bsizes2, blk_size_3=bsizes3)
2828 ALLOCATE (current_dest(nblks(1), nblks(2), nblks(3), 0:ngroups - 1))
2829 DO igroup = 0, ngroups - 1
2831 current_dest(:, :, :, igroup) = block_dest(:, :, :)
2832 CALL para_env%bcast(current_dest(:, :, :, igroup), source=igroup*group_size)
2836 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2837 n_batch = (nocc*ngroups)/batch_size
2838 IF (
modulo(nocc*ngroups, batch_size) .NE. 0) n_batch = n_batch + 1
2840 DO i_batch = 1, n_batch
2842 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2843 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2847 DO kblk = 1, nblks(3)
2848 DO jblk = 1, nblks(2)
2849 DO iblk = 1, nblks(1)
2850 DO igroup = 0, ngroups - 1
2851 IF (block_source(iblk, jblk, kblk) == -1) cycle
2854 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2857 tag = i_msg - (i_batch - 1)*batch_size
2860 ind(:) = [iblk, jblk, kblk]
2861 iatom = idx_to_at(ind(dim_at))
2862 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
2866 IF (para_env%mepos == block_source(iblk, jblk, kblk))
THEN
2867 CALL dbt_get_block(t3c_main, [iblk, jblk, kblk], blk, found)
2871 IF (block_source(iblk, jblk, kblk) == current_dest(iblk, jblk, kblk, igroup))
THEN
2872 IF (found)
CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(blk), blk)
2874 IF (para_env%mepos == block_source(iblk, jblk, kblk) .AND. found)
THEN
2875 ALLOCATE (send_buff(tag)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2876 send_buff(tag)%array(:, :, :) = blk(:, :, :)
2878 CALL para_env%isend(msgin=send_buff(tag)%array, &
2879 dest=current_dest(iblk, jblk, kblk, igroup), &
2880 request=send_req(is), tag=tag)
2883 IF (para_env%mepos == current_dest(iblk, jblk, kblk, igroup))
THEN
2884 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2886 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk, kblk), &
2887 request=recv_req(ir), tag=tag)
2891 IF (found)
DEALLOCATE (blk)
2900 DO i = 1, batch_size
2901 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
2906 DO kblk = 1, nblks(3)
2907 DO jblk = 1, nblks(2)
2908 DO iblk = 1, nblks(1)
2909 DO igroup = 0, ngroups - 1
2910 IF (block_source(iblk, jblk, kblk) == -1) cycle
2913 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2916 tag = i_msg - (i_batch - 1)*batch_size
2919 ind(:) = [iblk, jblk, kblk]
2920 iatom = idx_to_at(ind(dim_at))
2921 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
2924 IF (para_env%mepos == current_dest(iblk, jblk, kblk, igroup) .AND. &
2925 block_source(iblk, jblk, kblk) .NE. current_dest(iblk, jblk, kblk, igroup))
THEN
2927 ALLOCATE (blk(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2928 blk(:, :, :) = recv_buff(tag)%array(:, :, :)
2929 CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(blk), blk)
2938 DO i = 1, batch_size
2939 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
2941 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
2943 CALL dbt_finalize(t3c_sub)
2945 END SUBROUTINE copy_3c_to_subgroup
2957 SUBROUTINE gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
2958 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t, ks_t_sub
2959 INTEGER,
INTENT(IN) :: group_size
2960 INTEGER,
DIMENSION(:, :, :),
INTENT(IN) :: sparsity_pattern
2964 CHARACTER(len=*),
PARAMETER :: routinen =
'gather_ks_matrix'
2966 INTEGER :: b_img, dest, handle, i, i_spin, iatom, &
2967 igroup, ir, is, jatom, n_mess, natom, &
2968 nimg, nspins, source, tag
2970 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
2971 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2972 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2974 CALL timeset(routinen, handle)
2976 nimg =
SIZE(sparsity_pattern, 3)
2977 natom =
SIZE(sparsity_pattern, 2)
2978 nspins =
SIZE(ks_t, 1)
2982 DO i_spin = 1, nspins
2985 IF (sparsity_pattern(iatom, jatom, b_img) > -1) n_mess = n_mess + 1
2990 ALLOCATE (send_buff(n_mess), recv_buff(n_mess))
2991 ALLOCATE (send_req(n_mess), recv_req(n_mess))
2997 DO i_spin = 1, nspins
3000 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3005 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3006 CALL dbt_get_stored_coordinates(ks_t_sub(i_spin, b_img), [iatom, jatom], source)
3007 igroup = sparsity_pattern(iatom, jatom, b_img)
3008 source = source + igroup*group_size
3009 IF (para_env%mepos == source)
THEN
3010 CALL dbt_get_block(ks_t_sub(i_spin, b_img), [iatom, jatom], blk, found)
3011 IF (source == dest)
THEN
3012 IF (found)
CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3014 ALLOCATE (send_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3015 send_buff(n_mess)%array(:, :) = 0.0_dp
3017 send_buff(n_mess)%array(:, :) = blk(:, :)
3020 CALL para_env%isend(msgin=send_buff(n_mess)%array, dest=dest, &
3021 request=send_req(is), tag=tag)
3027 IF (para_env%mepos == dest .AND. source .NE. dest)
THEN
3028 ALLOCATE (recv_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3030 CALL para_env%irecv(msgout=recv_buff(n_mess)%array, source=source, &
3031 request=recv_req(ir), tag=tag)
3042 DO i_spin = 1, nspins
3045 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3048 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3049 IF (para_env%mepos == dest)
THEN
3050 IF (.NOT.
ASSOCIATED(recv_buff(n_mess)%array)) cycle
3051 ALLOCATE (blk(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3052 blk(:, :) = recv_buff(n_mess)%array(:, :)
3053 CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3062 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
3063 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
3065 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
3068 CALL timestop(handle)
3070 END SUBROUTINE gather_ks_matrix
3085 SUBROUTINE get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
3086 group_size, ngroups, para_env, para_env_sub, ri_data)
3087 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3088 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work, t_2c_ao_tmp, ks_t_split
3089 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t_sub
3090 INTEGER,
INTENT(IN) :: group_size, ngroups
3094 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_tensors'
3096 INTEGER :: handle, i, i_img, i_ri, i_spin, iproc, &
3097 j, natom, nblks, nimg, nspins
3098 INTEGER(int_8) :: nze
3099 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3101 INTEGER,
DIMENSION(2) :: pdims_2d
3102 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3103 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3106 TYPE(dbt_pgrid_type) :: pgrid_2d
3107 TYPE(dbt_type) :: work, work_sub
3109 CALL timeset(routinen, handle)
3113 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3115 natom =
SIZE(ri_data%bsizes_RI)
3116 nblks =
SIZE(ri_data%bsizes_RI_split)
3117 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3118 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3119 DO i_ri = 1, ri_data%ncell_RI
3120 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3121 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3126 bsizes_ri_ext, bsizes_ri_ext, &
3128 DEALLOCATE (dist1, dist2)
3131 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3133 DEALLOCATE (dist1, dist2)
3137 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3139 DEALLOCATE (dist1, dist2)
3140 CALL dbt_create(ks_t_split(1), ks_t_split(2))
3143 ri_data%bsizes_AO, ri_data%bsizes_AO, &
3145 DEALLOCATE (dist1, dist2)
3147 nspins =
SIZE(ks_t_sub, 1)
3148 nimg =
SIZE(ks_t_sub, 2)
3150 DO i_spin = 1, nspins
3151 CALL dbt_create(t_2c_ao_tmp(1), ks_t_sub(i_spin, i_img))
3158 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3160 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3162 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3164 DO i = 0, pdims_2d(1) - 1
3165 DO j = 0, pdims_2d(2) - 1
3166 dbcsr_pgrid(i, j) = iproc
3172 ALLOCATE (col_dist(natom), row_dist(natom))
3173 row_dist(:) = dist1(:)
3174 col_dist(:) = dist2(:)
3176 ALLOCATE (ri_blk_size(natom))
3177 ri_blk_size(:) = ri_data%bsizes_RI(:)
3180 row_dist=row_dist, col_dist=col_dist)
3181 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3182 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3185 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3186 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3190 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3191 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3192 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3193 CALL dbt_clear(work_sub)
3196 CALL dbt_destroy(work)
3197 CALL dbt_destroy(work_sub)
3198 CALL dbt_pgrid_destroy(pgrid_2d)
3200 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3201 CALL timestop(handle)
3203 END SUBROUTINE get_subgroup_2c_tensors
3218 SUBROUTINE get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
3219 group_size, ngroups, para_env, para_env_sub, ri_data)
3220 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_int, t_3c_work_2, t_3c_work_3
3221 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, t_3c_apc_sub
3222 INTEGER,
INTENT(IN) :: group_size, ngroups
3226 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_tensors'
3228 INTEGER :: batch_size, bfac, bo(2), handle, &
3229 handle2, i_blk, i_img, i_ri, i_spin, &
3230 ib, natom, nblks_ao, nblks_ri, nimg, &
3232 INTEGER(int_8) :: nze
3233 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3234 bsizes_stack, bsizes_tmp, dist1, &
3235 dist2, dist3, dist_stack, idx_to_at
3236 INTEGER,
DIMENSION(3) :: pdims
3238 TYPE(dbt_distribution_type) :: t_dist
3239 TYPE(dbt_pgrid_type) :: pgrid
3240 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3242 CALL timeset(routinen, handle)
3244 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3245 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3246 DO i_ri = 1, ri_data%ncell_RI
3247 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3253 natom =
SIZE(ri_data%bsizes_RI)
3254 nblks_ri = max(1, natom/bfac)
3255 ALLOCATE (bsizes_tmp(nblks_ri))
3256 DO i_blk = 1, nblks_ri
3257 bo =
get_limit(natom, nblks_ri, i_blk - 1)
3258 bsizes_tmp(i_blk) = sum(ri_data%bsizes_RI(bo(1):bo(2)))
3260 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3261 DO i_ri = 1, ri_data%ncell_RI
3262 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = bsizes_tmp(:)
3265 batch_size = ri_data%kp_stack_size
3266 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3267 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3268 DO ib = 1, batch_size
3269 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3273 natom =
SIZE(ri_data%bsizes_RI)
3275 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3276 tensor_dims=[
SIZE(bsizes_ri_ext_split), 1, batch_size*
SIZE(ri_data%bsizes_AO_split)])
3280 pgrid, bsizes_ri_ext_split, ri_data%bsizes_AO_split, &
3281 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(RI | AO AO)")
3282 nimg =
SIZE(t_3c_int)
3284 CALL dbt_create(t_3c_int(1), t_3c_int(i_img))
3288 ALLOCATE (dist_stack(batch_size*nblks_ao))
3289 DO ib = 1, batch_size
3290 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3293 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3294 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3295 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3296 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3297 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3298 CALL dbt_distribution_destroy(t_dist)
3299 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3303 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3304 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3305 DEALLOCATE (dist1, dist2, dist3)
3308 ri_data%pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3309 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3310 DEALLOCATE (dist1, dist2, dist3)
3313 CALL timeset(routinen//
"_ints", handle2)
3314 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
3316 CALL dbt_copy(ri_data%kp_t_3c_int(i_img), t_3c_int(i_img), move_data=.true.)
3319 ALLOCATE (ri_data%kp_t_3c_int(nimg))
3321 CALL dbt_create(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img))
3324 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, i_img), work_atom_block, order=[2, 1, 3])
3325 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, group_size, ngroups, para_env)
3326 CALL dbt_copy(work_atom_block_sub, t_3c_int(i_img), move_data=.true.)
3329 CALL timestop(handle2)
3330 CALL dbt_pgrid_destroy(pgrid)
3331 CALL dbt_destroy(work_atom_block)
3332 CALL dbt_destroy(work_atom_block_sub)
3336 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3337 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3341 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3342 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3343 DEALLOCATE (dist1, dist2, dist3)
3346 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3347 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3348 DEALLOCATE (dist1, dist2, dist3)
3352 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3353 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3356 ALLOCATE (dist_stack(batch_size*nblks_ao))
3357 DO ib = 1, batch_size
3358 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3361 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3362 CALL dbt_create(t_3c_work_2(1),
"work_2_stack", t_dist, [1], [2, 3], &
3363 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3364 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3365 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3366 CALL dbt_distribution_destroy(t_dist)
3367 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3370 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3372 nspins =
SIZE(t_3c_apc, 1)
3373 CALL timeset(routinen//
"_apc", handle2)
3375 DO i_spin = 1, nspins
3376 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3379 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3380 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, group_size, &
3381 ngroups, para_env, ri_data%iatom_to_subgroup, 1, idx_to_at)
3382 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3384 DO i_spin = 1, nspins
3385 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3388 CALL timestop(handle2)
3389 CALL dbt_pgrid_destroy(pgrid)
3390 CALL dbt_destroy(tmp)
3391 CALL dbt_destroy(work_atom_block)
3392 CALL dbt_destroy(work_atom_block_sub)
3394 CALL timestop(handle)
3396 END SUBROUTINE get_subgroup_3c_tensors
3418 SUBROUTINE get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
3419 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
3420 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
3421 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_inv, t_2c_bint, t_2c_metric
3422 TYPE(
dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3423 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work
3424 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
3426 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot, mat_der_pot_sub
3427 INTEGER,
INTENT(IN) :: group_size, ngroups
3431 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_derivs'
3433 INTEGER :: handle, i, i_img, i_ri, i_spin, i_xyz, &
3434 iatom, iproc, j, natom, nblks, nimg, &
3436 INTEGER(int_8) :: nze
3437 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3439 INTEGER,
DIMENSION(2) :: pdims_2d
3440 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3441 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3444 TYPE(dbt_pgrid_type) :: pgrid_2d
3445 TYPE(dbt_type) :: work, work_sub
3447 CALL timeset(routinen, handle)
3452 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3454 natom =
SIZE(ri_data%bsizes_RI)
3455 nblks =
SIZE(ri_data%bsizes_RI_split)
3456 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3457 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3458 DO i_ri = 1, ri_data%ncell_RI
3459 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3460 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3465 bsizes_ri_ext, bsizes_ri_ext, &
3467 DEALLOCATE (dist1, dist2)
3469 CALL dbt_create(t_2c_inv(1), t_2c_bint(1))
3470 CALL dbt_create(t_2c_inv(1), t_2c_metric(1))
3472 CALL dbt_create(t_2c_inv(1), t_2c_inv(iatom))
3473 CALL dbt_create(t_2c_inv(1), t_2c_bint(iatom))
3474 CALL dbt_create(t_2c_inv(1), t_2c_metric(iatom))
3476 CALL dbt_create(t_2c_inv(1), t_2c_work(1))
3477 CALL dbt_create(t_2c_inv(1), t_2c_work(2))
3478 CALL dbt_create(t_2c_inv(1), t_2c_work(3))
3479 CALL dbt_create(t_2c_inv(1), t_2c_work(4))
3482 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3484 DEALLOCATE (dist1, dist2)
3488 CALL copy_2c_to_subgroup(t_2c_inv(iatom), ri_data%t_2c_inv(1, iatom), group_size, ngroups, para_env)
3489 CALL copy_2c_to_subgroup(t_2c_bint(iatom), ri_data%t_2c_int(1, iatom), group_size, ngroups, para_env)
3490 CALL copy_2c_to_subgroup(t_2c_metric(iatom), ri_data%t_2c_pot(1, iatom), group_size, ngroups, para_env)
3496 CALL dbt_create(t_2c_inv(1), t_2c_der_metric_sub(iatom, i_xyz))
3497 CALL copy_2c_to_subgroup(t_2c_der_metric_sub(iatom, i_xyz), t_2c_der_metric(iatom, i_xyz), &
3498 group_size, ngroups, para_env)
3499 CALL dbt_destroy(t_2c_der_metric(iatom, i_xyz))
3505 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3507 DEALLOCATE (dist1, dist2)
3508 nspins =
SIZE(rho_ao_t, 1)
3509 nimg =
SIZE(rho_ao_t, 2)
3512 DO i_spin = 1, nspins
3513 IF (.NOT. (i_img == 1 .AND. i_spin == 1)) &
3514 CALL dbt_create(rho_ao_t_sub(1, 1), rho_ao_t_sub(i_spin, i_img))
3515 CALL copy_2c_to_subgroup(rho_ao_t_sub(i_spin, i_img), rho_ao_t(i_spin, i_img), &
3516 group_size, ngroups, para_env)
3517 CALL dbt_destroy(rho_ao_t(i_spin, i_img))
3523 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3525 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3527 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3529 DO i = 0, pdims_2d(1) - 1
3530 DO j = 0, pdims_2d(2) - 1
3531 dbcsr_pgrid(i, j) = iproc
3537 ALLOCATE (col_dist(natom), row_dist(natom))
3538 row_dist(:) = dist1(:)
3539 col_dist(:) = dist2(:)
3541 ALLOCATE (ri_blk_size(natom))
3542 ri_blk_size(:) = ri_data%bsizes_RI(:)
3545 row_dist=row_dist, col_dist=col_dist)
3546 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3547 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3551 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3552 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3556 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3557 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3558 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3559 CALL dbt_clear(work_sub)
3565 CALL dbcsr_create(mat_der_pot_sub(i_img, i_xyz), template=mat_2c_pot(1))
3566 CALL dbt_copy_matrix_to_tensor(mat_der_pot(i_img, i_xyz), work)
3571 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3572 CALL dbt_copy_tensor_to_matrix(work_sub, mat_der_pot_sub(i_img, i_xyz))
3573 CALL dbcsr_filter(mat_der_pot_sub(i_img, i_xyz), ri_data%filter_eps)
3574 CALL dbt_clear(work_sub)
3578 CALL dbt_destroy(work)
3579 CALL dbt_destroy(work_sub)
3580 CALL dbt_pgrid_destroy(pgrid_2d)
3582 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3584 CALL timestop(handle)
3586 END SUBROUTINE get_subgroup_2c_derivs
3606 SUBROUTINE get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_AO, t_3c_der_AO_sub, &
3607 t_3c_der_RI, t_3c_der_RI_sub, t_3c_apc, t_3c_apc_sub, &
3608 t_3c_der_stack, group_size, ngroups, para_env, para_env_sub, &
3610 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_work_2, t_3c_work_3
3611 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ao, t_3c_der_ao_sub, &
3612 t_3c_der_ri, t_3c_der_ri_sub, &
3613 t_3c_apc, t_3c_apc_sub
3614 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_der_stack
3615 INTEGER,
INTENT(IN) :: group_size, ngroups
3619 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_derivs'
3621 INTEGER :: batch_size, handle, i_img, i_ri, i_spin, &
3622 i_xyz, ib, nblks_ao, nblks_ri, nimg, &
3624 INTEGER(int_8) :: nze
3625 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3626 bsizes_stack, dist1, dist2, dist3, &
3627 dist_stack, idx_to_at
3629 TYPE(dbt_distribution_type) :: t_dist
3630 TYPE(dbt_pgrid_type) :: pgrid
3631 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3633 CALL timeset(routinen, handle)
3636 nblks_ri =
SIZE(ri_data%bsizes_RI)
3637 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3638 DO i_ri = 1, ri_data%ncell_RI
3639 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
3642 CALL dbt_get_info(ri_data%kp_t_3c_int(1), pdims=pdims)
3643 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3646 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3647 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3648 DEALLOCATE (dist1, dist2, dist3)
3651 ri_data%pgrid_2, bsizes_ri_ext, ri_data%bsizes_AO, &
3652 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3653 DEALLOCATE (dist1, dist2, dist3)
3654 CALL dbt_pgrid_destroy(pgrid)
3660 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ao_sub(i_img, i_xyz))
3664 CALL dbt_copy(t_3c_der_ao(i_img, i_xyz), work_atom_block, move_data=.true.)
3665 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3666 group_size, ngroups, para_env)
3667 CALL dbt_copy(work_atom_block_sub, t_3c_der_ao_sub(i_img, i_xyz), move_data=.true.)
3671 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ri_sub(i_img, i_xyz))
3675 CALL dbt_copy(t_3c_der_ri(i_img, i_xyz), work_atom_block, move_data=.true.)
3676 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3677 group_size, ngroups, para_env)
3678 CALL dbt_copy(work_atom_block_sub, t_3c_der_ri_sub(i_img, i_xyz), move_data=.true.)
3682 CALL dbt_destroy(t_3c_der_ri(i_img, i_xyz))
3683 CALL dbt_destroy(t_3c_der_ao(i_img, i_xyz))
3686 CALL dbt_destroy(work_atom_block_sub)
3687 CALL dbt_destroy(work_atom_block)
3690 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3691 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3692 DO i_ri = 1, ri_data%ncell_RI
3693 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3697 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3698 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3701 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3702 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3703 DEALLOCATE (dist1, dist2, dist3)
3706 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3707 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3708 DEALLOCATE (dist1, dist2, dist3)
3711 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3712 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3713 DEALLOCATE (dist1, dist2, dist3)
3715 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3717 nspins =
SIZE(t_3c_apc, 1)
3719 DO i_spin = 1, nspins
3720 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3723 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3724 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, group_size, &
3725 ngroups, para_env, ri_data%iatom_to_subgroup, 1, idx_to_at)
3726 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3728 DO i_spin = 1, nspins
3729 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3732 CALL dbt_destroy(tmp)
3733 CALL dbt_destroy(work_atom_block)
3734 CALL dbt_destroy(work_atom_block_sub)
3735 CALL dbt_pgrid_destroy(pgrid)
3738 batch_size = ri_data%kp_stack_size
3739 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3740 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3741 DO ib = 1, batch_size
3742 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3745 ALLOCATE (dist1(ri_data%ncell_RI*nblks_ri), dist2(nblks_ao), dist3(nblks_ao))
3746 CALL dbt_get_info(ri_data%kp_t_3c_int(1), proc_dist_1=dist1, proc_dist_2=dist2, &
3747 proc_dist_3=dist3, pdims=pdims)
3749 ALLOCATE (dist_stack(batch_size*nblks_ao))
3750 DO ib = 1, batch_size
3751 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3754 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3755 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3756 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3757 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3758 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3759 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3760 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(4))
3761 CALL dbt_distribution_destroy(t_dist)
3762 CALL dbt_pgrid_destroy(pgrid)
3763 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3766 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(1))
3767 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(2))
3768 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(3))
3769 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(4))
3770 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(5))
3771 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(6))
3774 ALLOCATE (dist1(nblks_ao), dist2(ri_data%ncell_RI*nblks_ri), dist3(nblks_ao))
3775 CALL dbt_get_info(t_3c_apc_sub(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, &
3776 proc_dist_3=dist3, pdims=pdims)
3778 ALLOCATE (dist_stack(batch_size*nblks_ao))
3779 DO ib = 1, batch_size
3780 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3783 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3784 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3785 CALL dbt_create(t_3c_work_2(1),
"work_3_stack", t_dist, [1], [2, 3], &
3786 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3787 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3788 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3789 CALL dbt_distribution_destroy(t_dist)
3790 CALL dbt_pgrid_destroy(pgrid)
3791 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3793 CALL timestop(handle)
3795 END SUBROUTINE get_subgroup_3c_derivs
3803 SUBROUTINE reorder_3c_ints(t_3c_ints, ri_data)
3804 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_ints
3807 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_ints'
3809 INTEGER :: handle, i_img,
idx, idx_empty, idx_full, &
3811 INTEGER(int_8) :: nze
3813 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
3815 CALL timeset(routinen, handle)
3818 ALLOCATE (t_3c_tmp(nimg))
3820 CALL dbt_create(t_3c_ints(i_img), t_3c_tmp(i_img))
3821 CALL dbt_copy(t_3c_ints(i_img), t_3c_tmp(i_img), move_data=.true.)
3826 ALLOCATE (ri_data%idx_to_img(nimg))
3828 idx_empty = nimg + 1
3833 idx_empty = idx_empty - 1
3834 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_empty), move_data=.true.)
3835 ri_data%idx_to_img(idx_empty) = i_img
3837 idx_full = idx_full + 1
3838 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_full), move_data=.true.)
3839 ri_data%idx_to_img(idx_full) = i_img
3841 CALL dbt_destroy(t_3c_tmp(i_img))
3845 ri_data%nimg_nze = idx_full
3847 ALLOCATE (ri_data%img_to_idx(nimg))
3849 ri_data%img_to_idx(ri_data%idx_to_img(
idx)) =
idx
3852 CALL timestop(handle)
3854 END SUBROUTINE reorder_3c_ints
3862 SUBROUTINE reorder_3c_derivs(t_3c_derivs, ri_data)
3863 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_derivs
3866 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_derivs'
3868 INTEGER :: handle, i_img, i_xyz,
idx, nimg
3869 INTEGER(int_8) :: nze
3871 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
3873 CALL timeset(routinen, handle)
3876 ALLOCATE (t_3c_tmp(nimg))
3878 CALL dbt_create(t_3c_derivs(1, 1), t_3c_tmp(i_img))
3883 CALL dbt_copy(t_3c_derivs(i_img, i_xyz), t_3c_tmp(i_img), move_data=.true.)
3886 idx = ri_data%img_to_idx(i_img)
3887 CALL dbt_copy(t_3c_tmp(i_img), t_3c_derivs(
idx, i_xyz), move_data=.true.)
3889 IF (nze > 0) ri_data%nimg_nze = max(
idx, ri_data%nimg_nze)
3894 CALL dbt_destroy(t_3c_tmp(i_img))
3897 CALL timestop(handle)
3899 END SUBROUTINE reorder_3c_derivs
3907 SUBROUTINE get_sparsity_pattern(pattern, ri_data, qs_env)
3908 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: pattern
3912 INTEGER :: iatom, j_img, jatom, mj_img, natom, nimg
3913 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bins
3914 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: tmp_pattern
3915 INTEGER,
DIMENSION(3) :: cell_j
3916 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
3917 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
3922 DIMENSION(:),
POINTER :: nl_iterator
3926 NULLIFY (nl_2c, nl_iterator, kpoints, cell_to_index, dft_control, index_to_cell, para_env)
3928 CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, para_env=para_env, natom=natom)
3929 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell, sab_nl=nl_2c)
3932 pattern(:, :, :) = 0
3939 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3940 IF (j_img > nimg .OR. j_img < 1) cycle
3942 mj_img = get_opp_index(j_img, qs_env)
3943 IF (mj_img > nimg .OR. mj_img < 1) cycle
3945 IF (ri_data%present_images(j_img) == 0) cycle
3947 pattern(iatom, jatom, j_img) = 1
3958 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3959 IF (j_img > nimg .OR. j_img < 1) cycle
3961 mj_img = get_opp_index(j_img, qs_env)
3962 IF (mj_img .LE. nimg .AND. mj_img > 0) cycle
3964 IF (ri_data%present_images(j_img) == 0) cycle
3966 pattern(iatom, jatom, j_img) = 1
3970 CALL para_env%sum(pattern)
3975 IF (pattern(iatom, iatom, j_img) .NE. 0)
THEN
3976 mj_img = get_opp_index(j_img, qs_env)
3977 IF (mj_img > nimg .OR. mj_img < 1) cycle
3978 pattern(iatom, iatom, mj_img) = 0
3985 ALLOCATE (bins(natom))
3988 ALLOCATE (tmp_pattern(natom, natom, nimg))
3989 tmp_pattern(:, :, :) = 0
3993 IF (pattern(iatom, jatom, j_img) == 0) cycle
3994 mj_img = get_opp_index(j_img, qs_env)
3997 IF (mj_img > nimg .OR. mj_img < 1)
THEN
3999 bins(iatom) = bins(iatom) + 1
4000 tmp_pattern(iatom, jatom, j_img) = 1
4003 IF (bins(iatom) > bins(jatom))
THEN
4004 bins(jatom) = bins(jatom) + 1
4005 tmp_pattern(jatom, iatom, mj_img) = 1
4007 bins(iatom) = bins(iatom) + 1
4008 tmp_pattern(iatom, jatom, j_img) = 1
4016 pattern(:, :, :) = tmp_pattern(:, :, :) - 1
4018 END SUBROUTINE get_sparsity_pattern
4028 SUBROUTINE get_sub_dist(sparsity_pattern, ngroups, ri_data)
4029 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: sparsity_pattern
4030 INTEGER,
INTENT(IN) :: ngroups
4033 INTEGER :: b_img, ctr, iat, iatom, igroup, jatom, &
4035 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: max_at_per_group
4037 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4039 natom =
SIZE(sparsity_pattern, 2)
4040 nimg =
SIZE(sparsity_pattern, 3)
4045 IF (.NOT.
ALLOCATED(ri_data%iatom_to_subgroup))
THEN
4046 ALLOCATE (ri_data%iatom_to_subgroup(natom), max_at_per_group(ngroups))
4048 NULLIFY (ri_data%iatom_to_subgroup(iatom)%array)
4049 ALLOCATE (ri_data%iatom_to_subgroup(iatom)%array(ngroups))
4050 ri_data%iatom_to_subgroup(iatom)%array(:) = .false.
4054 IF (ub*ngroups < natom) ub = ub + 1
4055 max_at_per_group(:) = max(1, ub)
4060 DO WHILE (
modulo(sum(max_at_per_group), natom) .NE. 0)
4061 igroup =
modulo(ctr, ngroups) + 1
4062 max_at_per_group(igroup) = max_at_per_group(igroup) + 1
4067 DO igroup = 1, ngroups
4068 DO iat = 1, max_at_per_group(igroup)
4069 iatom =
modulo(ctr, natom) + 1
4070 ri_data%iatom_to_subgroup(iatom)%array(igroup) = .true.
4076 ALLOCATE (bins(ngroups))
4081 IF (sparsity_pattern(iatom, jatom, b_img) == -1) cycle
4082 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4085 IF (any(ri_data%kp_cost > epsilon(0.0_dp)))
THEN
4086 cost = ri_data%kp_cost(iatom, jatom, b_img)
4088 cost = real(ri_data%bsizes_AO(iatom)*ri_data%bsizes_AO(jatom),
dp)
4090 bins(igroup + 1) = bins(igroup + 1) + cost
4091 sparsity_pattern(iatom, jatom, b_img) = igroup
4096 END SUBROUTINE get_sub_dist
4107 SUBROUTINE update_pattern_to_forces(force_pattern, scf_pattern, ngroups, ri_data, qs_env)
4108 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: force_pattern, scf_pattern
4109 INTEGER,
INTENT(IN) :: ngroups
4113 INTEGER :: b_img, iatom, igroup, jatom, mb_img, &
4115 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4117 natom =
SIZE(scf_pattern, 2)
4118 nimg =
SIZE(scf_pattern, 3)
4120 ALLOCATE (bins(ngroups))
4124 mb_img = get_opp_index(b_img, qs_env)
4128 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4131 IF (scf_pattern(iatom, jatom, b_img) > -1) cycle
4134 IF (mb_img > 0 .AND. mb_img .LE. nimg)
THEN
4135 IF (scf_pattern(jatom, iatom, mb_img) == -1) cycle
4136 bins(igroup + 1) = bins(igroup + 1) + ri_data%kp_cost(jatom, iatom, mb_img)
4137 force_pattern(iatom, jatom, b_img) = igroup
4143 END SUBROUTINE update_pattern_to_forces
4151 SUBROUTINE get_kp_and_ri_images(ri_data, qs_env)
4155 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_kp_and_ri_images'
4157 INTEGER :: cell_j(3), cell_k(3), handle, i_img, iatom, ikind, j_img, jatom, jcell, katom, &
4158 kcell, kp_index_lbounds(3), kp_index_ubounds(3), natom, ngroups, nimg, nkind, pcoord(3), &
4160 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_ri, &
4161 nri_per_atom, present_img, ri_cells
4162 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
4163 REAL(
dp) :: bump_fact, dij, dik, image_range, &
4164 ri_range, rij(3), rik(3)
4165 TYPE(dbt_type) :: t_dummy
4170 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4177 DIMENSION(:),
POINTER :: nl_iterator
4181 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4184 NULLIFY (qs_kind_set, dist_2d, nl_2c, nl_iterator, dft_control, &
4185 particle_set, kpoints, para_env, cell_to_index, hfx_section)
4187 CALL timeset(routinen, handle)
4189 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, &
4190 dft_control=dft_control, particle_set=particle_set, kpoints=kpoints, &
4191 para_env=para_env, natom=natom)
4192 nimg = dft_control%nimages
4194 kp_index_lbounds = lbound(cell_to_index)
4195 kp_index_ubounds = ubound(cell_to_index)
4200 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4207 CALL erfc_cutoff(ri_data%eps_schwarz, ri_data%hfx_pot%omega, ri_data%hfx_pot%cutoff_radius)
4211 ri_data%kp_RI_range = 0.0_dp
4212 ri_data%kp_image_range = 0.0_dp
4217 ri_data%kp_RI_range = max(ri_range, ri_data%kp_RI_range)
4221 CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, kind_radius=image_range)
4224 ri_data%kp_image_range = max(image_range, ri_data%kp_image_range)
4228 ri_data%kp_bump_rad = bump_fact*ri_data%kp_RI_range
4234 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
4236 ALLOCATE (present_img(nimg))
4245 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4246 IF (j_img > nimg .OR. j_img < 1) cycle
4248 IF (dij > ri_data%kp_image_range) cycle
4250 ri_data%nimg = max(j_img, ri_data%nimg)
4251 present_img(j_img) = 1
4256 CALL para_env%max(ri_data%nimg)
4257 IF (ri_data%nimg > nimg) &
4258 cpabort(
"Make sure the smallest exponent of the RI-HFX basis is larger than that of the ORB basis.")
4261 CALL para_env%sum(present_img)
4262 ALLOCATE (ri_data%present_images(ri_data%nimg))
4263 ri_data%present_images = 0
4264 DO i_img = 1, ri_data%nimg
4265 IF (present_img(i_img) > 0) ri_data%present_images(i_img) = 1
4269 ri_data%pgrid, ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4270 map1=[1, 2], map2=[3], name=
"(AO AO | RI)")
4272 CALL dbt_mp_environ_pgrid(ri_data%pgrid, pdims, pcoord)
4273 CALL mp_comm_t3c%create(ri_data%pgrid%mp_comm_2d, 3, pdims)
4275 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4276 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4277 CALL dbt_destroy(t_dummy)
4282 ri_data%ri_metric,
"HFX_3c_nl", qs_env, op_pos=2, sym_ij=.false., &
4285 ALLOCATE (ri_cells(nimg))
4288 ALLOCATE (nri_per_atom(natom))
4294 iatom=iatom, jatom=jatom, katom=katom)
4297 IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
4298 any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
4300 jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4301 IF (jcell > nimg .OR. jcell < 1) cycle
4303 IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
4304 any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
4306 kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
4307 IF (kcell > nimg .OR. kcell < 1) cycle
4309 IF (dik > ri_data%kp_RI_range) cycle
4312 IF (jcell == 1 .AND. iatom == jatom) nri_per_atom(iatom) = nri_per_atom(iatom) + ri_data%bsizes_RI(katom)
4316 CALL para_env%sum(ri_cells)
4317 CALL para_env%sum(nri_per_atom)
4319 ALLOCATE (ri_data%img_to_RI_cell(nimg))
4320 ri_data%ncell_RI = 0
4321 ri_data%img_to_RI_cell = 0
4323 IF (ri_cells(i_img) > 0)
THEN
4324 ri_data%ncell_RI = ri_data%ncell_RI + 1
4325 ri_data%img_to_RI_cell(i_img) = ri_data%ncell_RI
4329 ALLOCATE (ri_data%RI_cell_to_img(ri_data%ncell_RI))
4331 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
4335 IF (ri_data%unit_nr > 0)
THEN
4336 WRITE (ri_data%unit_nr, fmt=
"(/T3,A,I29)") &
4337 "KP-HFX_RI_INFO| Number of RI-KP parallel groups:", ngroups
4338 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F31.3,A)") &
4339 "KP-HFX_RI_INFO| RI basis extension radius:", ri_data%kp_RI_range*
angstrom,
" Ang"
4340 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F12.3,A, F6.3, A)") &
4341 "KP-HFX_RI_INFO| RI basis bump factor and bump radius:", bump_fact,
" /", &
4342 ri_data%kp_bump_rad*
angstrom,
" Ang"
4343 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I16,A)") &
4344 "KP-HFX_RI_INFO| The extended RI bases cover up to ", ri_data%ncell_RI,
" unit cells"
4345 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I18)") &
4346 "KP-HFX_RI_INFO| Average number of sgf in extended RI bases:", sum(nri_per_atom)/natom
4347 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F13.3,A)") &
4348 "KP-HFX_RI_INFO| Consider all image cells within a radius of ", ri_data%kp_image_range*
angstrom,
" Ang"
4349 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I27/)") &
4350 "KP-HFX_RI_INFO| Number of image cells considered: ", ri_data%nimg
4354 CALL timestop(handle)
4356 END SUBROUTINE get_kp_and_ri_images
4371 SUBROUTINE get_stack_tensors(res_stack, rho_stack, ints_stack, rho_template, ints_template, &
4372 stack_size, ri_data, qs_env)
4373 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: res_stack, rho_stack, ints_stack
4374 TYPE(dbt_type),
INTENT(INOUT) :: rho_template, ints_template
4375 INTEGER,
INTENT(IN) :: stack_size
4379 INTEGER :: is, nblks, nblks_3c(3), pdims_3d(3)
4380 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_stack, dist1, &
4381 dist2, dist3, dist_stack1, &
4382 dist_stack2, dist_stack3
4383 TYPE(dbt_distribution_type) :: t_dist
4384 TYPE(dbt_pgrid_type) :: pgrid
4391 nblks =
SIZE(ri_data%bsizes_AO_split)
4392 ALLOCATE (bsizes_stack(stack_size*nblks))
4393 DO is = 1, stack_size
4394 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
4397 ALLOCATE (dist1(nblks), dist2(nblks), dist_stack1(stack_size*nblks), dist_stack2(stack_size*nblks))
4398 CALL dbt_get_info(rho_template, proc_dist_1=dist1, proc_dist_2=dist2)
4399 DO is = 1, stack_size
4400 dist_stack1((is - 1)*nblks + 1:is*nblks) = dist1(:)
4401 dist_stack2((is - 1)*nblks + 1:is*nblks) = dist2(:)
4406 CALL dbt_distribution_new(t_dist, ri_data%pgrid_2d, dist_stack1, dist_stack2)
4407 CALL dbt_create(rho_stack(1),
"RHO_stack", t_dist, [1], [2], bsizes_stack, bsizes_stack)
4408 CALL dbt_distribution_destroy(t_dist)
4409 DEALLOCATE (dist1, dist2, dist_stack1, dist_stack2)
4412 CALL create_2c_tensor(rho_stack(2), dist1, dist2, ri_data%pgrid_2d, bsizes_stack, bsizes_stack, name=
"RHO_stack")
4413 DEALLOCATE (dist1, dist2)
4415 CALL dbt_get_info(ints_template, nblks_total=nblks_3c)
4416 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)))
4417 ALLOCATE (dist_stack3(stack_size*nblks_3c(3)), bsizes_ri_ext(nblks_3c(2)))
4418 CALL dbt_get_info(ints_template, proc_dist_1=dist1, proc_dist_2=dist2, &
4419 proc_dist_3=dist3, blk_size_2=bsizes_ri_ext)
4420 DO is = 1, stack_size
4421 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
4425 CALL dbt_distribution_new(t_dist, ri_data%pgrid_1, dist1, dist2, dist_stack3)
4426 CALL dbt_create(ints_stack(1),
"ints_stack", t_dist, [1, 2], [3], ri_data%bsizes_AO_split, &
4427 bsizes_ri_ext, bsizes_stack)
4428 CALL dbt_distribution_destroy(t_dist)
4429 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
4433 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid, tensor_dims=[nblks_3c(1), nblks_3c(2), stack_size*nblks_3c(3)])
4434 CALL create_3c_tensor(ints_stack(2), dist1, dist2, dist3, pgrid, ri_data%bsizes_AO_split, &
4435 bsizes_ri_ext, bsizes_stack, [1, 2], [3], name=
"ints_stack")
4436 DEALLOCATE (dist1, dist2, dist3)
4437 CALL dbt_pgrid_destroy(pgrid)
4440 CALL dbt_create(ints_stack(1), res_stack(1))
4441 CALL dbt_create(ints_stack(2), res_stack(2))
4443 END SUBROUTINE get_stack_tensors
4457 SUBROUTINE fill_3c_stack(t_3c_stack, t_3c_in, images, stack_dim, ri_data, filter_at, filter_dim, &
4458 idx_to_at, img_bounds)
4459 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_stack
4460 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_in
4461 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4462 INTEGER,
INTENT(IN) :: stack_dim
4464 INTEGER,
INTENT(IN),
OPTIONAL :: filter_at, filter_dim
4465 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: idx_to_at
4466 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2)
4468 INTEGER :: dest(3), i_img,
idx, ind(3), lb, nblks, &
4470 LOGICAL :: do_filter, found
4471 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4472 TYPE(dbt_iterator_type) :: iter
4477 nblks =
SIZE(ri_data%bsizes_AO_split)
4480 IF (
PRESENT(filter_at) .AND.
PRESENT(filter_dim) .AND.
PRESENT(idx_to_at)) do_filter = .true.
4485 IF (
PRESENT(img_bounds))
THEN
4487 ub = img_bounds(2) - 1
4493 IF (i_img == 0 .OR. i_img > nimg) cycle
4498 CALL dbt_iterator_start(iter, t_3c_in(i_img))
4499 DO WHILE (dbt_iterator_blocks_left(iter))
4500 CALL dbt_iterator_next_block(iter, ind)
4501 CALL dbt_get_block(t_3c_in(i_img), ind, blk, found)
4502 IF (.NOT. found) cycle
4505 IF (.NOT. idx_to_at(ind(filter_dim)) == filter_at) cycle
4508 IF (stack_dim == 1)
THEN
4509 dest = [(
idx - offset - 1)*nblks + ind(1), ind(2), ind(3)]
4510 ELSE IF (stack_dim == 2)
THEN
4511 dest = [ind(1), (
idx - offset - 1)*nblks + ind(2), ind(3)]
4513 dest = [ind(1), ind(2), (
idx - offset - 1)*nblks + ind(3)]
4516 CALL dbt_put_block(t_3c_stack, dest, shape(blk), blk)
4519 CALL dbt_iterator_stop(iter)
4522 CALL dbt_finalize(t_3c_stack)
4524 END SUBROUTINE fill_3c_stack
4536 SUBROUTINE fill_2c_stack(t_2c_stack, t_2c_in, images, stack_dim, ri_data, img_bounds, shift)
4537 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_stack
4538 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_in
4539 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4540 INTEGER,
INTENT(IN) :: stack_dim
4542 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2), shift
4544 INTEGER :: dest(2), i_img,
idx, ind(2), lb, &
4545 my_shift, nblks, nimg, offset, ub
4547 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
4548 TYPE(dbt_iterator_type) :: iter
4553 nblks =
SIZE(ri_data%bsizes_AO_split)
4558 IF (
PRESENT(img_bounds))
THEN
4560 ub = img_bounds(2) - 1
4565 IF (
PRESENT(shift)) my_shift = shift
4569 IF (i_img == 0 .OR. i_img > nimg) cycle
4573 CALL dbt_iterator_start(iter, t_2c_in(i_img))
4574 DO WHILE (dbt_iterator_blocks_left(iter))
4575 CALL dbt_iterator_next_block(iter, ind)
4576 CALL dbt_get_block(t_2c_in(i_img), ind, blk, found)
4577 IF (.NOT. found) cycle
4579 IF (stack_dim == 1)
THEN
4580 dest = [(
idx - offset - 1)*nblks + ind(1), (my_shift - 1)*nblks + ind(2)]
4582 dest = [(my_shift - 1)*nblks + ind(1), (
idx - offset - 1)*nblks + ind(2)]
4585 CALL dbt_put_block(t_2c_stack, dest, shape(blk), blk)
4588 CALL dbt_iterator_stop(iter)
4591 CALL dbt_finalize(t_2c_stack)
4593 END SUBROUTINE fill_2c_stack
4601 SUBROUTINE unstack_t_3c_apc(t_3c_apc, t_stacked, idx)
4602 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_apc, t_stacked
4603 INTEGER,
INTENT(IN) ::
idx
4605 INTEGER :: current_idx
4606 INTEGER,
DIMENSION(3) :: ind, nblks_3c
4608 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4609 TYPE(dbt_iterator_type) :: iter
4612 CALL dbt_get_info(t_3c_apc, nblks_total=nblks_3c)
4615 CALL dbt_iterator_start(iter, t_stacked)
4616 DO WHILE (dbt_iterator_blocks_left(iter))
4617 CALL dbt_iterator_next_block(iter, ind)
4620 current_idx = (ind(3) - 1)/nblks_3c(3) + 1
4621 IF (.NOT.
idx == current_idx) cycle
4623 CALL dbt_get_block(t_stacked, ind, blk, found)
4624 IF (.NOT. found) cycle
4626 CALL dbt_put_block(t_3c_apc, [ind(1), ind(2), ind(3) - (
idx - 1)*nblks_3c(3)], shape(blk), blk)
4629 CALL dbt_iterator_stop(iter)
4632 END SUBROUTINE unstack_t_3c_apc
4642 SUBROUTINE get_atom_3c_ints(t_3c_at, t_3c_ints, iatom, dim_at, idx_to_at)
4643 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_at, t_3c_ints
4644 INTEGER,
INTENT(IN) :: iatom, dim_at
4645 INTEGER,
DIMENSION(:),
INTENT(IN) :: idx_to_at
4647 INTEGER,
DIMENSION(3) :: ind
4649 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4650 TYPE(dbt_iterator_type) :: iter
4653 CALL dbt_iterator_start(iter, t_3c_ints)
4654 DO WHILE (dbt_iterator_blocks_left(iter))
4655 CALL dbt_iterator_next_block(iter, ind)
4656 IF (.NOT. idx_to_at(ind(dim_at)) == iatom) cycle
4658 CALL dbt_get_block(t_3c_ints, ind, blk, found)
4659 IF (.NOT. found) cycle
4661 CALL dbt_put_block(t_3c_at, ind, shape(blk), blk)
4664 CALL dbt_iterator_stop(iter)
4666 CALL dbt_finalize(t_3c_at)
4668 END SUBROUTINE get_atom_3c_ints
4679 SUBROUTINE precalc_derivatives(t_3c_der_RI, t_3c_der_AO, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
4680 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ri, t_3c_der_ao
4681 TYPE(
dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot
4682 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_2c_der_metric
4686 CHARACTER(LEN=*),
PARAMETER :: routinen =
'precalc_derivatives'
4688 INTEGER :: handle, handle2, i_img, i_mem, i_ri, &
4689 i_xyz, iatom, n_mem, natom, nblks_ri, &
4690 ncell_ri, nimg, nkind, nthreads
4691 INTEGER(int_8) :: nze
4692 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, dist_ao_1, &
4693 dist_ao_2, dist_ri, dist_ri_ext, dummy_end, dummy_start, end_blocks, start_blocks
4694 INTEGER,
DIMENSION(3) :: pcoord, pdims
4695 INTEGER,
DIMENSION(:),
POINTER :: col_bsize, row_bsize
4699 TYPE(
dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_metric
4700 TYPE(dbt_distribution_type) :: t_dist
4701 TYPE(dbt_pgrid_type) :: pgrid
4702 TYPE(dbt_type) :: t_3c_template
4703 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
4708 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4715 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4717 NULLIFY (qs_kind_set, dist_2d, nl_2c, particle_set, dft_control, para_env, row_bsize, col_bsize)
4719 CALL timeset(routinen, handle)
4721 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, natom=natom, &
4722 particle_set=particle_set, dft_control=dft_control, para_env=para_env)
4725 ncell_ri = ri_data%ncell_RI
4727 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4737 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[max(1, natom/(ri_data%n_mem*nthreads)), natom, natom])
4739 CALL create_3c_tensor(t_3c_template, dist_ao_1, dist_ao_2, dist_ri, pgrid, &
4740 ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4741 map1=[1, 2], map2=[3], name=
"tmp")
4742 CALL dbt_destroy(t_3c_template)
4745 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
4746 ALLOCATE (dist_ri_ext(natom*ncell_ri))
4747 ALLOCATE (bsizes_ri_ext(natom*ncell_ri))
4748 ALLOCATE (bsizes_ri_ext_split(nblks_ri*ncell_ri))
4749 DO i_ri = 1, ncell_ri
4750 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
4751 dist_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = dist_ri(:)
4752 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
4755 CALL dbt_distribution_new(t_dist, pgrid, dist_ao_1, dist_ao_2, dist_ri_ext)
4756 CALL dbt_create(t_3c_template,
"KP_3c_der", t_dist, [1, 2], [3], &
4757 ri_data%bsizes_AO, ri_data%bsizes_AO, bsizes_ri_ext)
4758 CALL dbt_distribution_destroy(t_dist)
4760 ALLOCATE (t_3c_der_ri_prv(nimg, 1, 3), t_3c_der_ao_prv(nimg, 1, 3))
4763 CALL dbt_create(t_3c_template, t_3c_der_ri_prv(i_img, 1, i_xyz))
4764 CALL dbt_create(t_3c_template, t_3c_der_ao_prv(i_img, 1, i_xyz))
4767 CALL dbt_destroy(t_3c_template)
4769 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
4770 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
4772 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4773 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4774 CALL dbt_pgrid_destroy(pgrid)
4777 "HFX_3c_nl", qs_env, op_pos=2, sym_jk=.false., own_dist=.true.)
4779 n_mem = ri_data%n_mem
4781 start_blocks, end_blocks)
4782 DEALLOCATE (dummy_start, dummy_end)
4784 CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, ri_data%pgrid_2, &
4785 bsizes_ri_ext_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
4786 map1=[1], map2=[2, 3], name=
"der (RI | AO AO)")
4789 CALL dbt_create(t_3c_template, t_3c_der_ri(i_img, i_xyz))
4790 CALL dbt_create(t_3c_template, t_3c_der_ao(i_img, i_xyz))
4796 nl_3c, basis_set_ao, basis_set_ao, basis_set_ri, &
4797 ri_data%ri_metric, der_eps=ri_data%eps_schwarz_forces, op_pos=2, &
4798 do_kpoints=.true., do_hfx_kpoints=.true., &
4799 bounds_k=[start_blocks(i_mem), end_blocks(i_mem)], &
4800 ri_range=ri_data%kp_RI_range, img_to_ri_cell=ri_data%img_to_RI_cell)
4802 CALL timeset(routinen//
"_cpy", handle2)
4809 CALL dbt_copy(t_3c_der_ao_prv(i_img, 1, i_xyz), t_3c_template, &
4810 order=[3, 2, 1], move_data=.true.)
4811 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4812 CALL dbt_copy(t_3c_template, t_3c_der_ao(i_img, i_xyz), &
4813 order=[1, 3, 2], move_data=.true., summation=.true.)
4819 CALL dbt_copy(t_3c_der_ri_prv(i_img, 1, i_xyz), t_3c_template, &
4820 order=[3, 2, 1], move_data=.true.)
4821 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4822 CALL dbt_copy(t_3c_template, t_3c_der_ri(i_img, i_xyz), &
4823 order=[1, 3, 2], move_data=.true., summation=.true.)
4827 CALL timestop(handle2)
4829 CALL dbt_destroy(t_3c_template)
4834 CALL dbt_destroy(t_3c_der_ri_prv(i_img, 1, i_xyz))
4835 CALL dbt_destroy(t_3c_der_ao_prv(i_img, 1, i_xyz))
4838 DEALLOCATE (t_3c_der_ri_prv, t_3c_der_ao_prv)
4841 CALL reorder_3c_derivs(t_3c_der_ri, ri_data)
4842 CALL reorder_3c_derivs(t_3c_der_ao, ri_data)
4844 CALL timeset(routinen//
"_2c", handle2)
4847 ALLOCATE (row_bsize(
SIZE(ri_data%bsizes_RI)))
4848 ALLOCATE (col_bsize(
SIZE(ri_data%bsizes_RI)))
4849 row_bsize(:) = ri_data%bsizes_RI
4850 col_bsize(:) = ri_data%bsizes_RI
4852 CALL dbcsr_create(dbcsr_template,
"2c_der", dbcsr_dist, dbcsr_type_no_symmetry, &
4853 row_bsize, col_bsize)
4855 DEALLOCATE (col_bsize, row_bsize)
4857 ALLOCATE (mat_der_metric(nimg, 3))
4860 CALL dbcsr_create(mat_der_pot(i_img, i_xyz), template=dbcsr_template)
4861 CALL dbcsr_create(mat_der_metric(i_img, i_xyz), template=dbcsr_template)
4868 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
4870 basis_set_ri, basis_set_ri, ri_data%hfx_pot, do_kpoints=.true.)
4875 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
4877 basis_set_ri, basis_set_ri, ri_data%ri_metric, do_kpoints=.true.)
4883 CALL dbt_create(ri_data%t_2c_inv(1, 1), t_2c_der_metric(iatom, i_xyz))
4884 CALL get_ext_2c_int(t_2c_der_metric(iatom, i_xyz), mat_der_metric(:, i_xyz), &
4885 iatom, iatom, 1, ri_data, qs_env)
4891 CALL timestop(handle2)
4893 CALL timestop(handle)
4895 END SUBROUTINE precalc_derivatives
4916 SUBROUTINE get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, img, pref, &
4917 ri_data, qs_env, work_virial, cell, particle_set, diag, offdiag)
4920 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_contr
4921 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_der
4922 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
4923 INTEGER,
INTENT(IN) :: img
4924 REAL(
dp),
INTENT(IN) :: pref
4927 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
4928 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
4930 POINTER :: particle_set
4931 LOGICAL,
INTENT(IN),
OPTIONAL ::
diag, offdiag
4933 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_2c_der_force'
4935 INTEGER :: handle, i_img, i_ri, i_xyz, iat, &
4936 iat_of_kind, ikind, j_img, j_ri, &
4937 j_xyz, jat, jat_of_kind, jkind, natom
4938 INTEGER,
DIMENSION(2) :: ind
4939 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
4940 LOGICAL :: found, my_diag, my_offdiag, use_virial
4941 REAL(
dp) :: new_force
4942 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :),
TARGET :: contr_blk, der_blk
4943 REAL(
dp),
DIMENSION(3) :: scoord
4944 TYPE(dbt_iterator_type) :: iter
4947 NULLIFY (kpoints, index_to_cell)
4952 CALL timeset(routinen, handle)
4954 use_virial = .false.
4955 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
4960 my_offdiag = .false.
4961 IF (
PRESENT(
diag)) my_offdiag = offdiag
4963 CALL get_qs_env(qs_env, kpoints=kpoints, natom=natom)
4972 CALL dbt_iterator_start(iter, t_2c_der(i_xyz))
4973 DO WHILE (dbt_iterator_blocks_left(iter))
4974 CALL dbt_iterator_next_block(iter, ind)
4977 IF ((my_diag .AND. .NOT. my_offdiag) .OR. (.NOT. my_diag .AND. my_offdiag))
THEN
4978 IF (my_diag .AND. (ind(1) .NE. ind(2))) cycle
4979 IF (my_offdiag .AND. (ind(1) == ind(2))) cycle
4982 CALL dbt_get_block(t_2c_der(i_xyz), ind, der_blk, found)
4984 CALL dbt_get_block(t_2c_contr, ind, contr_blk, found)
4990 new_force = pref*sum(der_blk(:, :)*contr_blk(:, :))
4992 i_ri = (ind(1) - 1)/natom + 1
4993 i_img = ri_data%RI_cell_to_img(i_ri)
4994 iat = ind(1) - (i_ri - 1)*natom
4995 iat_of_kind = atom_of_kind(iat)
4996 ikind = kind_of(iat)
4998 j_ri = (ind(2) - 1)/natom + 1
4999 j_img = ri_data%RI_cell_to_img(j_ri)
5000 jat = ind(2) - (j_ri - 1)*natom
5001 jat_of_kind = atom_of_kind(jat)
5002 jkind = kind_of(jat)
5006 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5009 IF (use_virial)
THEN
5012 scoord(:) = scoord(:) + real(index_to_cell(:, i_img),
dp)
5016 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5022 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5025 IF (use_virial)
THEN
5028 scoord(:) = scoord(:) + real(index_to_cell(:, j_img) + index_to_cell(:, img),
dp)
5032 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - new_force*scoord(j_xyz)
5036 DEALLOCATE (contr_blk)
5039 DEALLOCATE (der_blk)
5041 CALL dbt_iterator_stop(iter)
5045 CALL timestop(handle)
5071 idx_to_at_RI, idx_to_at_AO, i_images, lb_img, pref, &
5072 ri_data, qs_env, work_virial, cell, particle_set)
5075 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_contr
5076 TYPE(dbt_type),
DIMENSION(3),
INTENT(INOUT) :: t_3c_der_1, t_3c_der_2
5077 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of, idx_to_at_ri, &
5078 idx_to_at_ao, i_images
5079 INTEGER,
INTENT(IN) :: lb_img
5080 REAL(
dp),
INTENT(IN) :: pref
5083 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
5084 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
5086 POINTER :: particle_set
5088 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_force_from_3c_trace'
5090 INTEGER :: handle, i_ri, i_xyz, iat, iat_of_kind,
idx, ikind, j_xyz, jat, jat_of_kind, &
5091 jkind, kat, kat_of_kind, kkind, nblks_ao, nblks_ri, ri_img
5092 INTEGER,
DIMENSION(3) :: ind
5093 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
5094 LOGICAL :: found, found_1, found_2, use_virial
5095 REAL(
dp) :: new_force
5096 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :),
TARGET :: contr_blk, der_blk_1, der_blk_2, &
5098 REAL(
dp),
DIMENSION(3) :: scoord
5099 TYPE(dbt_iterator_type) :: iter
5102 NULLIFY (kpoints, index_to_cell)
5104 CALL timeset(routinen, handle)
5109 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
5110 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
5112 use_virial = .false.
5113 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
5120 CALL dbt_iterator_start(iter, t_3c_contr)
5121 DO WHILE (dbt_iterator_blocks_left(iter))
5122 CALL dbt_iterator_next_block(iter, ind)
5124 CALL dbt_get_block(t_3c_contr, ind, contr_blk, found)
5128 CALL dbt_get_block(t_3c_der_1(i_xyz), ind, der_blk_1, found_1)
5129 IF (.NOT. found_1)
THEN
5130 DEALLOCATE (der_blk_1)
5131 ALLOCATE (der_blk_1(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5132 der_blk_1(:, :, :) = 0.0_dp
5134 CALL dbt_get_block(t_3c_der_2(i_xyz), ind, der_blk_2, found_2)
5135 IF (.NOT. found_2)
THEN
5136 DEALLOCATE (der_blk_2)
5137 ALLOCATE (der_blk_2(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5138 der_blk_2(:, :, :) = 0.0_dp
5141 ALLOCATE (der_blk_3(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5142 der_blk_3(:, :, :) = -(der_blk_1(:, :, :) + der_blk_2(:, :, :))
5148 new_force = pref*sum(der_blk_1(:, :, :)*contr_blk(:, :, :))
5150 i_ri = (ind(1) - 1)/nblks_ri + 1
5151 ri_img = ri_data%RI_cell_to_img(i_ri)
5152 iat = idx_to_at_ri(ind(1) - (i_ri - 1)*nblks_ri)
5153 iat_of_kind = atom_of_kind(iat)
5154 ikind = kind_of(iat)
5157 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5160 IF (use_virial)
THEN
5163 scoord(:) = scoord(:) + real(index_to_cell(:, ri_img),
dp)
5167 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5172 new_force = pref*sum(der_blk_2(:, :, :)*contr_blk(:, :, :))
5173 jat = idx_to_at_ao(ind(2))
5174 jat_of_kind = atom_of_kind(jat)
5175 jkind = kind_of(jat)
5178 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5181 IF (use_virial)
THEN
5187 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5193 new_force = pref*sum(der_blk_3(:, :, :)*contr_blk(:, :, :))
5194 idx = (ind(3) - 1)/nblks_ao + 1
5195 kat = idx_to_at_ao(ind(3) - (
idx - 1)*nblks_ao)
5196 kat_of_kind = atom_of_kind(kat)
5197 kkind = kind_of(kat)
5200 force(kkind)%fock_4c(i_xyz, kat_of_kind) = force(kkind)%fock_4c(i_xyz, kat_of_kind) &
5203 IF (use_virial)
THEN
5205 scoord(:) = scoord(:) + real(index_to_cell(:, i_images(lb_img - 1 +
idx)),
dp)
5209 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5213 DEALLOCATE (der_blk_1, der_blk_2, der_blk_3)
5215 DEALLOCATE (contr_blk)
5218 CALL dbt_iterator_stop(iter)
5220 CALL timestop(handle)
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Types and set/get functions for auxiliary density matrix methods.
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public bussy2024
Handles all functions related to the CELL.
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
methods related to the blacs parallel environment
subroutine, public cp_blacs_env_release(blacs_env)
releases the given blacs_env
subroutine, public cp_blacs_env_create(blacs_env, para_env, blacs_grid_layout, blacs_repeatable, row_major, grid_2d)
allocates and initializes a type that represent a blacs context
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_distribution_release(dist)
...
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_finalize(matrix)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_release(matrix)
...
subroutine, public dbcsr_clear(matrix)
...
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, uplo_to_full)
used to replace the cholesky decomposition by the inverse
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Routines that link DBCSR and CP2K concepts together.
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_power(matrix, exponent, threshold, n_dependent, para_env, blacs_env, verbose, eigenvectors, eigenvalues)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
This is the start of a dbt_api, all publically needed functions are exported here....
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
subroutine, public distribution_2d_create(distribution_2d, blacs_env, local_rows_ptr, n_local_rows, local_cols_ptr, row_distribution_ptr, col_distribution_ptr, n_local_cols, n_row_distribution, n_col_distribution)
initializes the distribution_2d
subroutine, public distribution_2d_release(distribution_2d)
...
RI-methods for HFX and K-points. \auhtor Augustin Bussy (01.2023)
subroutine, public hfx_ri_update_forces_kp(qs_env, ri_data, nspins, hf_fraction, rho_ao, use_virial)
Update the K-points RI-HFX forces.
subroutine, public hfx_ri_update_ks_kp(qs_env, ri_data, ks_matrix, ehfx, rho_ao, geometry_did_change, nspins, hf_fraction)
Update the KS matrices for each real-space image.
subroutine, public get_force_from_3c_trace(force, t_3c_contr, t_3c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, deriv_dim)
This routines calculates the force contribution from a trace over 3D tensors, i.e....
subroutine, public get_idx_to_atom(idx_to_at, bsizes_split, bsizes_orig)
a small utility function that returns the atom corresponding to a block of a split tensor
subroutine, public hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_ri, t_2c_int_pot, t_3c_int, do_kpoints)
Calculate 2-center and 3-center integrals.
subroutine, public get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, do_ovlp)
Update the forces due to the derivative of the a 2-center product d/dR (Q|R)
Types and set/get functions for HFX.
Routines useful for iterative matrix calculations.
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
real(kind=dp), parameter, public cutoff_screen_factor
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Collection of simple mathematical functions and subroutines.
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Interface to the message passing library MPI.
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public angstrom
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
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.