Cabana_VerletList.hpp

Go to the documentation of this file.

/****************************************************************************
 * Copyright (c) 2018-2023 by the Cabana authors                            *
 * All rights reserved.                                                     *
 *                                                                          *
 * This file is part of the Cabana library. Cabana is distributed under a   *
 * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
 * the top-level directory.                                                 *
 *                                                                          *
 * SPDX-License-Identifier: BSD-3-Clause                                    *
 ****************************************************************************/

#ifndef CABANA_VERLETLIST_HPP
#define CABANA_VERLETLIST_HPP
 
#include <Cabana_LinkedCellList.hpp>
#include <Cabana_NeighborList.hpp>
#include <Cabana_Parallel.hpp>
 
#include <Kokkos_Core.hpp>
#include <Kokkos_Profiling_ScopedRegion.hpp>
 
#include <cassert>
 
namespace Cabana
{
//---------------------------------------------------------------------------//
// Verlet List Memory Layout Tag.
//---------------------------------------------------------------------------//
struct VerletLayoutCSR
{
};

struct VerletLayout2D
{
};
 
//---------------------------------------------------------------------------//
// Verlet List Data.
//---------------------------------------------------------------------------//
template <class MemorySpace, class LayoutTag>
struct VerletListData;

template <class MemorySpace>
struct VerletListData<MemorySpace, VerletLayoutCSR>
{
    using memory_space = MemorySpace;

    Kokkos::View<int*, memory_space> counts;

    Kokkos::View<int*, memory_space> offsets;

    Kokkos::View<int*, memory_space> neighbors;

    KOKKOS_INLINE_FUNCTION
    void addNeighbor( const int pid, const int nid ) const
    {
        neighbors( offsets( pid ) +
                   Kokkos::atomic_fetch_add( &counts( pid ), 1 ) ) = nid;
    }

    KOKKOS_INLINE_FUNCTION
    void setNeighbor( const int pid, const int nid, const int new_id ) const
    {
        neighbors( offsets( pid ) + nid ) = new_id;
    }
};

template <class MemorySpace>
struct VerletListData<MemorySpace, VerletLayout2D>
{
    using memory_space = MemorySpace;

    Kokkos::View<int*, memory_space> counts;

    Kokkos::View<int**, memory_space> neighbors;

    std::size_t max_n;

    KOKKOS_INLINE_FUNCTION
    void addNeighbor( const int pid, const int nid ) const
    {
        std::size_t count = Kokkos::atomic_fetch_add( &counts( pid ), 1 );
        if ( count < neighbors.extent( 1 ) )
            neighbors( pid, count ) = nid;
    }

    KOKKOS_INLINE_FUNCTION
    void setNeighbor( const int pid, const int nid, const int new_id ) const
    {
        neighbors( pid, nid ) = new_id;
    }
};
 
//---------------------------------------------------------------------------//
 
//---------------------------------------------------------------------------//
 
namespace Impl
{
 
//---------------------------------------------------------------------------//
// Verlet List Builder
//---------------------------------------------------------------------------//
template <class DeviceType, class RandomAccessPositionType, class RadiusType,
          class AlgorithmTag, class LayoutTag, class BuildOpTag,
          class ArrayType, std::size_t NumSpaceDim>
struct VerletListBuilder
{
    static constexpr std::size_t num_space_dim = NumSpaceDim;
 
    // Types.
    using device = DeviceType;
    using PositionValueType = typename RandomAccessPositionType::value_type;
    using memory_space = typename device::memory_space;
    using execution_space = typename device::execution_space;
 
    // List data.
    VerletListData<memory_space, LayoutTag> _data;
    // Background squared neighbor cutoff.
    PositionValueType rsqr;
    // Fixed or per-particle neighbor radius.
    RadiusType radius;
 
    // Positions.
    RandomAccessPositionType _position;
    std::size_t pid_begin, pid_end;
 
    // Binning Data.
    BinningData<memory_space> bin_data_1d;
    LinkedCellList<memory_space, PositionValueType, num_space_dim>
        linked_cell_list;
 
    // Check to count or refill.
    bool refill;
    bool count;
 
    // Maximum allocated neighbors per particle
    std::size_t alloc_n;
 
    // Constructor with a single cutoff radius.
    template <class PositionType>
    VerletListBuilder( PositionType positions, const std::size_t begin,
                       const std::size_t end,
                       const RadiusType neighborhood_radius,
                       const PositionValueType cell_size_ratio,
                       const ArrayType grid_min, const ArrayType grid_max,
                       const std::size_t max_neigh )
        : pid_begin( begin )
        , pid_end( end )
        , alloc_n( max_neigh )
    {
        init( positions, neighborhood_radius, cell_size_ratio, grid_min,
              grid_max );
        // This value is not currently used, but set to be consistent with the
        // variable cutoff case below.
        radius = neighborhood_radius;
    }
 
    // Constructor with a background radius (used for the LinkedCellList) and a
    // per-particle radius.
    template <class PositionType>
    VerletListBuilder( PositionType positions, const std::size_t begin,
                       const std::size_t end,
                       const PositionValueType background_radius,
                       const RadiusType neighborhood_radius,
                       const PositionValueType cell_size_ratio,
                       const ArrayType grid_min, const ArrayType grid_max,
                       const std::size_t max_neigh )
        : pid_begin( begin )
        , pid_end( end )
        , alloc_n( max_neigh )
    {
        assert( size( positions ) == size( neighborhood_radius ) );
        init( positions, background_radius, cell_size_ratio, grid_min,
              grid_max );
 
        // Store a shallow copy (not squared).
        // TODO: for cases where the radii never change, this could be better
        // optimized with a deep copy of the squared radius instead.
        radius = neighborhood_radius;
    }
 
    template <class PositionType>
    void init( PositionType positions,
               const PositionValueType neighborhood_radius,
               const PositionValueType cell_size_ratio,
               const ArrayType grid_min, const ArrayType grid_max )
    {
        count = true;
        refill = false;
 
        // Create the count view.
        _data.counts = Kokkos::View<int*, memory_space>( "num_neighbors",
                                                         size( positions ) );
 
        // Make a guess for the number of neighbors per particle for 2D lists.
        initCounts( LayoutTag() );
 
        // Shallow copy for random access read-only memory.
        _position = positions;
 
        double grid_size = cell_size_ratio * neighborhood_radius;
        ArrayType grid_delta;
        for ( std::size_t d = 0; d < num_space_dim; ++d )
            grid_delta[d] = grid_size;
 
        // Bin the particles in the grid. Don't actually sort them but make a
        // permutation vector. Note that we are binning all particles here and
        // not just the requested range. This is because all particles are
        // treated as candidates for neighbors.
        linked_cell_list =
            createLinkedCellList( _position, grid_delta, grid_min, grid_max,
                                  neighborhood_radius, cell_size_ratio );
        bin_data_1d = linked_cell_list.binningData();
 
        // We will use the square of the distance for neighbor determination.
        rsqr = neighborhood_radius * neighborhood_radius;
    }
 
    // Check if particle pair i-j is within cutoff, potentially with variable
    // radii.
    KOKKOS_INLINE_FUNCTION auto withinCutoff( [[maybe_unused]] const int i,
                                              const double dist_sqr ) const
    {
        // Square the radius on the fly if using a per-particle field to avoid a
        // deep copy.
        if constexpr ( is_slice<RadiusType>::value ||
                       Kokkos::is_view<RadiusType>::value )
            return dist_sqr <= radius( i ) * radius( i );
        // This value is already squared.
        else
            return dist_sqr <= rsqr;
    }
 
    // Check if potential neighbor j is NOT within cutoff, meaning particle i
    // should add instead for symmetry.
    KOKKOS_INLINE_FUNCTION auto
    neighborNotWithinCutoff( [[maybe_unused]] const int j,
                             [[maybe_unused]] const double dist_sqr ) const
    {
        // This neighbor needs to be added if they will not find this particle.
        if constexpr ( is_slice<RadiusType>::value ||
                       Kokkos::is_view<RadiusType>::value )
        {
            return dist_sqr >= radius( j ) * radius( j );
        }
        else
        {
            // For a fixed radius, this will never occur.
            return false;
        }
    }
 
    // Count neighbors, with consideration for self particle i and neighbor j.
    KOKKOS_INLINE_FUNCTION auto countNeighbor( const int i, const int j,
                                               const double dist_sqr ) const
    {
        int c = 0;
        // Always add self if within cutoff.
        if ( withinCutoff( i, dist_sqr ) )
        {
            c++;
            // Add neighbor if they will not find this particle.
            if ( neighborNotWithinCutoff( j, dist_sqr ) )
                c++;
        }
        return c;
    }
 
    // Add neighbors, with consideration for self particle i and neighbor j.
    KOKKOS_INLINE_FUNCTION void addNeighbor( const int i, const int j,
                                             const double dist_sqr ) const
    {
        // Always add self if within cutoff.
        if ( withinCutoff( i, dist_sqr ) )
        {
            _data.addNeighbor( i, j );
 
            // Add neighbor if they will not find this particle.
            if ( neighborNotWithinCutoff( j, dist_sqr ) )
                _data.addNeighbor( j, i );
        }
    }
 
    // Neighbor count team operator (only used for CSR lists).
    struct CountNeighborsTag
    {
    };
    using CountNeighborsPolicy =
        Kokkos::TeamPolicy<execution_space, CountNeighborsTag,
                           Kokkos::IndexType<int>,
                           Kokkos::Schedule<Kokkos::Dynamic>>;
 
    KOKKOS_INLINE_FUNCTION
    void
    operator()( const CountNeighborsTag&,
                const typename CountNeighborsPolicy::member_type& team ) const
    {
        // The league rank of the team is the cardinal cell index we are
        // working on.
        int cell = team.league_rank();
 
        // Get the stencil for this cell.
        Kokkos::Array<int, num_space_dim> min;
        Kokkos::Array<int, num_space_dim> max;
        linked_cell_list.getStencilCells( cell, min, max );
 
        // Operate on the particles in the bin.
        std::size_t b_offset = bin_data_1d.binOffset( cell );
        Kokkos::parallel_for(
            Kokkos::TeamThreadRange( team, 0, bin_data_1d.binSize( cell ) ),
            [&]( const int bi )
            {
                // Get the true particle id. The binned particle index is the
                // league rank of the team.
                std::size_t pid = linked_cell_list.permutation( bi + b_offset );
 
                if ( ( pid >= pid_begin ) && ( pid < pid_end ) )
                {
                    // Cache the particle coordinates.
                    Kokkos::Array<PositionValueType, num_space_dim> xp;
                    for ( std::size_t d = 0; d < num_space_dim; ++d )
                        xp[d] = _position( pid, d );
 
                    // Loop over the cell stencil.
                    int stencil_count =
                        linkedcell_stencil_count( team, min, max, xp, pid );
 
                    Kokkos::single( Kokkos::PerThread( team ), [&]()
                                    { _data.counts( pid ) = stencil_count; } );
                }
            } );
    }
 
    template <class TeamType, std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<3 == NSD, int>
    linkedcell_stencil_count( TeamType team, const Kokkos::Array<int, 3> min,
                              const Kokkos::Array<int, 3> max,
                              const Kokkos::Array<PositionValueType, 3> xp,
                              const std::size_t pid ) const
    {
        int stencil_count = 0;
        Kokkos::Array<int, 3> ijk;
        for ( int i = min[0]; i < max[0]; ++i )
            for ( int j = min[1]; j < max[1]; ++j )
                for ( int k = min[2]; k < max[2]; ++k )
                {
                    ijk = { i, j, k };
                    linkedcell_bin( team, ijk, xp, pid, stencil_count );
                }
        return stencil_count;
    }
 
    template <class TeamType, std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<2 == NSD, int>
    linkedcell_stencil_count( TeamType team, const Kokkos::Array<int, 2> min,
                              const Kokkos::Array<int, 2> max,
                              const Kokkos::Array<PositionValueType, 2> xp,
                              const std::size_t pid ) const
    {
        int stencil_count = 0;
        Kokkos::Array<int, 2> ij;
        for ( int i = min[0]; i < max[0]; ++i )
            for ( int j = min[1]; j < max[1]; ++j )
            {
                ij = { i, j };
                linkedcell_bin( team, ij, xp, pid, stencil_count );
            }
        return stencil_count;
    }
 
    template <class TeamType>
    KOKKOS_INLINE_FUNCTION void
    linkedcell_bin( TeamType team, const Kokkos::Array<int, num_space_dim> ijk,
                    const Kokkos::Array<PositionValueType, num_space_dim> xp,
                    const std::size_t pid, int& stencil_count ) const
    {
        // See if we should actually check this box for neighbors.
        if ( withinCutoff(
                 pid, linked_cell_list.cellStencil().grid.minDistanceToPoint(
                          xp, ijk ) ) )
        {
            std::size_t n_offset = linked_cell_list.binOffset( ijk );
            std::size_t num_n = linked_cell_list.binSize( ijk );
 
            // Check the particles in this bin to see if they are neighbors. If
            // they are add to the count for this bin.
            int cell_count = 0;
            neighbor_reduce( team, pid, xp, n_offset, num_n, cell_count,
                             BuildOpTag() );
            stencil_count += cell_count;
        }
    }
 
    // Neighbor count team vector loop (only used for CSR lists).
    KOKKOS_INLINE_FUNCTION void
    neighbor_reduce( const typename CountNeighborsPolicy::member_type& team,
                     const std::size_t pid,
                     const Kokkos::Array<double, num_space_dim> xp,
                     const int n_offset, const int num_n, int& cell_count,
                     TeamVectorOpTag ) const
    {
        Kokkos::parallel_reduce(
            Kokkos::ThreadVectorRange( team, num_n ),
            [&]( const int n, int& local_count )
            { neighbor_kernel( pid, xp, n_offset, n, local_count ); },
            cell_count );
    }
 
    // Neighbor count serial loop (only used for CSR lists).
    KOKKOS_INLINE_FUNCTION
    void neighbor_reduce( const typename CountNeighborsPolicy::member_type,
                          const std::size_t pid,
                          const Kokkos::Array<double, num_space_dim> xp,
                          const int n_offset, const int num_n, int& cell_count,
                          TeamOpTag ) const
    {
        for ( int n = 0; n < num_n; n++ )
            neighbor_kernel( pid, xp, n_offset, n, cell_count );
    }
 
    // Neighbor count kernel
    KOKKOS_INLINE_FUNCTION
    void neighbor_kernel( const int pid,
                          const Kokkos::Array<double, num_space_dim> xp,
                          const int n_offset, const int n,
                          int& local_count ) const
    {
        //  Get the true id of the candidate  neighbor.
        std::size_t nid = linked_cell_list.permutation( n_offset + n );
 
        // Cache the candidate neighbor particle coordinates.
        Kokkos::Array<double, num_space_dim> xn;
        for ( std::size_t d = 0; d < num_space_dim; ++d )
            xn[d] = _position( nid, d );
 
        // If this could be a valid neighbor, continue.
        if ( NeighborDiscriminator<AlgorithmTag>::isValid( pid, xp, nid, xn ) )
        {
            // Calculate the distance between the particle and its candidate
            // neighbor.
            PositionValueType dist_sqr = 0.0;
            for ( std::size_t d = 0; d < num_space_dim; ++d )
            {
                PositionValueType dx = xp[d] - xn[d];
                dist_sqr += dx * dx;
            }
 
            // If within the cutoff add to the count.
            local_count += countNeighbor( pid, nid, dist_sqr );
        }
    }
 
    // Process the CSR counts by computing offsets and allocating the neighbor
    // list.
    template <class KokkosMemorySpace>
    struct OffsetScanOp
    {
        using kokkos_mem_space = KokkosMemorySpace;
        Kokkos::View<int*, kokkos_mem_space> counts;
        Kokkos::View<int*, kokkos_mem_space> offsets;
        KOKKOS_INLINE_FUNCTION
        void operator()( const int i, int& update, const bool final_pass ) const
        {
            if ( final_pass )
                offsets( i ) = update;
            update += counts( i );
        }
    };
 
    void initCounts( VerletLayoutCSR ) {}
 
    void initCounts( VerletLayout2D )
    {
        if ( alloc_n > 0 )
        {
            count = false;
 
            _data.neighbors = Kokkos::View<int**, memory_space>(
                Kokkos::ViewAllocateWithoutInitializing( "neighbors" ),
                _data.counts.size(), alloc_n );
        }
    }
 
    void processCounts( VerletLayoutCSR )
    {
        // Allocate offsets.
        _data.offsets = Kokkos::View<int*, memory_space>(
            Kokkos::ViewAllocateWithoutInitializing( "neighbor_offsets" ),
            _data.counts.size() );
 
        // Calculate offsets from counts and the total number of counts.
        OffsetScanOp<memory_space> offset_op;
        offset_op.counts = _data.counts;
        offset_op.offsets = _data.offsets;
        int total_num_neighbor;
        Kokkos::RangePolicy<execution_space> range_policy(
            0, _data.counts.extent( 0 ) );
        Kokkos::parallel_scan( "Cabana::VerletListBuilder::offset_scan",
                               range_policy, offset_op, total_num_neighbor );
        Kokkos::fence();
 
        // Allocate the neighbor list.
        _data.neighbors = Kokkos::View<int*, memory_space>(
            Kokkos::ViewAllocateWithoutInitializing( "neighbors" ),
            total_num_neighbor );
 
        // Reset the counts. We count again when we fill.
        Kokkos::deep_copy( _data.counts, 0 );
    }
 
    // Process 2D counts by computing the maximum number of neighbors and
    // reallocating the 2D data structure if needed.
    void processCounts( VerletLayout2D )
    {
        // Calculate the maximum number of neighbors.
        auto counts = _data.counts;
        int max;
        Kokkos::Max<int> max_reduce( max );
        Kokkos::parallel_reduce(
            "Cabana::VerletListBuilder::reduce_max",
            Kokkos::RangePolicy<execution_space>( 0, _data.counts.size() ),
            KOKKOS_LAMBDA( const int i, int& value ) {
                if ( counts( i ) > value )
                    value = counts( i );
            },
            max_reduce );
        Kokkos::fence();
        _data.max_n = static_cast<std::size_t>( max );
 
        // Reallocate the neighbor list if previous size is exceeded.
        if ( count or _data.max_n > _data.neighbors.extent( 1 ) )
        {
            refill = true;
            Kokkos::deep_copy( _data.counts, 0 );
            _data.neighbors = Kokkos::View<int**, memory_space>(
                Kokkos::ViewAllocateWithoutInitializing( "neighbors" ),
                _data.counts.size(), _data.max_n );
        }
    }
 
    // Neighbor count team operator.
    struct FillNeighborsTag
    {
    };
    using FillNeighborsPolicy =
        Kokkos::TeamPolicy<execution_space, FillNeighborsTag,
                           Kokkos::IndexType<int>,
                           Kokkos::Schedule<Kokkos::Dynamic>>;
    KOKKOS_INLINE_FUNCTION
    void
    operator()( const FillNeighborsTag&,
                const typename FillNeighborsPolicy::member_type& team ) const
    {
        // The league rank of the team is the cardinal cell index we are
        // working on.
        int cell = team.league_rank();
 
        // Get the stencil for this cell.
        Kokkos::Array<int, num_space_dim> min;
        Kokkos::Array<int, num_space_dim> max;
        linked_cell_list.getStencilCells( cell, min, max );
 
        // Operate on the particles in the bin.
        std::size_t b_offset = bin_data_1d.binOffset( cell );
        Kokkos::parallel_for(
            Kokkos::TeamThreadRange( team, 0, bin_data_1d.binSize( cell ) ),
            [&]( const int bi )
            {
                // Get the true particle id. The binned particle index is the
                // league rank of the team.
                std::size_t pid = linked_cell_list.permutation( bi + b_offset );
 
                if ( ( pid >= pid_begin ) && ( pid < pid_end ) )
                {
                    // Cache the particle coordinates.
                    Kokkos::Array<PositionValueType, num_space_dim> xp;
                    for ( std::size_t d = 0; d < num_space_dim; ++d )
                        xp[d] = _position( pid, d );
 
                    linkedcell_stencil_fill( team, min, max, xp, pid );
                }
            } );
    }
 
    template <class TeamType, std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<3 == NSD, void>
    linkedcell_stencil_fill( TeamType team, const Kokkos::Array<int, 3> min,
                             const Kokkos::Array<int, 3> max,
                             const Kokkos::Array<PositionValueType, 3> xp,
                             const std::size_t pid ) const
    {
        Kokkos::Array<int, 3> ijk;
        for ( int i = min[0]; i < max[0]; ++i )
            for ( int j = min[1]; j < max[1]; ++j )
                for ( int k = min[2]; k < max[2]; ++k )
                {
                    ijk = { i, j, k };
                    linkedcell_bin( team, ijk, xp, pid );
                }
    }
 
    template <class TeamType, std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<2 == NSD, void>
    linkedcell_stencil_fill( TeamType team, const Kokkos::Array<int, 2> min,
                             const Kokkos::Array<int, 2> max,
                             const Kokkos::Array<PositionValueType, 2> xp,
                             const std::size_t pid ) const
    {
        Kokkos::Array<int, 2> ij;
        for ( int i = min[0]; i < max[0]; ++i )
            for ( int j = min[1]; j < max[1]; ++j )
            {
                ij = { i, j };
                linkedcell_bin( team, ij, xp, pid );
            }
    }
 
    template <class TeamType>
    KOKKOS_INLINE_FUNCTION void
    linkedcell_bin( TeamType team, const Kokkos::Array<int, num_space_dim> ijk,
                    const Kokkos::Array<PositionValueType, num_space_dim> xp,
                    const std::size_t pid ) const
    {
        // See if we should actually check this box for neighbors.
        if ( withinCutoff(
                 pid, linked_cell_list.cellStencil().grid.minDistanceToPoint(
                          xp, ijk ) ) )
        {
            // Check the particles in this bin to see if they are neighbors.
            std::size_t n_offset = linked_cell_list.binOffset( ijk );
            int num_n = linked_cell_list.binSize( ijk );
            neighbor_for( team, pid, xp, n_offset, num_n, BuildOpTag() );
        }
    }
 
    // Neighbor fill team vector loop.
    KOKKOS_INLINE_FUNCTION void
    neighbor_for( const typename FillNeighborsPolicy::member_type& team,
                  const std::size_t pid,
                  const Kokkos::Array<double, num_space_dim> xp,
                  const int n_offset, const int num_n, TeamVectorOpTag ) const
    {
        Kokkos::parallel_for( Kokkos::ThreadVectorRange( team, num_n ),
                              [&]( const int n )
                              { neighbor_kernel( pid, xp, n_offset, n ); } );
    }
 
    // Neighbor fill serial loop.
    KOKKOS_INLINE_FUNCTION
    void neighbor_for( const typename FillNeighborsPolicy::member_type team,
                       const std::size_t pid,
                       const Kokkos::Array<double, num_space_dim> xp,
                       const int n_offset, const int num_n, TeamOpTag ) const
    {
        for ( int n = 0; n < num_n; n++ )
            Kokkos::single( Kokkos::PerThread( team ), [&]()
                            { neighbor_kernel( pid, xp, n_offset, n ); } );
    }
 
    // Neighbor fill kernel.
    KOKKOS_INLINE_FUNCTION
    void neighbor_kernel( const int pid,
                          const Kokkos::Array<double, num_space_dim> xp,
                          const int n_offset, const int n ) const
    {
        //  Get the true id of the candidate neighbor.
        std::size_t nid = linked_cell_list.permutation( n_offset + n );
 
        // Cache the candidate neighbor particle coordinates.
        Kokkos::Array<double, num_space_dim> xn;
        for ( std::size_t d = 0; d < num_space_dim; ++d )
            xn[d] = _position( nid, d );
 
        // If this could be a valid neighbor, continue.
        if ( NeighborDiscriminator<AlgorithmTag>::isValid( pid, xp, nid, xn ) )
        {
            // Calculate the distance between the particle and its candidate
            // neighbor.
            PositionValueType dist_sqr = 0.0;
            for ( std::size_t d = 0; d < num_space_dim; ++d )
            {
                PositionValueType dx = xp[d] - xn[d];
                dist_sqr += dx * dx;
            }
 
            // If within the cutoff increment the neighbor count and add as a
            // neighbor at that index.
            addNeighbor( pid, nid, dist_sqr );
        }
    }
};
 
// Builder creation functions. This is only necessary to define the different
// random access types.
template <std::size_t NumSpaceDim, class DeviceType, class AlgorithmTag,
          class LayoutTag, class BuildOpTag, class ArrayType,
          class PositionType>
auto createVerletListBuilder(
    PositionType x, const std::size_t begin, const std::size_t end,
    const typename PositionType::value_type radius,
    const typename PositionType::value_type cell_size_ratio,
    const ArrayType grid_min, const ArrayType grid_max,
    const std::size_t max_neigh,
    typename std::enable_if<( is_slice<PositionType>::value ), int>::type* = 0 )
{
    using RandomAccessPositionType = typename PositionType::random_access_slice;
    return VerletListBuilder<DeviceType, RandomAccessPositionType,
                             typename PositionType::value_type, AlgorithmTag,
                             LayoutTag, BuildOpTag, ArrayType, NumSpaceDim>(
        x, begin, end, radius, cell_size_ratio, grid_min, grid_max, max_neigh );
}
 
template <std::size_t NumSpaceDim, class DeviceType, class AlgorithmTag,
          class LayoutTag, class BuildOpTag, class ArrayType,
          class PositionType>
auto createVerletListBuilder(
    PositionType x, const std::size_t begin, const std::size_t end,
    const typename PositionType::value_type radius,
    const typename PositionType::value_type cell_size_ratio,
    const ArrayType grid_min, const ArrayType grid_max,
    const std::size_t max_neigh,
    typename std::enable_if<( Kokkos::is_view<PositionType>::value ),
                            int>::type* = 0 )
{
    using memory_space = typename DeviceType::memory_space;
    using RandomAccessPositionType =
        Kokkos::View<typename PositionType::data_type, memory_space,
                     Kokkos::MemoryTraits<Kokkos::RandomAccess>>;
    return VerletListBuilder<DeviceType, RandomAccessPositionType,
                             typename PositionType::value_type, AlgorithmTag,
                             LayoutTag, BuildOpTag, ArrayType, NumSpaceDim>(
        x, begin, end, radius, cell_size_ratio, grid_min, grid_max, max_neigh );
}
 
template <std::size_t NumSpaceDim, class DeviceType, class AlgorithmTag,
          class LayoutTag, class BuildOpTag, class PositionType,
          class RadiusType, class ArrayType>
auto createVerletListBuilder(
    PositionType x, const std::size_t begin, const std::size_t end,
    const typename PositionType::value_type background_radius,
    const RadiusType radius,
    const typename PositionType::value_type cell_size_ratio,
    const ArrayType grid_min, const ArrayType grid_max,
    const std::size_t max_neigh,
    typename std::enable_if<( is_slice<PositionType>::value ), int>::type* = 0 )
{
    using RandomAccessPositionType = typename PositionType::random_access_slice;
    return VerletListBuilder<DeviceType, RandomAccessPositionType, RadiusType,
                             AlgorithmTag, LayoutTag, BuildOpTag, ArrayType,
                             NumSpaceDim>( x, begin, end, background_radius,
                                           radius, cell_size_ratio, grid_min,
                                           grid_max, max_neigh );
}
 
template <std::size_t NumSpaceDim, class DeviceType, class AlgorithmTag,
          class LayoutTag, class BuildOpTag, class PositionType,
          class RadiusType, class ArrayType>
auto createVerletListBuilder(
    PositionType x, const std::size_t begin, const std::size_t end,
    const typename PositionType::value_type background_radius,
    const RadiusType radius,
    const typename PositionType::value_type cell_size_ratio,
    const ArrayType grid_min, const ArrayType grid_max,
    const std::size_t max_neigh,
    typename std::enable_if<( Kokkos::is_view<PositionType>::value ),
                            int>::type* = 0 )
{
    using RandomAccessPositionType =
        Kokkos::View<typename PositionType::value_type**, DeviceType,
                     Kokkos::MemoryTraits<Kokkos::RandomAccess>>;
    return VerletListBuilder<DeviceType, RandomAccessPositionType, RadiusType,
                             AlgorithmTag, LayoutTag, BuildOpTag, ArrayType,
                             NumSpaceDim>( x, begin, end, background_radius,
                                           radius, cell_size_ratio, grid_min,
                                           grid_max, max_neigh );
}
 
//---------------------------------------------------------------------------//

} // end namespace Impl
 
//---------------------------------------------------------------------------//
template <class MemorySpace, class AlgorithmTag, class LayoutTag,
          class BuildTag = TeamVectorOpTag, std::size_t NumSpaceDim = 3>
class VerletList
{
  public:
    static_assert( Kokkos::is_memory_space<MemorySpace>::value, "" );

    static constexpr std::size_t num_space_dim = NumSpaceDim;

    using memory_space = MemorySpace;

    using execution_space = typename memory_space::execution_space;

    VerletListData<memory_space, LayoutTag> _data;

    VerletList() = default;

    template <class PositionType,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    VerletList(
        PositionType x, const std::size_t begin, const std::size_t end,
        const typename PositionType::value_type neighborhood_radius,
        const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
        const ArrayType<typename PositionType::value_type, num_space_dim,
                        Args...>
            grid_min,
        const ArrayType<typename PositionType::value_type, num_space_dim,
                        Args...>
            grid_max,
#else
        const ArrayType<typename PositionType::value_type, num_space_dim>
            grid_min,
        const ArrayType<typename PositionType::value_type, num_space_dim>
            grid_max,
#endif
        const std::size_t max_neigh = 0,
        typename std::enable_if<( is_slice<PositionType>::value ||
                                  Kokkos::is_view<PositionType>::value ),
                                int>::type* = 0 )
    {
        build( x, begin, end, neighborhood_radius, cell_size_ratio, grid_min,
               grid_max, max_neigh );
    }

    template <class PositionType,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    void
    build( PositionType x, const std::size_t begin, const std::size_t end,
           const typename PositionType::value_type neighborhood_radius,
           const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
           const ArrayType<typename PositionType::value_type, num_space_dim,
                           Args...>
               grid_min,
           const ArrayType<typename PositionType::value_type, num_space_dim,
                           Args...>
               grid_max,
#else
           const ArrayType<typename PositionType::value_type, num_space_dim>
               grid_min,
           const ArrayType<typename PositionType::value_type, num_space_dim>
               grid_max,
#endif
           const std::size_t max_neigh = 0,
           typename std::enable_if<( is_slice<PositionType>::value ||
                                     Kokkos::is_view<PositionType>::value ),
                                   int>::type* = 0 )
    {
        // Use the default execution space.
        build( execution_space{}, x, begin, end, neighborhood_radius,
               cell_size_ratio, grid_min, grid_max, max_neigh );
    }

    template <class PositionType, class RadiusType,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    VerletList(
        PositionType x, const std::size_t begin, const std::size_t end,
        const typename PositionType::value_type background_radius,
        RadiusType neighborhood_radius,
        const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
        const ArrayType<typename PositionType::value_type, num_space_dim,
                        Args...>
            grid_min,
        const ArrayType<typename PositionType::value_type, num_space_dim,
                        Args...>
            grid_max,
#else
        const ArrayType<typename PositionType::value_type, num_space_dim>
            grid_min,
        const ArrayType<typename PositionType::value_type, num_space_dim>
            grid_max,
#endif
        const std::size_t max_neigh = 0,
        typename std::enable_if<( is_slice<PositionType>::value ||
                                  Kokkos::is_view<PositionType>::value ),
                                int>::type* = 0 )
    {
        build( x, begin, end, background_radius, neighborhood_radius,
               cell_size_ratio, grid_min, grid_max, max_neigh );
    }

    template <class PositionType,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    VerletList(
        PositionType x,
        const typename PositionType::value_type neighborhood_radius,
        const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
        const ArrayType<typename PositionType::value_type, num_space_dim,
                        Args...>
            grid_min,
        const ArrayType<typename PositionType::value_type, num_space_dim,
                        Args...>
            grid_max,
#else
        const ArrayType<typename PositionType::value_type, num_space_dim>
            grid_min,
        const ArrayType<typename PositionType::value_type, num_space_dim>
            grid_max,
#endif
        const std::size_t max_neigh = 0,
        typename std::enable_if<( is_slice<PositionType>::value ||
                                  Kokkos::is_view<PositionType>::value ),
                                int>::type* = 0 )
    {
        build( x, neighborhood_radius, cell_size_ratio, grid_min, grid_max,
               max_neigh );
    }

    template <class PositionType>
    VerletList(
        PositionType x, const std::size_t begin, const std::size_t end,
        const typename PositionType::value_type neighborhood_radius,
        const typename PositionType::value_type cell_size_ratio,
        const typename PositionType::value_type grid_min[num_space_dim],
        const typename PositionType::value_type grid_max[num_space_dim],
        const std::size_t max_neigh = 0,
        typename std::enable_if<( is_slice<PositionType>::value ||
                                  Kokkos::is_view<PositionType>::value ),
                                int>::type* = 0 )
    {
        build( x, begin, end, neighborhood_radius, cell_size_ratio, grid_min,
               grid_max, max_neigh );
    }

    template <class PositionType>
    VerletList(
        PositionType x,
        const typename PositionType::value_type neighborhood_radius,
        const typename PositionType::value_type cell_size_ratio,
        const typename PositionType::value_type grid_min[num_space_dim],
        const typename PositionType::value_type grid_max[num_space_dim],
        const std::size_t max_neigh = 0,
        typename std::enable_if<( is_slice<PositionType>::value ||
                                  Kokkos::is_view<PositionType>::value ),
                                int>::type* = 0 )
    {
        build( x, neighborhood_radius, cell_size_ratio, grid_min, grid_max,
               max_neigh );
    }

    template <class PositionType, class RadiusType>
    VerletList(
        PositionType x, const std::size_t begin, const std::size_t end,
        const typename PositionType::value_type background_radius,
        RadiusType neighborhood_radius,
        const typename PositionType::value_type cell_size_ratio,
        const typename PositionType::value_type grid_min[num_space_dim],
        const typename PositionType::value_type grid_max[num_space_dim],
        const std::size_t max_neigh = 0,
        typename std::enable_if<( is_slice<PositionType>::value ||
                                  Kokkos::is_view<PositionType>::value ),
                                int>::type* = 0 )
    {
        build( x, begin, end, background_radius, neighborhood_radius,
               cell_size_ratio, grid_min, grid_max, max_neigh );
    }

    template <class PositionType, class RadiusType,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    void build( PositionType x, const std::size_t begin, const std::size_t end,
                const typename PositionType::value_type background_radius,
                RadiusType neighborhood_radius,
                const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
                const ArrayType<typename PositionType::value_type,
                                num_space_dim, Args...>
                    grid_min,
                const ArrayType<typename PositionType::value_type,
                                num_space_dim, Args...>
                    grid_max,
#else
                const ArrayType<typename PositionType::value_type,
                                num_space_dim>
                    grid_min,
                const ArrayType<typename PositionType::value_type,
                                num_space_dim>
                    grid_max,
#endif
                const std::size_t max_neigh = 0 )
    {
        build( execution_space{}, x, begin, end, background_radius,
               neighborhood_radius, cell_size_ratio, grid_min, grid_max,
               max_neigh );
    }

    template <class PositionType,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    void
    build( PositionType x,
           const typename PositionType::value_type neighborhood_radius,
           const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
           const ArrayType<typename PositionType::value_type, num_space_dim,
                           Args...>
               grid_min,
           const ArrayType<typename PositionType::value_type, num_space_dim,
                           Args...>
               grid_max,
#else
           const ArrayType<typename PositionType::value_type, num_space_dim>
               grid_min,
           const ArrayType<typename PositionType::value_type, num_space_dim>
               grid_max,
#endif
           const std::size_t max_neigh = 0,
           typename std::enable_if<( is_slice<PositionType>::value ||
                                     Kokkos::is_view<PositionType>::value ),
                                   int>::type* = 0 )
    {
        build( x, 0, size( x ), neighborhood_radius, cell_size_ratio, grid_min,
               grid_max, max_neigh );
    }

    template <class PositionType>
    void
    build( PositionType x, const std::size_t begin, const std::size_t end,
           const typename PositionType::value_type neighborhood_radius,
           const typename PositionType::value_type cell_size_ratio,
           const typename PositionType::value_type grid_min[num_space_dim],
           const typename PositionType::value_type grid_max[num_space_dim],
           const std::size_t max_neigh = 0,
           typename std::enable_if<( is_slice<PositionType>::value ||
                                     Kokkos::is_view<PositionType>::value ),
                                   int>::type* = 0 )
    {
        // Use the default execution space.
        build( execution_space{}, x, begin, end, neighborhood_radius,
               cell_size_ratio, grid_min, grid_max, max_neigh );
    }

    template <class PositionType, class ExecutionSpace>
    void
    build( ExecutionSpace exec_space, PositionType x, const std::size_t begin,
           const std::size_t end,
           const typename PositionType::value_type neighborhood_radius,
           const typename PositionType::value_type cell_size_ratio,
           const typename PositionType::value_type grid_min[num_space_dim],
           const typename PositionType::value_type grid_max[num_space_dim],
           const std::size_t max_neigh = 0,
           typename std::enable_if<( is_slice<PositionType>::value ||
                                     Kokkos::is_view<PositionType>::value ),
                                   int>::type* = 0 )
    {
        auto kokkos_min =
            copyArray<typename PositionType::value_type, num_space_dim>(
                grid_min );
        auto kokkos_max =
            copyArray<typename PositionType::value_type, num_space_dim>(
                grid_max );
        build( exec_space, x, begin, end, neighborhood_radius, cell_size_ratio,
               kokkos_min, kokkos_max, max_neigh );
    }

    template <class PositionType, class RadiusType>
    void build( PositionType x, const std::size_t begin, const std::size_t end,
                const typename PositionType::value_type background_radius,
                RadiusType neighborhood_radius,
                const typename PositionType::value_type cell_size_ratio,
                const typename PositionType::value_type grid_min[num_space_dim],
                const typename PositionType::value_type grid_max[num_space_dim],
                const std::size_t max_neigh = 0 )
    {
        // Use the default execution space.
        build( execution_space{}, x, begin, end, background_radius,
               neighborhood_radius, cell_size_ratio, grid_min, grid_max,
               max_neigh );
    }

    template <class PositionType>
    void
    build( PositionType x,
           const typename PositionType::value_type neighborhood_radius,
           const typename PositionType::value_type cell_size_ratio,
           const typename PositionType::value_type grid_min[num_space_dim],
           const typename PositionType::value_type grid_max[num_space_dim],
           const std::size_t max_neigh = 0,
           typename std::enable_if<( is_slice<PositionType>::value ||
                                     Kokkos::is_view<PositionType>::value ),
                                   int>::type* = 0 )
    {
        build( x, 0, size( x ), neighborhood_radius, cell_size_ratio, grid_min,
               grid_max, max_neigh );
    }

    template <class PositionType, class RadiusType, class ExecutionSpace>
    void build( ExecutionSpace exec_space, PositionType x,
                const std::size_t begin, const std::size_t end,
                const typename PositionType::value_type background_radius,
                RadiusType neighborhood_radius,
                const typename PositionType::value_type cell_size_ratio,
                const typename PositionType::value_type grid_min[num_space_dim],
                const typename PositionType::value_type grid_max[num_space_dim],
                const std::size_t max_neigh = 0 )
    {
        auto kokkos_min =
            copyArray<typename PositionType::value_type, num_space_dim>(
                grid_min );
        auto kokkos_max =
            copyArray<typename PositionType::value_type, num_space_dim>(
                grid_max );
        build( exec_space, x, begin, end, background_radius,
               neighborhood_radius, cell_size_ratio, kokkos_min, kokkos_max,
               max_neigh );
    }

    template <class PositionType, class ExecutionSpace,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    void
    build( ExecutionSpace, PositionType x, const std::size_t begin,
           const std::size_t end,
           const typename PositionType::value_type neighborhood_radius,
           const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
           const ArrayType<typename PositionType::value_type, num_space_dim,
                           Args...>
               grid_min,
           const ArrayType<typename PositionType::value_type, num_space_dim,
                           Args...>
               grid_max,
#else
           const ArrayType<typename PositionType::value_type, num_space_dim>
               grid_min,
           const ArrayType<typename PositionType::value_type, num_space_dim>
               grid_max,
#endif
           const std::size_t max_neigh = 0,
           typename std::enable_if<( is_slice<PositionType>::value ||
                                     Kokkos::is_view<PositionType>::value ),
                                   int>::type* = 0 )
    {
        Kokkos::Profiling::ScopedRegion region( "Cabana::VerletList::build" );
 
        static_assert( is_accessible_from<memory_space, ExecutionSpace>{}, "" );
 
        assert( end >= begin );
        assert( end <= size( x ) );
 
        using device_type = Kokkos::Device<ExecutionSpace, memory_space>;
        // Create a builder functor.
        auto builder =
            Impl::createVerletListBuilder<num_space_dim, device_type,
                                          AlgorithmTag, LayoutTag, BuildTag>(
                x, begin, end, neighborhood_radius, cell_size_ratio, grid_min,
                grid_max, max_neigh );
        buildImpl( builder );
    }

    template <class PositionType, class RadiusType, class ExecutionSpace,
              template <class, std::size_t, class...> class ArrayType,
              class... Args>
    void build( ExecutionSpace, PositionType x, const std::size_t begin,
                const std::size_t end,
                const typename PositionType::value_type background_radius,
                RadiusType neighborhood_radius,
                const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
                const ArrayType<typename PositionType::value_type,
                                num_space_dim, Args...>
                    grid_min,
                const ArrayType<typename PositionType::value_type,
                                num_space_dim, Args...>
                    grid_max,
#else
                const ArrayType<typename PositionType::value_type,
                                num_space_dim>
                    grid_min,
                const ArrayType<typename PositionType::value_type,
                                num_space_dim>
                    grid_max,
#endif
                const std::size_t max_neigh = 0 )
    {
        Kokkos::Profiling::ScopedRegion region( "Cabana::VerletList::build" );
 
        static_assert( is_accessible_from<memory_space, ExecutionSpace>{}, "" );
 
        assert( end >= begin );
        assert( end <= size( x ) );
 
        // Create a builder functor.
        using device_type = Kokkos::Device<ExecutionSpace, memory_space>;
        auto builder =
            Impl::createVerletListBuilder<num_space_dim, device_type,
                                          AlgorithmTag, LayoutTag, BuildTag>(
                x, begin, end, background_radius, neighborhood_radius,
                cell_size_ratio, grid_min, grid_max, max_neigh );
        buildImpl( builder );
    }

    template <class BuilderType>
    void buildImpl( BuilderType builder )
    {
        // For each particle in the range check each neighboring bin for
        // neighbor particles. Bins are at least the size of the
        // neighborhood radius so the bin in which the particle resides and
        // any surrounding bins are guaranteed to contain the neighboring
        // particles. For CSR lists, we count, then fill neighbors. For 2D
        // lists, we count and fill at the same time, unless the array size
        // is exceeded, at which point only counting is continued to
        // reallocate and refill.
        typename BuilderType::FillNeighborsPolicy fill_policy(
            builder.bin_data_1d.numBin(), Kokkos::AUTO, 4 );
        if ( builder.count )
        {
            typename BuilderType::CountNeighborsPolicy count_policy(
                builder.bin_data_1d.numBin(), Kokkos::AUTO, 4 );
            Kokkos::parallel_for( "Cabana::VerletList::count_neighbors",
                                  count_policy, builder );
        }
        else
        {
            builder.processCounts( LayoutTag() );
            Kokkos::parallel_for( "Cabana::VerletList::fill_neighbors",
                                  fill_policy, builder );
        }
        Kokkos::fence();
 
        // Process the counts by computing offsets and allocating the neighbor
        // list, if needed.
        builder.processCounts( LayoutTag() );
 
        // For each particle in the range fill (or refill) its part of the
        // neighbor list.
        if ( builder.count or builder.refill )
        {
            Kokkos::parallel_for( "Cabana::VerletList::fill_neighbors",
                                  fill_policy, builder );
            Kokkos::fence();
        }
 
        // Get the data from the builder.
        _data = builder._data;
    }

    KOKKOS_INLINE_FUNCTION
    void setNeighbor( const std::size_t particle_index,
                      const std::size_t neighbor_index,
                      const int new_index ) const
    {
        _data.setNeighbor( particle_index, neighbor_index, new_index );
    }
};
 
//---------------------------------------------------------------------------//
// VerletList creation functions.
//---------------------------------------------------------------------------//

template <class AlgorithmTag, class LayoutTag, class BuildTag,
          class PositionType,
          template <class, std::size_t, class...> class ArrayType,
          std::size_t NumSpaceDim, class... Args>
auto createVerletList(
    PositionType positions, const std::size_t begin, const std::size_t end,
    const typename PositionType::value_type radius,
    const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
    const ArrayType<typename PositionType::value_type, NumSpaceDim, Args...>
        grid_min,
    const ArrayType<typename PositionType::value_type, NumSpaceDim, Args...>
        grid_max,
#else
    const ArrayType<typename PositionType::value_type, NumSpaceDim> grid_min,
    const ArrayType<typename PositionType::value_type, NumSpaceDim> grid_max,
#endif
    const std::size_t max_neigh = 0 )
{
    using memory_space = typename PositionType::memory_space;
    return Cabana::VerletList<memory_space, AlgorithmTag, LayoutTag, BuildTag,
                              NumSpaceDim>( positions, begin, end, radius,
                                            cell_size_ratio, grid_min, grid_max,
                                            max_neigh );
}

template <class AlgorithmTag, class LayoutTag, class BuildTag,
          class PositionType,
          template <class, std::size_t, class...> class ArrayType,
          std::size_t NumSpaceDim, class... Args>
auto createVerletList(
    const PositionType& positions,
    const typename PositionType::value_type radius,
    const typename PositionType::value_type cell_size_ratio,
#ifdef KOKKOS_ENABLE_DEPRECATED_CODE_4
    const ArrayType<typename PositionType::value_type, NumSpaceDim, Args...>
        grid_min,
    const ArrayType<typename PositionType::value_type, NumSpaceDim, Args...>
        grid_max,
#else
    const ArrayType<typename PositionType::value_type, NumSpaceDim> grid_min,
    const ArrayType<typename PositionType::value_type, NumSpaceDim> grid_max,
#endif
    const std::size_t max_neigh = 0 )
{
    using memory_space = typename PositionType::memory_space;
    return Cabana::VerletList<memory_space, AlgorithmTag, LayoutTag, BuildTag,
                              NumSpaceDim>( positions, radius, cell_size_ratio,
                                            grid_min, grid_max, max_neigh );
}
 
//---------------------------------------------------------------------------//
// Neighbor list interface implementation.
//---------------------------------------------------------------------------//
template <class MemorySpace, class AlgorithmTag, class BuildTag,
          std::size_t Dim>
class NeighborList<
    VerletList<MemorySpace, AlgorithmTag, VerletLayoutCSR, BuildTag, Dim>>
{
  public:
    using memory_space = MemorySpace;
    using list_type =
        VerletList<MemorySpace, AlgorithmTag, VerletLayoutCSR, BuildTag, Dim>;

    KOKKOS_INLINE_FUNCTION
    static std::size_t totalNeighbor( const list_type& list )
    {
        // Size of the allocated memory gives total neighbors.
        return list._data.neighbors.extent( 0 );
    }

    KOKKOS_INLINE_FUNCTION
    static std::size_t maxNeighbor( const list_type& list )
    {
        std::size_t num_p = list._data.counts.size();
        return Impl::maxNeighbor( list, num_p );
    }

    KOKKOS_INLINE_FUNCTION
    static std::size_t numNeighbor( const list_type& list,
                                    const std::size_t particle_index )
    {
        return list._data.counts( particle_index );
    }

    KOKKOS_INLINE_FUNCTION
    static std::size_t getNeighbor( const list_type& list,
                                    const std::size_t particle_index,
                                    const std::size_t neighbor_index )
    {
        return list._data.neighbors( list._data.offsets( particle_index ) +
                                     neighbor_index );
    }
};
 
//---------------------------------------------------------------------------//
template <class MemorySpace, class AlgorithmTag, class BuildTag,
          std::size_t Dim>
class NeighborList<
    VerletList<MemorySpace, AlgorithmTag, VerletLayout2D, BuildTag, Dim>>
{
  public:
    using memory_space = MemorySpace;
    using list_type =
        VerletList<MemorySpace, AlgorithmTag, VerletLayout2D, BuildTag, Dim>;

    KOKKOS_INLINE_FUNCTION
    static std::size_t totalNeighbor( const list_type& list )
    {
        std::size_t num_p = list._data.counts.size();
        return Impl::totalNeighbor( list, num_p );
    }

    KOKKOS_INLINE_FUNCTION
    static std::size_t maxNeighbor( const list_type& list )
    {
        // Stored during neighbor search.
        return list._data.max_n;
    }

    KOKKOS_INLINE_FUNCTION
    static std::size_t numNeighbor( const list_type& list,
                                    const std::size_t particle_index )
    {
        return list._data.counts( particle_index );
    }

    KOKKOS_INLINE_FUNCTION
    static std::size_t getNeighbor( const list_type& list,
                                    const std::size_t particle_index,
                                    const std::size_t neighbor_index )
    {
        return list._data.neighbors( particle_index, neighbor_index );
    }
};
 
//---------------------------------------------------------------------------//
 
} // end namespace Cabana
 
#endif // end  CABANA_VERLETLIST_HPP