Application interface

neighborlist(;
    xpositions,
    #or
    positions,
    ypositions=nothing,
    cutoff,
    unitcell=nothing,
    parallel=true,
    show_progress=false,
    nbatches=(0,0),
)

Computes the list of pairs of particles closer to each other than cutoff.

• xpositions (or positions as a shortcut for single-set computations) provides the coordinates of the particles.
• If ypositions is also provided, returns the cross-pair neighbor list between the two sets.
• If unitcell is provided (as a vector of sides or a unit cell matrix), periodic boundary conditions are used.

Example

Compute the neighborlist within a set of Argon atoms, non-periodic:

julia> using CellListMap, PDBTools

julia> x = coor(read_pdb(CellListMap.argon_pdb_file));

julia> neighborlist(xpositions=x, cutoff=8.0, parallel=false)
857-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 20, 3.163779526466901)
 (1, 61, 4.08865164675529)
 (1, 67, 5.939772435456664)
 ⋮
 (78, 88, 7.0061163797598445)
 (88, 54, 7.933654063435482)

And now, considering the system periodic:

julia> using CellListMap, PDBTools

julia> x = coor(read_pdb(CellListMap.argon_pdb_file));

julia> neighborlist(xpositions=x, cutoff=8.0, unitcell=[21.0, 21.0, 21.0], parallel=false)
1143-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 7, 3.3638756414119397)
 (1, 20, 3.163779526466901)
 (1, 47, 6.243868666689442)
 ⋮
 (68, 38, 7.409287768713663)
 (68, 90, 7.8758006026725464)

For cross-pair neighbor lists between two sets:

julia> using CellListMap, PDBTools

julia> x = coor(read_pdb(CellListMap.argon_pdb_file, "index <= 50"));

julia> y = coor(read_pdb(CellListMap.argon_pdb_file, "index > 50"));

julia> CellListMap.neighborlist(xpositions=x, ypositions=y, cutoff=8.0, parallel=false)
439-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 11, 4.08865164675529)
 (1, 17, 5.939772435456664)
 (1, 30, 2.4572288423012236)
 ⋮
 (46, 48, 4.9269093987894745)
 (46, 1, 7.99947286297016)

InPlaceNeighborList(;
    xpositions::AbstractVecOrMat,
    #or
    positions::AbstractVecOrMat,
    ypositions::Union{AbstractVecOrMat,Nothing}=nothing,
    cutoff::T,
    unitcell::Union{AbstractVecOrMat,Nothing}=nothing,
    parallel::Bool=true,
    show_progress::Bool=false,
) where {T}

Function that initializes the InPlaceNeighborList structure, to be used for in-place computation of neighbor lists.

• xpositions (or positions as a shortcut for single-set computations) provides the coordinates of the particles.
• If ypositions is also provided, the neighbor list between the two sets is computed.
• If unitcell is provided, periodic boundary conditions will be used. The unitcell can be a vector of Orthorhombic box sides, or an actual unitcell matrix for general cells.
• If unitcell is not provided (value nothing), no periodic boundary conditions will be considered.

Examples

Here the neighborlist structure is constructed for the first time, and used to compute the neighbor lists with the mutating neighborlist! function:

julia> using CellListMap, StaticArrays

julia> x = rand(SVector{3,Float64}, 10^4);

julia> system = InPlaceNeighborList(xpositions=x, cutoff=0.1, unitcell=[1,1,1])
InPlaceNeighborList with types:
CellList{3, Float64}
Box{OrthorhombicCell, 3, Float64, Float64, 9}
Current list buffer length: 0

julia> neighborlist!(system)
210034-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 357, 0.09922225615002134)
 (1, 488, 0.043487074695938925)
 (1, 2209, 0.017779967072139684)
 ⋮
 (9596, 1653, 0.0897570322108541)
 (9596, 7927, 0.0898266280344037)

The coordinates of the system, its unitcell, or the cutoff can be changed with the update! function. If the number of pairs of the list does not change significantly, the new calculation is minimally allocating, or non-allocating at all, in particular if the computation is run without parallelization:

If the structure is used repeatedly for similar systems, the allocations will vanish, except for minor allocations used in the threading computation (if a non-parallel computation is executed, the allocations will vanish completely).

update!(
    sys::AbstractParticleSystem;
    xpositions = nothing,
    #or
    positions = nothing,
    ypositions = nothing,
    cutoff = nothing,
    unitcell = nothing,
    parallel = nothing,
)

Update one or more properties of sys in a single call. Only the keyword arguments that are provided (i.e. not nothing) are updated.

• xpositions (or positions as a shortcut for single-set systems): new coordinates for the first (or only) set of particles. Accepts the same types as the ParticleSystem constructor: a Vector{SVector}, a vector of vectors, or an (D, N) matrix. The internal storage is resized automatically if the number of particles changes.

• ypositions: new coordinates for the second set of particles (only valid for two-set systems). Same accepted types as xpositions.

• cutoff: new cutoff distance.

• unitcell: new unit cell. Must be of the same cell type (orthorhombic or triclinic) as the original system. Manual updating of non-periodic systems is not allowed.

• parallel: whether to use multi-threading (true or false).

Example

julia> using CellListMap, StaticArrays

julia> sys = ParticleSystem(
           xpositions = rand(SVector{3,Float64}, 100),
           unitcell = [1,1,1],
           cutoff = 0.1,
           output = 0.0,
       );

julia> new_x = rand(SVector{3,Float64}, 100);

julia> update!(sys; xpositions=new_x, cutoff=0.2, unitcell=[2,2,2], parallel=false);

update!(
    system::InPlaceNeighborList;
    xpositions = nothing,
    #or
    positions = nothing,
    ypositions = nothing,
    cutoff = nothing,
    unitcell = nothing,
)

Updates a InPlaceNeighborList system. Only the keyword arguments that are provided (i.e. not nothing) are updated. positions is accepted as an alias for xpositions.

Examples

For self-pairs computations

julia> x = rand(SVector{3,Float64}, 10^3);

julia> system = InPlaceNeighborList(xpositions=x, cutoff=0.1)
InPlaceNeighborList with types:
CellList{3, Float64}
Box{NonPeriodicCell, 3, Float64, Float64, 9}
Current list buffer length: 0

julia> neighborlist!(system);

julia> new_x = rand(SVector{3,Float64}, 10^3);

julia> update!(system; xpositions=new_x, cutoff=0.05)
InPlaceNeighborList with types:
CellList{3, Float64}
Box{NonPeriodicCell, 3, Float64, Float64, 9}
Current list buffer length: 1826

julia> neighborlist!(system)
224-element Vector{Tuple{Int64, Int64, Float64}}:
 (25, 486, 0.03897345036790646)
 ⋮
 (723, 533, 0.04795768478723409)
 (868, 920, 0.042087156715720137)

neighborlist(system::InPlaceNeighborList)

Computes the neighbor list in-place, given a InPlaceNeighborList system.

Example

In the following example, we compute the neighbor list of a set of random particles, and then we change the coordinates and recompute the neighbor list without reallocations (or with minimal allocations if run in parallel):

julia> using CellListMap, StaticArrays

julia> x = rand(SVector{3,Float64}, 10^4);

julia> system = InPlaceNeighborList(xpositions=x, cutoff=0.1, unitcell=[1,1,1], parallel=false);

julia> neighborlist!(system)
210034-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 2669, 0.04444346517920411)
 (1, 8475, 0.02554075837438248)
 ⋮
 (9463, 5955, 0.08698158178214915)
 (9463, 2308, 0.09482635540291776)

julia> x .= rand(SVector{3,Float64}, 10^4); # change coordinates

julia> @time neighborlist!(system; parallel=false)
  0.007978 seconds
209418-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 1253, 0.09420839394144173)
 (1, 2048, 0.01448095691254145)
 ⋮
 (9728, 6367, 0.08204145034963985)
 (9728, 2594, 0.06536277710826768)

NeighborPair{N,T,T2}

Structure that holds the information of a pair of particles that are neighbors within the cutoff distance.

Fields accessed by the user:

• i::Int: index of the first particle in the original array of coordinates.
• j::Int: index of the second particle in the original array of coordinates.
• x::SVector{N,T}: coordinates of the first particle (minimum-image adjusted).
• y::SVector{N,T}: coordinates of the second particle (minimum-image adjusted).
• d::T: Euclidean distance between the particles (computed lazily).
• d2::T2: squared Euclidean distance between the particles.

ParticleSystem(
    xpositions::AbstractVector{SVector{N,T}};
    ypositions::Union{AbstractVector{SVector{N,T}}, Nothing} = nothing,
    unitcell::Union{SVector, SMatrix, Nothing} = nothing,
    cutoff::Number,
    output,
    output_name::Symbol = :output,
    parallel::Bool = true,
    nbatches::Tuple{Int,Int} = (0, 0),
    lcell = 1,
    validate_coordinates = _validate_coordinates,
)

Type-stable constructor. Accepts a Vector{SVector{N,T}} for positions and an SVector (orthorhombic), SMatrix (triclinic), or nothing (non-periodic) for the unit cell. All type parameters are resolved at compile time.

ParticleSystem(;
    xpositions::Union{AbstractVector{<:AbstractVector},AbstractMatrix},
    #or
    xpositions::Union{AbstractVector{<:AbstractVector},AbstractMatrix},
    ypositions::Union{AbstractVector{<:AbstractVector},AbstractMatrix},
    # and
    unitcell::Union{Nothing,AbstractVecOrMat} = nothing,
    cutoff::Number,
    output,
    output_name::Symbol = :output,
    parallel::Bool = true,
    nbatches::Tuple{Int,Int} = (0, 0),
    validate_coordinates = _validate_coordinates,
)

Flexible keyword-only constructor. Accepts any supported coordinate type (vectors of vectors, matrices, or SVector arrays) and any AbstractVecOrMat for the unit cell. Converts inputs and delegates to the type-stable constructor.

Constructor of the ParticleSystem type given the positions of the particles.

• Positions can be provided as vectors of 2D or 3D vectors (preferentially static vectors from StaticArrays), or as (2,N) or (3,N) matrices (v0.8.28 is required for matrices).

• If only the xpositions array is provided, a single set of coordinates is considered, and the computation will be mapped for the N(N-1) pairs of this set.

• If the xpositions and ypositions arrays of coordinates are provided, the computation will be mapped to the N×M pairs of particles, being N and M the number of particles of each set of coordinates.

The unit cell (either a vector for Orthorhombic cells or a full unit cell matrix for Triclinic cells - where columns contain the lattice vectors), the cutoff used for the construction of the cell lists and the output variable of the calculations. If unitcell == nothing, the system is considered not-periodic, in which case artificial periodic boundaries will be built such that images are farther from each other than the cutoff.

output_name can be set to a symbol that best identifies the output variable. For instance, if output_name=:forces, the forces can be retrieved from the structure using the system.forces notation.

The parallel and nbatches flags control the parallelization scheme of computations. By default the parallelization is turned on and nbatches is set with heuristics that provide good efficiency in most cases.

After construction, use update!(system; xpositions=..., ypositions=..., cutoff=..., unitcell=..., parallel=...) to update any system properties before subsequent pairwise! calls.

The validate_coordinates function can be used to validate the coordinates before computations, and throw appropriate error messages. By default the validation checks if the coordinates are not missing or NaN. The function must have a single input parameter and be (x) -> nothing to skip any validation.

Example

In these examples, we compute the sum of the squared distances between the particles that are within the cutoff:

Single set of particles

julia> using CellListMap

julia> using PDBTools: read_pdb, coor

julia> positions = coor(read_pdb(CellListMap.argon_pdb_file));

julia> sys = ParticleSystem(
           xpositions = positions,
           unitcell = [21.0, 21.0, 21.0],
           cutoff = 8.0,
           output = 0.0,
        );

julia> pairwise!((pair,output) -> output += pair.d2, sys)
43774.54367600001

Two sets of particles

julia> using CellListMap, PDBTools

julia> xpositions = coor(read_pdb(CellListMap.argon_pdb_file))[1:50];

julia> ypositions = coor(read_pdb(CellListMap.argon_pdb_file))[51:100];

julia> sys = ParticleSystem(
           xpositions = xpositions,
           ypositions = ypositions,
           unitcell = [21.0, 21.0, 21.0],
           cutoff = 8.0,
           output = 0.0,
           parallel = false, # use true for parallelization
        );

julia> pairwise!((pair,output) -> output += pair.d2, sys)
21886.196785000004

ParticleSystemPositions{N,T}

Wrapper around a Vector{SVector{N,T}} that carries particle coordinates and an updated::Ref{Bool} flag. Any mutation through the standard array interface (indexing, push!, resize!, empty!, append!, in-place broadcast) sets the flag, triggering cell-list recomputation on the next pairwise! call.

Users interact with this type only through system.xpositions / system.ypositions, treating it as an ordinary vector. To replace the entire coordinate array or update other system properties, use update!.

AbstractParticleSystem{OutputName}

An abstract type representing a particle system that computes interactions between particles using cell lists. Can be used to control dispatch of functions that operate on particle systems, because all particle systems should be a subtype of this type.

Use the ParticleSystem constructor and interface for anything else than dispatch.

ParticleSystem1

Structure that carries the information necessary for pairwise! computations, for systems with one set of positions (thus, replacing the loops over N(N-1) pairs of particles of the set). Can be used to control dispatch.

Use the ParticleSystem constructor and interface for anything else than dispatch.

ParticleSystem2

Structure that carries the information necessary for pairwise! computations, for systems with two set of positions (thus, replacing the loops over N×M pairs of particles, being N and M the number of particles of each set). Can be used to control dispatch.

Use the ParticleSystem constructor and interface for anything else than dispatch.

pairwise!(
    f::Function, system::AbstractParticleSystem; 
    show_progress=true, reset=true,
)

Function that maps the f function into all pairs of particles of system that are found to be within the cutoff.

The function f receives a NeighborPair struct and the output:

function f(pair, output)
    # pair.i, pair.j: indices of the particles
    # pair.x, pair.y: coordinates (minimum-image adjusted)
    # pair.d: distance between particles
    # pair.d2: squared distance
    # update output
    return output
end

Thread-safety is taken care automatically in parallel executions.

pairwise is an alias to pairwise! for syntax consistency when the output variable is immutable.

If reset is set to false, the value of system.output will not be set to zero(typeof(system.output)) before the new accumulation.

Example

In this example we compute the sum of 1/(1+d) where d is the distance between particles of a set, for d < cutoff.

julia> sys = ParticleSystem(
           xpositions = rand(SVector{3,Float64},1000), 
           unitcell=[1,1,1], 
           cutoff = 0.1, 
           output = 0.0
           );

julia> pairwise!((pair, output) -> output += 1 / (1 + pair.d), sys)
1870.0274887950268

update!(
    sys::AbstractParticleSystem;
    xpositions = nothing,
    #or
    positions = nothing,
    ypositions = nothing,
    cutoff = nothing,
    unitcell = nothing,
    parallel = nothing,
)

Update one or more properties of sys in a single call. Only the keyword arguments that are provided (i.e. not nothing) are updated.

• xpositions (or positions as a shortcut for single-set systems): new coordinates for the first (or only) set of particles. Accepts the same types as the ParticleSystem constructor: a Vector{SVector}, a vector of vectors, or an (D, N) matrix. The internal storage is resized automatically if the number of particles changes.

• ypositions: new coordinates for the second set of particles (only valid for two-set systems). Same accepted types as xpositions.

• cutoff: new cutoff distance.

• unitcell: new unit cell. Must be of the same cell type (orthorhombic or triclinic) as the original system. Manual updating of non-periodic systems is not allowed.

• parallel: whether to use multi-threading (true or false).

Example

julia> using CellListMap, StaticArrays

julia> sys = ParticleSystem(
           xpositions = rand(SVector{3,Float64}, 100),
           unitcell = [1,1,1],
           cutoff = 0.1,
           output = 0.0,
       );

julia> new_x = rand(SVector{3,Float64}, 100);

julia> update!(sys; xpositions=new_x, cutoff=0.2, unitcell=[2,2,2], parallel=false);

resize_output!(sys::AbstractParticleSystem, n::Int)

Resizes the output array and the auxiliary output arrays used for multithreading, if the number of particles of the system changed.

The function will error if Base.resize! is not defined for the type of system.output. In this case, a Base.resize! method must be implemented by the user.

copy_output(x)

Defines how the output variable is copied. Identical to Base.copy(x) and implemented for the types in Union{Number, StaticArraysCore.FieldVector, StaticArraysCore.SVector}.

Other custom output types must have their copy_output method implemented.

Example

using CellListMap
# Custom data type
struct A x::Int end
# Custom output type (array of A)
output = [ A(0) for _ in 1:100 ]
# How to copy an array of `A`
CellListMap.copy_output(v::Vector{A}) = [ x for x in v ]

# Alternatively, in this case, one could have defined:
Base.copy(a::A) = a
CellListMap.copy_output(v::Vector{A}) = copy(v)

The user must guarantee that the copy is independent of the original array. For many custom types it is possible to define

CellListMap.copy_output(v::Vector{T}) where {T<:CustomType} = deepcopy(v)

reset_output(x)
reset_output!(x)

Function that defines how to reset (or zero) the output variable. For Union{Number, StaticArraysCore.FieldVector, StaticArraysCore.SVector} it is implemented as zero(x).

Other custom output types must have their reset_output! method implemented.

The function must return the variable itself. If it is immutable, a new instante of the variable must be created, with the reset value.

reset_output and reset_output! are aliases, and by convention reset_output! is preferred for mutable types.

Example

In this example, we define a reset_output function that will set to +Inf the minimum distance between particles (not always resetting means zeroing).

julia> using CellListMap

julia> struct MinimumDistance d::Float64 end

julia> CellListMap.reset_output(x::MinimumDistance) = MinimumDistance(+Inf)

julia> x = MinimumDistance(1.0)
MinimumDistance(1.0)

julia> CellListMap.reset_output(x)
MinimumDistance(Inf)

See the reducer help entry for a complete example of how to use reset_output.

reducer(x,y)
reducer!(x,y)

Defines how to reduce (combine, or merge) to variables computed in parallel to obtain a single instance of the variable with the reduced result.

reducer and reducer! are aliases, and reducer! is preferred, by convention for mutating functions.

The most common reducer is the sum, and this is how it is implemented for Union{Number, StaticArraysCore.FieldVector, StaticArraysCore.SVector}. For example, when computing energies, or forces, the total energy is the sum of the energies. The force on one particle is the sum of the forces between the particle and every other particle. Thus, the implemented reducer is the sum:

reducer(x,y) = +(x,y)

However, in many cases, reduction must be done differently. For instance, if the minimum distance between particles is to be computed, it is interesting to define a custom type and associated reducer. For example:

struct MinimumDistance d::Float64 end
reducer(x::MinimumDistance, y::MinimumDistance) = MinimumDistance(min(x.d, y.d))

The overloading of reducer allows the use of parallel computations for custom, complex data types, containing different types of variables, fields, or sizes.

The appropriate behavior of the reducer should be carefully inspected by the user to avoid spurious results.

Example

In this example we show how to obtain the minimum distance among argon atoms in a simulation box.

julia> using CellListMap, PDBTools

julia> positions = coor(read_pdb(CellListMap.argon_pdb_file));

julia> struct MinimumDistance d::Float64 end # Custom output type

julia> CellListMap.copy_output(d::MinimumDistance) = MinimumDistance(d.d) # Custom copy function for `Out`

julia> CellListMap.reset_output(d::MinimumDistance) = MinimumDistance(+Inf) # How to reset an array with elements of type `MinimumDistance`

julia> CellListMap.reducer(md1::MinimumDistance, md2::MinimumDistance) = MinimumDistance(min(md1.d, md2.d)) # Custom reduction function

julia> # Construct the system
       sys = ParticleSystem(;
           positions = positions,
           unitcell = [21,21,21],
           cutoff = 8.0,
           output = MinimumDistance(+Inf),
       );

julia> # Obtain the minimum distance between atoms:
       pairwise!((pair,output) -> pair.d < output.d ? MinimumDistance(pair.d) : output, sys)
MinimumDistance(2.1991993997816563)

reduce_output!(output, output_threaded)

Function that defines how to reduce the vector of output_threaded variables (one per batch) into the single output variable, after a parallel computation. The default implementation calls CellListMap.reducer(output, output_threaded[i]) for each batch sequentially.

For custom output types that are large collections (e.g. lists of pairs), overloading this function can avoid significant allocation overhead. The generic implementation performs nbatches incremental resize! calls on output, each copying the existing data to a new buffer. For output types where the final size is known only after computing all batches, this leads to O(N × nbatches) total memory traffic in the reduce step. A custom overload can compute the total size first, resize once, and then copy each batch directly.

Requirements for the overload

• The signature must be reduce_output!(output::MyType, output_threaded::Vector{<:MyType})
• The function must return output.

Example

The following shows a custom overload for a hypothetical PairList output type that holds a growing list of pairs. The overload resizes the output list once and copies each batch with copyto!, avoiding the O(N × nbatches) overhead of the generic implementation:

struct PairList
    n::Int
    list::Vector{Tuple{Int,Int,Float64}}
end

function CellListMap.reduce_output!(output::PairList, output_threaded::Vector{<:PairList})
    ntot = output.n + sum(nb.n for nb in output_threaded; init = 0)
    if length(output.list) < ntot
        resize!(output.list, ntot)
    end
    offset = output.n
    for nb in output_threaded
        output.list[offset + 1:offset + nb.n] .= @view(nb.list[1:nb.n])
        offset += nb.n
    end
    output.n = ntot
    return output
end

wrap_relative_to(x, xref, unit_cell_matrix::AbstractMatrix)

Wraps the coordinates of point x such that it is the minimum image relative to xref.

wrap_relative_to(x,xref,sides::AbstractVector)

Wraps the coordinates of point x such that it is the minimum image relative to xref, for an Orthorhombic cell of which only the sides are provided.

get_computing_box(sys::AbstractParticleSystem)

Retrieves the computing box of the system. The computing box is large enough to contain all coordinates of the particles, plus the cutoff.

get_nbatches(sys::AbstractParticleSystem)

Returns the number of batches for parallel computations of cell list construction and function mapping.