mirror of
https://github.com/cp2k/cp2k.git
synced 2026-07-28 14:15:19 -04:00
RIRS-GW: optimize OMP parallelization (#5399)
Co-authored-by: Claude <noreply@anthropic.com>
This commit is contained in:
Max Graml
GitHub
parent
fc8b99b17a
commit
2085e8ca3a
2 changed files with 501 additions and 392 deletions
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue