35 USE dbcsr_api,
ONLY: &
36 dbcsr_add, dbcsr_clear, dbcsr_copy, dbcsr_create, dbcsr_distribution_get, &
37 dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, dbcsr_dot, &
38 dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, &
39 dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
40 dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_put_block, dbcsr_release, &
41 dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric
43 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
44 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
45 dbt_distribution_destroy, dbt_distribution_new, dbt_distribution_type, dbt_filter, &
46 dbt_finalize, dbt_get_block, dbt_get_info, dbt_get_stored_coordinates, &
47 dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, &
48 dbt_iterator_type, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
49 dbt_pgrid_type, dbt_put_block, dbt_scale, dbt_type
109#include "./base/base_uses.f90"
118 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'hfx_ri_kp'
135 SUBROUTINE adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
136 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
137 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
138 TYPE(qs_environment_type),
POINTER :: qs_env
140 INTEGER :: i_img, i_RI, i_spin, iatom, natom, &
141 nblks_RI, nimg, nkind, nspins
142 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_RI_ext, dist1, dist2, dist3
143 TYPE(dft_control_type),
POINTER :: dft_control
144 TYPE(mp_para_env_type),
POINTER :: para_env
146 NULLIFY (dft_control, para_env)
153 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, para_env=para_env, nkind=nkind)
157 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
158 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
159 DO i_ri = 1, ri_data%ncell_RI
160 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
163 ALLOCATE (ri_data%t_3c_int_ctr_1(1, nimg))
165 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
166 ri_data%bsizes_AO_split, [1, 2], [3], name=
"(AO RI | AO)")
169 CALL dbt_create(ri_data%t_3c_int_ctr_1(1, 1), ri_data%t_3c_int_ctr_1(1, i_img))
171 DEALLOCATE (dist1, dist2, dist3)
173 ALLOCATE (ri_data%t_3c_int_ctr_2(1, 1))
175 ri_data%pgrid_1, ri_data%bsizes_AO_split, bsizes_ri_ext, &
176 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
177 DEALLOCATE (dist1, dist2, dist3)
180 DEALLOCATE (bsizes_ri_ext)
181 nblks_ri =
SIZE(ri_data%bsizes_RI)
182 ALLOCATE (bsizes_ri_ext(nblks_ri*ri_data%ncell_RI))
183 DO i_ri = 1, ri_data%ncell_RI
184 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
187 ALLOCATE (ri_data%t_2c_inv(1, natom), ri_data%t_2c_int(1, natom), ri_data%t_2c_pot(1, natom))
188 CALL create_2c_tensor(ri_data%t_2c_inv(1, 1), dist1, dist2, ri_data%pgrid_2d, &
189 bsizes_ri_ext, bsizes_ri_ext, &
191 DEALLOCATE (dist1, dist2)
192 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, 1))
193 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, 1))
195 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_inv(1, iatom))
196 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_int(1, iatom))
197 CALL dbt_create(ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, iatom))
200 ALLOCATE (ri_data%kp_cost(natom, natom, nimg))
201 ri_data%kp_cost = 0.0_dp
204 nspins = dft_control%nspins
205 ALLOCATE (ri_data%rho_ao_t(nspins, nimg), ri_data%ks_t(nspins, nimg))
206 CALL create_2c_tensor(ri_data%rho_ao_t(1, 1), dist1, dist2, ri_data%pgrid_2d, &
207 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
209 DEALLOCATE (dist1, dist2)
211 CALL dbt_create(dbcsr_template, ri_data%ks_t(1, 1))
213 IF (nspins == 2)
THEN
214 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(2, 1))
215 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(2, 1))
218 DO i_spin = 1, nspins
219 CALL dbt_create(ri_data%rho_ao_t(1, 1), ri_data%rho_ao_t(i_spin, i_img))
220 CALL dbt_create(ri_data%ks_t(1, 1), ri_data%ks_t(i_spin, i_img))
224 END SUBROUTINE adapt_ri_data_to_kp
232 SUBROUTINE hfx_ri_pre_scf_kp(dbcsr_template, ri_data, qs_env)
233 TYPE(dbcsr_type),
INTENT(INOUT) :: dbcsr_template
234 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
235 TYPE(qs_environment_type),
POINTER :: qs_env
237 CHARACTER(LEN=*),
PARAMETER :: routineN =
'hfx_ri_pre_scf_kp'
239 INTEGER :: handle, i_img, iatom, natom, nimg, nkind
240 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_op_pot, t_2c_op_RI
241 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_int
242 TYPE(dft_control_type),
POINTER :: dft_control
244 NULLIFY (dft_control)
246 CALL timeset(routinen, handle)
248 CALL get_qs_env(qs_env, dft_control=dft_control, natom=natom, nkind=nkind)
250 CALL cleanup_kp(ri_data)
253 IF (ri_data%flavor .NE.
ri_pmat) cpabort(
"K-points RI-HFX only with RHO flavor")
254 IF (ri_data%same_op) ri_data%same_op = .false.
255 IF (abs(ri_data%eps_pgf_orb - dft_control%qs_control%eps_pgf_orb) > 1.0e-16_dp) &
256 cpabort(
"RI%EPS_PGF_ORB and QS%EPS_PGF_ORB must be identical for RI-HFX k-points")
258 CALL get_kp_and_ri_images(ri_data, qs_env)
262 ALLOCATE (t_2c_op_pot(nimg), t_2c_op_ri(nimg))
263 ALLOCATE (t_3c_int(1, nimg))
267 CALL adapt_ri_data_to_kp(dbcsr_template, ri_data, qs_env)
272 CALL get_ext_2c_int(ri_data%t_2c_inv(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
276 CALL get_ext_2c_int(ri_data%t_2c_int(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, &
277 qs_env, off_diagonal=.true.)
278 CALL apply_bump(ri_data%t_2c_int(1, iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
281 CALL get_ext_2c_int(ri_data%t_2c_pot(1, iatom), t_2c_op_ri, iatom, iatom, 1, ri_data, qs_env, &
282 do_inverse=.true., skip_inverse=.true.)
283 CALL apply_bump(ri_data%t_2c_pot(1, iatom), iatom, ri_data, qs_env, from_left=.true., &
284 from_right=.true., debump=.true.)
289 CALL dbcsr_release(t_2c_op_ri(i_img))
292 ALLOCATE (ri_data%kp_mat_2c_pot(1, nimg))
294 CALL dbcsr_create(ri_data%kp_mat_2c_pot(1, i_img), template=t_2c_op_pot(i_img))
295 CALL dbcsr_copy(ri_data%kp_mat_2c_pot(1, i_img), t_2c_op_pot(i_img))
296 CALL dbcsr_release(t_2c_op_pot(i_img))
301 CALL precontract_3c_ints(t_3c_int, ri_data, qs_env)
304 CALL reorder_3c_ints(ri_data%t_3c_int_ctr_1(1, :), ri_data)
306 CALL timestop(handle)
308 END SUBROUTINE hfx_ri_pre_scf_kp
314 SUBROUTINE cleanup_kp(ri_data)
315 TYPE(hfx_ri_type),
INTENT(INOUT) :: ri_data
319 IF (
ALLOCATED(ri_data%kp_cost))
DEALLOCATE (ri_data%kp_cost)
320 IF (
ALLOCATED(ri_data%idx_to_img))
DEALLOCATE (ri_data%idx_to_img)
321 IF (
ALLOCATED(ri_data%img_to_idx))
DEALLOCATE (ri_data%img_to_idx)
322 IF (
ALLOCATED(ri_data%present_images))
DEALLOCATE (ri_data%present_images)
323 IF (
ALLOCATED(ri_data%img_to_RI_cell))
DEALLOCATE (ri_data%img_to_RI_cell)
324 IF (
ALLOCATED(ri_data%RI_cell_to_img))
DEALLOCATE (ri_data%RI_cell_to_img)
326 IF (
ALLOCATED(ri_data%kp_mat_2c_pot))
THEN
327 DO j = 1,
SIZE(ri_data%kp_mat_2c_pot, 2)
328 DO i = 1,
SIZE(ri_data%kp_mat_2c_pot, 1)
329 CALL dbcsr_release(ri_data%kp_mat_2c_pot(i, j))
332 DEALLOCATE (ri_data%kp_mat_2c_pot)
335 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
336 DO i = 1,
SIZE(ri_data%kp_t_3c_int)
337 CALL dbt_destroy(ri_data%kp_t_3c_int(i))
339 DEALLOCATE (ri_data%kp_t_3c_int)
342 IF (
ALLOCATED(ri_data%t_2c_inv))
THEN
343 DO j = 1,
SIZE(ri_data%t_2c_inv, 2)
344 DO i = 1,
SIZE(ri_data%t_2c_inv, 1)
345 CALL dbt_destroy(ri_data%t_2c_inv(i, j))
348 DEALLOCATE (ri_data%t_2c_inv)
351 IF (
ALLOCATED(ri_data%t_2c_int))
THEN
352 DO j = 1,
SIZE(ri_data%t_2c_int, 2)
353 DO i = 1,
SIZE(ri_data%t_2c_int, 1)
354 CALL dbt_destroy(ri_data%t_2c_int(i, j))
357 DEALLOCATE (ri_data%t_2c_int)
360 IF (
ALLOCATED(ri_data%t_2c_pot))
THEN
361 DO j = 1,
SIZE(ri_data%t_2c_pot, 2)
362 DO i = 1,
SIZE(ri_data%t_2c_pot, 1)
363 CALL dbt_destroy(ri_data%t_2c_pot(i, j))
366 DEALLOCATE (ri_data%t_2c_pot)
369 IF (
ALLOCATED(ri_data%t_3c_int_ctr_1))
THEN
370 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_1, 2)
371 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_1, 1)
372 CALL dbt_destroy(ri_data%t_3c_int_ctr_1(i, j))
375 DEALLOCATE (ri_data%t_3c_int_ctr_1)
378 IF (
ALLOCATED(ri_data%t_3c_int_ctr_2))
THEN
379 DO j = 1,
SIZE(ri_data%t_3c_int_ctr_2, 2)
380 DO i = 1,
SIZE(ri_data%t_3c_int_ctr_2, 1)
381 CALL dbt_destroy(ri_data%t_3c_int_ctr_2(i, j))
384 DEALLOCATE (ri_data%t_3c_int_ctr_2)
387 IF (
ALLOCATED(ri_data%rho_ao_t))
THEN
388 DO j = 1,
SIZE(ri_data%rho_ao_t, 2)
389 DO i = 1,
SIZE(ri_data%rho_ao_t, 1)
390 CALL dbt_destroy(ri_data%rho_ao_t(i, j))
393 DEALLOCATE (ri_data%rho_ao_t)
396 IF (
ALLOCATED(ri_data%ks_t))
THEN
397 DO j = 1,
SIZE(ri_data%ks_t, 2)
398 DO i = 1,
SIZE(ri_data%ks_t, 1)
399 CALL dbt_destroy(ri_data%ks_t(i, j))
402 DEALLOCATE (ri_data%ks_t)
405 END SUBROUTINE cleanup_kp
419 geometry_did_change, nspins, hf_fraction)
423 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: ks_matrix
424 REAL(kind=
dp),
INTENT(OUT) :: ehfx
425 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_ao
426 LOGICAL,
INTENT(IN) :: geometry_did_change
427 INTEGER,
INTENT(IN) :: nspins
428 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
430 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_ks_kp'
432 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_spin, iatom, &
433 iblk, igroup, jatom, mb_img, n_batch_nze, natom, ngroups, nimg, nimg_nze
434 INTEGER(int_8) :: nflop, nze
435 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges_at, batch_ranges_nze, &
437 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
438 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: sparsity_pattern
439 LOGICAL :: use_delta_p
440 REAL(
dp) :: etmp,
fac, occ, pfac, pref, t1, t2, t3, &
443 TYPE(dbcsr_type) :: ks_desymm, rho_desymm, tmp
444 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
445 TYPE(dbcsr_type),
POINTER :: dbcsr_template
446 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ks_t_split, t_2c_ao_tmp, t_2c_work, &
447 t_3c_int, t_3c_work_2, t_3c_work_3
448 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: ks_t, ks_t_sub, t_3c_apc, t_3c_apc_sub
452 NULLIFY (para_env, para_env_sub, blacs_env_sub, hfx_section, dbcsr_template)
454 CALL timeset(routinen, handle)
456 CALL get_qs_env(qs_env, para_env=para_env, natom=natom)
458 IF (nspins == 1)
THEN
459 fac = 0.5_dp*hf_fraction
461 fac = 1.0_dp*hf_fraction
464 IF (geometry_did_change)
THEN
465 CALL hfx_ri_pre_scf_kp(ks_matrix(1, 1)%matrix, ri_data, qs_env)
468 nimg_nze = ri_data%nimg_nze
494 IF (.NOT. use_delta_p) pfac = 0.0_dp
495 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, pfac, ri_data, qs_env)
497 ALLOCATE (ks_t(nspins, nimg))
499 DO i_spin = 1, nspins
500 CALL dbt_create(ri_data%ks_t(1, 1), ks_t(i_spin, i_img))
504 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
505 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
513 ALLOCATE (t_3c_apc(nspins, nimg))
515 DO i_spin = 1, nspins
516 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
519 CALL contract_pmat_3c(t_3c_apc, ri_data%rho_ao_t, ri_data, qs_env)
524 ri_data%kp_stack_size = batch_size
526 IF (mod(para_env%num_pe, ngroups) .NE. 0)
THEN
527 cpwarn(
"KP_NGROUPS must be an integer divisor of the total number of MPI ranks. It was set to 1.")
531 IF ((mod(ngroups, natom) .NE. 0) .AND. (mod(natom, ngroups) .NE. 0) .AND. geometry_did_change)
THEN
533 cpwarn(
"Better load balancing is reached if NGROUPS is a multiple/divisor of the number of atoms")
535 group_size = para_env%num_pe/ngroups
536 igroup = para_env%mepos/group_size
538 ALLOCATE (para_env_sub)
539 CALL para_env_sub%from_split(para_env, igroup)
543 ALLOCATE (sparsity_pattern(natom, natom, nimg))
544 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
545 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
548 ALLOCATE (mat_2c_pot(nimg), ks_t_sub(nspins, nimg), t_2c_ao_tmp(1), ks_t_split(2), t_2c_work(3))
549 CALL get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
550 group_size, ngroups, para_env, para_env_sub, ri_data)
552 ALLOCATE (t_3c_int(nimg), t_3c_apc_sub(nspins, nimg), t_3c_work_2(3), t_3c_work_3(3))
553 CALL get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
554 group_size, ngroups, para_env, para_env_sub, ri_data)
558 ALLOCATE (batch_ranges_at(natom + 1))
559 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
561 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
562 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
564 batch_ranges_at(iatom) = iblk
568 n_batch_nze = nimg_nze/batch_size
569 IF (
modulo(nimg_nze, batch_size) .NE. 0) n_batch_nze = n_batch_nze + 1
570 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
571 DO i_batch = 1, n_batch_nze
572 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
574 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
576 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
577 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
578 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
579 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
582 ri_data%kp_cost(:, :, :) = 0.0_dp
583 ALLOCATE (iapc_pairs(nimg, 2))
585 CALL dbt_batched_contract_init(ks_t_split(1))
586 CALL dbt_batched_contract_init(ks_t_split(2))
589 IF (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup) cycle
591 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp
597 CALL timeset(routinen//
"_2c", handle2)
598 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
599 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
600 dbcsr_template=dbcsr_template)
601 CALL dbt_copy(t_2c_work(1), t_2c_work(2), move_data=.true.)
602 CALL dbt_filter(t_2c_work(2), ri_data%filter_eps)
603 CALL timestop(handle2)
605 CALL dbt_batched_contract_init(t_2c_work(2))
606 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env)
607 CALL timeset(routinen//
"_3c", handle2)
610 DO i_batch = 1, n_batch_nze
611 CALL fill_3c_stack(t_3c_work_3(3), t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
612 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
613 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
614 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1), move_data=.true.)
616 CALL dbt_contract(1.0_dp, t_2c_work(2), t_3c_work_3(1), &
617 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
618 contract_1=[2], notcontract_1=[1], &
619 contract_2=[1], notcontract_2=[2, 3], &
620 filter_eps=ri_data%filter_eps, flop=nflop)
621 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
622 CALL dbt_copy(t_3c_work_3(2), t_3c_work_2(2), order=[2, 1, 3], move_data=.true.)
623 CALL dbt_copy(t_3c_work_3(3), t_3c_work_3(1))
627 DO i_spin = 1, nspins
628 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
629 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
630 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
633 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
634 CALL dbt_contract(-pref*
fac, t_3c_work_2(1), t_3c_work_2(2), &
635 1.0_dp, ks_t_split(i_spin), map_1=[1], map_2=[2], &
636 contract_1=[2, 3], notcontract_1=[1], &
637 contract_2=[2, 3], notcontract_2=[1], &
638 filter_eps=ri_data%filter_eps, &
639 move_data=i_spin == nspins, flop=nflop)
640 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
643 CALL timestop(handle2)
644 CALL dbt_batched_contract_finalize(t_2c_work(2))
647 ri_data%kp_cost(iatom, jatom, b_img) = t4 - t3
650 CALL dbt_batched_contract_finalize(ks_t_split(1))
651 CALL dbt_batched_contract_finalize(ks_t_split(2))
653 DO i_spin = 1, nspins
654 CALL dbt_copy(ks_t_split(i_spin), t_2c_ao_tmp(1), move_data=.true.)
655 CALL dbt_copy(t_2c_ao_tmp(1), ks_t_sub(i_spin, b_img), summation=.true.)
658 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
659 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
660 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
661 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
663 CALL para_env%sum(ri_data%dbcsr_nflop)
664 CALL para_env%sum(ri_data%kp_cost)
666 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
669 CALL gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
673 CALL dbt_copy(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img), move_data=.true.)
677 CALL dbt_destroy(t_2c_ao_tmp(1))
678 CALL dbt_destroy(ks_t_split(1))
679 CALL dbt_destroy(ks_t_split(2))
680 CALL dbt_destroy(t_2c_work(1))
681 CALL dbt_destroy(t_2c_work(2))
682 CALL dbt_destroy(t_3c_work_2(1))
683 CALL dbt_destroy(t_3c_work_2(2))
684 CALL dbt_destroy(t_3c_work_2(3))
685 CALL dbt_destroy(t_3c_work_3(1))
686 CALL dbt_destroy(t_3c_work_3(2))
687 CALL dbt_destroy(t_3c_work_3(3))
689 CALL dbt_destroy(t_3c_int(i_img))
690 CALL dbcsr_release(mat_2c_pot(i_img))
691 DO i_spin = 1, nspins
692 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
693 CALL dbt_destroy(ks_t_sub(i_spin, i_img))
696 IF (
ASSOCIATED(dbcsr_template))
THEN
697 CALL dbcsr_release(dbcsr_template)
698 DEALLOCATE (dbcsr_template)
703 CALL para_env_sub%free()
704 DEALLOCATE (para_env_sub)
709 CALL get_pmat_images(ri_data%rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
710 DO i_spin = 1, nspins
712 CALL dbt_copy(ks_t(i_spin, b_img), ri_data%ks_t(i_spin, b_img), summation=.true.)
715 mb_img = get_opp_index(b_img, qs_env)
716 IF (mb_img > 0 .AND. mb_img .LE. nimg)
THEN
717 CALL dbt_copy(ks_t(i_spin, mb_img), ri_data%ks_t(i_spin, b_img), order=[2, 1], summation=.true.)
722 DO i_spin = 1, nspins
723 CALL dbt_destroy(ks_t(i_spin, b_img))
728 CALL dbt_create(ri_data%ks_t(1, 1), t_2c_ao_tmp(1))
729 CALL dbcsr_create(tmp, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_symmetric)
730 CALL dbcsr_create(ks_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
731 CALL dbcsr_create(rho_desymm, template=ks_matrix(1, 1)%matrix, matrix_type=dbcsr_type_no_symmetry)
734 DO i_spin = 1, nspins
735 CALL dbt_filter(ri_data%ks_t(i_spin, i_img), ri_data%filter_eps)
736 CALL dbt_copy(ri_data%ks_t(i_spin, i_img), t_2c_ao_tmp(1))
737 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), ks_desymm)
738 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), tmp)
739 CALL dbcsr_add(ks_matrix(i_spin, i_img)%matrix, tmp, 1.0_dp, 1.0_dp)
741 CALL dbt_copy(ri_data%rho_ao_t(i_spin, i_img), t_2c_ao_tmp(1))
742 CALL dbt_copy_tensor_to_matrix(t_2c_ao_tmp(1), rho_desymm)
744 CALL dbcsr_dot(ks_desymm, rho_desymm, etmp)
745 ehfx = ehfx + 0.5_dp*etmp
747 IF (.NOT. use_delta_p)
CALL dbt_clear(ri_data%ks_t(i_spin, i_img))
750 CALL dbcsr_release(rho_desymm)
751 CALL dbcsr_release(ks_desymm)
752 CALL dbcsr_release(tmp)
753 CALL dbt_destroy(t_2c_ao_tmp(1))
755 CALL timestop(handle)
774 INTEGER,
INTENT(IN) :: nspins
775 REAL(kind=
dp),
INTENT(IN) :: hf_fraction
776 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: rho_ao
777 LOGICAL,
INTENT(IN),
OPTIONAL :: use_virial
779 CHARACTER(LEN=*),
PARAMETER :: routinen =
'hfx_ri_update_forces_kp'
781 INTEGER :: b_img, batch_size, group_size, handle, handle2, i_batch, i_img, i_loop, i_spin, &
782 i_xyz, iatom, iblk, igroup, j_xyz, jatom, k_xyz, n_batch, natom, ngroups, nimg, nimg_nze
783 INTEGER(int_8) :: nflop, nze
784 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, batch_ranges_at, &
785 batch_ranges_nze, dist1, dist2, &
786 i_images, idx_to_at_ao, idx_to_at_ri, &
788 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: iapc_pairs
789 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: force_pattern, sparsity_pattern
790 INTEGER,
DIMENSION(2, 1) :: bounds_iat, bounds_jat
791 LOGICAL :: use_virial_prv
792 REAL(
dp) ::
fac, occ, pref, t1, t2
793 REAL(
dp),
DIMENSION(3, 3) :: work_virial
797 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:) :: mat_2c_pot
798 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_pot, mat_der_pot_sub
799 TYPE(dbcsr_type),
POINTER :: dbcsr_template
800 TYPE(dbt_type) :: t_2c_r, t_2c_r_split
801 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_2c_bint, t_2c_binv, t_2c_der_pot, &
802 t_2c_inv, t_2c_metric, t_2c_work, &
803 t_3c_der_stack, t_3c_work_2, &
805 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
806 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, &
814 NULLIFY (para_env, para_env_sub, hfx_section, blacs_env_sub, dbcsr_template, force, atomic_kind_set, &
815 virial, particle_set, cell)
817 CALL timeset(routinen, handle)
819 use_virial_prv = .false.
820 IF (
PRESENT(use_virial)) use_virial_prv = use_virial
822 IF (nspins == 1)
THEN
823 fac = 0.5_dp*hf_fraction
825 fac = 1.0_dp*hf_fraction
828 CALL get_qs_env(qs_env, natom=natom, para_env=para_env, force=force, cell=cell, virial=virial, &
829 atomic_kind_set=atomic_kind_set, particle_set=particle_set)
832 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
833 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
835 ALLOCATE (idx_to_at_ri(
SIZE(ri_data%bsizes_RI_split)))
836 CALL get_idx_to_atom(idx_to_at_ri, ri_data%bsizes_RI_split, ri_data%bsizes_RI)
839 ALLOCATE (t_3c_der_ri(nimg, 3), t_3c_der_ao(nimg, 3), mat_der_pot(nimg, 3), t_2c_der_metric(natom, 3))
843 CALL precalc_derivatives(t_3c_der_ri, t_3c_der_ao, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
846 ALLOCATE (rho_ao_t(nspins, nimg))
848 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
850 DEALLOCATE (dist1, dist2)
851 IF (nspins == 2)
CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(2, 1))
853 DO i_spin = 1, nspins
854 CALL dbt_create(rho_ao_t(1, 1), rho_ao_t(i_spin, i_img))
857 CALL get_pmat_images(rho_ao_t, rho_ao, 0.0_dp, ri_data, qs_env)
860 ALLOCATE (t_3c_apc(nspins, nimg))
862 DO i_spin = 1, nspins
863 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_apc(i_spin, i_img))
866 CALL contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
871 group_size = para_env%num_pe/ngroups
872 igroup = para_env%mepos/group_size
874 ALLOCATE (para_env_sub)
875 CALL para_env_sub%from_split(para_env, igroup)
879 ALLOCATE (sparsity_pattern(natom, natom, nimg))
880 CALL get_sparsity_pattern(sparsity_pattern, ri_data, qs_env)
881 CALL get_sub_dist(sparsity_pattern, ngroups, ri_data)
884 ALLOCATE (t_2c_inv(natom), mat_2c_pot(nimg), rho_ao_t_sub(nspins, nimg), t_2c_work(5), &
885 t_2c_der_metric_sub(natom, 3), mat_der_pot_sub(nimg, 3), t_2c_bint(natom), &
886 t_2c_metric(natom), t_2c_binv(natom))
887 CALL get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
888 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
889 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
890 CALL dbt_create(t_2c_work(1), t_2c_r)
891 CALL dbt_create(t_2c_work(5), t_2c_r_split)
893 ALLOCATE (t_2c_der_pot(3))
895 CALL dbt_create(t_2c_r, t_2c_der_pot(i_xyz))
899 ALLOCATE (t_3c_work_2(3), t_3c_work_3(4), t_3c_der_stack(6), t_3c_der_ao_sub(nimg, 3), &
900 t_3c_der_ri_sub(nimg, 3), t_3c_apc_sub(nspins, nimg))
901 CALL get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_ao, t_3c_der_ao_sub, &
902 t_3c_der_ri, t_3c_der_ri_sub, t_3c_apc, t_3c_apc_sub, t_3c_der_stack, &
903 group_size, ngroups, para_env, para_env_sub, ri_data)
906 ALLOCATE (batch_ranges_at(natom + 1))
907 batch_ranges_at(natom + 1) =
SIZE(ri_data%bsizes_AO_split) + 1
909 DO iblk = 1,
SIZE(ri_data%bsizes_AO_split)
910 IF (idx_to_at_ao(iblk) == iatom + 1)
THEN
912 batch_ranges_at(iatom) = iblk
916 CALL dbt_batched_contract_init(t_3c_work_3(1), batch_range_2=batch_ranges_at)
917 CALL dbt_batched_contract_init(t_3c_work_3(2), batch_range_2=batch_ranges_at)
918 CALL dbt_batched_contract_init(t_3c_work_3(3), batch_range_2=batch_ranges_at)
919 CALL dbt_batched_contract_init(t_3c_work_2(1), batch_range_1=batch_ranges_at)
920 CALL dbt_batched_contract_init(t_3c_work_2(2), batch_range_1=batch_ranges_at)
923 nimg_nze = ri_data%nimg_nze
924 batch_size = ri_data%kp_stack_size
925 n_batch = nimg_nze/batch_size
926 IF (
modulo(nimg_nze, batch_size) .NE. 0) n_batch = n_batch + 1
927 ALLOCATE (batch_ranges_nze(n_batch + 1))
928 DO i_batch = 1, n_batch
929 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
931 batch_ranges_nze(n_batch + 1) = nimg_nze + 1
936 CALL dbt_create(t_2c_inv(iatom), t_2c_binv(iatom))
937 CALL dbt_copy(t_2c_inv(iatom), t_2c_binv(iatom))
938 CALL apply_bump(t_2c_binv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.false.)
939 CALL apply_bump(t_2c_inv(iatom), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
944 ALLOCATE (iapc_pairs(nimg, 2), i_images(nimg))
945 ALLOCATE (force_pattern(natom, natom, nimg))
946 force_pattern(:, :, :) = -1
955 IF (i_loop == 1 .AND. (.NOT. sparsity_pattern(iatom, jatom, b_img) == igroup)) cycle
956 IF (i_loop == 2 .AND. (.NOT. force_pattern(iatom, jatom, b_img) == igroup)) cycle
959 CALL timeset(routinen//
"_2c_1", handle2)
960 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
961 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
962 dbcsr_template=dbcsr_template)
963 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_inv(jatom), &
964 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
965 contract_1=[2], notcontract_1=[1], &
966 contract_2=[1], notcontract_2=[2], &
967 filter_eps=ri_data%filter_eps, flop=nflop)
968 CALL dbt_copy(t_2c_work(2), t_2c_work(5), move_data=.true.)
969 CALL dbt_filter(t_2c_work(5), ri_data%filter_eps)
970 CALL timestop(handle2)
972 CALL timeset(routinen//
"_3c", handle2)
973 bounds_iat(:, 1) = [sum(ri_data%bsizes_AO(1:iatom - 1)) + 1, sum(ri_data%bsizes_AO(1:iatom))]
974 bounds_jat(:, 1) = [sum(ri_data%bsizes_AO(1:jatom - 1)) + 1, sum(ri_data%bsizes_AO(1:jatom))]
975 CALL dbt_clear(t_2c_r_split)
977 DO i_spin = 1, nspins
978 CALL dbt_batched_contract_init(rho_ao_t_sub(i_spin, b_img))
981 CALL get_iapc_pairs(iapc_pairs, b_img, ri_data, qs_env, i_images)
982 DO i_batch = 1, n_batch
986 CALL dbt_clear(t_3c_der_stack(i_xyz))
987 CALL fill_3c_stack(t_3c_der_stack(i_xyz), t_3c_der_ri_sub(:, i_xyz), &
988 iapc_pairs(:, 1), 3, ri_data, filter_at=jatom, &
989 filter_dim=2, idx_to_at=idx_to_at_ao, &
990 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
992 CALL dbt_clear(t_3c_der_stack(3 + i_xyz))
993 CALL fill_3c_stack(t_3c_der_stack(3 + i_xyz), t_3c_der_ao_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)])
999 DO i_spin = 1, nspins
1001 CALL dbt_clear(t_3c_work_2(3))
1002 CALL fill_3c_stack(t_3c_work_2(3), t_3c_apc_sub(i_spin, :), iapc_pairs(:, 2), 3, &
1003 ri_data, filter_at=iatom, filter_dim=1, idx_to_at=idx_to_at_ao, &
1004 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1007 CALL dbt_copy(t_3c_work_2(3), t_3c_work_2(1), move_data=.true.)
1011 CALL dbt_contract(1.0_dp, rho_ao_t_sub(i_spin, b_img), t_3c_work_2(1), &
1012 0.0_dp, t_3c_work_2(2), map_1=[1], map_2=[2, 3], &
1013 contract_1=[1], notcontract_1=[2], &
1014 contract_2=[1], notcontract_2=[2, 3], &
1015 bounds_1=bounds_iat, bounds_2=bounds_jat, &
1016 filter_eps=ri_data%filter_eps, flop=nflop)
1017 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1023 CALL dbt_copy(t_3c_work_2(2), t_3c_work_3(1), order=[2, 1, 3], move_data=.true.)
1024 CALL dbt_batched_contract_init(t_2c_work(5))
1025 CALL dbt_contract(1.0_dp, t_2c_work(5), t_3c_work_3(1), &
1026 0.0_dp, t_3c_work_3(2), map_1=[1], map_2=[2, 3], &
1027 contract_1=[1], notcontract_1=[2], &
1028 contract_2=[1], notcontract_2=[2, 3], &
1029 filter_eps=ri_data%filter_eps, flop=nflop)
1030 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1031 CALL dbt_batched_contract_finalize(t_2c_work(5))
1034 CALL dbt_copy(t_3c_work_3(2), t_3c_work_3(4), move_data=.true.)
1035 IF (use_virial_prv)
THEN
1037 t_3c_der_stack(4:6), atom_of_kind, kind_of, &
1038 idx_to_at_ri, idx_to_at_ao, i_images, &
1039 batch_ranges_nze(i_batch), 2.0_dp*pref, &
1040 ri_data, qs_env, work_virial, cell, particle_set)
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, &
1048 CALL dbt_clear(t_3c_work_3(4))
1052 IF (i_loop == 2) cycle
1055 CALL fill_3c_stack(t_3c_work_3(4), ri_data%kp_t_3c_int, iapc_pairs(:, 1), 3, ri_data, &
1056 filter_at=jatom, filter_dim=2, idx_to_at=idx_to_at_ao, &
1057 img_bounds=[batch_ranges_nze(i_batch), batch_ranges_nze(i_batch + 1)])
1058 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(3), move_data=.true.)
1060 CALL dbt_batched_contract_init(t_2c_r_split)
1061 CALL dbt_contract(1.0_dp, t_3c_work_3(1), t_3c_work_3(3), &
1062 1.0_dp, t_2c_r_split, map_1=[1], map_2=[2], &
1063 contract_1=[2, 3], notcontract_1=[1], &
1064 contract_2=[2, 3], notcontract_2=[1], &
1065 filter_eps=ri_data%filter_eps, flop=nflop)
1066 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1067 CALL dbt_batched_contract_finalize(t_2c_r_split)
1068 CALL dbt_copy(t_3c_work_3(4), t_3c_work_3(1))
1071 DO i_spin = 1, nspins
1072 CALL dbt_batched_contract_finalize(rho_ao_t_sub(i_spin, b_img))
1074 CALL timestop(handle2)
1076 IF (i_loop == 2) cycle
1078 IF (iatom == jatom .AND. b_img == 1) pref = 0.5_dp*pref
1080 CALL timeset(routinen//
"_2c_2", handle2)
1082 CALL dbt_copy(t_2c_r_split, t_2c_r, move_data=.true.)
1084 CALL get_ext_2c_int(t_2c_work(1), mat_2c_pot, iatom, jatom, b_img, ri_data, qs_env, &
1085 blacs_env_ext=blacs_env_sub, para_env_ext=para_env_sub, &
1086 dbcsr_template=dbcsr_template)
1099 CALL get_ext_2c_int(t_2c_der_pot(i_xyz), mat_der_pot_sub(:, i_xyz), iatom, jatom, &
1100 b_img, ri_data, qs_env, blacs_env_ext=blacs_env_sub, &
1101 para_env_ext=para_env_sub, dbcsr_template=dbcsr_template)
1104 IF (use_virial_prv)
THEN
1105 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1106 b_img, pref, ri_data, qs_env, work_virial, cell, particle_set)
1108 CALL get_2c_der_force(force, t_2c_r, t_2c_der_pot, atom_of_kind, kind_of, &
1109 b_img, pref, ri_data, qs_env)
1113 CALL dbt_clear(t_2c_der_pot(i_xyz))
1117 CALL dbt_contract(1.0_dp, t_2c_metric(iatom), t_2c_r, &
1118 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1119 contract_1=[2], notcontract_1=[1], &
1120 contract_2=[1], notcontract_2=[2], &
1121 filter_eps=ri_data%filter_eps, flop=nflop)
1122 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1123 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_work(1), &
1124 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1125 contract_1=[2], notcontract_1=[1], &
1126 contract_2=[2], notcontract_2=[1], &
1127 filter_eps=ri_data%filter_eps, flop=nflop)
1128 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1132 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(iatom), &
1133 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1134 contract_1=[2], notcontract_1=[1], &
1135 contract_2=[1], notcontract_2=[2], &
1136 filter_eps=ri_data%filter_eps, flop=nflop)
1137 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1139 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(3), &
1140 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1141 contract_1=[1], notcontract_1=[2], &
1142 contract_2=[1], notcontract_2=[2], &
1143 filter_eps=ri_data%filter_eps, flop=nflop)
1144 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1146 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1147 CALL get_2c_bump_forces(force, t_2c_work(4), iatom, atom_of_kind, kind_of, pref, &
1148 ri_data, qs_env, work_virial)
1151 CALL dbt_contract(1.0_dp, t_2c_binv(iatom), t_2c_work(2), &
1152 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1153 contract_1=[1], notcontract_1=[2], &
1154 contract_2=[1], notcontract_2=[2], &
1155 filter_eps=ri_data%filter_eps, flop=nflop)
1156 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1158 IF (use_virial_prv)
THEN
1159 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1160 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1161 diag=.true., offdiag=.false.)
1163 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1164 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1168 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1169 CALL apply_bump(t_2c_work(2), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1171 IF (use_virial_prv)
THEN
1172 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1173 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1174 diag=.false., offdiag=.true.)
1176 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(iatom, :), atom_of_kind, &
1177 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1182 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(iatom), &
1183 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1184 contract_1=[2], notcontract_1=[1], &
1185 contract_2=[1], notcontract_2=[2], &
1186 filter_eps=ri_data%filter_eps, flop=nflop)
1187 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1189 CALL dbt_contract(1.0_dp, t_2c_bint(iatom), t_2c_work(4), &
1190 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1191 contract_1=[1], notcontract_1=[2], &
1192 contract_2=[1], notcontract_2=[2], &
1193 filter_eps=ri_data%filter_eps, flop=nflop)
1194 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1196 CALL get_2c_bump_forces(force, t_2c_work(2), iatom, atom_of_kind, kind_of, -pref, &
1197 ri_data, qs_env, work_virial)
1200 CALL dbt_contract(1.0_dp, t_2c_work(1), t_2c_r, &
1201 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1202 contract_1=[1], notcontract_1=[2], &
1203 contract_2=[1], notcontract_2=[2], &
1204 filter_eps=ri_data%filter_eps, flop=nflop)
1205 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1207 CALL dbt_contract(1.0_dp, t_2c_work(2), t_2c_metric(jatom), &
1208 0.0_dp, t_2c_work(3), map_1=[1], map_2=[2], &
1209 contract_1=[2], notcontract_1=[1], &
1210 contract_2=[1], notcontract_2=[2], &
1211 filter_eps=ri_data%filter_eps, flop=nflop)
1212 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1216 CALL dbt_contract(1.0_dp, t_2c_work(3), t_2c_binv(jatom), &
1217 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1218 contract_1=[2], notcontract_1=[1], &
1219 contract_2=[1], notcontract_2=[2], &
1220 filter_eps=ri_data%filter_eps, flop=nflop)
1221 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1223 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(3), &
1224 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1225 contract_1=[1], notcontract_1=[2], &
1226 contract_2=[1], notcontract_2=[2], &
1227 filter_eps=ri_data%filter_eps, flop=nflop)
1228 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1230 CALL dbt_copy(t_2c_work(2), t_2c_work(4), summation=.true.)
1231 CALL get_2c_bump_forces(force, t_2c_work(4), jatom, atom_of_kind, kind_of, pref, &
1232 ri_data, qs_env, work_virial)
1235 CALL dbt_contract(1.0_dp, t_2c_binv(jatom), t_2c_work(2), &
1236 0.0_dp, t_2c_work(4), map_1=[1], map_2=[2], &
1237 contract_1=[1], notcontract_1=[2], &
1238 contract_2=[1], notcontract_2=[2], &
1239 filter_eps=ri_data%filter_eps, flop=nflop)
1240 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1242 IF (use_virial_prv)
THEN
1243 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1244 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1245 diag=.true., offdiag=.false.)
1247 CALL get_2c_der_force(force, t_2c_work(4), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1248 kind_of, 1, -pref, ri_data, qs_env,
diag=.true., offdiag=.false.)
1252 CALL dbt_copy(t_2c_work(4), t_2c_work(2))
1253 CALL apply_bump(t_2c_work(2), jatom, ri_data, qs_env, from_left=.true., from_right=.true.)
1255 IF (use_virial_prv)
THEN
1256 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1257 kind_of, 1, -pref, ri_data, qs_env, work_virial, cell, particle_set, &
1258 diag=.false., offdiag=.true.)
1260 CALL get_2c_der_force(force, t_2c_work(2), t_2c_der_metric_sub(jatom, :), atom_of_kind, &
1261 kind_of, 1, -pref, ri_data, qs_env,
diag=.false., offdiag=.true.)
1266 CALL dbt_contract(1.0_dp, t_2c_work(4), t_2c_bint(jatom), &
1267 0.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1268 contract_1=[2], notcontract_1=[1], &
1269 contract_2=[1], notcontract_2=[2], &
1270 filter_eps=ri_data%filter_eps, flop=nflop)
1271 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1273 CALL dbt_contract(1.0_dp, t_2c_bint(jatom), t_2c_work(4), &
1274 1.0_dp, t_2c_work(2), map_1=[1], map_2=[2], &
1275 contract_1=[1], notcontract_1=[2], &
1276 contract_2=[1], notcontract_2=[2], &
1277 filter_eps=ri_data%filter_eps, flop=nflop)
1278 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1280 CALL get_2c_bump_forces(force, t_2c_work(2), jatom, atom_of_kind, kind_of, -pref, &
1281 ri_data, qs_env, work_virial)
1283 CALL timestop(handle2)
1288 IF (i_loop == 1)
THEN
1289 CALL update_pattern_to_forces(force_pattern, sparsity_pattern, ngroups, ri_data, qs_env)
1293 CALL dbt_batched_contract_finalize(t_3c_work_3(1))
1294 CALL dbt_batched_contract_finalize(t_3c_work_3(2))
1295 CALL dbt_batched_contract_finalize(t_3c_work_3(3))
1296 CALL dbt_batched_contract_finalize(t_3c_work_2(1))
1297 CALL dbt_batched_contract_finalize(t_3c_work_2(2))
1299 IF (use_virial_prv)
THEN
1303 virial%pv_fock_4c(i_xyz, j_xyz) = virial%pv_fock_4c(i_xyz, j_xyz) &
1304 + work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1312 CALL para_env_sub%free()
1313 DEALLOCATE (para_env_sub)
1315 CALL para_env%sync()
1317 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
1320 IF (
ASSOCIATED(dbcsr_template))
THEN
1321 CALL dbcsr_release(dbcsr_template)
1322 DEALLOCATE (dbcsr_template)
1324 CALL dbt_destroy(t_2c_r)
1325 CALL dbt_destroy(t_2c_r_split)
1326 CALL dbt_destroy(t_2c_work(1))
1327 CALL dbt_destroy(t_2c_work(2))
1328 CALL dbt_destroy(t_2c_work(3))
1329 CALL dbt_destroy(t_2c_work(4))
1330 CALL dbt_destroy(t_2c_work(5))
1331 CALL dbt_destroy(t_3c_work_2(1))
1332 CALL dbt_destroy(t_3c_work_2(2))
1333 CALL dbt_destroy(t_3c_work_2(3))
1334 CALL dbt_destroy(t_3c_work_3(1))
1335 CALL dbt_destroy(t_3c_work_3(2))
1336 CALL dbt_destroy(t_3c_work_3(3))
1337 CALL dbt_destroy(t_3c_work_3(4))
1338 CALL dbt_destroy(t_3c_der_stack(1))
1339 CALL dbt_destroy(t_3c_der_stack(2))
1340 CALL dbt_destroy(t_3c_der_stack(3))
1341 CALL dbt_destroy(t_3c_der_stack(4))
1342 CALL dbt_destroy(t_3c_der_stack(5))
1343 CALL dbt_destroy(t_3c_der_stack(6))
1345 CALL dbt_destroy(t_2c_der_pot(i_xyz))
1348 CALL dbt_destroy(t_2c_inv(iatom))
1349 CALL dbt_destroy(t_2c_binv(iatom))
1350 CALL dbt_destroy(t_2c_bint(iatom))
1351 CALL dbt_destroy(t_2c_metric(iatom))
1353 CALL dbt_destroy(t_2c_der_metric_sub(iatom, i_xyz))
1357 CALL dbcsr_release(mat_2c_pot(i_img))
1358 DO i_spin = 1, nspins
1359 CALL dbt_destroy(rho_ao_t_sub(i_spin, i_img))
1360 CALL dbt_destroy(t_3c_apc_sub(i_spin, i_img))
1365 CALL dbt_destroy(t_3c_der_ri_sub(i_img, i_xyz))
1366 CALL dbt_destroy(t_3c_der_ao_sub(i_img, i_xyz))
1367 CALL dbcsr_release(mat_der_pot_sub(i_img, i_xyz))
1371 CALL timestop(handle)
1386 SUBROUTINE apply_bump(t_2c_inout, atom_i, ri_data, qs_env, from_left, from_right, debump)
1387 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_inout
1388 INTEGER,
INTENT(IN) :: atom_i
1391 LOGICAL,
INTENT(IN),
OPTIONAL :: from_left, from_right, debump
1393 INTEGER :: i_img, i_ri, iatom, ind(2), j_img, j_ri, &
1394 jatom, natom, nblks(2), nimg, nkind
1395 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1396 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1397 LOGICAL :: found, my_debump, my_left, my_right
1398 REAL(
dp) :: bval, r0, r1, ri(3), rj(3), rref(3), &
1400 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1402 TYPE(dbt_iterator_type) :: iter
1405 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1407 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1409 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1410 kpoints=kpoints, particle_set=particle_set)
1411 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1414 IF (
PRESENT(debump)) my_debump = debump
1417 IF (
PRESENT(from_left)) my_left = from_left
1420 IF (
PRESENT(from_right)) my_right = from_right
1421 cpassert(my_left .OR. my_right)
1423 CALL dbt_get_info(t_2c_inout, nblks_total=nblks)
1424 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1425 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1430 r1 = ri_data%kp_RI_range
1431 r0 = ri_data%kp_bump_rad
1432 rref =
pbc(particle_set(atom_i)%r, cell)
1437 CALL dbt_iterator_start(iter, t_2c_inout)
1438 DO WHILE (dbt_iterator_blocks_left(iter))
1439 CALL dbt_iterator_next_block(iter, ind)
1440 CALL dbt_get_block(t_2c_inout, ind, blk, found)
1441 IF (.NOT. found) cycle
1443 i_ri = (ind(1) - 1)/natom + 1
1444 i_img = ri_data%RI_cell_to_img(i_ri)
1445 iatom = ind(1) - (i_ri - 1)*natom
1448 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
1450 j_ri = (ind(2) - 1)/natom + 1
1451 j_img = ri_data%RI_cell_to_img(j_ri)
1452 jatom = ind(2) - (j_ri - 1)*natom
1455 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1457 IF (.NOT. my_debump)
THEN
1458 IF (my_left) blk(:, :) = blk(:, :)*bump(norm2(ri - rref), r0, r1)
1459 IF (my_right) blk(:, :) = blk(:, :)*bump(norm2(rj - rref), r0, r1)
1463 bval = bump(norm2(ri - rref), r0, r1)
1464 IF (my_left .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1465 bval = bump(norm2(rj - rref), r0, r1)
1466 IF (my_right .AND. bval > epsilon(1.0_dp)) blk(:, :) = blk(:, :)/bval
1469 CALL dbt_put_block(t_2c_inout, ind, shape(blk), blk)
1473 CALL dbt_iterator_stop(iter)
1475 CALL dbt_filter(t_2c_inout, ri_data%filter_eps)
1477 END SUBROUTINE apply_bump
1491 SUBROUTINE get_2c_bump_forces(force, t_2c_in, atom_i, atom_of_kind, kind_of, pref, ri_data, &
1492 qs_env, work_virial)
1494 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_in
1495 INTEGER,
INTENT(IN) :: atom_i
1496 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
1497 REAL(
dp),
INTENT(IN) :: pref
1500 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT) :: work_virial
1502 INTEGER :: i, i_img, i_ri, i_xyz, iat_of_kind, iatom, ikind, ind(2), j_img, j_ri, j_xyz, &
1503 jat_of_kind, jatom, jkind, natom, nblks(2), nimg, nkind
1504 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1505 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1507 REAL(
dp) :: new_force, r0, r1, ri(3), rj(3), &
1508 rref(3), scoord(3), x
1509 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
1511 TYPE(dbt_iterator_type) :: iter
1514 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1516 NULLIFY (qs_kind_set, particle_set, kpoints, index_to_cell, cell_to_index, cell)
1518 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
1519 kpoints=kpoints, particle_set=particle_set)
1520 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1522 CALL dbt_get_info(t_2c_in, nblks_total=nblks)
1523 cpassert(nblks(1) == ri_data%ncell_RI*natom)
1524 cpassert(nblks(2) == ri_data%ncell_RI*natom)
1529 r1 = ri_data%kp_RI_range
1530 r0 = ri_data%kp_bump_rad
1531 rref =
pbc(particle_set(atom_i)%r, cell)
1533 iat_of_kind = atom_of_kind(atom_i)
1534 ikind = kind_of(atom_i)
1540 CALL dbt_iterator_start(iter, t_2c_in)
1541 DO WHILE (dbt_iterator_blocks_left(iter))
1542 CALL dbt_iterator_next_block(iter, ind)
1543 IF (ind(1) .NE. ind(2)) cycle
1545 CALL dbt_get_block(t_2c_in, ind, blk, found)
1546 IF (.NOT. found) cycle
1549 j_ri = (ind(2) - 1)/natom + 1
1550 j_img = ri_data%RI_cell_to_img(j_ri)
1551 jatom = ind(2) - (j_ri - 1)*natom
1552 jat_of_kind = atom_of_kind(jatom)
1553 jkind = kind_of(jatom)
1556 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
1557 x = norm2(rj - rref)
1558 IF (x < r0 .OR. x > r1) cycle
1561 DO i = 1,
SIZE(blk, 1)
1562 new_force = new_force + blk(i, i)
1564 new_force = pref*new_force*dbump(x, r0, r1)
1570 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) + &
1571 new_force*(rj(i_xyz) - rref(i_xyz))/x
1577 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1578 + new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1583 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) - &
1584 new_force*(rj(i_xyz) - rref(i_xyz))/x
1590 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) &
1591 - new_force*scoord(j_xyz)*(rj(i_xyz) - rref(i_xyz))/x
1597 CALL dbt_iterator_stop(iter)
1600 END SUBROUTINE get_2c_bump_forces
1609 FUNCTION bump(x, r0, r1)
RESULT(b)
1610 REAL(
dp),
INTENT(IN) :: x, r0, r1
1618 r = (x - r0)/(r1 - r0)
1619 b = -6.0_dp*r**5 + 15.0_dp*r**4 - 10.0_dp*r**3 + 1.0_dp
1620 IF (x .GE. r1) b = 0.0_dp
1621 IF (x .LE. r0) b = 1.0_dp
1632 FUNCTION dbump(x, r0, r1)
RESULT(b)
1633 REAL(
dp),
INTENT(IN) :: x, r0, r1
1638 r = (x - r0)/(r1 - r0)
1639 b = (-30.0_dp*r**4 + 60.0_dp*r**3 - 30.0_dp*r**2)/(r1 - r0)
1640 IF (x .GE. r1) b = 0.0_dp
1641 IF (x .LE. r0) b = 0.0_dp
1652 FUNCTION get_apc_index_from_ib(i_index, b_index, qs_env)
RESULT(apc_index)
1653 INTEGER,
INTENT(IN) :: i_index, b_index
1655 INTEGER :: apc_index
1657 INTEGER,
DIMENSION(3) :: cell_apc
1658 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1659 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1663 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1666 cell_apc(:) = index_to_cell(:, i_index) + index_to_cell(:, b_index)
1668 IF (any([cell_apc(1), cell_apc(2), cell_apc(3)] < lbound(cell_to_index)) .OR. &
1669 any([cell_apc(1), cell_apc(2), cell_apc(3)] > ubound(cell_to_index)))
THEN
1673 apc_index = cell_to_index(cell_apc(1), cell_apc(2), cell_apc(3))
1676 END FUNCTION get_apc_index_from_ib
1685 FUNCTION get_apc_index(a_index, c_index, qs_env)
RESULT(i_index)
1686 INTEGER,
INTENT(IN) :: a_index, c_index
1690 INTEGER,
DIMENSION(3) :: cell_i
1691 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1692 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1696 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1698 cell_i(:) = index_to_cell(:, a_index) + index_to_cell(:, c_index)
1700 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1701 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1705 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1708 END FUNCTION get_apc_index
1717 FUNCTION get_i_index(apc_index, b_index, qs_env)
RESULT(i_index)
1718 INTEGER,
INTENT(IN) :: apc_index, b_index
1722 INTEGER,
DIMENSION(3) :: cell_i
1723 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1724 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1728 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1730 cell_i(:) = index_to_cell(:, apc_index) - index_to_cell(:, b_index)
1732 IF (any([cell_i(1), cell_i(2), cell_i(3)] < lbound(cell_to_index)) .OR. &
1733 any([cell_i(1), cell_i(2), cell_i(3)] > ubound(cell_to_index)))
THEN
1737 i_index = cell_to_index(cell_i(1), cell_i(2), cell_i(3))
1740 END FUNCTION get_i_index
1751 SUBROUTINE get_ac_pairs(ac_pairs, apc_index, ri_data, qs_env)
1752 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: ac_pairs
1753 INTEGER,
INTENT(IN) :: apc_index
1757 INTEGER :: a_index, actual_img, c_index, nimg
1759 nimg =
SIZE(ac_pairs, 1)
1764 DO a_index = 1, nimg
1765 actual_img = ri_data%idx_to_img(a_index)
1767 c_index = get_i_index(apc_index, actual_img, qs_env)
1768 ac_pairs(a_index, 1) = a_index
1769 ac_pairs(a_index, 2) = c_index
1773 END SUBROUTINE get_ac_pairs
1785 SUBROUTINE get_iapc_pairs(iapc_pairs, b_index, ri_data, qs_env, actual_i_img)
1786 INTEGER,
DIMENSION(:, :),
INTENT(INOUT) :: iapc_pairs
1787 INTEGER,
INTENT(IN) :: b_index
1790 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: actual_i_img
1792 INTEGER :: actual_img, apc_index, i_index, nimg
1794 nimg =
SIZE(iapc_pairs, 1)
1795 IF (
PRESENT(actual_i_img)) actual_i_img(:) = 0
1797 iapc_pairs(:, :) = 0
1800 DO i_index = 1, nimg
1801 actual_img = ri_data%idx_to_img(i_index)
1802 apc_index = get_apc_index_from_ib(actual_img, b_index, qs_env)
1803 IF (apc_index == 0) cycle
1804 iapc_pairs(i_index, 1) = i_index
1805 iapc_pairs(i_index, 2) = apc_index
1806 IF (
PRESENT(actual_i_img)) actual_i_img(i_index) = actual_img
1809 END SUBROUTINE get_iapc_pairs
1818 FUNCTION get_opp_index(a_index, qs_env)
RESULT(opp_index)
1819 INTEGER,
INTENT(IN) :: a_index
1821 INTEGER :: opp_index
1823 INTEGER,
DIMENSION(3) :: opp_cell
1824 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1825 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1828 NULLIFY (kpoints, cell_to_index, index_to_cell)
1831 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
1833 opp_cell(:) = -index_to_cell(:, a_index)
1835 IF (any([opp_cell(1), opp_cell(2), opp_cell(3)] < lbound(cell_to_index)) .OR. &
1836 any([opp_cell(1), opp_cell(2), opp_cell(3)] > ubound(cell_to_index)))
THEN
1840 opp_index = cell_to_index(opp_cell(1), opp_cell(2), opp_cell(3))
1843 END FUNCTION get_opp_index
1854 SUBROUTINE get_pmat_images(rho_ao_t, rho_ao, scale_prev_p, ri_data, qs_env)
1855 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t
1856 TYPE(dbcsr_p_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao
1857 REAL(
dp),
INTENT(IN) :: scale_prev_p
1861 INTEGER :: cell_j(3), i_img, i_spin, iatom, icol, &
1862 irow, j_img, jatom, mi_img, mj_img, &
1864 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1867 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock, pblock_desymm
1868 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks, rho_desymm
1869 TYPE(dbt_type) :: tmp
1873 DIMENSION(:),
POINTER :: nl_iterator
1875 POINTER :: sab_nl, sab_nl_nosym
1878 NULLIFY (rho_desymm, kpoints, sab_nl_nosym, scf_env, matrix_ks, dft_control, &
1879 sab_nl, nl_iterator, cell_to_index, pblock, pblock_desymm)
1881 CALL get_qs_env(qs_env, kpoints=kpoints, scf_env=scf_env, matrix_ks_kp=matrix_ks, dft_control=dft_control)
1882 CALL get_kpoint_info(kpoints, sab_nl_nosym=sab_nl_nosym, cell_to_index=cell_to_index, sab_nl=sab_nl)
1884 IF (dft_control%do_admm)
THEN
1885 CALL get_admm_env(qs_env%admm_env, matrix_ks_aux_fit_kp=matrix_ks)
1888 nspins =
SIZE(matrix_ks, 1)
1891 ALLOCATE (rho_desymm(nspins, nimg))
1893 DO i_spin = 1, nspins
1894 ALLOCATE (rho_desymm(i_spin, i_img)%matrix)
1895 CALL dbcsr_create(rho_desymm(i_spin, i_img)%matrix, template=matrix_ks(i_spin, i_img)%matrix, &
1896 matrix_type=dbcsr_type_no_symmetry)
1900 CALL dbt_create(rho_desymm(1, 1)%matrix, tmp)
1907 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
1908 IF (j_img > nimg .OR. j_img < 1) cycle
1911 IF (iatom == jatom)
fac = 0.5_dp
1912 mj_img = get_opp_index(j_img, qs_env)
1914 IF (mj_img == 0)
fac = 1.0_dp
1918 IF (iatom > jatom)
THEN
1924 DO i_spin = 1, nspins
1925 CALL dbcsr_get_block_p(rho_ao(i_spin, j_img)%matrix, irow, icol, pblock, found)
1926 IF (.NOT. found) cycle
1929 CALL dbcsr_get_block_p(rho_desymm(i_spin, j_img)%matrix, iatom, jatom, pblock_desymm, found)
1930 IF (.NOT. found) cycle
1932 IF (iatom > jatom)
THEN
1933 pblock_desymm(:, :) =
fac*transpose(pblock(:, :))
1935 pblock_desymm(:, :) =
fac*pblock(:, :)
1942 DO i_spin = 1, nspins
1943 CALL dbt_scale(rho_ao_t(i_spin, i_img), scale_prev_p)
1945 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, i_img)%matrix, tmp)
1946 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), summation=.true., move_data=.true.)
1949 mi_img = get_opp_index(i_img, qs_env)
1950 IF (mi_img > 0 .AND. mi_img .LE. nimg)
THEN
1951 CALL dbt_copy_matrix_to_tensor(rho_desymm(i_spin, mi_img)%matrix, tmp)
1952 CALL dbt_copy(tmp, rho_ao_t(i_spin, i_img), order=[2, 1], summation=.true., move_data=.true.)
1954 CALL dbt_filter(rho_ao_t(i_spin, i_img), ri_data%filter_eps)
1959 DO i_spin = 1, nspins
1960 CALL dbcsr_release(rho_desymm(i_spin, i_img)%matrix)
1961 DEALLOCATE (rho_desymm(i_spin, i_img)%matrix)
1965 CALL dbt_destroy(tmp)
1966 DEALLOCATE (rho_desymm)
1968 END SUBROUTINE get_pmat_images
1987 SUBROUTINE get_ext_2c_int(t_2c_pot, mat_orig, atom_i, atom_j, img_b, ri_data, qs_env, do_inverse, &
1988 para_env_ext, blacs_env_ext, dbcsr_template, off_diagonal, skip_inverse)
1989 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_pot
1990 TYPE(dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_orig
1991 INTEGER,
INTENT(IN) :: atom_i, atom_j, img_b
1994 LOGICAL,
INTENT(IN),
OPTIONAL :: do_inverse
1997 TYPE(dbcsr_type),
OPTIONAL,
POINTER :: dbcsr_template
1998 LOGICAL,
INTENT(IN),
OPTIONAL :: off_diagonal, skip_inverse
2000 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_ext_2c_int'
2002 INTEGER :: blk, group, handle, handle2, i_img, i_ri, iatom, iblk, ikind, img_tot, j_img, &
2003 j_ri, jatom, jblk, jkind, n_dependent, natom, nblks_ri, nimg, nkind
2004 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2
2005 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: present_atoms_i, present_atoms_j
2006 INTEGER,
DIMENSION(3) :: cell_b, cell_i, cell_j, cell_tot
2007 INTEGER,
DIMENSION(:),
POINTER :: col_dist, col_dist_ext, ri_blk_size_ext, &
2008 row_dist, row_dist_ext
2009 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell, pgrid
2010 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
2011 LOGICAL :: do_inverse_prv, found, my_offd, &
2012 skip_inverse_prv, use_template
2013 REAL(
dp) :: bfac, dij, r0, r1, threshold
2014 REAL(
dp),
DIMENSION(3) :: ri, rij, rj, rref, scoord
2015 REAL(
dp),
DIMENSION(:, :),
POINTER :: pblock
2018 TYPE(dbcsr_distribution_type) :: dbcsr_dist, dbcsr_dist_ext
2019 TYPE(dbcsr_iterator_type) :: dbcsr_iter
2020 TYPE(dbcsr_type) :: work, work_tight, work_tight_inv
2021 TYPE(dbt_type) :: t_2c_tmp
2024 DIMENSION(:),
TARGET :: basis_set_ri
2028 DIMENSION(:),
POINTER :: nl_iterator
2032 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2034 NULLIFY (qs_kind_set, nl_2c, nl_iterator, cell, kpoints, cell_to_index, index_to_cell, dist_2d, &
2035 para_env, pblock, blacs_env, particle_set, col_dist, row_dist, pgrid, &
2036 col_dist_ext, row_dist_ext)
2038 CALL timeset(routinen, handle)
2043 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, cell=cell, &
2044 kpoints=kpoints, para_env=para_env, blacs_env=blacs_env, particle_set=particle_set)
2045 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
2047 do_inverse_prv = .false.
2048 IF (
PRESENT(do_inverse)) do_inverse_prv = do_inverse
2049 IF (do_inverse_prv)
THEN
2050 cpassert(atom_i == atom_j)
2053 skip_inverse_prv = .false.
2054 IF (
PRESENT(skip_inverse)) skip_inverse_prv = skip_inverse
2057 IF (
PRESENT(off_diagonal)) my_offd = off_diagonal
2059 IF (
PRESENT(para_env_ext)) para_env => para_env_ext
2060 IF (
PRESENT(blacs_env_ext)) blacs_env => blacs_env_ext
2062 nimg =
SIZE(mat_orig)
2064 CALL timeset(routinen//
"_nl_iter", handle2)
2067 ALLOCATE (dist1(natom), dist2(natom))
2069 dist1(iatom) = mod(iatom, blacs_env%num_pe(1))
2070 dist2(iatom) = mod(iatom, blacs_env%num_pe(2))
2074 ALLOCATE (basis_set_ri(nkind))
2078 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
2080 ALLOCATE (present_atoms_i(natom, nimg), present_atoms_j(natom, nimg))
2086 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom, r=rij, cell=cell_j, &
2087 ikind=ikind, jkind=jkind)
2091 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
2092 IF (j_img > nimg .OR. j_img < 1) cycle
2094 IF (iatom == atom_i .AND. dij .LE. ri_data%kp_RI_range) present_atoms_i(jatom, j_img) = 1
2095 IF (iatom == atom_j .AND. dij .LE. ri_data%kp_RI_range) present_atoms_j(jatom, j_img) = 1
2100 CALL timestop(handle2)
2102 CALL para_env%sum(present_atoms_i)
2103 CALL para_env%sum(present_atoms_j)
2107 use_template = .false.
2108 IF (
PRESENT(dbcsr_template))
THEN
2109 IF (
ASSOCIATED(dbcsr_template)) use_template = .true.
2112 IF (use_template)
THEN
2113 CALL dbcsr_create(work, template=dbcsr_template)
2115 CALL dbcsr_get_info(mat_orig(1), distribution=dbcsr_dist)
2116 CALL dbcsr_distribution_get(dbcsr_dist, row_dist=row_dist, col_dist=col_dist, group=group, pgrid=pgrid)
2117 ALLOCATE (row_dist_ext(ri_data%ncell_RI*natom), col_dist_ext(ri_data%ncell_RI*natom))
2118 ALLOCATE (ri_blk_size_ext(ri_data%ncell_RI*natom))
2119 DO i_ri = 1, ri_data%ncell_RI
2120 row_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = row_dist(:)
2121 col_dist_ext((i_ri - 1)*natom + 1:i_ri*natom) = col_dist(:)
2122 ri_blk_size_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2125 CALL dbcsr_distribution_new(dbcsr_dist_ext, group=group, pgrid=pgrid, &
2126 row_dist=row_dist_ext, col_dist=col_dist_ext)
2127 CALL dbcsr_create(work, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2128 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2129 CALL dbcsr_distribution_release(dbcsr_dist_ext)
2130 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2132 IF (
PRESENT(dbcsr_template))
THEN
2133 ALLOCATE (dbcsr_template)
2134 CALL dbcsr_create(dbcsr_template, template=work)
2138 cell_b(:) = index_to_cell(:, img_b)
2140 i_ri = ri_data%img_to_RI_cell(i_img)
2141 IF (i_ri == 0) cycle
2142 cell_i(:) = index_to_cell(:, i_img)
2144 j_ri = ri_data%img_to_RI_cell(j_img)
2145 IF (j_ri == 0) cycle
2146 cell_j(:) = index_to_cell(:, j_img)
2147 cell_tot = cell_j - cell_i + cell_b
2149 IF (any([cell_tot(1), cell_tot(2), cell_tot(3)] < lbound(cell_to_index)) .OR. &
2150 any([cell_tot(1), cell_tot(2), cell_tot(3)] > ubound(cell_to_index))) cycle
2151 img_tot = cell_to_index(cell_tot(1), cell_tot(2), cell_tot(3))
2152 IF (img_tot > nimg .OR. img_tot < 1) cycle
2154 CALL dbcsr_iterator_start(dbcsr_iter, mat_orig(img_tot))
2155 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
2156 CALL dbcsr_iterator_next_block(dbcsr_iter, row=iatom, column=jatom, blk=blk)
2157 IF (present_atoms_i(iatom, i_img) == 0) cycle
2158 IF (present_atoms_j(jatom, j_img) == 0) cycle
2159 IF (my_offd .AND. (i_ri - 1)*natom + iatom == (j_ri - 1)*natom + jatom) cycle
2161 CALL dbcsr_get_block_p(mat_orig(img_tot), iatom, jatom, pblock, found)
2162 IF (.NOT. found) cycle
2164 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2167 CALL dbcsr_iterator_stop(dbcsr_iter)
2171 CALL dbcsr_finalize(work)
2173 IF (do_inverse_prv)
THEN
2175 r1 = ri_data%kp_RI_range
2176 r0 = ri_data%kp_bump_rad
2179 nblks_ri = sum(present_atoms_i)
2180 ALLOCATE (col_dist_ext(nblks_ri), row_dist_ext(nblks_ri), ri_blk_size_ext(nblks_ri))
2183 i_ri = ri_data%img_to_RI_cell(i_img)
2184 IF (i_ri == 0) cycle
2186 IF (present_atoms_i(iatom, i_img) == 0) cycle
2188 col_dist_ext(iblk) = col_dist(iatom)
2189 row_dist_ext(iblk) = row_dist(iatom)
2190 ri_blk_size_ext(iblk) = ri_data%bsizes_RI(iatom)
2194 CALL dbcsr_distribution_new(dbcsr_dist_ext, group=group, pgrid=pgrid, &
2195 row_dist=row_dist_ext, col_dist=col_dist_ext)
2196 CALL dbcsr_create(work_tight, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2197 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2198 CALL dbcsr_create(work_tight_inv, dist=dbcsr_dist_ext, name=
"RI_ext", matrix_type=dbcsr_type_no_symmetry, &
2199 row_blk_size=ri_blk_size_ext, col_blk_size=ri_blk_size_ext)
2200 CALL dbcsr_distribution_release(dbcsr_dist_ext)
2201 DEALLOCATE (col_dist_ext, row_dist_ext, ri_blk_size_ext)
2205 rref =
pbc(particle_set(atom_i)%r, cell)
2209 i_ri = ri_data%img_to_RI_cell(i_img)
2210 IF (i_ri == 0) cycle
2212 IF (present_atoms_i(iatom, i_img) == 0) cycle
2216 CALL scaled_to_real(ri, scoord(:) + index_to_cell(:, i_img), cell)
2220 j_ri = ri_data%img_to_RI_cell(j_img)
2221 IF (j_ri == 0) cycle
2223 IF (present_atoms_j(jatom, j_img) == 0) cycle
2227 CALL scaled_to_real(rj, scoord(:) + index_to_cell(:, j_img), cell)
2229 CALL dbcsr_get_block_p(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock, found)
2230 IF (.NOT. found) cycle
2233 IF (iblk .NE. jblk) bfac = bump(norm2(ri - rref), r0, r1)*bump(norm2(rj - rref), r0, r1)
2234 CALL dbcsr_put_block(work_tight, iblk, jblk, bfac*pblock(:, :))
2239 CALL dbcsr_finalize(work_tight)
2240 CALL dbcsr_clear(work)
2242 IF (.NOT. skip_inverse_prv)
THEN
2243 SELECT CASE (ri_data%t2c_method)
2245 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
2246 CALL invert_hotelling(work_tight_inv, work_tight, threshold=threshold, silent=.false.)
2248 CALL dbcsr_copy(work_tight_inv, work_tight)
2251 upper_to_full=.true.)
2253 CALL dbcsr_copy(work_tight_inv, work_tight)
2254 CALL cp_dbcsr_power(work_tight_inv, -1.0_dp, ri_data%eps_eigval, n_dependent, &
2255 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
2258 CALL dbcsr_copy(work_tight_inv, work_tight)
2266 i_ri = ri_data%img_to_RI_cell(i_img)
2267 IF (i_ri == 0) cycle
2269 IF (present_atoms_i(iatom, i_img) == 0) cycle
2274 j_ri = ri_data%img_to_RI_cell(j_img)
2275 IF (j_ri == 0) cycle
2277 IF (present_atoms_j(jatom, j_img) == 0) cycle
2280 CALL dbcsr_get_block_p(work_tight_inv, iblk, jblk, pblock, found)
2281 IF (.NOT. found) cycle
2283 CALL dbcsr_put_block(work, (i_ri - 1)*natom + iatom, (j_ri - 1)*natom + jatom, pblock)
2288 CALL dbcsr_finalize(work)
2290 CALL dbcsr_release(work_tight)
2291 CALL dbcsr_release(work_tight_inv)
2294 CALL dbt_create(work, t_2c_tmp)
2295 CALL dbt_copy_matrix_to_tensor(work, t_2c_tmp)
2296 CALL dbt_copy(t_2c_tmp, t_2c_pot, move_data=.true.)
2297 CALL dbt_filter(t_2c_pot, ri_data%filter_eps)
2299 CALL dbt_destroy(t_2c_tmp)
2300 CALL dbcsr_release(work)
2302 CALL timestop(handle)
2304 END SUBROUTINE get_ext_2c_int
2314 SUBROUTINE contract_pmat_3c(t_3c_apc, rho_ao_t, ri_data, qs_env)
2315 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, rho_ao_t
2319 CHARACTER(len=*),
PARAMETER :: routinen =
'contract_pmat_3c'
2321 INTEGER :: apc_img, batch_size, handle, i_batch, &
2322 i_img, i_spin, j_batch, n_batch_img, &
2323 n_batch_nze, nimg, nimg_nze, nspins
2324 INTEGER(int_8) :: nflop, nze
2325 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges_img, batch_ranges_nze, &
2327 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: ac_pairs
2328 REAL(
dp) :: occ, t1, t2
2329 TYPE(dbt_type) :: t_3c_tmp
2330 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: ints_stack, res_stack, rho_stack
2333 CALL timeset(routinen, handle)
2335 CALL get_qs_env(qs_env, dft_control=dft_control)
2338 nimg_nze = ri_data%nimg_nze
2339 nspins = dft_control%nspins
2341 CALL dbt_create(t_3c_apc(1, 1), t_3c_tmp)
2343 batch_size = nimg/ri_data%n_mem
2346 n_batch_img = nimg/batch_size
2347 IF (
modulo(nimg, batch_size) .NE. 0) n_batch_img = n_batch_img + 1
2348 ALLOCATE (batch_ranges_img(n_batch_img + 1))
2349 DO i_batch = 1, n_batch_img
2350 batch_ranges_img(i_batch) = (i_batch - 1)*batch_size + 1
2352 batch_ranges_img(n_batch_img + 1) = nimg + 1
2355 n_batch_nze = nimg_nze/batch_size
2356 IF (
modulo(nimg_nze, batch_size) .NE. 0) n_batch_nze = n_batch_nze + 1
2357 ALLOCATE (batch_ranges_nze(n_batch_nze + 1))
2358 DO i_batch = 1, n_batch_nze
2359 batch_ranges_nze(i_batch) = (i_batch - 1)*batch_size + 1
2361 batch_ranges_nze(n_batch_nze + 1) = nimg_nze + 1
2364 ALLOCATE (rho_stack(2), ints_stack(2), res_stack(2))
2365 CALL get_stack_tensors(res_stack, rho_stack, ints_stack, rho_ao_t(1, 1), &
2366 ri_data%t_3c_int_ctr_1(1, 1), batch_size, ri_data, qs_env)
2368 ALLOCATE (ac_pairs(nimg, 2), int_indices(nimg_nze))
2369 DO i_img = 1, nimg_nze
2370 int_indices(i_img) = i_img
2374 DO j_batch = 1, n_batch_nze
2376 CALL fill_3c_stack(ints_stack(1), ri_data%t_3c_int_ctr_1(1, :), int_indices, 3, ri_data, &
2377 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)])
2378 CALL dbt_copy(ints_stack(1), ints_stack(2), move_data=.true.)
2380 DO i_spin = 1, nspins
2381 DO i_batch = 1, n_batch_img
2383 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2384 CALL get_ac_pairs(ac_pairs, apc_img, ri_data, qs_env)
2385 CALL fill_2c_stack(rho_stack(1), rho_ao_t(i_spin, :), ac_pairs(:, 2), 1, ri_data, &
2386 img_bounds=[batch_ranges_nze(j_batch), batch_ranges_nze(j_batch + 1)], &
2387 shift=apc_img - batch_ranges_img(i_batch) + 1)
2392 CALL dbt_copy(rho_stack(1), rho_stack(2), move_data=.true.)
2395 CALL dbt_batched_contract_init(rho_stack(2))
2396 CALL dbt_contract(1.0_dp, ints_stack(2), rho_stack(2), &
2397 0.0_dp, res_stack(2), map_1=[1, 2], map_2=[3], &
2398 contract_1=[3], notcontract_1=[1, 2], &
2399 contract_2=[1], notcontract_2=[2], &
2400 filter_eps=ri_data%filter_eps, flop=nflop)
2401 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2402 CALL dbt_batched_contract_finalize(rho_stack(2))
2403 CALL dbt_copy(res_stack(2), res_stack(1), move_data=.true.)
2405 DO apc_img = batch_ranges_img(i_batch), batch_ranges_img(i_batch + 1) - 1
2407 CALL unstack_t_3c_apc(t_3c_tmp, res_stack(1), apc_img - batch_ranges_img(i_batch) + 1)
2408 CALL dbt_copy(t_3c_tmp, t_3c_apc(i_spin, apc_img), summation=.true., move_data=.true.)
2414 DEALLOCATE (batch_ranges_img)
2415 DEALLOCATE (batch_ranges_nze)
2417 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
2419 CALL dbt_destroy(rho_stack(1))
2420 CALL dbt_destroy(rho_stack(2))
2421 CALL dbt_destroy(ints_stack(1))
2422 CALL dbt_destroy(ints_stack(2))
2423 CALL dbt_destroy(res_stack(1))
2424 CALL dbt_destroy(res_stack(2))
2425 CALL dbt_destroy(t_3c_tmp)
2427 CALL timestop(handle)
2429 END SUBROUTINE contract_pmat_3c
2437 SUBROUTINE precontract_3c_ints(t_3c_int, ri_data, qs_env)
2438 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_int
2442 CHARACTER(len=*),
PARAMETER :: routinen =
'precontract_3c_ints'
2444 INTEGER :: batch_size, handle, i_batch, i_img, &
2445 i_ri, iatom, is, n_batch, natom, &
2446 nblks, nblks_3c(3), nimg
2447 INTEGER(int_8) :: nflop
2448 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: batch_ranges, bsizes_ri_ext, bsizes_ri_ext_split, &
2449 bsizes_stack, dist1, dist2, dist3, dist_stack3, idx_to_at_ao, int_indices
2450 TYPE(dbt_distribution_type) :: t_dist
2451 TYPE(dbt_type) :: t_2c_ri_tmp(2), t_3c_tmp(3)
2453 CALL timeset(routinen, handle)
2458 ALLOCATE (int_indices(nimg))
2460 int_indices(i_img) = i_img
2463 ALLOCATE (idx_to_at_ao(
SIZE(ri_data%bsizes_AO_split)))
2464 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
2466 nblks =
SIZE(ri_data%bsizes_RI_split)
2467 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
2468 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
2469 DO i_ri = 1, ri_data%ncell_RI
2470 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
2471 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
2474 bsizes_ri_ext, bsizes_ri_ext, &
2476 DEALLOCATE (dist1, dist2)
2478 bsizes_ri_ext_split, bsizes_ri_ext_split, &
2480 DEALLOCATE (dist1, dist2)
2483 batch_size = nimg/ri_data%n_mem
2484 n_batch = nimg/batch_size
2485 IF (
modulo(nimg, batch_size) .NE. 0) n_batch = n_batch + 1
2486 ALLOCATE (batch_ranges(n_batch + 1))
2487 DO i_batch = 1, n_batch
2488 batch_ranges(i_batch) = (i_batch - 1)*batch_size + 1
2490 batch_ranges(n_batch + 1) = nimg + 1
2492 nblks =
SIZE(ri_data%bsizes_AO_split)
2493 ALLOCATE (bsizes_stack(batch_size*nblks))
2494 DO is = 1, batch_size
2495 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
2498 CALL dbt_get_info(t_3c_int(1, 1), nblks_total=nblks_3c)
2499 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)), dist_stack3(batch_size*nblks_3c(3)))
2500 CALL dbt_get_info(t_3c_int(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, proc_dist_3=dist3)
2501 DO is = 1, batch_size
2502 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
2505 CALL dbt_distribution_new(t_dist, ri_data%pgrid, dist1, dist2, dist_stack3)
2506 CALL dbt_create(t_3c_tmp(1),
"ints_stack", t_dist, [1], [2, 3], bsizes_ri_ext_split, &
2507 ri_data%bsizes_AO_split, bsizes_stack)
2508 CALL dbt_distribution_destroy(t_dist)
2509 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
2511 CALL dbt_create(t_3c_tmp(1), t_3c_tmp(2))
2512 CALL dbt_create(t_3c_int(1, 1), t_3c_tmp(3))
2515 CALL dbt_copy(ri_data%t_2c_inv(1, iatom), t_2c_ri_tmp(1))
2516 CALL apply_bump(t_2c_ri_tmp(1), iatom, ri_data, qs_env, from_left=.true., from_right=.true.)
2517 CALL dbt_copy(t_2c_ri_tmp(1), t_2c_ri_tmp(2), move_data=.true.)
2519 CALL dbt_batched_contract_init(t_2c_ri_tmp(2))
2520 DO i_batch = 1, n_batch
2522 CALL fill_3c_stack(t_3c_tmp(1), t_3c_int(1, :), int_indices, 3, ri_data, &
2523 img_bounds=[batch_ranges(i_batch), batch_ranges(i_batch + 1)], &
2524 filter_at=iatom, filter_dim=2, idx_to_at=idx_to_at_ao)
2526 CALL dbt_contract(1.0_dp, t_2c_ri_tmp(2), t_3c_tmp(1), &
2527 0.0_dp, t_3c_tmp(2), map_1=[1], map_2=[2, 3], &
2528 contract_1=[2], notcontract_1=[1], &
2529 contract_2=[1], notcontract_2=[2, 3], &
2530 filter_eps=ri_data%filter_eps, flop=nflop)
2531 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2533 DO i_img = batch_ranges(i_batch), batch_ranges(i_batch + 1) - 1
2534 CALL unstack_t_3c_apc(t_3c_tmp(3), t_3c_tmp(2), i_img - batch_ranges(i_batch) + 1)
2535 CALL dbt_copy(t_3c_tmp(3), ri_data%t_3c_int_ctr_1(1, i_img), summation=.true., &
2536 order=[2, 1, 3], move_data=.true.)
2538 CALL dbt_clear(t_3c_tmp(1))
2540 CALL dbt_batched_contract_finalize(t_2c_ri_tmp(2))
2543 CALL dbt_destroy(t_2c_ri_tmp(1))
2544 CALL dbt_destroy(t_2c_ri_tmp(2))
2545 CALL dbt_destroy(t_3c_tmp(1))
2546 CALL dbt_destroy(t_3c_tmp(2))
2547 CALL dbt_destroy(t_3c_tmp(3))
2550 CALL dbt_destroy(t_3c_int(1, i_img))
2553 CALL timestop(handle)
2555 END SUBROUTINE precontract_3c_ints
2566 SUBROUTINE copy_2c_to_subgroup(t2c_sub, t2c_main, group_size, ngroups, para_env)
2567 TYPE(dbt_type),
INTENT(INOUT) :: t2c_sub, t2c_main
2568 INTEGER,
INTENT(IN) :: group_size, ngroups
2571 INTEGER :: batch_size, i, i_batch, i_msg, iblk, &
2572 igroup, iproc, ir, is, jblk, n_batch, &
2574 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2
2575 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: block_dest, block_source
2576 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: current_dest
2577 INTEGER,
DIMENSION(2) :: ind, nblks
2579 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
2580 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2581 TYPE(dbt_iterator_type) :: iter
2582 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2588 CALL dbt_get_info(t2c_main, nblks_total=nblks)
2591 ALLOCATE (block_source(nblks(1), nblks(2)))
2595 CALL dbt_iterator_start(iter, t2c_main)
2596 DO WHILE (dbt_iterator_blocks_left(iter))
2597 CALL dbt_iterator_next_block(iter, ind)
2598 CALL dbt_get_block(t2c_main, ind, blk, found)
2599 IF (.NOT. found) cycle
2601 block_source(ind(1), ind(2)) = para_env%mepos
2606 CALL dbt_iterator_stop(iter)
2609 CALL para_env%sum(nocc)
2610 CALL para_env%sum(block_source)
2611 block_source = block_source + para_env%num_pe - 1
2612 IF (nocc == 0)
RETURN
2615 igroup = para_env%mepos/group_size
2616 ALLOCATE (block_dest(nblks(1), nblks(2)))
2618 DO jblk = 1, nblks(2)
2619 DO iblk = 1, nblks(1)
2620 IF (block_source(iblk, jblk) == -1) cycle
2622 CALL dbt_get_stored_coordinates(t2c_sub, [iblk, jblk], iproc)
2623 block_dest(iblk, jblk) = igroup*group_size + iproc
2627 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)))
2628 CALL dbt_get_info(t2c_main, blk_size_1=bsizes1, blk_size_2=bsizes2)
2630 ALLOCATE (current_dest(nblks(1), nblks(2), 0:ngroups - 1))
2631 DO igroup = 0, ngroups - 1
2633 current_dest(:, :, igroup) = block_dest(:, :)
2634 CALL para_env%bcast(current_dest(:, :, igroup), source=igroup*group_size)
2638 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2639 n_batch = (nocc*ngroups)/batch_size
2640 IF (
modulo(nocc*ngroups, batch_size) .NE. 0) n_batch = n_batch + 1
2642 DO i_batch = 1, n_batch
2644 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2645 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2649 DO jblk = 1, nblks(2)
2650 DO iblk = 1, nblks(1)
2651 DO igroup = 0, ngroups - 1
2652 IF (block_source(iblk, jblk) == -1) cycle
2655 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2658 tag = i_msg - (i_batch - 1)*batch_size
2661 IF (para_env%mepos == block_source(iblk, jblk))
THEN
2662 CALL dbt_get_block(t2c_main, [iblk, jblk], blk, found)
2666 IF (block_source(iblk, jblk) == current_dest(iblk, jblk, igroup))
THEN
2667 IF (found)
CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2669 IF (para_env%mepos == block_source(iblk, jblk) .AND. found)
THEN
2670 ALLOCATE (send_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2671 send_buff(tag)%array(:, :) = blk(:, :)
2673 CALL para_env%isend(msgin=send_buff(tag)%array, dest=current_dest(iblk, jblk, igroup), &
2674 request=send_req(is), tag=tag)
2677 IF (para_env%mepos == current_dest(iblk, jblk, igroup))
THEN
2678 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk)))
2680 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk), &
2681 request=recv_req(ir), tag=tag)
2685 IF (found)
DEALLOCATE (blk)
2693 DO i = 1, batch_size
2694 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
2699 DO jblk = 1, nblks(2)
2700 DO iblk = 1, nblks(1)
2701 DO igroup = 0, ngroups - 1
2702 IF (block_source(iblk, jblk) == -1) cycle
2705 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2708 tag = i_msg - (i_batch - 1)*batch_size
2710 IF (para_env%mepos == current_dest(iblk, jblk, igroup) .AND. &
2711 block_source(iblk, jblk) .NE. current_dest(iblk, jblk, igroup))
THEN
2713 ALLOCATE (blk(bsizes1(iblk), bsizes2(jblk)))
2714 blk(:, :) = recv_buff(tag)%array(:, :)
2715 CALL dbt_put_block(t2c_sub, [iblk, jblk], shape(blk), blk)
2723 DO i = 1, batch_size
2724 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
2726 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
2728 CALL dbt_finalize(t2c_sub)
2730 END SUBROUTINE copy_2c_to_subgroup
2744 SUBROUTINE copy_3c_to_subgroup(t3c_sub, t3c_main, group_size, ngroups, para_env, iatom_to_subgroup, &
2746 TYPE(dbt_type),
INTENT(INOUT) :: t3c_sub, t3c_main
2747 INTEGER,
INTENT(IN) :: group_size, ngroups
2750 INTENT(INOUT),
OPTIONAL :: iatom_to_subgroup
2751 INTEGER,
INTENT(IN),
OPTIONAL :: dim_at
2752 INTEGER,
DIMENSION(:),
OPTIONAL :: idx_to_at
2754 INTEGER :: batch_size, i, i_batch, i_msg, iatom, &
2755 iblk, igroup, iproc, ir, is, jblk, &
2756 kblk, n_batch, nocc, tag
2757 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes1, bsizes2, bsizes3
2758 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: block_dest, block_source
2759 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :, :) :: current_dest
2760 INTEGER,
DIMENSION(3) :: ind, nblks
2761 LOGICAL :: filter_at, found
2762 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
2763 TYPE(
cp_3d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2764 TYPE(dbt_iterator_type) :: iter
2765 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2771 CALL dbt_get_info(t3c_main, nblks_total=nblks)
2775 IF (
PRESENT(iatom_to_subgroup) .AND.
PRESENT(dim_at) .AND.
PRESENT(idx_to_at))
THEN
2777 cpassert(nblks(dim_at) ==
SIZE(idx_to_at))
2781 ALLOCATE (block_source(nblks(1), nblks(2), nblks(3)))
2785 CALL dbt_iterator_start(iter, t3c_main)
2786 DO WHILE (dbt_iterator_blocks_left(iter))
2787 CALL dbt_iterator_next_block(iter, ind)
2788 CALL dbt_get_block(t3c_main, ind, blk, found)
2789 IF (.NOT. found) cycle
2791 block_source(ind(1), ind(2), ind(3)) = para_env%mepos
2796 CALL dbt_iterator_stop(iter)
2799 CALL para_env%sum(nocc)
2800 CALL para_env%sum(block_source)
2801 block_source = block_source + para_env%num_pe - 1
2802 IF (nocc == 0)
RETURN
2805 igroup = para_env%mepos/group_size
2806 ALLOCATE (block_dest(nblks(1), nblks(2), nblks(3)))
2808 DO kblk = 1, nblks(3)
2809 DO jblk = 1, nblks(2)
2810 DO iblk = 1, nblks(1)
2811 IF (block_source(iblk, jblk, kblk) == -1) cycle
2813 CALL dbt_get_stored_coordinates(t3c_sub, [iblk, jblk, kblk], iproc)
2814 block_dest(iblk, jblk, kblk) = igroup*group_size + iproc
2819 ALLOCATE (bsizes1(nblks(1)), bsizes2(nblks(2)), bsizes3(nblks(3)))
2820 CALL dbt_get_info(t3c_main, blk_size_1=bsizes1, blk_size_2=bsizes2, blk_size_3=bsizes3)
2822 ALLOCATE (current_dest(nblks(1), nblks(2), nblks(3), 0:ngroups - 1))
2823 DO igroup = 0, ngroups - 1
2825 current_dest(:, :, :, igroup) = block_dest(:, :, :)
2826 CALL para_env%bcast(current_dest(:, :, :, igroup), source=igroup*group_size)
2830 batch_size = min(para_env%get_tag_ub(), 128000, nocc*ngroups)
2831 n_batch = (nocc*ngroups)/batch_size
2832 IF (
modulo(nocc*ngroups, batch_size) .NE. 0) n_batch = n_batch + 1
2834 DO i_batch = 1, n_batch
2836 ALLOCATE (send_buff(batch_size), recv_buff(batch_size))
2837 ALLOCATE (send_req(batch_size), recv_req(batch_size))
2841 DO kblk = 1, nblks(3)
2842 DO jblk = 1, nblks(2)
2843 DO iblk = 1, nblks(1)
2844 DO igroup = 0, ngroups - 1
2845 IF (block_source(iblk, jblk, kblk) == -1) cycle
2848 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2851 tag = i_msg - (i_batch - 1)*batch_size
2854 ind(:) = [iblk, jblk, kblk]
2855 iatom = idx_to_at(ind(dim_at))
2856 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
2860 IF (para_env%mepos == block_source(iblk, jblk, kblk))
THEN
2861 CALL dbt_get_block(t3c_main, [iblk, jblk, kblk], blk, found)
2865 IF (block_source(iblk, jblk, kblk) == current_dest(iblk, jblk, kblk, igroup))
THEN
2866 IF (found)
CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(blk), blk)
2868 IF (para_env%mepos == block_source(iblk, jblk, kblk) .AND. found)
THEN
2869 ALLOCATE (send_buff(tag)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2870 send_buff(tag)%array(:, :, :) = blk(:, :, :)
2872 CALL para_env%isend(msgin=send_buff(tag)%array, &
2873 dest=current_dest(iblk, jblk, kblk, igroup), &
2874 request=send_req(is), tag=tag)
2877 IF (para_env%mepos == current_dest(iblk, jblk, kblk, igroup))
THEN
2878 ALLOCATE (recv_buff(tag)%array(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2880 CALL para_env%irecv(msgout=recv_buff(tag)%array, source=block_source(iblk, jblk, kblk), &
2881 request=recv_req(ir), tag=tag)
2885 IF (found)
DEALLOCATE (blk)
2894 DO i = 1, batch_size
2895 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
2900 DO kblk = 1, nblks(3)
2901 DO jblk = 1, nblks(2)
2902 DO iblk = 1, nblks(1)
2903 DO igroup = 0, ngroups - 1
2904 IF (block_source(iblk, jblk, kblk) == -1) cycle
2907 IF (i_msg < (i_batch - 1)*batch_size + 1 .OR. i_msg > i_batch*batch_size) cycle
2910 tag = i_msg - (i_batch - 1)*batch_size
2913 ind(:) = [iblk, jblk, kblk]
2914 iatom = idx_to_at(ind(dim_at))
2915 IF (.NOT. iatom_to_subgroup(iatom)%array(igroup + 1)) cycle
2918 IF (para_env%mepos == current_dest(iblk, jblk, kblk, igroup) .AND. &
2919 block_source(iblk, jblk, kblk) .NE. current_dest(iblk, jblk, kblk, igroup))
THEN
2921 ALLOCATE (blk(bsizes1(iblk), bsizes2(jblk), bsizes3(kblk)))
2922 blk(:, :, :) = recv_buff(tag)%array(:, :, :)
2923 CALL dbt_put_block(t3c_sub, [iblk, jblk, kblk], shape(blk), blk)
2932 DO i = 1, batch_size
2933 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
2935 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
2937 CALL dbt_finalize(t3c_sub)
2939 END SUBROUTINE copy_3c_to_subgroup
2951 SUBROUTINE gather_ks_matrix(ks_t, ks_t_sub, group_size, sparsity_pattern, para_env, ri_data)
2952 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t, ks_t_sub
2953 INTEGER,
INTENT(IN) :: group_size
2954 INTEGER,
DIMENSION(:, :, :),
INTENT(IN) :: sparsity_pattern
2958 CHARACTER(len=*),
PARAMETER :: routinen =
'gather_ks_matrix'
2960 INTEGER :: b_img, dest, handle, i, i_spin, iatom, &
2961 igroup, ir, is, jatom, n_mess, natom, &
2962 nimg, nspins, source, tag
2964 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
2965 TYPE(
cp_2d_r_p_type),
ALLOCATABLE,
DIMENSION(:) :: recv_buff, send_buff
2966 TYPE(
mp_request_type),
ALLOCATABLE,
DIMENSION(:) :: recv_req, send_req
2968 CALL timeset(routinen, handle)
2970 nimg =
SIZE(sparsity_pattern, 3)
2971 natom =
SIZE(sparsity_pattern, 2)
2972 nspins =
SIZE(ks_t, 1)
2976 DO i_spin = 1, nspins
2979 IF (sparsity_pattern(iatom, jatom, b_img) > -1) n_mess = n_mess + 1
2984 ALLOCATE (send_buff(n_mess), recv_buff(n_mess))
2985 ALLOCATE (send_req(n_mess), recv_req(n_mess))
2991 DO i_spin = 1, nspins
2994 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
2999 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3000 CALL dbt_get_stored_coordinates(ks_t_sub(i_spin, b_img), [iatom, jatom], source)
3001 igroup = sparsity_pattern(iatom, jatom, b_img)
3002 source = source + igroup*group_size
3003 IF (para_env%mepos == source)
THEN
3004 CALL dbt_get_block(ks_t_sub(i_spin, b_img), [iatom, jatom], blk, found)
3005 IF (source == dest)
THEN
3006 IF (found)
CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3008 ALLOCATE (send_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3009 send_buff(n_mess)%array(:, :) = 0.0_dp
3011 send_buff(n_mess)%array(:, :) = blk(:, :)
3014 CALL para_env%isend(msgin=send_buff(n_mess)%array, dest=dest, &
3015 request=send_req(is), tag=tag)
3021 IF (para_env%mepos == dest .AND. source .NE. dest)
THEN
3022 ALLOCATE (recv_buff(n_mess)%array(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3024 CALL para_env%irecv(msgout=recv_buff(n_mess)%array, source=source, &
3025 request=recv_req(ir), tag=tag)
3036 DO i_spin = 1, nspins
3039 IF (sparsity_pattern(iatom, jatom, b_img) < 0) cycle
3042 CALL dbt_get_stored_coordinates(ks_t(i_spin, b_img), [iatom, jatom], dest)
3043 IF (para_env%mepos == dest)
THEN
3044 IF (.NOT.
ASSOCIATED(recv_buff(n_mess)%array)) cycle
3045 ALLOCATE (blk(ri_data%bsizes_AO(iatom), ri_data%bsizes_AO(jatom)))
3046 blk(:, :) = recv_buff(n_mess)%array(:, :)
3047 CALL dbt_put_block(ks_t(i_spin, b_img), [iatom, jatom], shape(blk), blk)
3056 IF (
ASSOCIATED(send_buff(i)%array))
DEALLOCATE (send_buff(i)%array)
3057 IF (
ASSOCIATED(recv_buff(i)%array))
DEALLOCATE (recv_buff(i)%array)
3059 DEALLOCATE (send_buff, recv_buff, send_req, recv_req)
3062 CALL timestop(handle)
3064 END SUBROUTINE gather_ks_matrix
3079 SUBROUTINE get_subgroup_2c_tensors(mat_2c_pot, t_2c_work, t_2c_ao_tmp, ks_t_split, ks_t_sub, &
3080 group_size, ngroups, para_env, para_env_sub, ri_data)
3081 TYPE(dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3082 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work, t_2c_ao_tmp, ks_t_split
3083 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: ks_t_sub
3084 INTEGER,
INTENT(IN) :: group_size, ngroups
3088 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_tensors'
3090 INTEGER :: handle, i, i_img, i_ri, i_spin, iproc, &
3091 j, natom, nblks, nimg, nspins
3092 INTEGER(int_8) :: nze
3093 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3095 INTEGER,
DIMENSION(2) :: pdims_2d
3096 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3097 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3099 TYPE(dbcsr_distribution_type) :: dbcsr_dist_sub
3100 TYPE(dbt_pgrid_type) :: pgrid_2d
3101 TYPE(dbt_type) :: work, work_sub
3103 CALL timeset(routinen, handle)
3107 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3109 natom =
SIZE(ri_data%bsizes_RI)
3110 nblks =
SIZE(ri_data%bsizes_RI_split)
3111 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3112 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3113 DO i_ri = 1, ri_data%ncell_RI
3114 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3115 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3120 bsizes_ri_ext, bsizes_ri_ext, &
3122 DEALLOCATE (dist1, dist2)
3125 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3127 DEALLOCATE (dist1, dist2)
3131 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3133 DEALLOCATE (dist1, dist2)
3134 CALL dbt_create(ks_t_split(1), ks_t_split(2))
3137 ri_data%bsizes_AO, ri_data%bsizes_AO, &
3139 DEALLOCATE (dist1, dist2)
3141 nspins =
SIZE(ks_t_sub, 1)
3142 nimg =
SIZE(ks_t_sub, 2)
3144 DO i_spin = 1, nspins
3145 CALL dbt_create(t_2c_ao_tmp(1), ks_t_sub(i_spin, i_img))
3152 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3154 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3156 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3158 DO i = 0, pdims_2d(1) - 1
3159 DO j = 0, pdims_2d(2) - 1
3160 dbcsr_pgrid(i, j) = iproc
3166 ALLOCATE (col_dist(natom), row_dist(natom))
3167 row_dist(:) = dist1(:)
3168 col_dist(:) = dist2(:)
3170 ALLOCATE (ri_blk_size(natom))
3171 ri_blk_size(:) = ri_data%bsizes_RI(:)
3173 CALL dbcsr_distribution_new(dbcsr_dist_sub, group=para_env_sub%get_handle(), pgrid=dbcsr_pgrid, &
3174 row_dist=row_dist, col_dist=col_dist)
3175 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3176 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3179 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3180 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3184 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3185 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3186 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3187 CALL dbt_clear(work_sub)
3190 CALL dbt_destroy(work)
3191 CALL dbt_destroy(work_sub)
3192 CALL dbt_pgrid_destroy(pgrid_2d)
3193 CALL dbcsr_distribution_release(dbcsr_dist_sub)
3194 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3195 CALL timestop(handle)
3197 END SUBROUTINE get_subgroup_2c_tensors
3212 SUBROUTINE get_subgroup_3c_tensors(t_3c_int, t_3c_work_2, t_3c_work_3, t_3c_apc, t_3c_apc_sub, &
3213 group_size, ngroups, para_env, para_env_sub, ri_data)
3214 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_int, t_3c_work_2, t_3c_work_3
3215 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_apc, t_3c_apc_sub
3216 INTEGER,
INTENT(IN) :: group_size, ngroups
3220 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_tensors'
3222 INTEGER :: batch_size, bfac, bo(2), handle, &
3223 handle2, i_blk, i_img, i_ri, i_spin, &
3224 ib, natom, nblks_ao, nblks_ri, nimg, &
3226 INTEGER(int_8) :: nze
3227 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3228 bsizes_stack, bsizes_tmp, dist1, &
3229 dist2, dist3, dist_stack, idx_to_at
3230 INTEGER,
DIMENSION(3) :: pdims
3232 TYPE(dbt_distribution_type) :: t_dist
3233 TYPE(dbt_pgrid_type) :: pgrid
3234 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3236 CALL timeset(routinen, handle)
3238 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3239 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3240 DO i_ri = 1, ri_data%ncell_RI
3241 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3247 natom =
SIZE(ri_data%bsizes_RI)
3248 nblks_ri = max(1, natom/bfac)
3249 ALLOCATE (bsizes_tmp(nblks_ri))
3250 DO i_blk = 1, nblks_ri
3251 bo =
get_limit(natom, nblks_ri, i_blk - 1)
3252 bsizes_tmp(i_blk) = sum(ri_data%bsizes_RI(bo(1):bo(2)))
3254 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3255 DO i_ri = 1, ri_data%ncell_RI
3256 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = bsizes_tmp(:)
3259 batch_size = ri_data%kp_stack_size
3260 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3261 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3262 DO ib = 1, batch_size
3263 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3267 natom =
SIZE(ri_data%bsizes_RI)
3269 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3270 tensor_dims=[
SIZE(bsizes_ri_ext_split), 1, batch_size*
SIZE(ri_data%bsizes_AO_split)])
3274 pgrid, bsizes_ri_ext_split, ri_data%bsizes_AO_split, &
3275 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(RI | AO AO)")
3276 nimg =
SIZE(t_3c_int)
3278 CALL dbt_create(t_3c_int(1), t_3c_int(i_img))
3282 ALLOCATE (dist_stack(batch_size*nblks_ao))
3283 DO ib = 1, batch_size
3284 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3287 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3288 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3289 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3290 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3291 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3292 CALL dbt_distribution_destroy(t_dist)
3293 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3297 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3298 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3299 DEALLOCATE (dist1, dist2, dist3)
3302 ri_data%pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3303 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3304 DEALLOCATE (dist1, dist2, dist3)
3307 CALL timeset(routinen//
"_ints", handle2)
3308 IF (
ALLOCATED(ri_data%kp_t_3c_int))
THEN
3310 CALL dbt_copy(ri_data%kp_t_3c_int(i_img), t_3c_int(i_img), move_data=.true.)
3313 ALLOCATE (ri_data%kp_t_3c_int(nimg))
3315 CALL dbt_create(t_3c_int(i_img), ri_data%kp_t_3c_int(i_img))
3318 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, i_img), work_atom_block, order=[2, 1, 3])
3319 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, group_size, ngroups, para_env)
3320 CALL dbt_copy(work_atom_block_sub, t_3c_int(i_img), move_data=.true.)
3323 CALL timestop(handle2)
3324 CALL dbt_pgrid_destroy(pgrid)
3325 CALL dbt_destroy(work_atom_block)
3326 CALL dbt_destroy(work_atom_block_sub)
3330 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3331 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3335 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3336 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3337 DEALLOCATE (dist1, dist2, dist3)
3340 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3341 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3342 DEALLOCATE (dist1, dist2, dist3)
3346 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3347 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3350 ALLOCATE (dist_stack(batch_size*nblks_ao))
3351 DO ib = 1, batch_size
3352 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3355 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3356 CALL dbt_create(t_3c_work_2(1),
"work_2_stack", t_dist, [1], [2, 3], &
3357 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3358 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3359 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3360 CALL dbt_distribution_destroy(t_dist)
3361 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3364 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3366 nspins =
SIZE(t_3c_apc, 1)
3367 CALL timeset(routinen//
"_apc", handle2)
3369 DO i_spin = 1, nspins
3370 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3373 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3374 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, group_size, &
3375 ngroups, para_env, ri_data%iatom_to_subgroup, 1, idx_to_at)
3376 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3378 DO i_spin = 1, nspins
3379 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3382 CALL timestop(handle2)
3383 CALL dbt_pgrid_destroy(pgrid)
3384 CALL dbt_destroy(tmp)
3385 CALL dbt_destroy(work_atom_block)
3386 CALL dbt_destroy(work_atom_block_sub)
3388 CALL timestop(handle)
3390 END SUBROUTINE get_subgroup_3c_tensors
3412 SUBROUTINE get_subgroup_2c_derivs(t_2c_inv, t_2c_bint, t_2c_metric, mat_2c_pot, t_2c_work, rho_ao_t, &
3413 rho_ao_t_sub, t_2c_der_metric, t_2c_der_metric_sub, mat_der_pot, &
3414 mat_der_pot_sub, group_size, ngroups, para_env, para_env_sub, ri_data)
3415 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_inv, t_2c_bint, t_2c_metric
3416 TYPE(dbcsr_type),
DIMENSION(:),
INTENT(INOUT) :: mat_2c_pot
3417 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_work
3418 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: rho_ao_t, rho_ao_t_sub, t_2c_der_metric, &
3420 TYPE(dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot, mat_der_pot_sub
3421 INTEGER,
INTENT(IN) :: group_size, ngroups
3425 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_2c_derivs'
3427 INTEGER :: handle, i, i_img, i_ri, i_spin, i_xyz, &
3428 iatom, iproc, j, natom, nblks, nimg, &
3430 INTEGER(int_8) :: nze
3431 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3433 INTEGER,
DIMENSION(2) :: pdims_2d
3434 INTEGER,
DIMENSION(:),
POINTER :: col_dist, ri_blk_size, row_dist
3435 INTEGER,
DIMENSION(:, :),
POINTER :: dbcsr_pgrid
3437 TYPE(dbcsr_distribution_type) :: dbcsr_dist_sub
3438 TYPE(dbt_pgrid_type) :: pgrid_2d
3439 TYPE(dbt_type) :: work, work_sub
3441 CALL timeset(routinen, handle)
3446 CALL dbt_pgrid_create(para_env_sub, pdims_2d, pgrid_2d)
3448 natom =
SIZE(ri_data%bsizes_RI)
3449 nblks =
SIZE(ri_data%bsizes_RI_split)
3450 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*natom))
3451 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks))
3452 DO i_ri = 1, ri_data%ncell_RI
3453 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
3454 bsizes_ri_ext_split((i_ri - 1)*nblks + 1:i_ri*nblks) = ri_data%bsizes_RI_split(:)
3459 bsizes_ri_ext, bsizes_ri_ext, &
3461 DEALLOCATE (dist1, dist2)
3463 CALL dbt_create(t_2c_inv(1), t_2c_bint(1))
3464 CALL dbt_create(t_2c_inv(1), t_2c_metric(1))
3466 CALL dbt_create(t_2c_inv(1), t_2c_inv(iatom))
3467 CALL dbt_create(t_2c_inv(1), t_2c_bint(iatom))
3468 CALL dbt_create(t_2c_inv(1), t_2c_metric(iatom))
3470 CALL dbt_create(t_2c_inv(1), t_2c_work(1))
3471 CALL dbt_create(t_2c_inv(1), t_2c_work(2))
3472 CALL dbt_create(t_2c_inv(1), t_2c_work(3))
3473 CALL dbt_create(t_2c_inv(1), t_2c_work(4))
3476 bsizes_ri_ext_split, bsizes_ri_ext_split, &
3478 DEALLOCATE (dist1, dist2)
3482 CALL copy_2c_to_subgroup(t_2c_inv(iatom), ri_data%t_2c_inv(1, iatom), group_size, ngroups, para_env)
3483 CALL copy_2c_to_subgroup(t_2c_bint(iatom), ri_data%t_2c_int(1, iatom), group_size, ngroups, para_env)
3484 CALL copy_2c_to_subgroup(t_2c_metric(iatom), ri_data%t_2c_pot(1, iatom), group_size, ngroups, para_env)
3490 CALL dbt_create(t_2c_inv(1), t_2c_der_metric_sub(iatom, i_xyz))
3491 CALL copy_2c_to_subgroup(t_2c_der_metric_sub(iatom, i_xyz), t_2c_der_metric(iatom, i_xyz), &
3492 group_size, ngroups, para_env)
3493 CALL dbt_destroy(t_2c_der_metric(iatom, i_xyz))
3499 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3501 DEALLOCATE (dist1, dist2)
3502 nspins =
SIZE(rho_ao_t, 1)
3503 nimg =
SIZE(rho_ao_t, 2)
3506 DO i_spin = 1, nspins
3507 IF (.NOT. (i_img == 1 .AND. i_spin == 1)) &
3508 CALL dbt_create(rho_ao_t_sub(1, 1), rho_ao_t_sub(i_spin, i_img))
3509 CALL copy_2c_to_subgroup(rho_ao_t_sub(i_spin, i_img), rho_ao_t(i_spin, i_img), &
3510 group_size, ngroups, para_env)
3511 CALL dbt_destroy(rho_ao_t(i_spin, i_img))
3517 ri_data%bsizes_RI, ri_data%bsizes_RI, &
3519 CALL dbt_create(ri_data%kp_mat_2c_pot(1, 1), work)
3521 ALLOCATE (dbcsr_pgrid(0:pdims_2d(1) - 1, 0:pdims_2d(2) - 1))
3523 DO i = 0, pdims_2d(1) - 1
3524 DO j = 0, pdims_2d(2) - 1
3525 dbcsr_pgrid(i, j) = iproc
3531 ALLOCATE (col_dist(natom), row_dist(natom))
3532 row_dist(:) = dist1(:)
3533 col_dist(:) = dist2(:)
3535 ALLOCATE (ri_blk_size(natom))
3536 ri_blk_size(:) = ri_data%bsizes_RI(:)
3538 CALL dbcsr_distribution_new(dbcsr_dist_sub, group=para_env_sub%get_handle(), pgrid=dbcsr_pgrid, &
3539 row_dist=row_dist, col_dist=col_dist)
3540 CALL dbcsr_create(mat_2c_pot(1), dist=dbcsr_dist_sub, name=
"sub", matrix_type=dbcsr_type_no_symmetry, &
3541 row_blk_size=ri_blk_size, col_blk_size=ri_blk_size)
3545 IF (i_img > 1)
CALL dbcsr_create(mat_2c_pot(i_img), template=mat_2c_pot(1))
3546 CALL dbt_copy_matrix_to_tensor(ri_data%kp_mat_2c_pot(1, i_img), work)
3550 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3551 CALL dbt_copy_tensor_to_matrix(work_sub, mat_2c_pot(i_img))
3552 CALL dbcsr_filter(mat_2c_pot(i_img), ri_data%filter_eps)
3553 CALL dbt_clear(work_sub)
3559 CALL dbcsr_create(mat_der_pot_sub(i_img, i_xyz), template=mat_2c_pot(1))
3560 CALL dbt_copy_matrix_to_tensor(mat_der_pot(i_img, i_xyz), work)
3561 CALL dbcsr_release(mat_der_pot(i_img, i_xyz))
3565 CALL copy_2c_to_subgroup(work_sub, work, group_size, ngroups, para_env)
3566 CALL dbt_copy_tensor_to_matrix(work_sub, mat_der_pot_sub(i_img, i_xyz))
3567 CALL dbcsr_filter(mat_der_pot_sub(i_img, i_xyz), ri_data%filter_eps)
3568 CALL dbt_clear(work_sub)
3572 CALL dbt_destroy(work)
3573 CALL dbt_destroy(work_sub)
3574 CALL dbt_pgrid_destroy(pgrid_2d)
3575 CALL dbcsr_distribution_release(dbcsr_dist_sub)
3576 DEALLOCATE (col_dist, row_dist, ri_blk_size, dbcsr_pgrid)
3578 CALL timestop(handle)
3580 END SUBROUTINE get_subgroup_2c_derivs
3600 SUBROUTINE get_subgroup_3c_derivs(t_3c_work_2, t_3c_work_3, t_3c_der_AO, t_3c_der_AO_sub, &
3601 t_3c_der_RI, t_3c_der_RI_sub, t_3c_apc, t_3c_apc_sub, &
3602 t_3c_der_stack, group_size, ngroups, para_env, para_env_sub, &
3604 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_work_2, t_3c_work_3
3605 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ao, t_3c_der_ao_sub, &
3606 t_3c_der_ri, t_3c_der_ri_sub, &
3607 t_3c_apc, t_3c_apc_sub
3608 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_der_stack
3609 INTEGER,
INTENT(IN) :: group_size, ngroups
3613 CHARACTER(len=*),
PARAMETER :: routinen =
'get_subgroup_3c_derivs'
3615 INTEGER :: batch_size, handle, i_img, i_ri, i_spin, &
3616 i_xyz, ib, nblks_ao, nblks_ri, nimg, &
3618 INTEGER(int_8) :: nze
3619 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, &
3620 bsizes_stack, dist1, dist2, dist3, &
3621 dist_stack, idx_to_at
3623 TYPE(dbt_distribution_type) :: t_dist
3624 TYPE(dbt_pgrid_type) :: pgrid
3625 TYPE(dbt_type) :: tmp, work_atom_block, work_atom_block_sub
3627 CALL timeset(routinen, handle)
3630 nblks_ri =
SIZE(ri_data%bsizes_RI)
3631 ALLOCATE (bsizes_ri_ext(ri_data%ncell_RI*nblks_ri))
3632 DO i_ri = 1, ri_data%ncell_RI
3633 bsizes_ri_ext((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI(:)
3636 CALL dbt_get_info(ri_data%kp_t_3c_int(1), pdims=pdims)
3637 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3640 pgrid, bsizes_ri_ext, ri_data%bsizes_AO, &
3641 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3642 DEALLOCATE (dist1, dist2, dist3)
3645 ri_data%pgrid_2, bsizes_ri_ext, ri_data%bsizes_AO, &
3646 ri_data%bsizes_AO, [1], [2, 3], name=
"(RI | AO AO)")
3647 DEALLOCATE (dist1, dist2, dist3)
3648 CALL dbt_pgrid_destroy(pgrid)
3654 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ao_sub(i_img, i_xyz))
3658 CALL dbt_copy(t_3c_der_ao(i_img, i_xyz), work_atom_block, move_data=.true.)
3659 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3660 group_size, ngroups, para_env)
3661 CALL dbt_copy(work_atom_block_sub, t_3c_der_ao_sub(i_img, i_xyz), move_data=.true.)
3665 CALL dbt_create(ri_data%kp_t_3c_int(1), t_3c_der_ri_sub(i_img, i_xyz))
3669 CALL dbt_copy(t_3c_der_ri(i_img, i_xyz), work_atom_block, move_data=.true.)
3670 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, &
3671 group_size, ngroups, para_env)
3672 CALL dbt_copy(work_atom_block_sub, t_3c_der_ri_sub(i_img, i_xyz), move_data=.true.)
3676 CALL dbt_destroy(t_3c_der_ri(i_img, i_xyz))
3677 CALL dbt_destroy(t_3c_der_ao(i_img, i_xyz))
3680 CALL dbt_destroy(work_atom_block_sub)
3681 CALL dbt_destroy(work_atom_block)
3684 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
3685 ALLOCATE (bsizes_ri_ext_split(ri_data%ncell_RI*nblks_ri))
3686 DO i_ri = 1, ri_data%ncell_RI
3687 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
3691 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid, &
3692 tensor_dims=[1,
SIZE(bsizes_ri_ext_split), batch_size*
SIZE(ri_data%bsizes_AO_split)])
3695 pgrid, ri_data%bsizes_AO, bsizes_ri_ext, &
3696 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3697 DEALLOCATE (dist1, dist2, dist3)
3700 ri_data%pgrid_1, ri_data%bsizes_AO, bsizes_ri_ext, &
3701 ri_data%bsizes_AO, [1], [2, 3], name=
"(AO RI | AO)")
3702 DEALLOCATE (dist1, dist2, dist3)
3705 pgrid, ri_data%bsizes_AO_split, bsizes_ri_ext_split, &
3706 ri_data%bsizes_AO_split, [1], [2, 3], name=
"(AO RI | AO)")
3707 DEALLOCATE (dist1, dist2, dist3)
3709 ALLOCATE (idx_to_at(
SIZE(ri_data%bsizes_AO)))
3711 nspins =
SIZE(t_3c_apc, 1)
3713 DO i_spin = 1, nspins
3714 CALL dbt_create(tmp, t_3c_apc_sub(i_spin, i_img))
3717 CALL dbt_copy(t_3c_apc(i_spin, i_img), work_atom_block, move_data=.true.)
3718 CALL copy_3c_to_subgroup(work_atom_block_sub, work_atom_block, group_size, &
3719 ngroups, para_env, ri_data%iatom_to_subgroup, 1, idx_to_at)
3720 CALL dbt_copy(work_atom_block_sub, t_3c_apc_sub(i_spin, i_img), move_data=.true.)
3722 DO i_spin = 1, nspins
3723 CALL dbt_destroy(t_3c_apc(i_spin, i_img))
3726 CALL dbt_destroy(tmp)
3727 CALL dbt_destroy(work_atom_block)
3728 CALL dbt_destroy(work_atom_block_sub)
3729 CALL dbt_pgrid_destroy(pgrid)
3732 batch_size = ri_data%kp_stack_size
3733 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
3734 ALLOCATE (bsizes_stack(batch_size*nblks_ao))
3735 DO ib = 1, batch_size
3736 bsizes_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = ri_data%bsizes_AO_split(:)
3739 ALLOCATE (dist1(ri_data%ncell_RI*nblks_ri), dist2(nblks_ao), dist3(nblks_ao))
3740 CALL dbt_get_info(ri_data%kp_t_3c_int(1), proc_dist_1=dist1, proc_dist_2=dist2, &
3741 proc_dist_3=dist3, pdims=pdims)
3743 ALLOCATE (dist_stack(batch_size*nblks_ao))
3744 DO ib = 1, batch_size
3745 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3748 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3749 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3750 CALL dbt_create(t_3c_work_3(1),
"work_3_stack", t_dist, [1], [2, 3], &
3751 bsizes_ri_ext_split, ri_data%bsizes_AO_split, bsizes_stack)
3752 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(2))
3753 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(3))
3754 CALL dbt_create(t_3c_work_3(1), t_3c_work_3(4))
3755 CALL dbt_distribution_destroy(t_dist)
3756 CALL dbt_pgrid_destroy(pgrid)
3757 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3760 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(1))
3761 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(2))
3762 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(3))
3763 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(4))
3764 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(5))
3765 CALL dbt_create(t_3c_work_3(1), t_3c_der_stack(6))
3768 ALLOCATE (dist1(nblks_ao), dist2(ri_data%ncell_RI*nblks_ri), dist3(nblks_ao))
3769 CALL dbt_get_info(t_3c_apc_sub(1, 1), proc_dist_1=dist1, proc_dist_2=dist2, &
3770 proc_dist_3=dist3, pdims=pdims)
3772 ALLOCATE (dist_stack(batch_size*nblks_ao))
3773 DO ib = 1, batch_size
3774 dist_stack((ib - 1)*nblks_ao + 1:ib*nblks_ao) = dist3(:)
3777 CALL dbt_pgrid_create(para_env_sub, pdims, pgrid)
3778 CALL dbt_distribution_new(t_dist, pgrid, dist1, dist2, dist_stack)
3779 CALL dbt_create(t_3c_work_2(1),
"work_3_stack", t_dist, [1], [2, 3], &
3780 ri_data%bsizes_AO_split, bsizes_ri_ext_split, bsizes_stack)
3781 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(2))
3782 CALL dbt_create(t_3c_work_2(1), t_3c_work_2(3))
3783 CALL dbt_distribution_destroy(t_dist)
3784 CALL dbt_pgrid_destroy(pgrid)
3785 DEALLOCATE (dist1, dist2, dist3, dist_stack)
3787 CALL timestop(handle)
3789 END SUBROUTINE get_subgroup_3c_derivs
3797 SUBROUTINE reorder_3c_ints(t_3c_ints, ri_data)
3798 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_ints
3801 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_ints'
3803 INTEGER :: handle, i_img,
idx, idx_empty, idx_full, &
3805 INTEGER(int_8) :: nze
3807 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
3809 CALL timeset(routinen, handle)
3812 ALLOCATE (t_3c_tmp(nimg))
3814 CALL dbt_create(t_3c_ints(i_img), t_3c_tmp(i_img))
3815 CALL dbt_copy(t_3c_ints(i_img), t_3c_tmp(i_img), move_data=.true.)
3820 ALLOCATE (ri_data%idx_to_img(nimg))
3822 idx_empty = nimg + 1
3827 idx_empty = idx_empty - 1
3828 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_empty), move_data=.true.)
3829 ri_data%idx_to_img(idx_empty) = i_img
3831 idx_full = idx_full + 1
3832 CALL dbt_copy(t_3c_tmp(i_img), t_3c_ints(idx_full), move_data=.true.)
3833 ri_data%idx_to_img(idx_full) = i_img
3835 CALL dbt_destroy(t_3c_tmp(i_img))
3839 ri_data%nimg_nze = idx_full
3841 ALLOCATE (ri_data%img_to_idx(nimg))
3843 ri_data%img_to_idx(ri_data%idx_to_img(
idx)) =
idx
3846 CALL timestop(handle)
3848 END SUBROUTINE reorder_3c_ints
3856 SUBROUTINE reorder_3c_derivs(t_3c_derivs, ri_data)
3857 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_derivs
3860 CHARACTER(LEN=*),
PARAMETER :: routinen =
'reorder_3c_derivs'
3862 INTEGER :: handle, i_img, i_xyz,
idx, nimg
3863 INTEGER(int_8) :: nze
3865 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_3c_tmp
3867 CALL timeset(routinen, handle)
3870 ALLOCATE (t_3c_tmp(nimg))
3872 CALL dbt_create(t_3c_derivs(1, 1), t_3c_tmp(i_img))
3877 CALL dbt_copy(t_3c_derivs(i_img, i_xyz), t_3c_tmp(i_img), move_data=.true.)
3880 idx = ri_data%img_to_idx(i_img)
3881 CALL dbt_copy(t_3c_tmp(i_img), t_3c_derivs(
idx, i_xyz), move_data=.true.)
3883 IF (nze > 0) ri_data%nimg_nze = max(
idx, ri_data%nimg_nze)
3888 CALL dbt_destroy(t_3c_tmp(i_img))
3891 CALL timestop(handle)
3893 END SUBROUTINE reorder_3c_derivs
3901 SUBROUTINE get_sparsity_pattern(pattern, ri_data, qs_env)
3902 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: pattern
3906 INTEGER :: iatom, j_img, jatom, mj_img, natom, nimg
3907 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bins
3908 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :) :: tmp_pattern
3909 INTEGER,
DIMENSION(3) :: cell_j
3910 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
3911 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
3916 DIMENSION(:),
POINTER :: nl_iterator
3920 NULLIFY (nl_2c, nl_iterator, kpoints, cell_to_index, dft_control, index_to_cell, para_env)
3922 CALL get_qs_env(qs_env, kpoints=kpoints, dft_control=dft_control, para_env=para_env, natom=natom)
3923 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index, index_to_cell=index_to_cell, sab_nl=nl_2c)
3926 pattern(:, :, :) = 0
3933 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3934 IF (j_img > nimg .OR. j_img < 1) cycle
3936 mj_img = get_opp_index(j_img, qs_env)
3937 IF (mj_img > nimg .OR. mj_img < 1) cycle
3939 IF (ri_data%present_images(j_img) == 0) cycle
3941 pattern(iatom, jatom, j_img) = 1
3952 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
3953 IF (j_img > nimg .OR. j_img < 1) cycle
3955 mj_img = get_opp_index(j_img, qs_env)
3956 IF (mj_img .LE. nimg .AND. mj_img > 0) cycle
3958 IF (ri_data%present_images(j_img) == 0) cycle
3960 pattern(iatom, jatom, j_img) = 1
3964 CALL para_env%sum(pattern)
3969 IF (pattern(iatom, iatom, j_img) .NE. 0)
THEN
3970 mj_img = get_opp_index(j_img, qs_env)
3971 IF (mj_img > nimg .OR. mj_img < 1) cycle
3972 pattern(iatom, iatom, mj_img) = 0
3979 ALLOCATE (bins(natom))
3982 ALLOCATE (tmp_pattern(natom, natom, nimg))
3983 tmp_pattern(:, :, :) = 0
3987 IF (pattern(iatom, jatom, j_img) == 0) cycle
3988 mj_img = get_opp_index(j_img, qs_env)
3991 IF (mj_img > nimg .OR. mj_img < 1)
THEN
3993 bins(iatom) = bins(iatom) + 1
3994 tmp_pattern(iatom, jatom, j_img) = 1
3997 IF (bins(iatom) > bins(jatom))
THEN
3998 bins(jatom) = bins(jatom) + 1
3999 tmp_pattern(jatom, iatom, mj_img) = 1
4001 bins(iatom) = bins(iatom) + 1
4002 tmp_pattern(iatom, jatom, j_img) = 1
4010 pattern(:, :, :) = tmp_pattern(:, :, :) - 1
4012 END SUBROUTINE get_sparsity_pattern
4022 SUBROUTINE get_sub_dist(sparsity_pattern, ngroups, ri_data)
4023 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: sparsity_pattern
4024 INTEGER,
INTENT(IN) :: ngroups
4027 INTEGER :: b_img, ctr, iat, iatom, igroup, jatom, &
4029 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: max_at_per_group
4031 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4033 natom =
SIZE(sparsity_pattern, 2)
4034 nimg =
SIZE(sparsity_pattern, 3)
4039 IF (.NOT.
ALLOCATED(ri_data%iatom_to_subgroup))
THEN
4040 ALLOCATE (ri_data%iatom_to_subgroup(natom), max_at_per_group(ngroups))
4042 NULLIFY (ri_data%iatom_to_subgroup(iatom)%array)
4043 ALLOCATE (ri_data%iatom_to_subgroup(iatom)%array(ngroups))
4044 ri_data%iatom_to_subgroup(iatom)%array(:) = .false.
4048 IF (ub*ngroups < natom) ub = ub + 1
4049 max_at_per_group(:) = max(1, ub)
4054 DO WHILE (
modulo(sum(max_at_per_group), natom) .NE. 0)
4055 igroup =
modulo(ctr, ngroups) + 1
4056 max_at_per_group(igroup) = max_at_per_group(igroup) + 1
4061 DO igroup = 1, ngroups
4062 DO iat = 1, max_at_per_group(igroup)
4063 iatom =
modulo(ctr, natom) + 1
4064 ri_data%iatom_to_subgroup(iatom)%array(igroup) = .true.
4070 ALLOCATE (bins(ngroups))
4075 IF (sparsity_pattern(iatom, jatom, b_img) == -1) cycle
4076 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4079 IF (any(ri_data%kp_cost > epsilon(0.0_dp)))
THEN
4080 cost = ri_data%kp_cost(iatom, jatom, b_img)
4082 cost = real(ri_data%bsizes_AO(iatom)*ri_data%bsizes_AO(jatom),
dp)
4084 bins(igroup + 1) = bins(igroup + 1) + cost
4085 sparsity_pattern(iatom, jatom, b_img) = igroup
4090 END SUBROUTINE get_sub_dist
4101 SUBROUTINE update_pattern_to_forces(force_pattern, scf_pattern, ngroups, ri_data, qs_env)
4102 INTEGER,
DIMENSION(:, :, :),
INTENT(INOUT) :: force_pattern, scf_pattern
4103 INTEGER,
INTENT(IN) :: ngroups
4107 INTEGER :: b_img, iatom, igroup, jatom, mb_img, &
4109 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: bins
4111 natom =
SIZE(scf_pattern, 2)
4112 nimg =
SIZE(scf_pattern, 3)
4114 ALLOCATE (bins(ngroups))
4118 mb_img = get_opp_index(b_img, qs_env)
4122 igroup = minloc(bins, 1, mask=ri_data%iatom_to_subgroup(iatom)%array) - 1
4125 IF (scf_pattern(iatom, jatom, b_img) > -1) cycle
4128 IF (mb_img > 0 .AND. mb_img .LE. nimg)
THEN
4129 IF (scf_pattern(jatom, iatom, mb_img) == -1) cycle
4130 bins(igroup + 1) = bins(igroup + 1) + ri_data%kp_cost(jatom, iatom, mb_img)
4131 force_pattern(iatom, jatom, b_img) = igroup
4137 END SUBROUTINE update_pattern_to_forces
4145 SUBROUTINE get_kp_and_ri_images(ri_data, qs_env)
4149 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_kp_and_ri_images'
4151 INTEGER :: cell_j(3), cell_k(3), handle, i_img, iatom, ikind, j_img, jatom, jcell, katom, &
4152 kcell, kp_index_lbounds(3), kp_index_ubounds(3), natom, ngroups, nimg, nkind, pcoord(3), &
4154 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_ri, &
4155 nri_per_atom, present_img, ri_cells
4156 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
4157 REAL(
dp) :: bump_fact, dij, dik, image_range, &
4158 ri_range, rij(3), rik(3)
4159 TYPE(dbt_type) :: t_dummy
4164 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4171 DIMENSION(:),
POINTER :: nl_iterator
4175 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4178 NULLIFY (qs_kind_set, dist_2d, nl_2c, nl_iterator, dft_control, &
4179 particle_set, kpoints, para_env, cell_to_index, hfx_section)
4181 CALL timeset(routinen, handle)
4183 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, &
4184 dft_control=dft_control, particle_set=particle_set, kpoints=kpoints, &
4185 para_env=para_env, natom=natom)
4186 nimg = dft_control%nimages
4188 kp_index_lbounds = lbound(cell_to_index)
4189 kp_index_ubounds = ubound(cell_to_index)
4194 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4201 CALL erfc_cutoff(ri_data%eps_schwarz, ri_data%hfx_pot%omega, ri_data%hfx_pot%cutoff_radius)
4205 ri_data%kp_RI_range = 0.0_dp
4206 ri_data%kp_image_range = 0.0_dp
4211 ri_data%kp_RI_range = max(ri_range, ri_data%kp_RI_range)
4215 CALL get_gto_basis_set(basis_set_ri(ikind)%gto_basis_set, kind_radius=image_range)
4218 ri_data%kp_image_range = max(image_range, ri_data%kp_image_range)
4222 ri_data%kp_bump_rad = bump_fact*ri_data%kp_RI_range
4228 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
4230 ALLOCATE (present_img(nimg))
4239 j_img = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4240 IF (j_img > nimg .OR. j_img < 1) cycle
4242 IF (dij > ri_data%kp_image_range) cycle
4244 ri_data%nimg = max(j_img, ri_data%nimg)
4245 present_img(j_img) = 1
4250 CALL para_env%max(ri_data%nimg)
4251 IF (ri_data%nimg > nimg) &
4252 cpabort(
"Make sure the smallest exponent of the RI-HFX basis is larger than that of the ORB basis.")
4255 CALL para_env%sum(present_img)
4256 ALLOCATE (ri_data%present_images(ri_data%nimg))
4257 ri_data%present_images = 0
4258 DO i_img = 1, ri_data%nimg
4259 IF (present_img(i_img) > 0) ri_data%present_images(i_img) = 1
4263 ri_data%pgrid, ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4264 map1=[1, 2], map2=[3], name=
"(AO AO | RI)")
4266 CALL dbt_mp_environ_pgrid(ri_data%pgrid, pdims, pcoord)
4267 CALL mp_comm_t3c%create(ri_data%pgrid%mp_comm_2d, 3, pdims)
4269 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4270 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4271 CALL dbt_destroy(t_dummy)
4276 ri_data%ri_metric,
"HFX_3c_nl", qs_env, op_pos=2, sym_ij=.false., &
4279 ALLOCATE (ri_cells(nimg))
4282 ALLOCATE (nri_per_atom(natom))
4288 iatom=iatom, jatom=jatom, katom=katom)
4291 IF (any([cell_j(1), cell_j(2), cell_j(3)] < kp_index_lbounds) .OR. &
4292 any([cell_j(1), cell_j(2), cell_j(3)] > kp_index_ubounds)) cycle
4294 jcell = cell_to_index(cell_j(1), cell_j(2), cell_j(3))
4295 IF (jcell > nimg .OR. jcell < 1) cycle
4297 IF (any([cell_k(1), cell_k(2), cell_k(3)] < kp_index_lbounds) .OR. &
4298 any([cell_k(1), cell_k(2), cell_k(3)] > kp_index_ubounds)) cycle
4300 kcell = cell_to_index(cell_k(1), cell_k(2), cell_k(3))
4301 IF (kcell > nimg .OR. kcell < 1) cycle
4303 IF (dik > ri_data%kp_RI_range) cycle
4306 IF (jcell == 1 .AND. iatom == jatom) nri_per_atom(iatom) = nri_per_atom(iatom) + ri_data%bsizes_RI(katom)
4310 CALL para_env%sum(ri_cells)
4311 CALL para_env%sum(nri_per_atom)
4313 ALLOCATE (ri_data%img_to_RI_cell(nimg))
4314 ri_data%ncell_RI = 0
4315 ri_data%img_to_RI_cell = 0
4317 IF (ri_cells(i_img) > 0)
THEN
4318 ri_data%ncell_RI = ri_data%ncell_RI + 1
4319 ri_data%img_to_RI_cell(i_img) = ri_data%ncell_RI
4323 ALLOCATE (ri_data%RI_cell_to_img(ri_data%ncell_RI))
4325 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
4329 IF (ri_data%unit_nr > 0)
THEN
4330 WRITE (ri_data%unit_nr, fmt=
"(/T3,A,I29)") &
4331 "KP-HFX_RI_INFO| Number of RI-KP parallel groups:", ngroups
4332 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F31.3,A)") &
4333 "KP-HFX_RI_INFO| RI basis extension radius:", ri_data%kp_RI_range*
angstrom,
" Ang"
4334 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F12.3,A, F6.3, A)") &
4335 "KP-HFX_RI_INFO| RI basis bump factor and bump radius:", bump_fact,
" /", &
4336 ri_data%kp_bump_rad*
angstrom,
" Ang"
4337 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I16,A)") &
4338 "KP-HFX_RI_INFO| The extended RI bases cover up to ", ri_data%ncell_RI,
" unit cells"
4339 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I18)") &
4340 "KP-HFX_RI_INFO| Average number of sgf in extended RI bases:", sum(nri_per_atom)/natom
4341 WRITE (ri_data%unit_nr, fmt=
"(T3,A,F13.3,A)") &
4342 "KP-HFX_RI_INFO| Consider all image cells within a radius of ", ri_data%kp_image_range*
angstrom,
" Ang"
4343 WRITE (ri_data%unit_nr, fmt=
"(T3,A,I27/)") &
4344 "KP-HFX_RI_INFO| Number of image cells considered: ", ri_data%nimg
4348 CALL timestop(handle)
4350 END SUBROUTINE get_kp_and_ri_images
4365 SUBROUTINE get_stack_tensors(res_stack, rho_stack, ints_stack, rho_template, ints_template, &
4366 stack_size, ri_data, qs_env)
4367 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: res_stack, rho_stack, ints_stack
4368 TYPE(dbt_type),
INTENT(INOUT) :: rho_template, ints_template
4369 INTEGER,
INTENT(IN) :: stack_size
4373 INTEGER :: is, nblks, nblks_3c(3), pdims_3d(3)
4374 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_stack, dist1, &
4375 dist2, dist3, dist_stack1, &
4376 dist_stack2, dist_stack3
4377 TYPE(dbt_distribution_type) :: t_dist
4378 TYPE(dbt_pgrid_type) :: pgrid
4385 nblks =
SIZE(ri_data%bsizes_AO_split)
4386 ALLOCATE (bsizes_stack(stack_size*nblks))
4387 DO is = 1, stack_size
4388 bsizes_stack((is - 1)*nblks + 1:is*nblks) = ri_data%bsizes_AO_split(:)
4391 ALLOCATE (dist1(nblks), dist2(nblks), dist_stack1(stack_size*nblks), dist_stack2(stack_size*nblks))
4392 CALL dbt_get_info(rho_template, proc_dist_1=dist1, proc_dist_2=dist2)
4393 DO is = 1, stack_size
4394 dist_stack1((is - 1)*nblks + 1:is*nblks) = dist1(:)
4395 dist_stack2((is - 1)*nblks + 1:is*nblks) = dist2(:)
4400 CALL dbt_distribution_new(t_dist, ri_data%pgrid_2d, dist_stack1, dist_stack2)
4401 CALL dbt_create(rho_stack(1),
"RHO_stack", t_dist, [1], [2], bsizes_stack, bsizes_stack)
4402 CALL dbt_distribution_destroy(t_dist)
4403 DEALLOCATE (dist1, dist2, dist_stack1, dist_stack2)
4406 CALL create_2c_tensor(rho_stack(2), dist1, dist2, ri_data%pgrid_2d, bsizes_stack, bsizes_stack, name=
"RHO_stack")
4407 DEALLOCATE (dist1, dist2)
4409 CALL dbt_get_info(ints_template, nblks_total=nblks_3c)
4410 ALLOCATE (dist1(nblks_3c(1)), dist2(nblks_3c(2)), dist3(nblks_3c(3)))
4411 ALLOCATE (dist_stack3(stack_size*nblks_3c(3)), bsizes_ri_ext(nblks_3c(2)))
4412 CALL dbt_get_info(ints_template, proc_dist_1=dist1, proc_dist_2=dist2, &
4413 proc_dist_3=dist3, blk_size_2=bsizes_ri_ext)
4414 DO is = 1, stack_size
4415 dist_stack3((is - 1)*nblks_3c(3) + 1:is*nblks_3c(3)) = dist3(:)
4419 CALL dbt_distribution_new(t_dist, ri_data%pgrid_1, dist1, dist2, dist_stack3)
4420 CALL dbt_create(ints_stack(1),
"ints_stack", t_dist, [1, 2], [3], ri_data%bsizes_AO_split, &
4421 bsizes_ri_ext, bsizes_stack)
4422 CALL dbt_distribution_destroy(t_dist)
4423 DEALLOCATE (dist1, dist2, dist3, dist_stack3)
4427 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid, tensor_dims=[nblks_3c(1), nblks_3c(2), stack_size*nblks_3c(3)])
4428 CALL create_3c_tensor(ints_stack(2), dist1, dist2, dist3, pgrid, ri_data%bsizes_AO_split, &
4429 bsizes_ri_ext, bsizes_stack, [1, 2], [3], name=
"ints_stack")
4430 DEALLOCATE (dist1, dist2, dist3)
4431 CALL dbt_pgrid_destroy(pgrid)
4434 CALL dbt_create(ints_stack(1), res_stack(1))
4435 CALL dbt_create(ints_stack(2), res_stack(2))
4437 END SUBROUTINE get_stack_tensors
4451 SUBROUTINE fill_3c_stack(t_3c_stack, t_3c_in, images, stack_dim, ri_data, filter_at, filter_dim, &
4452 idx_to_at, img_bounds)
4453 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_stack
4454 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_in
4455 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4456 INTEGER,
INTENT(IN) :: stack_dim
4458 INTEGER,
INTENT(IN),
OPTIONAL :: filter_at, filter_dim
4459 INTEGER,
DIMENSION(:),
INTENT(INOUT),
OPTIONAL :: idx_to_at
4460 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2)
4462 INTEGER :: dest(3), i_img,
idx, ind(3), lb, nblks, &
4464 LOGICAL :: do_filter, found
4465 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4466 TYPE(dbt_iterator_type) :: iter
4471 nblks =
SIZE(ri_data%bsizes_AO_split)
4474 IF (
PRESENT(filter_at) .AND.
PRESENT(filter_dim) .AND.
PRESENT(idx_to_at)) do_filter = .true.
4479 IF (
PRESENT(img_bounds))
THEN
4481 ub = img_bounds(2) - 1
4487 IF (i_img == 0 .OR. i_img > nimg) cycle
4492 CALL dbt_iterator_start(iter, t_3c_in(i_img))
4493 DO WHILE (dbt_iterator_blocks_left(iter))
4494 CALL dbt_iterator_next_block(iter, ind)
4495 CALL dbt_get_block(t_3c_in(i_img), ind, blk, found)
4496 IF (.NOT. found) cycle
4499 IF (.NOT. idx_to_at(ind(filter_dim)) == filter_at) cycle
4502 IF (stack_dim == 1)
THEN
4503 dest = [(
idx - offset - 1)*nblks + ind(1), ind(2), ind(3)]
4504 ELSE IF (stack_dim == 2)
THEN
4505 dest = [ind(1), (
idx - offset - 1)*nblks + ind(2), ind(3)]
4507 dest = [ind(1), ind(2), (
idx - offset - 1)*nblks + ind(3)]
4510 CALL dbt_put_block(t_3c_stack, dest, shape(blk), blk)
4513 CALL dbt_iterator_stop(iter)
4516 CALL dbt_finalize(t_3c_stack)
4518 END SUBROUTINE fill_3c_stack
4530 SUBROUTINE fill_2c_stack(t_2c_stack, t_2c_in, images, stack_dim, ri_data, img_bounds, shift)
4531 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_stack
4532 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_in
4533 INTEGER,
DIMENSION(:),
INTENT(INOUT) :: images
4534 INTEGER,
INTENT(IN) :: stack_dim
4536 INTEGER,
INTENT(IN),
OPTIONAL :: img_bounds(2), shift
4538 INTEGER :: dest(2), i_img,
idx, ind(2), lb, &
4539 my_shift, nblks, nimg, offset, ub
4541 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: blk
4542 TYPE(dbt_iterator_type) :: iter
4547 nblks =
SIZE(ri_data%bsizes_AO_split)
4552 IF (
PRESENT(img_bounds))
THEN
4554 ub = img_bounds(2) - 1
4559 IF (
PRESENT(shift)) my_shift = shift
4563 IF (i_img == 0 .OR. i_img > nimg) cycle
4567 CALL dbt_iterator_start(iter, t_2c_in(i_img))
4568 DO WHILE (dbt_iterator_blocks_left(iter))
4569 CALL dbt_iterator_next_block(iter, ind)
4570 CALL dbt_get_block(t_2c_in(i_img), ind, blk, found)
4571 IF (.NOT. found) cycle
4573 IF (stack_dim == 1)
THEN
4574 dest = [(
idx - offset - 1)*nblks + ind(1), (my_shift - 1)*nblks + ind(2)]
4576 dest = [(my_shift - 1)*nblks + ind(1), (
idx - offset - 1)*nblks + ind(2)]
4579 CALL dbt_put_block(t_2c_stack, dest, shape(blk), blk)
4582 CALL dbt_iterator_stop(iter)
4585 CALL dbt_finalize(t_2c_stack)
4587 END SUBROUTINE fill_2c_stack
4595 SUBROUTINE unstack_t_3c_apc(t_3c_apc, t_stacked, idx)
4596 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_apc, t_stacked
4597 INTEGER,
INTENT(IN) ::
idx
4599 INTEGER :: current_idx
4600 INTEGER,
DIMENSION(3) :: ind, nblks_3c
4602 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4603 TYPE(dbt_iterator_type) :: iter
4606 CALL dbt_get_info(t_3c_apc, nblks_total=nblks_3c)
4609 CALL dbt_iterator_start(iter, t_stacked)
4610 DO WHILE (dbt_iterator_blocks_left(iter))
4611 CALL dbt_iterator_next_block(iter, ind)
4614 current_idx = (ind(3) - 1)/nblks_3c(3) + 1
4615 IF (.NOT.
idx == current_idx) cycle
4617 CALL dbt_get_block(t_stacked, ind, blk, found)
4618 IF (.NOT. found) cycle
4620 CALL dbt_put_block(t_3c_apc, [ind(1), ind(2), ind(3) - (
idx - 1)*nblks_3c(3)], shape(blk), blk)
4623 CALL dbt_iterator_stop(iter)
4626 END SUBROUTINE unstack_t_3c_apc
4636 SUBROUTINE get_atom_3c_ints(t_3c_at, t_3c_ints, iatom, dim_at, idx_to_at)
4637 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_at, t_3c_ints
4638 INTEGER,
INTENT(IN) :: iatom, dim_at
4639 INTEGER,
DIMENSION(:),
INTENT(IN) :: idx_to_at
4641 INTEGER,
DIMENSION(3) :: ind
4643 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: blk
4644 TYPE(dbt_iterator_type) :: iter
4647 CALL dbt_iterator_start(iter, t_3c_ints)
4648 DO WHILE (dbt_iterator_blocks_left(iter))
4649 CALL dbt_iterator_next_block(iter, ind)
4650 IF (.NOT. idx_to_at(ind(dim_at)) == iatom) cycle
4652 CALL dbt_get_block(t_3c_ints, ind, blk, found)
4653 IF (.NOT. found) cycle
4655 CALL dbt_put_block(t_3c_at, ind, shape(blk), blk)
4658 CALL dbt_iterator_stop(iter)
4660 CALL dbt_finalize(t_3c_at)
4662 END SUBROUTINE get_atom_3c_ints
4673 SUBROUTINE precalc_derivatives(t_3c_der_RI, t_3c_der_AO, mat_der_pot, t_2c_der_metric, ri_data, qs_env)
4674 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_der_ri, t_3c_der_ao
4675 TYPE(dbcsr_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_der_pot
4676 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_2c_der_metric
4680 CHARACTER(LEN=*),
PARAMETER :: routinen =
'precalc_derivatives'
4682 INTEGER :: handle, handle2, i_img, i_mem, i_ri, &
4683 i_xyz, iatom, n_mem, natom, nblks_ri, &
4684 ncell_ri, nimg, nkind, nthreads
4685 INTEGER(int_8) :: nze
4686 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: bsizes_ri_ext, bsizes_ri_ext_split, dist_ao_1, &
4687 dist_ao_2, dist_ri, dist_ri_ext, dummy_end, dummy_start, end_blocks, start_blocks
4688 INTEGER,
DIMENSION(3) :: pcoord, pdims
4689 INTEGER,
DIMENSION(:),
POINTER :: col_bsize, row_bsize
4691 TYPE(dbcsr_distribution_type) :: dbcsr_dist
4692 TYPE(dbcsr_type) :: dbcsr_template
4693 TYPE(dbcsr_type),
ALLOCATABLE,
DIMENSION(:, :) :: mat_der_metric
4694 TYPE(dbt_distribution_type) :: t_dist
4695 TYPE(dbt_pgrid_type) :: pgrid
4696 TYPE(dbt_type) :: t_3c_template
4697 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
4702 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri
4709 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
4711 NULLIFY (qs_kind_set, dist_2d, nl_2c, particle_set, dft_control, para_env, row_bsize, col_bsize)
4713 CALL timeset(routinen, handle)
4715 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, natom=natom, &
4716 particle_set=particle_set, dft_control=dft_control, para_env=para_env)
4719 ncell_ri = ri_data%ncell_RI
4721 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
4731 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[max(1, natom/(ri_data%n_mem*nthreads)), natom, natom])
4733 CALL create_3c_tensor(t_3c_template, dist_ao_1, dist_ao_2, dist_ri, pgrid, &
4734 ri_data%bsizes_AO, ri_data%bsizes_AO, ri_data%bsizes_RI, &
4735 map1=[1, 2], map2=[3], name=
"tmp")
4736 CALL dbt_destroy(t_3c_template)
4739 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
4740 ALLOCATE (dist_ri_ext(natom*ncell_ri))
4741 ALLOCATE (bsizes_ri_ext(natom*ncell_ri))
4742 ALLOCATE (bsizes_ri_ext_split(nblks_ri*ncell_ri))
4743 DO i_ri = 1, ncell_ri
4744 bsizes_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = ri_data%bsizes_RI(:)
4745 dist_ri_ext((i_ri - 1)*natom + 1:i_ri*natom) = dist_ri(:)
4746 bsizes_ri_ext_split((i_ri - 1)*nblks_ri + 1:i_ri*nblks_ri) = ri_data%bsizes_RI_split(:)
4749 CALL dbt_distribution_new(t_dist, pgrid, dist_ao_1, dist_ao_2, dist_ri_ext)
4750 CALL dbt_create(t_3c_template,
"KP_3c_der", t_dist, [1, 2], [3], &
4751 ri_data%bsizes_AO, ri_data%bsizes_AO, bsizes_ri_ext)
4752 CALL dbt_distribution_destroy(t_dist)
4754 ALLOCATE (t_3c_der_ri_prv(nimg, 1, 3), t_3c_der_ao_prv(nimg, 1, 3))
4757 CALL dbt_create(t_3c_template, t_3c_der_ri_prv(i_img, 1, i_xyz))
4758 CALL dbt_create(t_3c_template, t_3c_der_ao_prv(i_img, 1, i_xyz))
4761 CALL dbt_destroy(t_3c_template)
4763 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
4764 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
4766 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
4767 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
4768 CALL dbt_pgrid_destroy(pgrid)
4771 "HFX_3c_nl", qs_env, op_pos=2, sym_jk=.false., own_dist=.true.)
4773 n_mem = ri_data%n_mem
4775 start_blocks, end_blocks)
4776 DEALLOCATE (dummy_start, dummy_end)
4778 CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, ri_data%pgrid_2, &
4779 bsizes_ri_ext_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
4780 map1=[1], map2=[2, 3], name=
"der (RI | AO AO)")
4783 CALL dbt_create(t_3c_template, t_3c_der_ri(i_img, i_xyz))
4784 CALL dbt_create(t_3c_template, t_3c_der_ao(i_img, i_xyz))
4790 nl_3c, basis_set_ao, basis_set_ao, basis_set_ri, &
4791 ri_data%ri_metric, der_eps=ri_data%eps_schwarz_forces, op_pos=2, &
4792 do_kpoints=.true., do_hfx_kpoints=.true., &
4793 bounds_k=[start_blocks(i_mem), end_blocks(i_mem)], &
4794 ri_range=ri_data%kp_RI_range, img_to_ri_cell=ri_data%img_to_RI_cell)
4796 CALL timeset(routinen//
"_cpy", handle2)
4803 CALL dbt_copy(t_3c_der_ao_prv(i_img, 1, i_xyz), t_3c_template, &
4804 order=[3, 2, 1], move_data=.true.)
4805 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4806 CALL dbt_copy(t_3c_template, t_3c_der_ao(i_img, i_xyz), &
4807 order=[1, 3, 2], move_data=.true., summation=.true.)
4813 CALL dbt_copy(t_3c_der_ri_prv(i_img, 1, i_xyz), t_3c_template, &
4814 order=[3, 2, 1], move_data=.true.)
4815 CALL dbt_filter(t_3c_template, ri_data%filter_eps)
4816 CALL dbt_copy(t_3c_template, t_3c_der_ri(i_img, i_xyz), &
4817 order=[1, 3, 2], move_data=.true., summation=.true.)
4821 CALL timestop(handle2)
4823 CALL dbt_destroy(t_3c_template)
4828 CALL dbt_destroy(t_3c_der_ri_prv(i_img, 1, i_xyz))
4829 CALL dbt_destroy(t_3c_der_ao_prv(i_img, 1, i_xyz))
4832 DEALLOCATE (t_3c_der_ri_prv, t_3c_der_ao_prv)
4835 CALL reorder_3c_derivs(t_3c_der_ri, ri_data)
4836 CALL reorder_3c_derivs(t_3c_der_ao, ri_data)
4838 CALL timeset(routinen//
"_2c", handle2)
4841 ALLOCATE (row_bsize(
SIZE(ri_data%bsizes_RI)))
4842 ALLOCATE (col_bsize(
SIZE(ri_data%bsizes_RI)))
4843 row_bsize(:) = ri_data%bsizes_RI
4844 col_bsize(:) = ri_data%bsizes_RI
4846 CALL dbcsr_create(dbcsr_template,
"2c_der", dbcsr_dist, dbcsr_type_no_symmetry, &
4847 row_bsize, col_bsize)
4848 CALL dbcsr_distribution_release(dbcsr_dist)
4849 DEALLOCATE (col_bsize, row_bsize)
4851 ALLOCATE (mat_der_metric(nimg, 3))
4854 CALL dbcsr_create(mat_der_pot(i_img, i_xyz), template=dbcsr_template)
4855 CALL dbcsr_create(mat_der_metric(i_img, i_xyz), template=dbcsr_template)
4858 CALL dbcsr_release(dbcsr_template)
4862 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
4864 basis_set_ri, basis_set_ri, ri_data%hfx_pot, do_kpoints=.true.)
4869 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
4871 basis_set_ri, basis_set_ri, ri_data%ri_metric, do_kpoints=.true.)
4877 CALL dbt_create(ri_data%t_2c_inv(1, 1), t_2c_der_metric(iatom, i_xyz))
4878 CALL get_ext_2c_int(t_2c_der_metric(iatom, i_xyz), mat_der_metric(:, i_xyz), &
4879 iatom, iatom, 1, ri_data, qs_env)
4882 CALL dbcsr_release(mat_der_metric(i_img, i_xyz))
4885 CALL timestop(handle2)
4887 CALL timestop(handle)
4889 END SUBROUTINE precalc_derivatives
4910 SUBROUTINE get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, img, pref, &
4911 ri_data, qs_env, work_virial, cell, particle_set, diag, offdiag)
4914 TYPE(dbt_type),
INTENT(INOUT) :: t_2c_contr
4915 TYPE(dbt_type),
DIMENSION(:),
INTENT(INOUT) :: t_2c_der
4916 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
4917 INTEGER,
INTENT(IN) :: img
4918 REAL(
dp),
INTENT(IN) :: pref
4921 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
4922 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
4924 POINTER :: particle_set
4925 LOGICAL,
INTENT(IN),
OPTIONAL ::
diag, offdiag
4927 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_2c_der_force'
4929 INTEGER :: handle, i_img, i_ri, i_xyz, iat, &
4930 iat_of_kind, ikind, j_img, j_ri, &
4931 j_xyz, jat, jat_of_kind, jkind, natom
4932 INTEGER,
DIMENSION(2) :: ind
4933 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
4934 LOGICAL :: found, my_diag, my_offdiag, use_virial
4935 REAL(
dp) :: new_force
4936 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :),
TARGET :: contr_blk, der_blk
4937 REAL(
dp),
DIMENSION(3) :: scoord
4938 TYPE(dbt_iterator_type) :: iter
4941 NULLIFY (kpoints, index_to_cell)
4946 CALL timeset(routinen, handle)
4948 use_virial = .false.
4949 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
4954 my_offdiag = .false.
4955 IF (
PRESENT(
diag)) my_offdiag = offdiag
4957 CALL get_qs_env(qs_env, kpoints=kpoints, natom=natom)
4966 CALL dbt_iterator_start(iter, t_2c_der(i_xyz))
4967 DO WHILE (dbt_iterator_blocks_left(iter))
4968 CALL dbt_iterator_next_block(iter, ind)
4971 IF ((my_diag .AND. .NOT. my_offdiag) .OR. (.NOT. my_diag .AND. my_offdiag))
THEN
4972 IF (my_diag .AND. (ind(1) .NE. ind(2))) cycle
4973 IF (my_offdiag .AND. (ind(1) == ind(2))) cycle
4976 CALL dbt_get_block(t_2c_der(i_xyz), ind, der_blk, found)
4978 CALL dbt_get_block(t_2c_contr, ind, contr_blk, found)
4984 new_force = pref*sum(der_blk(:, :)*contr_blk(:, :))
4986 i_ri = (ind(1) - 1)/natom + 1
4987 i_img = ri_data%RI_cell_to_img(i_ri)
4988 iat = ind(1) - (i_ri - 1)*natom
4989 iat_of_kind = atom_of_kind(iat)
4990 ikind = kind_of(iat)
4992 j_ri = (ind(2) - 1)/natom + 1
4993 j_img = ri_data%RI_cell_to_img(j_ri)
4994 jat = ind(2) - (j_ri - 1)*natom
4995 jat_of_kind = atom_of_kind(jat)
4996 jkind = kind_of(jat)
5000 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5003 IF (use_virial)
THEN
5006 scoord(:) = scoord(:) + real(index_to_cell(:, i_img),
dp)
5010 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5016 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5019 IF (use_virial)
THEN
5022 scoord(:) = scoord(:) + real(index_to_cell(:, j_img) + index_to_cell(:, img),
dp)
5026 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) - new_force*scoord(j_xyz)
5030 DEALLOCATE (contr_blk)
5033 DEALLOCATE (der_blk)
5035 CALL dbt_iterator_stop(iter)
5039 CALL timestop(handle)
5065 idx_to_at_RI, idx_to_at_AO, i_images, lb_img, pref, &
5066 ri_data, qs_env, work_virial, cell, particle_set)
5069 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_contr
5070 TYPE(dbt_type),
DIMENSION(3),
INTENT(INOUT) :: t_3c_der_1, t_3c_der_2
5071 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of, idx_to_at_ri, &
5072 idx_to_at_ao, i_images
5073 INTEGER,
INTENT(IN) :: lb_img
5074 REAL(
dp),
INTENT(IN) :: pref
5077 REAL(
dp),
DIMENSION(3, 3),
INTENT(INOUT),
OPTIONAL :: work_virial
5078 TYPE(
cell_type),
OPTIONAL,
POINTER :: cell
5080 POINTER :: particle_set
5082 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_force_from_3c_trace'
5084 INTEGER :: handle, i_ri, i_xyz, iat, iat_of_kind,
idx, ikind, j_xyz, jat, jat_of_kind, &
5085 jkind, kat, kat_of_kind, kkind, nblks_ao, nblks_ri, ri_img
5086 INTEGER,
DIMENSION(3) :: ind
5087 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
5088 LOGICAL :: found, found_1, found_2, use_virial
5089 REAL(
dp) :: new_force
5090 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :, :),
TARGET :: contr_blk, der_blk_1, der_blk_2, &
5092 REAL(
dp),
DIMENSION(3) :: scoord
5093 TYPE(dbt_iterator_type) :: iter
5096 NULLIFY (kpoints, index_to_cell)
5098 CALL timeset(routinen, handle)
5103 nblks_ri =
SIZE(ri_data%bsizes_RI_split)
5104 nblks_ao =
SIZE(ri_data%bsizes_AO_split)
5106 use_virial = .false.
5107 IF (
PRESENT(work_virial) .AND.
PRESENT(cell) .AND.
PRESENT(particle_set)) use_virial = .true.
5114 CALL dbt_iterator_start(iter, t_3c_contr)
5115 DO WHILE (dbt_iterator_blocks_left(iter))
5116 CALL dbt_iterator_next_block(iter, ind)
5118 CALL dbt_get_block(t_3c_contr, ind, contr_blk, found)
5122 CALL dbt_get_block(t_3c_der_1(i_xyz), ind, der_blk_1, found_1)
5123 IF (.NOT. found_1)
THEN
5124 DEALLOCATE (der_blk_1)
5125 ALLOCATE (der_blk_1(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5126 der_blk_1(:, :, :) = 0.0_dp
5128 CALL dbt_get_block(t_3c_der_2(i_xyz), ind, der_blk_2, found_2)
5129 IF (.NOT. found_2)
THEN
5130 DEALLOCATE (der_blk_2)
5131 ALLOCATE (der_blk_2(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5132 der_blk_2(:, :, :) = 0.0_dp
5135 ALLOCATE (der_blk_3(
SIZE(contr_blk, 1),
SIZE(contr_blk, 2),
SIZE(contr_blk, 3)))
5136 der_blk_3(:, :, :) = -(der_blk_1(:, :, :) + der_blk_2(:, :, :))
5142 new_force = pref*sum(der_blk_1(:, :, :)*contr_blk(:, :, :))
5144 i_ri = (ind(1) - 1)/nblks_ri + 1
5145 ri_img = ri_data%RI_cell_to_img(i_ri)
5146 iat = idx_to_at_ri(ind(1) - (i_ri - 1)*nblks_ri)
5147 iat_of_kind = atom_of_kind(iat)
5148 ikind = kind_of(iat)
5151 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
5154 IF (use_virial)
THEN
5157 scoord(:) = scoord(:) + real(index_to_cell(:, ri_img),
dp)
5161 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5166 new_force = pref*sum(der_blk_2(:, :, :)*contr_blk(:, :, :))
5167 jat = idx_to_at_ao(ind(2))
5168 jat_of_kind = atom_of_kind(jat)
5169 jkind = kind_of(jat)
5172 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
5175 IF (use_virial)
THEN
5181 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5187 new_force = pref*sum(der_blk_3(:, :, :)*contr_blk(:, :, :))
5188 idx = (ind(3) - 1)/nblks_ao + 1
5189 kat = idx_to_at_ao(ind(3) - (
idx - 1)*nblks_ao)
5190 kat_of_kind = atom_of_kind(kat)
5191 kkind = kind_of(kat)
5194 force(kkind)%fock_4c(i_xyz, kat_of_kind) = force(kkind)%fock_4c(i_xyz, kat_of_kind) &
5197 IF (use_virial)
THEN
5199 scoord(:) = scoord(:) + real(index_to_cell(:, i_images(lb_img - 1 +
idx)),
dp)
5203 work_virial(i_xyz, j_xyz) = work_virial(i_xyz, j_xyz) + new_force*scoord(j_xyz)
5207 DEALLOCATE (der_blk_1, der_blk_2, der_blk_3)
5209 DEALLOCATE (contr_blk)
5212 CALL dbt_iterator_stop(iter)
5214 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)
...
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...
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, upper_to_full)
used to replace the cholesky decomposition by the inverse
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_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, 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, 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.