diff --git a/src/gw_integrals.F b/src/gw_integrals.F index ffa69da351..ec0ad15614 100644 --- a/src/gw_integrals.F +++ b/src/gw_integrals.F @@ -20,7 +20,6 @@ MODULE gw_integrals get_cell,& pbc USE cp_array_utils, ONLY: cp_2d_r_p_type - USE cp_control_types, ONLY: dft_control_type USE cp_files, ONLY: close_file,& open_file USE gamma, ONLY: init_md_ftable @@ -41,7 +40,6 @@ MODULE gw_integrals USE particle_types, ONLY: particle_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type - USE qs_kind_types, ONLY: qs_kind_type USE t_c_g0, ONLY: get_lmax_init,& init @@ -54,12 +52,393 @@ MODULE gw_integrals CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_integrals' - PUBLIC :: build_3c_integral_block + PUBLIC :: build_3c_integral_block, build_3c_integral_block_ctx, & + gw_3c_context_type, gw_3c_context_create, gw_3c_context_release, & + gw_3c_workspace_type, gw_3c_workspace_create, gw_3c_workspace_release + +! ************************************************************************************************** +!> \brief Shared read-only context for repeated 3-center integral block builds: screening +!> parameters, basis maxima, contracted sphi tables, and the one-time gamma / +!> truncated-Coulomb table initializations. Create and release OUTSIDE any OMP parallel +!> region; creation is MPI-collective when the potential is truncated. +! ************************************************************************************************** + TYPE gw_3c_context_type + TYPE(libint_potential_type) :: potential_parameter + INTEGER :: op_ij = do_potential_id, & + op_jk = do_potential_id + REAL(KIND=dp) :: dr_ij = 0.0_dp, dr_jk = 0.0_dp, & + dr_ik = 0.0_dp + INTEGER :: maxli = 0, maxlj = 0, maxlk = 0, & + max_am = 0, m_max = 0 + INTEGER :: max_ncoi = 0, max_ncoj = 0, max_ncok = 0 + INTEGER :: max_nsgfi = 0, max_nsgfj = 0, & + max_nsgfk = 0, max_nset = 0, natom = 0 + TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi => NULL(), tspj => NULL(), & + spk => NULL() + TYPE(gto_basis_set_p_type), DIMENSION(:), ALLOCATABLE :: basis_i, basis_j, basis_k + INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set => NULL() + TYPE(cell_type), POINTER :: cell => NULL() + REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp + END TYPE gw_3c_context_type + +! ************************************************************************************************** +!> \brief Per-thread workspace for 3-center integral block builds: libint object + contraction +!> buffers. Each thread creates its own (inside the parallel region is fine). +! ************************************************************************************************** + TYPE gw_3c_workspace_type + TYPE(cp_libint_t) :: lib + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpp_buffer, ccp_buffer + END TYPE gw_3c_workspace_type CONTAINS ! ************************************************************************************************** -!> \brief ... +!> \brief Builds the shared 3c-integral context: screening radii, basis maxima, contracted sphi +!> tables, and the one-time gamma / truncated-Coulomb table initializations. +!> \param ctx ... +!> \param qs_env ... +!> \param potential_parameter ... +!> \param basis_j ... +!> \param basis_k ... +!> \param basis_i ... +! ************************************************************************************************** + SUBROUTINE gw_3c_context_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i) + + TYPE(gw_3c_context_type), INTENT(OUT) :: ctx + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(libint_potential_type), INTENT(IN) :: potential_parameter + TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_j, basis_k, basis_i + + CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_3c_context_create' + + INTEGER :: egfi, handle, ibasis, ilist, imax, iset, & + jset, kset, nbasis, ncoi, sgfi, unit_id + INTEGER, DIMENSION(:), POINTER :: npgfi, npgfj, npgfk, nsgfi, nsgfj, nsgfk + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(gto_basis_set_type), POINTER :: basis_set + TYPE(mp_para_env_type), POINTER :: para_env + + CALL timeset(routineN, handle) + + ctx%potential_parameter = potential_parameter + ctx%op_ij = potential_parameter%potential_type + ctx%op_jk = do_potential_id + + IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_ij == do_potential_short) THEN + ctx%dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor + ctx%dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor + ELSEIF (ctx%op_ij == do_potential_coulomb) THEN + ctx%dr_ij = 1000000.0_dp + ctx%dr_ik = 1000000.0_dp + END IF + + NULLIFY (atomic_kind_set, para_env) + CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, cell=ctx%cell, natom=ctx%natom, & + para_env=para_env, particle_set=ctx%particle_set) + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=ctx%kind_of) + CALL get_cell(cell=ctx%cell, h=ctx%hmat) + + ctx%basis_i = basis_i + ctx%basis_j = basis_j + ctx%basis_k = basis_k + + ! max l per basis for libint; max nset/nco/nsgf for the LIBXSMM contraction buffers + nbasis = SIZE(basis_i) + DO ibasis = 1, nbasis + CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, & + nset=iset, nsgf_set=nsgfi, npgf=npgfi) + ctx%maxli = MAX(ctx%maxli, imax) + ctx%max_nset = MAX(ctx%max_nset, iset) + ctx%max_nsgfi = MAX(ctx%max_nsgfi, MAXVAL(nsgfi)) + ctx%max_ncoi = MAX(ctx%max_ncoi, MAXVAL(npgfi)*ncoset(ctx%maxli)) + END DO + DO ibasis = 1, nbasis + CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, & + nset=jset, nsgf_set=nsgfj, npgf=npgfj) + ctx%maxlj = MAX(ctx%maxlj, imax) + ctx%max_nset = MAX(ctx%max_nset, jset) + ctx%max_nsgfj = MAX(ctx%max_nsgfj, MAXVAL(nsgfj)) + ctx%max_ncoj = MAX(ctx%max_ncoj, MAXVAL(npgfj)*ncoset(ctx%maxlj)) + END DO + DO ibasis = 1, nbasis + CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, & + nset=kset, nsgf_set=nsgfk, npgf=npgfk) + ctx%maxlk = MAX(ctx%maxlk, imax) + ctx%max_nset = MAX(ctx%max_nset, kset) + ctx%max_nsgfk = MAX(ctx%max_nsgfk, MAXVAL(nsgfk)) + ctx%max_ncok = MAX(ctx%max_ncok, MAXVAL(npgfk)*ncoset(ctx%maxlk)) + END DO + ctx%m_max = ctx%maxli + ctx%maxlj + ctx%maxlk + ctx%max_am = MAX(ctx%maxli, ctx%maxlj, ctx%maxlk) + + ! contiguous (and for j transposed) sphi copies, shared read-only across threads + ALLOCATE (ctx%spi(ctx%max_nset, nbasis), ctx%tspj(ctx%max_nset, nbasis), & + ctx%spk(ctx%max_nset, nbasis)) + DO ibasis = 1, nbasis + DO iset = 1, ctx%max_nset + NULLIFY (ctx%spi(iset, ibasis)%array) + NULLIFY (ctx%tspj(iset, ibasis)%array) + NULLIFY (ctx%spk(iset, ibasis)%array) + END DO + END DO + DO ilist = 1, 3 + DO ibasis = 1, nbasis + IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set + IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set + IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set + DO iset = 1, basis_set%nset + ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset)) + sgfi = basis_set%first_sgf(1, iset) + egfi = sgfi + basis_set%nsgf_set(iset) - 1 + IF (ilist == 1) THEN + ALLOCATE (ctx%spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset))) + ctx%spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi) + ELSE IF (ilist == 2) THEN + ALLOCATE (ctx%tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi)) + ctx%tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi)) + ELSE + ALLOCATE (ctx%spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset))) + ctx%spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi) + END IF + END DO + END DO + END DO + + ! one-time table inits; the truncated-Coulomb init reads a file + bcasts => MPI-collective, + ! must happen here and never inside the per-block path or an OMP region + IF (ctx%op_ij == do_potential_truncated .OR. ctx%op_jk == do_potential_truncated) THEN + IF (ctx%m_max > get_lmax_init()) THEN + IF (para_env%mepos == 0) THEN + CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename) + END IF + CALL init(ctx%m_max, unit_id, para_env%mepos, para_env) + IF (para_env%mepos == 0) THEN + CALL close_file(unit_id) + END IF + END IF + END IF + CALL init_md_ftable(nmax=ctx%m_max) + + CALL timestop(handle) + + END SUBROUTINE gw_3c_context_create + +! ************************************************************************************************** +!> \brief Releases the shared 3c-integral context. +!> \param ctx ... +! ************************************************************************************************** + SUBROUTINE gw_3c_context_release(ctx) + + TYPE(gw_3c_context_type), INTENT(INOUT) :: ctx + + INTEGER :: ibasis, iset + + DO iset = 1, SIZE(ctx%spi, 1) + DO ibasis = 1, SIZE(ctx%spi, 2) + IF (ASSOCIATED(ctx%spi(iset, ibasis)%array)) DEALLOCATE (ctx%spi(iset, ibasis)%array) + IF (ASSOCIATED(ctx%tspj(iset, ibasis)%array)) DEALLOCATE (ctx%tspj(iset, ibasis)%array) + IF (ASSOCIATED(ctx%spk(iset, ibasis)%array)) DEALLOCATE (ctx%spk(iset, ibasis)%array) + END DO + END DO + DEALLOCATE (ctx%spi, ctx%tspj, ctx%spk) + NULLIFY (ctx%spi, ctx%tspj, ctx%spk, ctx%particle_set, ctx%cell) + IF (ALLOCATED(ctx%kind_of)) DEALLOCATE (ctx%kind_of) + IF (ALLOCATED(ctx%basis_i)) DEALLOCATE (ctx%basis_i) + IF (ALLOCATED(ctx%basis_j)) DEALLOCATE (ctx%basis_j) + IF (ALLOCATED(ctx%basis_k)) DEALLOCATE (ctx%basis_k) + + END SUBROUTINE gw_3c_context_release + +! ************************************************************************************************** +!> \brief Creates a per-thread 3c workspace: libint object + LIBXSMM contraction buffers. +!> \param ws ... +!> \param ctx ... +! ************************************************************************************************** + SUBROUTINE gw_3c_workspace_create(ws, ctx) + + TYPE(gw_3c_workspace_type), INTENT(OUT) :: ws + TYPE(gw_3c_context_type), INTENT(IN) :: ctx + + CALL cp_libint_init_3eri(ws%lib, ctx%max_am) + CALL cp_libint_set_contrdepth(ws%lib, 1) + ALLOCATE (ws%cpp_buffer(ctx%max_nsgfj*ctx%max_ncok), & + ws%ccp_buffer(ctx%max_nsgfj*ctx%max_nsgfk*ctx%max_ncoi)) + + END SUBROUTINE gw_3c_workspace_create + +! ************************************************************************************************** +!> \brief Releases a per-thread 3c workspace. +!> \param ws ... +! ************************************************************************************************** + SUBROUTINE gw_3c_workspace_release(ws) + + TYPE(gw_3c_workspace_type), INTENT(INOUT) :: ws + + CALL cp_libint_cleanup_3eri(ws%lib) + DEALLOCATE (ws%cpp_buffer, ws%ccp_buffer) + + END SUBROUTINE gw_3c_workspace_release + +! ************************************************************************************************** +!> \brief Computes the 3c integral block (mu(atom_j) nu(atom_k) | P(atom_i)) for ONE atom triple, +!> accumulating into the caller-zeroed int_3c at the given block offsets. +!> Thread-safe: reads the frozen ctx + read-only module tables, mutates only its arguments +!> and the per-thread ws. +!> \param int_3c pre-zeroed target; the triple's contribution is accumulated in place +!> \param ctx shared context from gw_3c_context_create +!> \param ws per-thread workspace from gw_3c_workspace_create +!> \param atom_j ... +!> \param atom_k ... +!> \param atom_i ... +!> \param cell_j ... +!> \param cell_k ... +!> \param cell_i ... +!> \param j_offset block offset of atom_j's first sgf in int_3c dim 1 (default 0) +!> \param k_offset block offset of atom_k's first sgf in int_3c dim 2 (default 0) +!> \param i_offset block offset of atom_i's first RI sgf in int_3c dim 3 (default 0) +!> \param screened .TRUE. if the kind-radius screens killed the whole triple (int_3c untouched) +! ************************************************************************************************** + SUBROUTINE build_3c_integral_block_ctx(int_3c, ctx, ws, atom_j, atom_k, atom_i, & + cell_j, cell_k, cell_i, & + j_offset, k_offset, i_offset, screened) + + REAL(KIND=dp), DIMENSION(:, :, :) :: int_3c + TYPE(gw_3c_context_type), INTENT(IN) :: ctx + TYPE(gw_3c_workspace_type), INTENT(INOUT) :: ws + INTEGER, INTENT(IN) :: atom_j, atom_k, atom_i + INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: cell_j, cell_k, cell_i + INTEGER, INTENT(IN), OPTIONAL :: j_offset, k_offset, i_offset + LOGICAL, INTENT(OUT), OPTIONAL :: screened + + INTEGER :: block_end_i, block_end_j, block_end_k, block_start_i, block_start_j, & + block_start_k, ikind, iset, jkind, jset, kkind, kset, my_i_offset, my_j_offset, & + my_k_offset, ncoi, ncoj, ncok, nseti, nsetj, nsetk, sgfi, sgfj, sgfk + INTEGER, DIMENSION(3) :: my_cell_i, my_cell_j, my_cell_k + INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, & + lmin_k, npgfi, npgfj, npgfk, nsgfi, & + nsgfj, nsgfk + INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k + REAL(KIND=dp) :: dij, dik, djk, kind_radius_i, & + kind_radius_j, kind_radius_k, sijk_ext + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr + REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk + REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k + REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, zeti, zetj, zetk + + IF (PRESENT(screened)) screened = .FALSE. + + my_cell_i(1:3) = 0 + IF (PRESENT(cell_i)) my_cell_i(1:3) = cell_i(1:3) + my_cell_j(1:3) = 0 + IF (PRESENT(cell_j)) my_cell_j(1:3) = cell_j(1:3) + my_cell_k(1:3) = 0 + IF (PRESENT(cell_k)) my_cell_k(1:3) = cell_k(1:3) + my_i_offset = 0 + IF (PRESENT(i_offset)) my_i_offset = i_offset + my_j_offset = 0 + IF (PRESENT(j_offset)) my_j_offset = j_offset + my_k_offset = 0 + IF (PRESENT(k_offset)) my_k_offset = k_offset + + ri = pbc(ctx%particle_set(atom_i)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_i, dp)) + rj = pbc(ctx%particle_set(atom_j)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_j, dp)) + rk = pbc(ctx%particle_set(atom_k)%r(1:3), ctx%cell) + MATMUL(ctx%hmat, REAL(my_cell_k, dp)) + + rjk(1:3) = rk(1:3) - rj(1:3) + rij(1:3) = rj(1:3) - ri(1:3) + rik(1:3) = rk(1:3) - ri(1:3) + + djk = NORM2(rjk) + dij = NORM2(rij) + dik = NORM2(rik) + + ikind = ctx%kind_of(atom_i) + jkind = ctx%kind_of(atom_j) + kkind = ctx%kind_of(atom_k) + + CALL get_gto_basis_set(ctx%basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, & + lmax=lmax_i, lmin=lmin_i, npgf=npgfi, nset=nseti, & + nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, & + zet=zeti, kind_radius=kind_radius_i) + CALL get_gto_basis_set(ctx%basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, & + lmax=lmax_j, lmin=lmin_j, npgf=npgfj, nset=nsetj, & + nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, & + zet=zetj, kind_radius=kind_radius_j) + CALL get_gto_basis_set(ctx%basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, & + lmax=lmax_k, lmin=lmin_k, npgf=npgfk, nset=nsetk, & + nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, & + zet=zetk, kind_radius=kind_radius_k) + + IF (kind_radius_j + kind_radius_i + ctx%dr_ij < dij .OR. & + kind_radius_j + kind_radius_k + ctx%dr_jk < djk .OR. & + kind_radius_k + kind_radius_i + ctx%dr_ik < dik) THEN + IF (PRESENT(screened)) screened = .TRUE. + RETURN + END IF + + DO iset = 1, nseti + DO jset = 1, nsetj + IF (set_radius_j(jset) + set_radius_i(iset) + ctx%dr_ij < dij) CYCLE + DO kset = 1, nsetk + IF (set_radius_j(jset) + set_radius_k(kset) + ctx%dr_jk < djk) CYCLE + IF (set_radius_k(kset) + set_radius_i(iset) + ctx%dr_ik < dik) CYCLE + + ncoi = npgfi(iset)*ncoset(lmax_i(iset)) + ncoj = npgfj(jset)*ncoset(lmax_j(jset)) + ncok = npgfk(kset)*ncoset(lmax_k(kset)) + + sgfi = first_sgf_i(1, iset) + sgfj = first_sgf_j(1, jset) + sgfk = first_sgf_k(1, kset) + + IF (ncoj*ncok*ncoi <= 0) CYCLE + ALLOCATE (sijk(ncoj, ncok, ncoi)) + sijk(:, :, :) = 0.0_dp + + CALL eri_3center(sijk, & + lmin_j(jset), lmax_j(jset), npgfj(jset), zetj(:, jset), & + rpgf_j(:, jset), rj, & + lmin_k(kset), lmax_k(kset), npgfk(kset), zetk(:, kset), & + rpgf_k(:, kset), rk, & + lmin_i(iset), lmax_i(iset), npgfi(iset), zeti(:, iset), & + rpgf_i(:, iset), ri, & + djk, dij, dik, ws%lib, ctx%potential_parameter, & + int_abc_ext=sijk_ext) + + ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset))) + CALL abc_contract_xsmm(sijk_contr, sijk, ctx%tspj(jset, jkind)%array, & + ctx%spk(kset, kkind)%array, ctx%spi(iset, ikind)%array, & + ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), & + nsgfi(iset), ws%cpp_buffer, ws%ccp_buffer) + DEALLOCATE (sijk) + + block_start_j = sgfj + my_j_offset + block_end_j = sgfj + nsgfj(jset) - 1 + my_j_offset + block_start_k = sgfk + my_k_offset + block_end_k = sgfk + nsgfk(kset) - 1 + my_k_offset + block_start_i = sgfi + my_i_offset + block_end_i = sgfi + nsgfi(iset) - 1 + my_i_offset + + int_3c(block_start_j:block_end_j, & + block_start_k:block_end_k, & + block_start_i:block_end_i) = & + int_3c(block_start_j:block_end_j, & + block_start_k:block_end_k, & + block_start_i:block_end_i) + & + sijk_contr(:, :, :) + DEALLOCATE (sijk_contr) + + END DO + END DO + END DO + + END SUBROUTINE build_3c_integral_block_ctx + +! ************************************************************************************************** +!> \brief Standalone wrapper: creates a transient ctx + ws, builds the 3c integral block(s) for the +!> requested atom triple(s) via build_3c_integral_block_ctx, then releases them. Reuse a +!> persistent ctx with build_3c_integral_block_ctx directly in OMP / repeated-call paths. !> \param int_3c ... !> \param qs_env ... !> \param potential_parameter ... @@ -94,176 +473,24 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'build_3c_integral_block' - INTEGER :: at_i, at_j, at_k, block_end_i, block_end_j, block_end_k, block_start_i, & - block_start_j, block_start_k, egfi, handle, i, i_offset, ibasis, ikind, ilist, imax, is, & - iset, j_offset, jkind, js, jset, k_offset, kkind, ks, kset, m_max, max_ncoi, max_ncoj, & - max_ncok, max_nset, max_nsgfi, max_nsgfj, max_nsgfk, maxli, maxlj, maxlk, natom, nbasis, & - ncoi, ncoj, ncok, nseti, nsetj, nsetk, op_ij, op_jk, sgfi, sgfj, sgfk, unit_id - INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of - INTEGER, DIMENSION(3) :: my_cell_i, my_cell_j, my_cell_k - INTEGER, DIMENSION(:), POINTER :: lmax_i, lmax_j, lmax_k, lmin_i, lmin_j, & - lmin_k, npgfi, npgfj, npgfk, nsgfi, & - nsgfj, nsgfk - INTEGER, DIMENSION(:, :), POINTER :: first_sgf_i, first_sgf_j, first_sgf_k - REAL(KIND=dp) :: dij, dik, djk, dr_ij, dr_ik, dr_jk, & - kind_radius_i, kind_radius_j, & - kind_radius_k, sijk_ext - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ccp_buffer, cpp_buffer, & - max_contraction_i, max_contraction_j, & - max_contraction_k - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: sijk, sijk_contr - REAL(KIND=dp), DIMENSION(3) :: ri, rij, rik, rj, rjk, rk - REAL(KIND=dp), DIMENSION(3, 3) :: hmat - REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_i, set_radius_j, set_radius_k - REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgf_i, rpgf_j, rpgf_k, sphi_i, sphi_j, & - sphi_k, zeti, zetj, zetk - TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(cell_type), POINTER :: cell - TYPE(cp_2d_r_p_type), DIMENSION(:, :), POINTER :: spi, spk, tspj - TYPE(cp_libint_t) :: lib - TYPE(dft_control_type), POINTER :: dft_control - TYPE(gto_basis_set_type), POINTER :: basis_set - TYPE(mp_para_env_type), POINTER :: para_env - TYPE(particle_type), DIMENSION(:), POINTER :: particle_set - TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + INTEGER :: at_i, at_j, at_k, handle, my_i_offset, & + my_j_offset, my_k_offset + TYPE(gw_3c_context_type) :: ctx + TYPE(gw_3c_workspace_type) :: ws CALL timeset(routineN, handle) - op_ij = potential_parameter%potential_type - op_jk = do_potential_id + CALL gw_3c_context_create(ctx, qs_env, potential_parameter, basis_j, basis_k, basis_i) + CALL gw_3c_workspace_create(ws, ctx) - dr_ij = 0.0_dp; dr_jk = 0.0_dp; dr_ik = 0.0_dp - - IF (op_ij == do_potential_truncated .OR. op_ij == do_potential_short) THEN - dr_ij = potential_parameter%cutoff_radius*cutoff_screen_factor - dr_ik = potential_parameter%cutoff_radius*cutoff_screen_factor - ELSEIF (op_ij == do_potential_coulomb) THEN - dr_ij = 1000000.0_dp - dr_ik = 1000000.0_dp - END IF - - NULLIFY (qs_kind_set, atomic_kind_set) - - ! get stuff - CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, & - natom=natom, dft_control=dft_control, para_env=para_env, & - particle_set=particle_set) - CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) - CALL get_cell(cell=cell, h=hmat) - - !Need the max l for each basis for libint and max nset, nco and nsgf for LIBXSMM contraction - nbasis = SIZE(basis_i) - max_nsgfi = 0 - max_ncoi = 0 - max_nset = 0 - maxli = 0 - DO ibasis = 1, nbasis - CALL get_gto_basis_set(gto_basis_set=basis_i(ibasis)%gto_basis_set, maxl=imax, & - nset=iset, nsgf_set=nsgfi, npgf=npgfi) - maxli = MAX(maxli, imax) - max_nset = MAX(max_nset, iset) - max_nsgfi = MAX(max_nsgfi, MAXVAL(nsgfi)) - max_ncoi = MAX(max_ncoi, MAXVAL(npgfi)*ncoset(maxli)) - END DO - max_nsgfj = 0 - max_ncoj = 0 - maxlj = 0 - DO ibasis = 1, nbasis - CALL get_gto_basis_set(gto_basis_set=basis_j(ibasis)%gto_basis_set, maxl=imax, & - nset=jset, nsgf_set=nsgfj, npgf=npgfj) - maxlj = MAX(maxlj, imax) - max_nset = MAX(max_nset, jset) - max_nsgfj = MAX(max_nsgfj, MAXVAL(nsgfj)) - max_ncoj = MAX(max_ncoj, MAXVAL(npgfj)*ncoset(maxlj)) - END DO - max_nsgfk = 0 - max_ncok = 0 - maxlk = 0 - DO ibasis = 1, nbasis - CALL get_gto_basis_set(gto_basis_set=basis_k(ibasis)%gto_basis_set, maxl=imax, & - nset=kset, nsgf_set=nsgfk, npgf=npgfk) - maxlk = MAX(maxlk, imax) - max_nset = MAX(max_nset, kset) - max_nsgfk = MAX(max_nsgfk, MAXVAL(nsgfk)) - max_ncok = MAX(max_ncok, MAXVAL(npgfk)*ncoset(maxlk)) - END DO - m_max = maxli + maxlj + maxlk - - !To minimize expensive memory opsand generally optimize contraction, pre-allocate - !contiguous sphi arrays (and transposed in the cas of sphi_i) - - NULLIFY (tspj, spi, spk) - ALLOCATE (spi(max_nset, nbasis), tspj(max_nset, nbasis), spk(max_nset, nbasis)) - - DO ibasis = 1, nbasis - DO iset = 1, max_nset - NULLIFY (spi(iset, ibasis)%array) - NULLIFY (tspj(iset, ibasis)%array) - - NULLIFY (spk(iset, ibasis)%array) - END DO - END DO - - DO ilist = 1, 3 - DO ibasis = 1, nbasis - IF (ilist == 1) basis_set => basis_i(ibasis)%gto_basis_set - IF (ilist == 2) basis_set => basis_j(ibasis)%gto_basis_set - IF (ilist == 3) basis_set => basis_k(ibasis)%gto_basis_set - - DO iset = 1, basis_set%nset - - ncoi = basis_set%npgf(iset)*ncoset(basis_set%lmax(iset)) - sgfi = basis_set%first_sgf(1, iset) - egfi = sgfi + basis_set%nsgf_set(iset) - 1 - - IF (ilist == 1) THEN - ALLOCATE (spi(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset))) - spi(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi) - - ELSE IF (ilist == 2) THEN - ALLOCATE (tspj(iset, ibasis)%array(basis_set%nsgf_set(iset), ncoi)) - tspj(iset, ibasis)%array(:, :) = TRANSPOSE(basis_set%sphi(1:ncoi, sgfi:egfi)) - - ELSE - ALLOCATE (spk(iset, ibasis)%array(ncoi, basis_set%nsgf_set(iset))) - spk(iset, ibasis)%array(:, :) = basis_set%sphi(1:ncoi, sgfi:egfi) - END IF - - END DO !iset - END DO !ibasis - END DO !ilist - - !Init the truncated Coulomb operator - IF (op_ij == do_potential_truncated .OR. op_jk == do_potential_truncated) THEN - - IF (m_max > get_lmax_init()) THEN - IF (para_env%mepos == 0) THEN - CALL open_file(unit_number=unit_id, file_name=potential_parameter%filename) - END IF - CALL init(m_max, unit_id, para_env%mepos, para_env) - IF (para_env%mepos == 0) THEN - CALL close_file(unit_id) - END IF - END IF - END IF - - CALL init_md_ftable(nmax=m_max) - - CALL cp_libint_init_3eri(lib, MAX(maxli, maxlj, maxlk)) - CALL cp_libint_set_contrdepth(lib, 1) - - !pre-allocate contraction buffers - ALLOCATE (cpp_buffer(max_nsgfj*max_ncok), ccp_buffer(max_nsgfj*max_nsgfk*max_ncoi)) int_3c(:, :, :) = 0.0_dp ! loop over all RI atoms - DO at_i = 1, natom - + DO at_i = 1, ctx%natom ! loop over all AO atoms - DO at_j = 1, natom - + DO at_j = 1, ctx%natom ! loop over all AO atoms - DO at_k = 1, natom + DO at_k = 1, ctx%natom IF (PRESENT(atom_i)) THEN IF (at_i /= atom_i) CYCLE @@ -275,174 +502,36 @@ CONTAINS IF (at_k /= atom_k) CYCLE END IF - my_cell_i(1:3) = 0 - IF (PRESENT(cell_i)) my_cell_i(1:3) = cell_i(1:3) - my_cell_j(1:3) = 0 - IF (PRESENT(cell_j)) my_cell_j(1:3) = cell_j(1:3) - my_cell_k(1:3) = 0 - IF (PRESENT(cell_k)) my_cell_k(1:3) = cell_k(1:3) + IF (PRESENT(atom_j)) THEN + my_j_offset = 0 + ELSE + CPASSERT(PRESENT(j_bf_start_from_atom)) + my_j_offset = j_bf_start_from_atom(at_j) - 1 + END IF + IF (PRESENT(atom_k)) THEN + my_k_offset = 0 + ELSE + CPASSERT(PRESENT(k_bf_start_from_atom)) + my_k_offset = k_bf_start_from_atom(at_k) - 1 + END IF + IF (PRESENT(atom_i)) THEN + my_i_offset = 0 + ELSE + CPASSERT(PRESENT(i_bf_start_from_atom)) + my_i_offset = i_bf_start_from_atom(at_i) - 1 + END IF - ri = pbc(particle_set(at_i)%r(1:3), cell) + MATMUL(hmat, REAL(my_cell_i, dp)) - rj = pbc(particle_set(at_j)%r(1:3), cell) + MATMUL(hmat, REAL(my_cell_j, dp)) - rk = pbc(particle_set(at_k)%r(1:3), cell) + MATMUL(hmat, REAL(my_cell_k, dp)) - - rjk(1:3) = rk(1:3) - rj(1:3) - rij(1:3) = rj(1:3) - ri(1:3) - rik(1:3) = rk(1:3) - ri(1:3) - - djk = NORM2(rjk) - dij = NORM2(rij) - dik = NORM2(rik) - - ikind = kind_of(at_i) - jkind = kind_of(at_j) - kkind = kind_of(at_k) - - CALL get_gto_basis_set(basis_i(ikind)%gto_basis_set, first_sgf=first_sgf_i, & - lmax=lmax_i, lmin=lmin_i, npgf=npgfi, nset=nseti, & - nsgf_set=nsgfi, pgf_radius=rpgf_i, set_radius=set_radius_i, & - sphi=sphi_i, zet=zeti, kind_radius=kind_radius_i) - - CALL get_gto_basis_set(basis_j(jkind)%gto_basis_set, first_sgf=first_sgf_j, & - lmax=lmax_j, lmin=lmin_j, npgf=npgfj, nset=nsetj, & - nsgf_set=nsgfj, pgf_radius=rpgf_j, set_radius=set_radius_j, & - sphi=sphi_j, zet=zetj, kind_radius=kind_radius_j) - - CALL get_gto_basis_set(basis_k(kkind)%gto_basis_set, first_sgf=first_sgf_k, & - lmax=lmax_k, lmin=lmin_k, npgf=npgfk, nset=nsetk, & - nsgf_set=nsgfk, pgf_radius=rpgf_k, set_radius=set_radius_k, & - sphi=sphi_k, zet=zetk, kind_radius=kind_radius_k) - - IF (kind_radius_j + kind_radius_i + dr_ij < dij) CYCLE - IF (kind_radius_j + kind_radius_k + dr_jk < djk) CYCLE - IF (kind_radius_k + kind_radius_i + dr_ik < dik) CYCLE - - ALLOCATE (max_contraction_i(nseti)) - max_contraction_i = 0.0_dp - DO iset = 1, nseti - sgfi = first_sgf_i(1, iset) - max_contraction_i(iset) = MAXVAL([(SUM(ABS(sphi_i(:, i))), i=sgfi, & - sgfi + nsgfi(iset) - 1)]) - END DO - - ALLOCATE (max_contraction_j(nsetj)) - max_contraction_j = 0.0_dp - DO jset = 1, nsetj - sgfj = first_sgf_j(1, jset) - max_contraction_j(jset) = MAXVAL([(SUM(ABS(sphi_j(:, i))), i=sgfj, & - sgfj + nsgfj(jset) - 1)]) - END DO - - ALLOCATE (max_contraction_k(nsetk)) - max_contraction_k = 0.0_dp - DO kset = 1, nsetk - sgfk = first_sgf_k(1, kset) - max_contraction_k(kset) = MAXVAL([(SUM(ABS(sphi_k(:, i))), i=sgfk, & - sgfk + nsgfk(kset) - 1)]) - END DO - - DO iset = 1, nseti - - DO jset = 1, nsetj - - IF (set_radius_j(jset) + set_radius_i(iset) + dr_ij < dij) CYCLE - - DO kset = 1, nsetk - - IF (set_radius_j(jset) + set_radius_k(kset) + dr_jk < djk) CYCLE - IF (set_radius_k(kset) + set_radius_i(iset) + dr_ik < dik) CYCLE - - ncoi = npgfi(iset)*ncoset(lmax_i(iset)) - ncoj = npgfj(jset)*ncoset(lmax_j(jset)) - ncok = npgfk(kset)*ncoset(lmax_k(kset)) - - sgfi = first_sgf_i(1, iset) - sgfj = first_sgf_j(1, jset) - sgfk = first_sgf_k(1, kset) - - IF (ncoj*ncok*ncoi <= 0) CYCLE - ALLOCATE (sijk(ncoj, ncok, ncoi)) - sijk(:, :, :) = 0.0_dp - - is = iset - js = jset - ks = kset - - CALL eri_3center(sijk, & - lmin_j(js), lmax_j(js), npgfj(js), zetj(:, js), & - rpgf_j(:, js), rj, & - lmin_k(ks), lmax_k(ks), npgfk(ks), zetk(:, ks), & - rpgf_k(:, ks), rk, & - lmin_i(is), lmax_i(is), npgfi(is), zeti(:, is), & - rpgf_i(:, is), ri, & - djk, dij, dik, lib, potential_parameter, & - int_abc_ext=sijk_ext) - - ALLOCATE (sijk_contr(nsgfj(jset), nsgfk(kset), nsgfi(iset))) - CALL abc_contract_xsmm(sijk_contr, sijk, tspj(jset, jkind)%array, & - spk(kset, kkind)%array, spi(iset, ikind)%array, & - ncoj, ncok, ncoi, nsgfj(jset), nsgfk(kset), & - nsgfi(iset), cpp_buffer, ccp_buffer) - DEALLOCATE (sijk) - - IF (PRESENT(atom_j)) THEN - j_offset = 0 - ELSE - CPASSERT(PRESENT(j_bf_start_from_atom)) - j_offset = j_bf_start_from_atom(at_j) - 1 - END IF - IF (PRESENT(atom_k)) THEN - k_offset = 0 - ELSE - CPASSERT(PRESENT(k_bf_start_from_atom)) - k_offset = k_bf_start_from_atom(at_k) - 1 - END IF - IF (PRESENT(atom_i)) THEN - i_offset = 0 - ELSE - CPASSERT(PRESENT(i_bf_start_from_atom)) - i_offset = i_bf_start_from_atom(at_i) - 1 - END IF - - block_start_j = sgfj + j_offset - block_end_j = sgfj + nsgfj(jset) - 1 + j_offset - block_start_k = sgfk + k_offset - block_end_k = sgfk + nsgfk(kset) - 1 + k_offset - block_start_i = sgfi + i_offset - block_end_i = sgfi + nsgfi(iset) - 1 + i_offset - - int_3c(block_start_j:block_end_j, & - block_start_k:block_end_k, & - block_start_i:block_end_i) = & - int_3c(block_start_j:block_end_j, & - block_start_k:block_end_k, & - block_start_i:block_end_i) + & - sijk_contr(:, :, :) - DEALLOCATE (sijk_contr) - - END DO - - END DO - - END DO - - DEALLOCATE (max_contraction_i, max_contraction_j, max_contraction_k) + CALL build_3c_integral_block_ctx(int_3c, ctx, ws, at_j, at_k, at_i, & + cell_j=cell_j, cell_k=cell_k, cell_i=cell_i, & + j_offset=my_j_offset, k_offset=my_k_offset, & + i_offset=my_i_offset) END DO ! atom_k (AO) END DO ! atom_j (AO) END DO ! atom_i (RI) - CALL cp_libint_cleanup_3eri(lib) - - DO iset = 1, max_nset - DO ibasis = 1, nbasis - IF (ASSOCIATED(spi(iset, ibasis)%array)) DEALLOCATE (spi(iset, ibasis)%array) - IF (ASSOCIATED(tspj(iset, ibasis)%array)) DEALLOCATE (tspj(iset, ibasis)%array) - - IF (ASSOCIATED(spk(iset, ibasis)%array)) DEALLOCATE (spk(iset, ibasis)%array) - END DO - END DO - DEALLOCATE (spi, tspj, spk) + CALL gw_3c_workspace_release(ws) + CALL gw_3c_context_release(ctx) CALL timestop(handle) diff --git a/src/gw_large_cell_gamma_ri_rs.F b/src/gw_large_cell_gamma_ri_rs.F index 7ef7f0267c..ea5bcaaa55 100644 --- a/src/gw_large_cell_gamma_ri_rs.F +++ b/src/gw_large_cell_gamma_ri_rs.F @@ -34,7 +34,13 @@ MODULE gw_large_cell_Gamma_ri_rs cp_logger_type USE cp_output_handling, ONLY: cp_p_file,& cp_print_key_should_output - USE gw_integrals, ONLY: build_3c_integral_block + USE gw_integrals, ONLY: build_3c_integral_block_ctx,& + gw_3c_context_create,& + gw_3c_context_release,& + gw_3c_context_type,& + gw_3c_workspace_create,& + gw_3c_workspace_release,& + gw_3c_workspace_type USE gw_large_cell_gamma, ONLY: G_occ_vir,& compute_QP_energies,& delete_unnecessary_files,& @@ -627,15 +633,14 @@ CONTAINS ri_blk_sizes, row_dist_grid REAL(KIND=dp) :: cutoff_ri, dist, t1 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: d_vec_local - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_local, int_2d_local, & - phi_global, phi_local, rho_pair_local, & - Z_blk, Z_global_temp - REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_local + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: D_local, d_lp_local, phi_global, & + phi_local, Z_blk, Z_global_temp REAL(KIND=dp), DIMENSION(3) :: pos_P REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_ptr TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_distribution_type) :: dist_phi, dist_Z TYPE(dbcsr_iterator_type) :: iter + TYPE(gw_3c_context_type) :: ctx_3c TYPE(mp_para_env_type), POINTER :: para_env TYPE(particle_type), DIMENSION(:), POINTER :: particle_set TYPE(section_vals_type), POINTER :: input @@ -697,6 +702,11 @@ CONTAINS END DO max_loc_ri = MAXVAL(ri_blk_sizes) + ! shared 3c-integral context for all atoms of this rank (collective: t_c_g0 init) + CALL gw_3c_context_create(ctx_3c, qs_env, bs_env%ri_metric, & + basis_j=bs_env%basis_set_AO, basis_k=bs_env%basis_set_AO, & + basis_i=bs_env%basis_set_RI) + ! Global mapping array: 4 bytes * N_grid_total (~4MB per million grid points) ALLOCATE (map_global_to_local(n_grid_total)) @@ -814,14 +824,11 @@ CONTAINS ! D. Build Local RHS Matrix (d_lp_local) ! --------------------------------------------------------------------- ALLOCATE (d_lp_local(n_local_grid, n_loc_ri)) - ALLOCATE (rho_pair_local(n_local_grid, max_ao_size*max_ao_size)) - ALLOCATE (int_3c_local(max_ao_size, max_ao_size, n_loc_ri)) - ALLOCATE (int_2d_local(max_ao_size*max_ao_size, n_loc_ri)) d_lp_local = 0.0_dp - CALL compute_d_lp(qs_env, bs_env, phi_local, d_lp_local, n_local_grid, & - n_loc_ri, atom_P, rho_pair_local, int_3c_local, int_2d_local, max_ao_size) + CALL compute_d_lp(bs_env, ctx_3c, phi_local, d_lp_local, n_local_grid, & + n_loc_ri, atom_P, max_ao_size) !$OMP PARALLEL DO DEFAULT(NONE) & !$OMP SHARED(n_loc_ri, n_local_grid, d_lp_local, d_vec_local) & @@ -887,7 +894,6 @@ CONTAINS DEALLOCATE (Z_global_temp, Z_blk) DEALLOCATE (D_local, d_vec_local, d_lp_local) - DEALLOCATE (rho_pair_local, int_3c_local, int_2d_local) DEALLOCATE (local_grid_idx, phi_local) @@ -896,6 +902,8 @@ CONTAINS ! Clean up the massive global array once the loop finishes DEALLOCATE (phi_global) + CALL gw_3c_context_release(ctx_3c) + CALL dbcsr_finalize(mat_Z_lP) IF (bs_env%unit_nr > 0) THEN @@ -924,95 +932,107 @@ CONTAINS END SUBROUTINE compute_coeff_Z_lP ! ************************************************************************************************** -!> \brief Computes the dense localized Right-Hand Side (RHS) matrix d_lp -!> \param qs_env ... +!> \brief Computes the dense localized RHS d_lp(l,P) = Σ_{μν} Φ_μ(r_l)·Φ_ν(r_l)·(μν|P) for one +!> RI atom P, OMP-threaded over (atom_j, atom_k) AO-pair blocks: per thread, build the 3c +!> block, then contract grid-chunked pair densities into d_lp via OMP REDUCTION(+:d_lp) +!> (cross-thread accumulation order differs from serial at eps level). !> \param bs_env ... -!> \param phi_val ... -!> \param d_lp ... -!> \param n_grid_total ... -!> \param n_loc_ri ... +!> \param ctx shared 3c-integral context (gw_3c_context_create) +!> \param phi_val Φ_μ(r_l) on the local-sphere grid (n_grid_total × n_ao) +!> \param d_lp output (n_grid_total × n_loc_ri), zeroed by the caller, accumulated here +!> \param n_grid_total number of local-sphere grid rows +!> \param n_loc_ri number of RI functions of atom_P !> \param atom_P ... -!> \param rho_pair ... -!> \param int_3c ... -!> \param int_2d ... !> \param max_ao_size ... ! ************************************************************************************************** - SUBROUTINE compute_d_lp(qs_env, bs_env, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, & - rho_pair, int_3c, int_2d, max_ao_size) + SUBROUTINE compute_d_lp(bs_env, ctx, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, & + max_ao_size) - TYPE(qs_environment_type), POINTER :: qs_env TYPE(post_scf_bandstructure_type), POINTER :: bs_env + TYPE(gw_3c_context_type), INTENT(IN) :: ctx REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: phi_val REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: d_lp - INTEGER, INTENT(IN) :: n_grid_total, n_loc_ri, atom_P - REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: rho_pair - REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: int_3c - REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: int_2d - INTEGER, INTENT(IN) :: max_ao_size + INTEGER, INTENT(IN) :: n_grid_total, n_loc_ri, atom_P, & + max_ao_size CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_d_lp' + INTEGER, PARAMETER :: grid_chunk = 1024 - INTEGER :: atom_j, atom_k, handle, j, jk_idx, & - jsize, jstart, k, ksize, kstart, l, ri + INTEGER :: atom_j, atom_k, c, handle, j, jk_idx, & + jsize, jstart, k, ksize, kstart, l, & + l0, ri + LOGICAL :: screened + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: int_2d_prv, rho_chunk + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: int_3c_prv + TYPE(gw_3c_workspace_type) :: ws CALL timeset(routineN, handle) - DO atom_j = 1, SIZE(bs_env%i_ao_start_from_atom) - jstart = bs_env%i_ao_start_from_atom(atom_j) - jsize = bs_env%i_ao_end_from_atom(atom_j) - jstart + 1 +!$OMP PARALLEL DEFAULT(NONE) & +!$OMP SHARED(bs_env, ctx, phi_val, d_lp, n_grid_total, n_loc_ri, atom_P, max_ao_size) & +!$OMP PRIVATE(atom_j, atom_k, c, j, jk_idx, jsize, jstart, k, ksize, kstart, l, l0, ri, & +!$OMP screened, int_2d_prv, rho_chunk, int_3c_prv, ws) + CALL gw_3c_workspace_create(ws, ctx) + ALLOCATE (int_3c_prv(max_ao_size, max_ao_size, n_loc_ri)) + ALLOCATE (int_2d_prv(max_ao_size*max_ao_size, n_loc_ri)) + ALLOCATE (rho_chunk(grid_chunk, max_ao_size*max_ao_size)) + +!$OMP DO COLLAPSE(2) SCHEDULE(DYNAMIC) REDUCTION(+:d_lp) + DO atom_j = 1, SIZE(bs_env%i_ao_start_from_atom) DO atom_k = 1, SIZE(bs_env%i_ao_start_from_atom) + jstart = bs_env%i_ao_start_from_atom(atom_j) + jsize = bs_env%i_ao_end_from_atom(atom_j) - jstart + 1 kstart = bs_env%i_ao_start_from_atom(atom_k) ksize = bs_env%i_ao_end_from_atom(atom_k) - kstart + 1 - int_3c(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp + int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri) = 0.0_dp ! Compute B_{μν,P} = (μν|P) - CALL build_3c_integral_block(int_3c(1:jsize, 1:ksize, 1:n_loc_ri), & - qs_env, bs_env%ri_metric, & - basis_j=bs_env%basis_set_AO, & - basis_k=bs_env%basis_set_AO, & - basis_i=bs_env%basis_set_RI, & - atom_k=atom_k, atom_j=atom_j, atom_i=atom_P) + CALL build_3c_integral_block_ctx(int_3c_prv(1:jsize, 1:ksize, 1:n_loc_ri), & + ctx, ws, atom_j=atom_j, atom_k=atom_k, atom_i=atom_P, & + screened=screened) - ! Build Pair Density: ρ(l, μν) = Φ_μ(r_l) \times Φ_ν(r_l) - !$OMP PARALLEL DO DEFAULT(NONE) & - !$OMP SHARED(ksize, jsize, n_grid_total, rho_pair, phi_val, jstart, kstart) & - !$OMP PRIVATE(k, j, l, jk_idx) & - !$OMP COLLAPSE(2) - DO k = 1, ksize - DO j = 1, jsize - jk_idx = (k - 1)*jsize + j - DO l = 1, n_grid_total - rho_pair(l, jk_idx) = phi_val(l, jstart + j - 1)*phi_val(l, kstart + k - 1) - END DO - END DO - END DO - !$OMP END PARALLEL DO + ! screened pair => 3c block is identically zero, nothing to contract + IF (screened) CYCLE ! Flatten 3D B_{μν, P} tensor to 2D B_{(μν), P} matrix for BLAS - !$OMP PARALLEL DO DEFAULT(NONE) & - !$OMP SHARED(n_loc_ri, ksize, jsize, int_2d, int_3c) & - !$OMP PRIVATE(ri, k, j, jk_idx) & - !$OMP COLLAPSE(2) DO ri = 1, n_loc_ri DO k = 1, ksize DO j = 1, jsize jk_idx = (k - 1)*jsize + j - int_2d(jk_idx, ri) = int_3c(j, k, ri) + int_2d_prv(jk_idx, ri) = int_3c_prv(j, k, ri) END DO END DO END DO - !$OMP END PARALLEL DO - ! Perform Contraction: d_{l, P} += ρ(l, μν) \times B_{(μν), P} - CALL dgemm("N", "N", n_grid_total, n_loc_ri, jsize*ksize, & - 1.0_dp, rho_pair, n_grid_total, & - int_2d, max_ao_size*max_ao_size, & - 1.0_dp, d_lp(1, 1), n_grid_total) + ! Pair density ρ(l, μν) = Φ_μ(r_l) Φ_ν(r_l) in grid chunks, contracted on the fly: + ! d_{l,P} += ρ(l, μν) B_{(μν),P} (dgemm runs serially inside the parallel region) + DO l0 = 1, n_grid_total, grid_chunk + c = MIN(grid_chunk, n_grid_total - l0 + 1) + DO k = 1, ksize + DO j = 1, jsize + jk_idx = (k - 1)*jsize + j + DO l = 1, c + rho_chunk(l, jk_idx) = phi_val(l0 + l - 1, jstart + j - 1)* & + phi_val(l0 + l - 1, kstart + k - 1) + END DO + END DO + END DO + CALL dgemm("N", "N", c, n_loc_ri, jsize*ksize, & + 1.0_dp, rho_chunk, grid_chunk, & + int_2d_prv, max_ao_size*max_ao_size, & + 1.0_dp, d_lp(l0, 1), n_grid_total) + END DO END DO END DO +!$OMP END DO + + DEALLOCATE (int_3c_prv, int_2d_prv, rho_chunk) + CALL gw_3c_workspace_release(ws) + +!$OMP END PARALLEL CALL timestop(handle)