Application interface

neighborlist(x, cutoff; unitcell=nothing, parallel=true, show_progress=false)

Computes the list of pairs of particles in x which are closer to each other than cutoff. If the keyword parameter unitcell is provided (as a vector of sides or a general unit cell matrix, periodic boundary conditions are considered).

Example

Compute the neighborlist between within a set Argon atoms, considering the system non-periodic (do not provide a unitcell):

julia> using CellListMap, PDBTools

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

julia> neighborlist(x, 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(x, 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)

neighborlist(
    x, y, cutoff; 
    unitcell=nothing, 
    parallel=true, 
    show_progress=false, 
    nbatches=(0,0)
)

Computes the list of pairs of particles of x which are closer than r to the particles of y.

Examples

Compute the neighborlist between two sets of Argon atoms, considering the system non-periodic (do not provide a unitcell):

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(x, y, 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)

Now, considering the system periodic:

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(x, y, 8.0; unitcell = [21.0, 21.0, 21.0], parallel=false)
584-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 13, 7.0177634180502215)
 (1, 24, 7.97689645513632)
 (1, 29, 3.177029085967527)
 ⋮
 (18, 10, 6.9654396670725)
 (18, 37, 6.222988130894417)

InPlaceNeighborList(;
    x::AbstractVecOrMat,
    y::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.

• If only x is provided, the neighbor list of the set is computed.
• If x and y are provided, the neighbor list between the 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 unicell is not provide (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(x=x, cutoff=0.1, unitcell=[1,1,1]) 
InPlaceNeighborList with types: 
CellList{3, Float64}
Box{OrthorhombicCell, 3, Float64, Float64, 9}
Current list buffer size: 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):

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

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

julia> @time neighborlist!(system);
  0.008004 seconds (228 allocations: 16.728 MiB)

julia> update!(system, rand(SVector{3,Float64}, 10^4); cutoff = 0.1, unitcell = [1,1,1]);

julia> @time neighborlist!(system);
  0.024811 seconds (167 allocations: 7.887 MiB)

julia> update!(system, rand(SVector{3,Float64}, 10^4); cutoff = 0.1, unitcell = [1,1,1]);

julia> @time neighborlist!(system);
  0.005213 seconds (164 allocations: 1.439 MiB)

julia> update!(system, rand(SVector{3,Float64}, 10^4); cutoff = 0.1, unitcell = [1,1,1]);

julia> @time neighborlist!(system);
  0.005276 seconds (162 allocations: 15.359 KiB)

update!(system::InPlaceNeighborList, x::AbstractVecOrMat; cutoff=nothing, unitcell=nothing)
update!(system::InPlaceNeighborList, x::AbstractVecOrMat, y::AbstractVecOrMat; cutoff=nothing, unitcell=nothing)

Updates a InPlaceNeighborList system, by updating the coordinates, cutoff, and unitcell.

Examples

For self-pairs computations

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

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

julia> neighborlist!(system);

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

julia> update!(system, new_x; cutoff = 0.05)
InPlaceNeighborList with types: 
CellList{3, Float64}
Box{NonPeriodicCell, 3, Float64, Float64, 9}
Current list buffer size: 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(x=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 (see https://m3g.github.io/CellListMap.jl/stable/parallelization/#Number-of-batches)). By default the parallelization is turned on and nbatches is set with heuristics that may provide good efficiency in most cases.

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 internal updated flag. When coordinates are mutated through the supported interface, the flag is set automatically, so that cell lists are recomputed on the next call to pairwise!.

A ParticleSystemPositions can be constructed from:

• A vector of vectors (e.g. Vector{Vector{Float64}}).
• A vector of SVectors (e.g. Vector{SVector{3,Float64}}).
• An (D, M) matrix, where D is the dimension and M the number of particles.

Mutating interface

The following functions mutate the positions and flag the array as updated, triggering recomputation of the cell lists on the next pairwise! call:

Read-only interface

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_cutoff!(system, cutoff)

Function to update the cutoff` of the system.

This function can be used to update the system geometry in iterative schemes.

Example

Here we initialize a particle system with a cutoff of 8.0 and then update the cutoff to 10.0.

julia> using CellListMap, PDBTools

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

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

julia> update_cutoff!(sys, 10.0)
ParticleSystem1{default_output_name} of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 21.0 0.0 0.0; 0.0 21.0 0.0; 0.0 0.0 21.0 ]
      cutoff = 10.0
      number of computing cells on each dimension = [5, 5, 5]
      computing cell sizes = [10.5, 10.5, 10.5] (lcell: 1)
      Total number of cells = 125
    CellListMap.CellList{3, Float64}
      100 real particles.
      8 cells with real particles.
      800 particles in computing box, including images.
    Parallelization auxiliary data set for 8 batch(es).
    Type of output variable (default_output_name): Float64

update_unitcell!(system, unitcell)

Function to update the unit cell of the system. The unicell must be of the same type (OrthorhombicCell, TriclinicCell) of the original system (changing the type of unit cell requires reconstructing the system).

The unitcell can be a N×N matrix or a vector of dimension N, where N is the dimension of the system (2D or 3D).

This function can be used to update the system geometry in iterative schemes, where the size of the simulation box changes during the simulation.

Example

julia> using CellListMap, StaticArrays, PDBTools

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

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

julia> update_unitcell!(sys, [30.0, 30.0, 30.0])
ParticleSystem1{default_output_name} of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 30.0 0.0 0.0; 0.0 30.0 0.0; 0.0 0.0 30.0 ]
      cutoff = 8.0
      number of computing cells on each dimension = [6, 6, 6]
      computing cell sizes = [10.0, 10.0, 10.0] (lcell: 1)
      Total number of cells = 216
    CellListMap.CellList{3, Float64}
      100 real particles.
      8 cells with real particles.
      800 particles in computing box, including images.
    Parallelization auxiliary data set for 8 batch(es).
    Type of output variable (default_output_name): Float64

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.