Skip to content

Commit

Permalink
Low-scaling RPA/GW: Reduce memory bottleneck of 3-center integral cal…
Browse files Browse the repository at this point in the history 
…culation
  • Loading branch information
@pseewald
pseewald committed Feb 24, 2020
1 parent 24abbe5 commit b5fa45f
Show file tree
Hide file tree
Showing 6 changed files with 108 additions and 51 deletions.
10 changes: 5 additions & 5 deletions src/mp2_gpw.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp
TYPE(dbcsr_t_type) :: t_3c_M
TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O, t_3c_overl_int
TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_O
TYPE(dbcsr_type), POINTER :: mo_coeff_all, mo_coeff_all_beta, mo_coeff_gw, mo_coeff_gw_beta, &
mo_coeff_o, mo_coeff_o_beta, mo_coeff_v, mo_coeff_v_beta
TYPE(dft_control_type), POINTER :: dft_control
Expand Down Expand Up @@ -402,7 +402,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
do_bse, gd_B_all, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, &
gw_corr_lev_occ, gw_corr_lev_virt, &
do_im_time, do_kpoints_cubic_RPA, kpoints, &
t_3c_overl_int, t_3c_M, t_3c_O, &
t_3c_M, t_3c_O, &
mp2_env%ri_metric, &
gd_B_occ_bse, gd_B_virt_bse, &
BIb_C_beta, BIb_C_gw_beta, gd_B_virtual_beta, &
Expand All @@ -422,7 +422,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
starts_array_mc_block, ends_array_mc_block, &
gw_corr_lev_occ, gw_corr_lev_virt, &
do_im_time, do_kpoints_cubic_RPA, kpoints, &
t_3c_overl_int, t_3c_M, t_3c_O, &
t_3c_M, t_3c_O, &
mp2_env%ri_metric, gd_B_occ_bse, gd_B_virt_bse)
END IF
Expand Down Expand Up @@ -549,7 +549,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
unit_nr, my_do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, &
mat_munu, &
mat_P_global, &
t_3c_overl_int, t_3c_M, t_3c_O, &
t_3c_M, t_3c_O, &
starts_array_mc, ends_array_mc, &
starts_array_mc_block, ends_array_mc_block, &
mp2_env%mp2_gpw%eps_filter, BIb_C_beta, homo_beta, Eigenval_beta, &
Expand All @@ -569,7 +569,7 @@ SUBROUTINE mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T
unit_nr, my_do_ri_sos_laplace_mp2, my_do_gw, do_im_time, do_bse, matrix_s, &
mat_munu, &
mat_P_global, &
t_3c_overl_int, t_3c_M, t_3c_O, &
t_3c_M, t_3c_O, &
starts_array_mc, ends_array_mc, &
starts_array_mc_block, ends_array_mc_block, &
mp2_env%mp2_gpw%eps_filter)
Expand Down
76 changes: 49 additions & 27 deletions src/mp2_integrals.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ MODULE mp2_integrals
dbcsr_release_p, dbcsr_scalar, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
dbcsr_type_real_8
USE dbcsr_tensor_api, ONLY: &
dbcsr_t_contract, dbcsr_t_copy, dbcsr_t_create, dbcsr_t_destroy, &
dbcsr_t_clear, dbcsr_t_contract, dbcsr_t_copy, dbcsr_t_create, dbcsr_t_destroy, &
dbcsr_t_distribution_destroy, dbcsr_t_distribution_new, dbcsr_t_distribution_type, &
dbcsr_t_filter, dbcsr_t_get_block, dbcsr_t_get_stored_coordinates, &
dbcsr_t_mp_environ_pgrid, dbcsr_t_pgrid_create, dbcsr_t_pgrid_create_expert, &
Expand Down Expand Up @@ -153,7 +153,6 @@ MODULE mp2_integrals
!> \param do_im_time ...
!> \param do_kpoints_cubic_RPA ...
!> \param kpoints ...
!> \param t_3c_overl_int ...
!> \param t_3c_M ...
!> \param t_3c_O ...
!> \param ri_metric ...
Expand Down Expand Up @@ -181,7 +180,7 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
starts_array_mc_block, ends_array_mc_block, &
gw_corr_lev_occ, gw_corr_lev_virt, &
do_im_time, do_kpoints_cubic_RPA, kpoints, &
t_3c_overl_int, t_3c_M, t_3c_O, &
t_3c_M, t_3c_O, &
ri_metric, gd_B_occ_bse, gd_B_virt_bse, BIb_C_beta, BIb_C_gw_beta, &
gd_B_virtual_beta, homo_beta, mo_coeff_o_beta, mo_coeff_v_beta, &
mo_coeff_all_beta, mo_coeff_gw_beta)
Expand Down Expand Up @@ -219,7 +218,6 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
INTEGER, INTENT(IN) :: gw_corr_lev_occ, gw_corr_lev_virt
LOGICAL, INTENT(IN) :: do_im_time, do_kpoints_cubic_RPA
TYPE(kpoint_type), POINTER :: kpoints
TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_overl_int
TYPE(dbcsr_t_type), INTENT(OUT) :: t_3c_M
TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :), &
INTENT(OUT) :: t_3c_O
Expand All @@ -235,19 +233,19 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
CHARACTER(LEN=*), PARAMETER :: routineN = 'mp2_ri_gpw_compute_in', &
routineP = moduleN//':'//routineN

INTEGER :: cut_memory, eri_method, gw_corr_lev_total, handle, handle2, handle4, i, &
INTEGER :: cm, cut_memory, eri_method, gw_corr_lev_total, handle, handle2, handle4, i, &
i_counter, itmp(2), j, LLL, max_row_col_local, max_row_col_local_beta, &
max_row_col_local_gw, max_row_col_local_occ_bse, max_row_col_local_virt_bse, min_bsize, &
mp_comm_t3c_2, my_B_all_end, my_B_all_size, my_B_all_start, my_B_occ_bse_end, &
my_B_occ_bse_size, my_B_occ_bse_start, my_B_size, my_B_size_beta, my_B_virt_bse_end, &
my_B_virt_bse_size, my_B_virt_bse_start, my_B_virtual_end, my_B_virtual_end_beta, &
my_B_virtual_start, my_B_virtual_start_beta, my_group_L_end, my_group_L_size, &
my_group_L_start, natom, ngroup
INTEGER :: ngroup_im_time, ngroup_im_time_P, nimg, nkind, num_small_eigen, potential_type, &
ri_metric_type, virtual, virtual_beta
INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI, sizes_AO, &
sizes_AO_split, sizes_RI, &
sizes_RI_split, sub_proc_map
my_group_L_start, natom
INTEGER :: ngroup, ngroup_im_time, ngroup_im_time_P, nimg, nkind, num_small_eigen, &
potential_type, ri_metric_type, virtual, virtual_beta
INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI, &
ends_array_mc_block_int, ends_array_mc_int, sizes_AO, sizes_AO_split, sizes_RI, &
sizes_RI_split, starts_array_mc_block_int, starts_array_mc_int, sub_proc_map
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: local_col_row_info, local_col_row_info_beta, &
local_col_row_info_gw, local_col_row_info_occ_bse, local_col_row_info_virt_bse
INTEGER, DIMENSION(3) :: pcoord, pdims, pdims_t3c, periodic
Expand All @@ -267,6 +265,7 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_munu_local_L
TYPE(dbcsr_t_pgrid_type) :: pgrid_t3c_M, pgrid_t3c_overl
TYPE(dbcsr_t_type) :: t_3c_overl_int_template
TYPE(dbcsr_t_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_overl_int
TYPE(dbcsr_type) :: matrix_bse_ab, matrix_bse_anu, matrix_bse_ij, matrix_bse_inu, &
matrix_ia_jb, matrix_ia_jb_beta, matrix_ia_jnu, matrix_ia_jnu_beta, matrix_in_jm, &
matrix_in_jm_beta, matrix_in_jnu, matrix_in_jnu_beta
Expand Down Expand Up @@ -760,8 +759,15 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_ao)
CALL create_contraction_batches(sizes_AO, cut_memory, starts_array_mc_int, ends_array_mc_int, &
starts_array_mc_block_int, ends_array_mc_block_int)
DEALLOCATE (starts_array_mc_int, ends_array_mc_int)
CALL create_3c_tensor(t_3c_overl_int_template, dist_RI, dist_AO_1, dist_AO_2, pgrid_t3c_overl, &
sizes_RI, sizes_AO, sizes_AO, map1=[1, 2], map2=[3], name="O (RI AO | AO)")
sizes_RI, sizes_AO, sizes_AO, map1=[1, 2], map2=[3], &
starts_array_block_2=starts_array_mc_block_int, &
ends_array_block_2=ends_array_mc_block_int, &
name="O (RI AO | AO)")
CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set)
CALL dbcsr_t_mp_environ_pgrid(pgrid_t3c_overl, pdims, pcoord)
Expand Down Expand Up @@ -841,27 +847,43 @@ SUBROUTINE mp2_ri_gpw_compute_in(BIb_C, BIb_C_gw, BIb_C_bse_ij, BIb_C_bse_ab, gd
ENDDO
ENDDO
CALL build_3c_integrals(t_3c_overl_int, &
qs_env%mp2_env%mp2_gpw%eps_filter, &
qs_env, &
nl_3c, &
basis_i=basis_set_ri_aux, &
basis_j=basis_set_ao, basis_k=basis_set_ao, &
potential_parameter=ri_metric, &
do_kpoints=do_kpoints_cubic_RPA)
! build integrals in batches and copy to optimized format
! note: integrals are stored in terms of atomic blocks. To avoid a memory bottleneck,
! integrals are calculated in batches and copied to optimized format with subatomic blocks
DO cm = 1, cut_memory
CALL build_3c_integrals(t_3c_overl_int, &
qs_env%mp2_env%mp2_gpw%eps_filter, &
qs_env, &
nl_3c, &
basis_i=basis_set_ri_aux, &
basis_j=basis_set_ao, basis_k=basis_set_ao, &
potential_parameter=ri_metric, &
do_kpoints=do_kpoints_cubic_RPA, &
bounds_j=[starts_array_mc_block_int(cm), ends_array_mc_block_int(cm)])
CALL timeset(routineN//"_copy_3c", handle4)
! copy integral tensor t_3c_overl_int to t_3c_O tensor optimized for contraction
DO i = 1, SIZE(t_3c_overl_int, 1)
DO j = 1, SIZE(t_3c_overl_int, 2)
CALL dbcsr_t_copy(t_3c_overl_int(i, j), t_3c_O(i, j), order=[1, 3, 2], &
summation=.TRUE., move_data=.TRUE.)
CALL dbcsr_t_clear(t_3c_overl_int(i, j))
CALL dbcsr_t_filter(t_3c_O(i, j), qs_env%mp2_env%mp2_gpw%eps_filter)
ENDDO
ENDDO
CALL timestop(handle4)
ENDDO
DEALLOCATE (basis_set_ri_aux, basis_set_ao)
! copy integral tensor t_3c_overl_int to t_3c_O tensor optimized for contraction
CALL timeset(routineN//"_copy_3c", handle4)
DO i = 1, SIZE(t_3c_overl_int, 1)
DO j = 1, SIZE(t_3c_overl_int, 2)
CALL dbcsr_t_copy(t_3c_overl_int(i, j), t_3c_O(i, j), order=[1, 3, 2], move_data=.FALSE.)
CALL dbcsr_t_filter(t_3c_O(i, j), qs_env%mp2_env%mp2_gpw%eps_filter)
CALL dbcsr_t_destroy(t_3c_overl_int(i, j))
ENDDO
ENDDO
CALL timestop(handle4)
DEALLOCATE (t_3c_overl_int)
DEALLOCATE (basis_set_ri_aux, basis_set_ao)
CALL neighbor_list_3c_destroy(nl_3c)
Expand Down
45 changes: 43 additions & 2 deletions src/qs_tensors.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -390,9 +390,31 @@ SUBROUTINE neighbor_list_3c_iterator_create(iterator, ijk_nl)
iterator%iter_level = 0
iterator%ijk_nl = ijk_nl

iterator%bounds_i = 0
iterator%bounds_j = 0
iterator%bounds_k = 0

CALL timestop(handle)
END SUBROUTINE

! **************************************************************************************************
!> \brief impose atomic upper and lower bounds
!> \param iterator ...
!> \param bounds_i ...
!> \param bounds_j ...
!> \param bounds_k ...
! **************************************************************************************************
SUBROUTINE nl_3c_iter_set_bounds(iterator, bounds_i, bounds_j, bounds_k)
TYPE(neighbor_list_3c_iterator_type), &
INTENT(INOUT) :: iterator
INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k

IF (PRESENT(bounds_i)) iterator%bounds_i = bounds_i
IF (PRESENT(bounds_j)) iterator%bounds_j = bounds_j
IF (PRESENT(bounds_k)) iterator%bounds_k = bounds_k

END SUBROUTINE

! **************************************************************************************************
!> \brief Destroy 3c-nl iterator
!> \param iterator ...
Expand Down Expand Up @@ -458,7 +480,19 @@ RECURSIVE FUNCTION neighbor_list_3c_iterate(iterator) RESULT(iter_stat)
CPASSERT(jatom_1 == jatom_2)
jatom = jatom_1

skip_this = .TRUE.
skip_this = .FALSE.
IF ((iterator%bounds_i(1) > 0 .AND. iatom < iterator%bounds_i(1)) &
.OR. (iterator%bounds_i(2) > 0 .AND. iatom > iterator%bounds_i(2))) skip_this = .TRUE.
IF ((iterator%bounds_j(1) > 0 .AND. jatom < iterator%bounds_j(1)) &
.OR. (iterator%bounds_j(2) > 0 .AND. jatom > iterator%bounds_j(2))) skip_this = .TRUE.
IF ((iterator%bounds_k(1) > 0 .AND. katom < iterator%bounds_k(1)) &
.OR. (iterator%bounds_k(2) > 0 .AND. katom > iterator%bounds_k(2))) skip_this = .TRUE.

IF (skip_this) THEN
iter_stat = neighbor_list_3c_iterate(iterator)
RETURN
ENDIF

SELECT CASE (iterator%ijk_nl%sym)
CASE (symmetric_none)
skip_this = .FALSE.
Expand Down Expand Up @@ -778,12 +812,16 @@ SUBROUTINE alloc_block_3c(t3c, nl_3c, basis_i, basis_j, basis_k, qs_env, potenti
!> \param do_kpoints ...
!> this routine requires that libint has been static initialised somewhere else
!> \param desymmetrize ...
!> \param bounds_i ...
!> \param bounds_j ...
!> \param bounds_k ...
! **************************************************************************************************
SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
nl_3c, basis_i, basis_j, basis_k, &
potential_parameter, &
int_eps, &
op_pos, do_kpoints, desymmetrize)
op_pos, do_kpoints, desymmetrize, &
bounds_i, bounds_j, bounds_k)
TYPE(dbcsr_t_type), DIMENSION(:, :), INTENT(INOUT) :: t3c
REAL(KIND=dp), INTENT(IN) :: filter_eps
TYPE(qs_environment_type), POINTER :: qs_env
Expand All @@ -793,6 +831,7 @@ SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
REAL(KIND=dp), INTENT(IN), OPTIONAL :: int_eps
INTEGER, INTENT(IN), OPTIONAL :: op_pos
LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints, desymmetrize
INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: bounds_i, bounds_j, bounds_k

CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integrals', &
routineP = moduleN//':'//routineN
Expand Down Expand Up @@ -936,6 +975,8 @@ SUBROUTINE build_3c_integrals(t3c, filter_eps, qs_env, &
ENDIF

CALL neighbor_list_3c_iterator_create(nl_3c_iter, nl_3c)
CALL nl_3c_iter_set_bounds(nl_3c_iter, bounds_i, bounds_j, bounds_k)

DO WHILE (neighbor_list_3c_iterate(nl_3c_iter) == 0)
CALL get_3c_iterator_info(nl_3c_iter, ikind=ikind, jkind=jkind, kkind=kkind, &
iatom=iatom, jatom=jatom, katom=katom, &
Expand Down
1 change: 1 addition & 0 deletions src/qs_tensors_types.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@ MODULE qs_tensors_types
TYPE(neighbor_list_iterator_p_type), DIMENSION(:), POINTER :: iter_jk => NULL()
INTEGER :: iter_level
TYPE(neighbor_list_3c_type) :: ijk_nl
INTEGER, DIMENSION(2) :: bounds_i, bounds_j, bounds_k
END TYPE

CONTAINS
Expand Down

0 comments on commit b5fa45f

Please sign in to comment.