Help entries

struct AuxThreaded{N, T}

Internal function or structure - interface may change.

Extended help

• particles_per_batch::Int64

• idxs::Vector{UnitRange{Int64}}: Default: Vector{UnitRange{Int}}(undef, 0)

• lists::Array{CellList{N, T}, 1} where {N, T}: Default: Vector{CellList{N, T}}(undef, 0)

Auxiliary structure to carry threaded lists and ranges of particles to be considered by each thread on parallel construction.

AuxThreaded(cl::CellListPair{N,T}) where {N,T}

Constructor for the AuxThreaded type for lists of disjoint particle sets, to be passed to UpdateCellList! for in-place update of cell lists.

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(3) for i in 1:50_000 ];

julia> y = [ 250*rand(3) for i in 1:10_000 ];

julia> cl = CellList(x,y,box);

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> cl = UpdateCellList!(x,box,cl,aux)
CellList{3, Float64}
  100000 real particles.
  31190 cells with real particles.
  1134378 particles in computing box, including images.

AuxThreaded(cl::CellList{N,T}) where {N,T}

Constructor for the AuxThreaded type, to be passed to UpdateCellList! for in-place update of cell lists.

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(3) for _ in 1:100_000 ];

julia> cl = CellList(x,box);

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> cl = UpdateCellList!(x,box,cl,aux)
CellList{3, Float64}
  100000 real particles.
  31190 cells with real particles.
  1134378 particles in computing box, including images.

struct Box{UnitCellType, N, T, TSQ, M, TR}

Internal function or structure - interface may change.

Extended help

• input_unit_cell::CellListMap.UnitCell

• aligned_unit_cell::CellListMap.UnitCell

• rotation::StaticArraysCore.SMatrix{N, N, TR, M} where {N, M, TR}

• inv_rotation::StaticArraysCore.SMatrix{N, N, TR, M} where {N, M, TR}

• lcell::Int64

• nc::StaticArraysCore.SVector{N, Int64} where N

• cutoff::Any

• cutoff_sqr::Any

• computing_box::Tuple{StaticArraysCore.SVector{N, T}, StaticArraysCore.SVector{N, T}} where {N, T}

• cell_size::StaticArraysCore.SVector

Structure that contains some data required to compute the linked cells. To be initialized with the box size and cutoff.

Examples

julia> using CellListMap

julia> sides = [250,250,250];

julia> cutoff = 10;

julia> box = Box(sides,cutoff)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [250.0 0.0 0.0; 0.0 250.0 0.0; 0.0 0.0 250.0]
  cutoff: 10.0
  number of computing cells on each dimension: [27, 27, 27]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 19683

julia> box = Box([ 10  0  0 
                    0 10  5
                    0  0 10 ], 1)
Box{TriclinicCell, 3, Float64, 9}
unit cell matrix: [10.0 0.0 0.0; 0.0 10.0 5.0; 0.0 0.0 10.0]
cutoff: 1.0
number of computing cells on each dimension: [12, 17, 12]
computing cell sizes: [1.0, 1.0, 1.0] (lcell: 1)
Total number of cells: 2448

Box(unitcell::Limits, cutoff; lcell::Int=1)

This constructor receives the output of limits(x) or limits(x,y) where x and y are the coordinates of the particles involved, and constructs a Box with size larger than the maximum coordinates ranges of all particles plus twice the cutoff. This is used to emulate pairwise interactions in non-periodic boxes. The output box is an NonPeriodicCell box type, which internally is treated as Orthorhombic with boundaries that guarantee that particles do not see images of each other.

Examples

julia> x = [ [100,100,100] .* rand(3) for i in 1:100_000 ];

julia> box = Box(limits(x),10)
Box{NonPeriodicCell, 3}
  unit cell matrix = [ 110.0, 0.0, 0.0; 0.0, 110.0, 0.0; 0.0, 0.0, 110.0 ]
  cutoff = 10.0
  number of computing cells on each dimension = [12, 12, 12]
  computing cell sizes = [11.0, 11.0, 11.0] (lcell: 1)
  Total number of cells = 1728

julia> y = [ [150,150,50] .* rand(3) for i in 1:100_000 ];

julia> box = Box(limits(x,y),10)
Box{NonPeriodicCell, 3}
  unit cell matrix = [ 160.0, 0.0, 0.0; 0.0, 160.0, 0.0; 0.0, 0.0, 110.0 ]
  cutoff = 10.0
  number of computing cells on each dimension = [17, 17, 12]
  computing cell sizes = [10.67, 10.67, 11.0] (lcell: 1)
  Total number of cells = 3468

Box(unit_cell_matrix::AbstractMatrix, cutoff, lcell::Int=1, UnitCellType=TriclinicCell)

Construct box structure given the cell matrix of lattice vectors. This constructor will always return a TriclinicCell box type, unless the UnitCellType parameter is set manually to OrthorhombicCell

Example

julia> unit_cell = [ 100   50    0 
                       0  120    0
                       0    0  130 ];

julia> box = Box(unit_cell,10)
Box{TriclinicCell, 3, Float64, 9}
  unit cell matrix: [100.0 50.0 0.0; 0.0 120.0 0.0; 0.0 0.0 130.0]
  cutoff: 10.0
  number of computing cells on each dimension: [17, 14, 15]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 3570

Box(sides::AbstractVector, cutoff, lcell::Int=1, UnitCellType=OrthorhombicCell)

For orthorhombic unit cells, Box can be initialized with a vector of the length of each side.

Example

julia> box = Box([120,150,100],10)
Box{OrthorhombicCell, 3, Float64, 9}
  unit cell matrix: [120.0 0.0 0.0; 0.0 150.0 0.0; 0.0 0.0 100.0]
  cutoff: 10.0
  number of computing cells on each dimension: [14, 17, 12]
  computing cell sizes: [10.0, 10.0, 10.0] (lcell: 1)
  Total number of cells: 2856

struct Cell{N, T}

Internal function or structure - interface may change.

Extended help

• linear_index::Int64

• cartesian_index::CartesianIndex

• center::StaticArraysCore.SVector

• contains_real::Bool

• n_particles::Int64

• particles::Array{CellListMap.ParticleWithIndex{N, T}, 1} where {N, T}

This structure contains the cell linear index and the information about if this cell is in the border of the box (such that its neighboring cells need to be wrapped)

struct CellList{N, T}

Internal function or structure - interface may change.

Extended help

• n_real_particles::Int64: Number of real particles.

• number_of_cells::Int64: Number of cells.

• n_particles::Int64: mutable number of particles in the computing box.

• n_cells_with_real_particles::Int64: mutable number of cells with real particles.

• n_cells_with_particles::Int64: mutable number of cells with particles, real or images.

• cell_indices::Vector{Int64}: Auxiliary array that contains the indexes in list of the cells with particles, real or images.

• cell_indices_real::Vector{Int64}: Auxiliary array that contains the indexes in the cells with real particles.

• cells::Array{CellListMap.Cell{N, T}, 1} where {N, T}: Vector containing cell lists of cells with particles.

• nbatches::CellListMap.NumberOfBatches: Number of batches for the parallel calculations.

• projected_particles::Array{Array{CellListMap.ProjectedParticle{N, T}, 1}, 1} where {N, T}: Auxiliar array to store projected particles.

Structure that contains the cell lists information.

CellList(x::AbstractMatrix, y::AbstractMatrix, box::Box{UnitCellType,N,T}; kargs...) where {UnitCellType,N,T}

Reinterprets the matrices x and y as vectors of static vectors and calls the equivalent function with the reinterprted input. The first dimension of the matrices must be the dimension of the points (2 or 3).

CellList(x::AbstractMatrix, box::Box{UnitCellType,N,T}; kargs...) where {UnitCellType,N,T}

Reinterprets the matrix x as vectors of static vectors and calls the equivalent function with the reinterprted input. The first dimension of the matrix must be the dimension of the points (2 or 3).

CellList(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box{UnitCellType,N,T};
    parallel::Bool=true,
    nbatches::Tuple{Int,Int}=(0,0),
    autoswap::Bool=true
) where {UnitCellType,N,T}

Function that will initialize a CellListPair structure from scratch, given two vectors of particle coordinates and a Box, which contain the size of the system, cutoff, etc. By default, the cell list will be constructed for smallest vector, but this is not always the optimal choice. Using autoswap=false the cell list is constructed for the second (y)

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> y = [ 250*rand(SVector{3,Float64}) for i in 1:10000 ];

julia> cl = CellList(x,y,box)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   10000 particles in the reference vector.
   961 cells with real particles of target vector.

julia> cl = CellList(x,y,box,autoswap=false)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   1000 particles in the reference vector.
   7389 cells with real particles of target vector.

CellList(
    x::AbstractVector{AbstractVector},
    box::Box{UnitCellType,N,T};
    parallel::Bool=true,
    nbatches::Tuple{Int,Int}=(0,0)
) where {UnitCellType,N,T}

Function that will initialize a CellList structure from scratch, given a vector or particle coordinates (a vector of vectors, typically of static vectors) and a Box, which contain the size ofthe system, cutoff, etc. Except for small systems, the number of parallel batches is equal to the number of threads, but it can be tunned for optimal performance in some cases.

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> cl = CellList(x,box)
CellList{3, Float64}
  100000 real particles.
  15600 cells with real particles.
  126276 particles in computing box, including images.

struct CellListPair{V, N, T, Swap}

Internal function or structure - interface may change.

Extended help

• ref::Any

• target::CellList

Structure that will cointain the cell lists of two independent sets of particles for cross-computation of interactions

mutable struct InPlaceNeighborList{B, C, A, NB<:CellListMap.NeighborList}

Internal function or structure - interface may change.

Structure that containst the system information for neighborlist computations. All fields are internal.

Extended help

• box::Any

• cl::Any

• aux::Any

• nb::CellListMap.NeighborList

• nb_threaded::Vector{NB} where NB<:CellListMap.NeighborList

• parallel::Bool

• show_progress::Bool

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)

struct Limits{N, T}

Internal function or structure - interface may change.

Extended help

• limits::StaticArraysCore.SVector

Structure that contains the maximum lengths on each direction, to dispatch on the construction of boxes without periodic boundary conditions.

struct NumberOfBatches

Internal function or structure - interface may change.

Extended help

• build_cell_lists::Int64

• map_computation::Int64

Structure to define the number of batches used in the parallel splitting of the calculations of the cell list construction and of the map_pairwise computation. It is initialized with a standard heuristic that returns at most the number of threads, but may return a smaller number if the system is small. The two parameters can be tunned for optimal performance of each step of the calculation (cell list construction and mapping of interactions). The construction of the cell lists require a larger number of particles for threading to be effective, Thus by default the system size that allows multi-threading is greater for this part of the calculation.

struct ParticleWithIndex{N, T}

Internal function or structure - interface may change.

Extended help

• index::Int64

• real::Bool

• coordinates::StaticArraysCore.SVector

Copies particle coordinates and associated index, to build contiguous particle lists in memory when building the cell lists. This strategy duplicates the particle coordinates data, but is probably worth the effort.

struct ProjectedParticle{N, T}

Internal function or structure - interface may change.

Extended help

• index::Int64

• xproj::Any

• coordinates::StaticArraysCore.SVector

Auxiliary structure to contain projected particles. Types of scalars are chosen such that with a SVector{3,Float64} the complete struct has 32bytes.

Internal function or structure - interface may change.

Structures to control dispatch on swapped vs. not swapped cell list pairs.

UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box,
    cl:CellListPair,
    parallel=true
)

Function that will update a previously allocated CellListPair structure, given new updated particle positions, for example. This method will allocate new aux threaded auxiliary arrays. For a non-allocating version, see the UpdateCellList!(x,y,box,cl,aux) method.

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> y = [ 250*rand(SVector{3,Float64}) for i in 1:10000 ];

julia> cl = CellList(x,y,box);

julia> cl = UpdateCellList!(x,y,box,cl); # update lists

UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,
    cl:CellList,
    parallel=true
)

Function that will update a previously allocated CellList structure, given new updated particle positions. This function will allocate new threaded auxiliary arrays in parallel calculations. To preallocate these auxiliary arrays, use the UpdateCellList!(x,box,cl,aux) method instead.

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> cl = CellList(x,box);

julia> box = Box([260,260,260],10);

julia> x = [ 260*rand(SVector{3,Float64}) for i in 1:1000 ];

julia> cl = UpdateCellList!(x,box,cl); # update lists

UpdateCellList!(
    x::AbstractMatrix,
    y::AbstractMatrix,
    box::Box,
    cl_pair::CellListPair,
    aux::Union{Nothing,AuxThreaded};
    parallel::Bool=true
) where {UnitCellType,N}

Reinterprets the matrices x and y as vectors of static vectors and calls the equivalent function with the reinterprted input. The first dimension of the matrices must be the dimension of the points (2 or 3).

UpdateCellList!(
    x::AbstractMatrix,
    y::AbstractMatrix,
    box::Box{UnitCellType,N},
    cl_pair::CellListPair;
    parallel::Bool=true
) where {UnitCellType,N}

Reinterprets the matrices x and y as vectors of static vectors and calls the equivalent function with the reinterprted input. The first dimension of the matrices must be the dimension of the points (2 or 3).

UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    y::AbstractVector{<:AbstractVector},
    box::Box,
    cl_pair::CellListPair,
    aux::Union{Nothing,AuxThreaded};
    parallel::Bool=true
)

This function will update the cl_pair structure that contains the cell lists for disjoint sets of particles. It receives the preallocated aux structure to avoid reallocating auxiliary arrays necessary for the threaded construct of the lists.

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(3) for i in 1:50_000 ];

julia> y = [ 250*rand(3) for i in 1:10_000 ];

julia> cl = CellList(x,y,box)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   50000 particles in the reference vector.
   7381 cells with real particles of target vector.

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> x = [ 250*rand(3) for i in 1:50_000 ];

julia> y = [ 250*rand(3) for i in 1:10_000 ];

julia> cl = UpdateCellList!(x,y,box,cl,aux)
CellList{3, Float64}
  10000 real particles.
  7358 cells with real particles.
  12591 particles in computing box, including images.

To illustrate the expected ammount of allocations, which are a consequence of thread spawning only:

julia> using BenchmarkTools

julia> @btime UpdateCellList!($x,$y,$box,$cl,$aux)
  715.661 μs (41 allocations: 3.88 KiB)
CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
   50000 particles in the reference vector.
   7414 cells with real particles of target vector.
   
julia> @btime UpdateCellList!($x,$y,$box,$cl,$aux,parallel=false)
   13.042 ms (0 allocations: 0 bytes)
 CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64}
    50000 particles in the reference vector.
    15031 cells with real particles of target vector.
 

UpdateCellList!(
    x::AbstractMatrix,
    box::Box,
    cl::CellList{N,T},
    aux::Union{Nothing,AuxThreaded{N,T}};
    parallel::Bool=true
) where {N,T}

Reinterprets the matrix x as vectors of static vectors and calls the equivalent function with the reinterprted input. The first dimension of the matrix must be the dimension of the points (2 or 3).

function UpdateCellList!(
    x::AbstractMatrix,
    box::Box,
    cl::CellList{N,T};
    parallel::Bool=true
) where {N,T}

Reinterprets the matrix x as vectors of static vectors and calls the equivalent function with the reinterprted input. The first dimension of the matrix must be the dimension of the points (2 or 3).

UpdateCellList!(
    x::AbstractVector{<:AbstractVector},
    box::Box,
    cl::CellList{N,T},
    aux::Union{Nothing,AuxThreaded{N,T}};
    parallel::Bool=true
) where {N,T}

Function that updates the cell list cl new coordinates x and possibly a new box box, and receives a preallocated aux structure of auxiliary vectors for threaded cell list construction. Given a preallocated aux vector, allocations in this function should be minimal, only associated with the spawning threads, or to expansion of the cell lists if the number of cells or number of particles increased.

Example

julia> box = Box([250,250,250],10);

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> cl = CellList(x,box);

julia> aux = CellListMap.AuxThreaded(cl)
CellListMap.AuxThreaded{3, Float64}
 Auxiliary arrays for nthreads = 8

julia> x = [ 250*rand(SVector{3,Float64}) for i in 1:100000 ];

julia> UpdateCellList!(x,box,cl,aux)
CellList{3, Float64}
  100000 real particles.
  15599 cells with real particles.
  125699 particles in computing box, including images.

To illustrate the expected ammount of allocations, which are a consequence of thread spawning only:

julia> using BenchmarkTools

julia> @btime UpdateCellList!($x,$box,$cl,$aux)
  16.384 ms (41 allocations: 3.88 KiB)
CellList{3, Float64}
  100000 real particles.
  15599 cells with real particles.
  125699 particles in computing box, including images.

julia> @btime UpdateCellList!($x,$box,$cl,$aux,parallel=false)
  20.882 ms (0 allocations: 0 bytes)
CellList{3, Float64}
  100000 real particles.
  15603 cells with real particles.
  125896 particles in computing box, including images.

_promote_types(cell,cutoff)

Internal function or structure - interface may change.

Extended help

Promotes the types of the unit cell matrix (or sides) and cutoff to floats if one or both were input as integers.

add_particle_to_celllist!(
    ip,
    x::SVector{N,T},
    box,
    cl::CellList{N,T};
    real_particle::Bool=true
) where {N,T}

Internal function or structure - interface may change.

Extended help

Adds one particle to the cell lists, updating all necessary arrays.

add_particles!(x,box,ishift,cl::CellList{N,T}) where {N,T}

Internal function or structure - interface may change.

Extended help

Add all particles in vector x to the cell list cl. ishift is the shift in particle index, meaning that particle i of vector x corresponds to the particle with original index i+ishift. The shift is used to construct cell lists from fractions of the original set of particles in parallel list construction.

align_cell(m::StaticMatrix)
align_cell!(m::AbstractMatrix)

Internal function or structure - interface may change.

Extended help

These functions rotate the unit cell matrix such that the largest lattice vector is oriented along the x-axis and, for 3D cells, also that the the plane formed by the largest and second largest lattice vectors is oriented perpendicular to the z-axis.

append_particles!(cell1::Cell,cell2::Cell)

Internal function or structure - interface may change.

Extended help

Add the particles of cell2 to cell1, updating the cell data and, if necessary, resizing (increasing) the particles array of cell1

cell_cartesian_indices(nc::SVector{N,Int}, i1D) where {N}

Internal function or structure - interface may change.

Extended help

Given the linear index of the cell in the cell list, returns the cartesian indices of the cell (for arbitrary dimension N).

cell_center(c::CartesianIndex{N},box::Box{UnitCellType,N,T}) where {UnitCellType,N,T}

Internal function or structure - interface may change.

Extended help

Computes the geometric center of a computing cell, to be used in the projection of points. Returns a SVector{N,T}

cell_limits(m::AbstractMatrix)

Internal function or structure - interface may change.

For 2D and 3D matrices, returns the maximum and minimum coordinates of all vertices.

cell_linear_index(nc::SVector{N,Int}, indices) where N

Internal function or structure - interface may change.

Extended help

Returns the index of the cell, in the 1D representation, from its cartesian coordinates.

cell_matrix_from_sides(sides::AbstractVector)

Internal function or structure - interface may change.

Extended help

Function that returns the Orthorhombic unit cell matrix given a sides vector. This function is type-unstable if the input is not static.

Example

julia> CellListMap.cell_matrix_from_sides([1,1,1])
3×3 SMatrix{3, 3, Int64, 9} with indices SOneTo(3)×SOneTo(3):
 1  0  0
 0  1  0
 0  0  1

cell_vertices(m::AbstractMatrix)

Internal function or structure - interface may change.

Function that returns the vertices of a unit cell in 2D or 3D, given the unit cell matrix.

check_unit_cell(box::Box)

Internal function or structure - interface may change.

Extended help

Checks if the unit cell satisfies the conditions for using the minimum-image convention.

copydata!(cell1::Cell,cell2::Cell)

Internal function or structure - interface may change.

Extended help

Copies the data from cell2 to cell1, meaning that particles are copied element-wise from cell2 to cell1, with the particles array of cell1 being resized (increased) if necessary.

current_and_neighbor_cells(box::Box{UnitCellType,N}) where {UnitCellType,N}

Internal function or structure - interface may change.

Extended help

Returns an iterator over all neighbor cells, including the center one.

dot(x::AbstractVector{T1},y::AbstractVector{T2}) where {T1,T2}

Internal function or structure - interface may change.

Extended help

LinearAlgebra.dot is slower than this for standard arrays (likely more accurate, but that is not relevant here).

draw_cell(m::AbstractMatrix; aspect_ratio=:auto)

Internal function or structure - interface may change.

Draw the unit cell in a 2D or 3D plot. Requires using Plots.

draw_cell_vertices(m::AbstractMatrix)

Internal function or structure - interface may change.

Function that returns the vertices of a unit cell matrix in 2D or 3D, as a vector of static vectors, in a proper order for ploting the cell (the first vertex, in the origin, is repeated at the end of the list, to close the figure)

draw_computing_cell(x,box::Box{UnitCellType,2}) where UnitCellType
draw_computing_cell(cl::CellList,box::Box{UnitCellType,2},x) where UnitCellType

Internal function or structure - interface may change.

This function creates a plot of the computing cell, in two dimensions.

draw_computing_cell(x,box::Box{UnitCellType,3}) where UnitCellType

Internal function or structure - interface may change.

This function creates a plot of the computing cell, in three dimensions.

fastmod1(x)

Internal function or structure - interface may change.

Computes mod(x,1), quickly, using x - floor(x). Maybe irrelevant.

in_computing_box(x::SVector{N},box::Box) where N

Internal function or structure - interface may change.

Extended help

Function that evaluates if a particle is inside the computing bounding box, defined by the maximum and minimum unit aligned cell coordinates.

limits(x,y)

Returns the lengths of a orthorhombic box that encompasses all the particles defined in x and y, to used to set a box without effective periodic boundary conditions.

limits(x)

Returns the lengths of a orthorhombic box that encompasses all the particles defined in x, to be used to set a box without effective periodic boundary conditions.

map_naive!(f::Function, output, x::AbstractVector, box::Box)
map_naive!(f::Function, output, x::AbstractVector, y::AbstractVector, box::Box)

Internal function or structure - interface may change.

Extended help

Function that uses the naive pairwise mapping algorithm, for testing.

map_pairwise(args...;kargs...) = map_pairwise!(args...;kargs...)

is an alias for map_pairwise! which is defined for two reasons: first, if the output of the funciton is immutable, it may be clearer to call this version, from a coding perspective. Second, the python interface through juliacall does not accept the bang as a valid character.

map_pairwise!(
    f::Function,
    output,
    box::Box,
    cl::CellList
    ;parallel::Bool=true,
    show_progress::Bool=false
)

This function will run over every pair of particles which are closer than box.cutoff and compute the Euclidean distance between the particles, considering the periodic boundary conditions given in the Box structure. If the distance is smaller than the cutoff, a function f of the coordinates of the two particles will be computed.

The function f receives six arguments as input:

f(x,y,i,j,d2,output)

Which are the coordinates of one particle, the coordinates of the second particle, the index of the first particle, the index of the second particle, the squared distance between them, and the output variable. It has also to return the same output variable. Thus, f may or not mutate output, but in either case it must return it. With that, it is possible to compute an average property of the distance of the particles or, for example, build a histogram. The squared distance d2 is computed internally for comparison with the cutoff, and is passed to the f because many times it is used for the desired computation.

Example

Computing the mean absolute difference in x position between random particles, remembering the number of pairs of n particles is n(n-1)/2. The function does not use the indices or the distance, such that we remove them from the parameters by using a closure.

julia> n = 100_000;

julia> box = Box([250,250,250],10);

julia> x = [ SVector{3,Float64}(sides .* rand(3)) for i in 1:n ];

julia> cl = CellList(x,box);

julia> f(x,y,sum_dx) = sum_dx + abs(x[1] - y[1])

julia> normalization = N / (N*(N-1)/2) # (number of particles) / (number of pairs)

julia> avg_dx = normalization * map_parwise!((x,y,i,j,d2,sum_dx) -> f(x,y,sum_dx), 0.0, box, cl)

map_pairwise!(f::Function,output,box::Box,cl::CellListPair)

The same but to evaluate some function between pairs of the particles of the vectors.

merge_cell_lists!(cl::CellList,aux::CellList)

Internal function or structure - interface may change.

Extended help

Merges an auxiliary aux cell list to cl, and returns the modified cl. Used to merge cell lists computed in parallel threads.

nbatches(cl)

Returns the number of batches for parallel processing that will be used in the pairwise function mappings associated to cell list cl. It returns the cl.nbatches.map_computation value. This function is important because it must be used to set the number of copies of custom preallocated output arrays.

A second argument can be provided, which may be :map or :build, in which case the function returns either the number of batches used for pairwise mapping or for the construction of the cell lists. Since this second value is internal and does not affect the interface, it can be usually ignored.

Example

julia> x = rand(3,1000); box = Box([1,1,1],0.1);

julia> cl = CellList(x,box,nbatches=(2,16));

julia> nbatches(cl)
16

julia> nbatches(cl,:map)
16

julia> nbatches(cl,:build)
2

neighbor_cells(box::Box{UnitCellType,N}) where {UnitCellType,N}

Internal function or structure - interface may change.

Extended help

Function that returns the iterator of the cartesian indices of all neighboring cells of a cell where the computing cell index is box.lcell.

neighbor_cells_forward(box::Box{UnitCellType,N}) where UnitCellType

Internal function or structure - interface may change.

Extended help

Function that returns the iterator of the cartesian indices of the cells that must be evaluated (forward, i. e. to avoid repeated interactions) if the cells have sides of length box.cell_size. N can be 2 or 3, for two- or three-dimensional systems.

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

Computes the list of pairs of particles of x which are closer than r to the particles of y. The autoswap option will swap x and y to try to optimize the cost of the construction of the cell list.

Example

julia> x = [ rand(3) for i in 1:10_000 ];

julia> y = [ rand(3) for i in 1:1_000 ];

julia> CellListMap.neighborlist(x,y,0.05)
5006-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 269, 0.04770884036497686)
 (25, 892, 0.03850515231540869)
 ⋮
 (9952, 749, 0.048875643578313456)
 (9984, 620, 0.04101242499363183)

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

julia> using CellListMap

julia> x = [ rand(3) for i in 1:10_000 ];

julia> neighborlist(x,0.05)
24848-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 1055, 0.022977369806392412)
 (1, 5086, 0.026650609138167428)
 ⋮
 (9989, 3379, 0.0467653507446483)
 (9989, 5935, 0.02432728985151653)

norm(v::AbstractVector{T}) where T

Internal function or structure - interface may change.

Extended help

norm_sqr from LinearAlgebra is not documented and is slower than this for standard arrays. Thus we define our own norm(x) = norm_sqr(x).

norm_sqr(v::AbstractVector{T}) where T

Internal function or structure - interface may change.

Extended help

norm_sqr from LinearAlgebra is not documented and is slower than this for standard arrays.

normalize(v::AbstractVector{T}) where T

Internal function or structure - interface may change.

Extended help

norm_sqr from LinearAlgebra is not documented and is slower than this for standard arrays. Thus we define our own normalize(v)

particle_cell(x::SVector{N,T}, box::Box) where {N,T}

Internal function or structure - interface may change.

Extended help

Returns the coordinates of the computing cell to which a particle belongs, given its coordinates and the cell_size vector. The computing box is always Orthorhombic, and the first computing box with positive coordinates has indexes Box.lcell + 1.

particles_per_cell(cl)

Internal function or structure - interface may change.

Returns the average number of real particles per computing cell.

partition!(by, x::AbstractVector)

Internal function or structure - interface may change.

Extended help

Function that reorders x vector by putting in the first positions the elements with values satisfying by(el). Returns the number of elements that satisfy the condition.

pathological_coordinates(N)

Internal function or structure - interface may change.

Extended help

Function to generate some coordinates with pathological properties, for testing. Returns x, y, sides and cutoff.

project_particles!(projected_particles,cellⱼ,cellᵢ,Δc,Δc_norm,box)

Internal function or structure - interface may change.

Extended help

Projects all particles of the cell cellⱼ into unnitary vector Δc with direction connecting the centers of cellⱼ and cellᵢ. Modifies projected_particles, and returns a view of `projected particles, where only the particles for which the projection on the direction of the cell centers still allows the particle to be within the cutoff distance of any point of the other cell.

reduce(output, output_threaded)

Most common reduction function, which sums the elements of the output. Here, output_threaded is a vector containing nbatches(cl) copies of the output variable (a scalar or an array). Custom reduction functions must replace this one if the reduction operation is not a simple sum. The output_threaded array is, by default, created automatically by copying the given output variable nbatches(cl) times.

Examples

Scalar reduction:

julia> output = 0.; output_threaded = [ 1, 2 ];

julia> CellListMap.reduce(output,output_threaded)
3

Array reduction:

julia> output = [0,0]; output_threaded = [ [1,1], [2,2] ];

julia> CellListMap.reduce(output,output_threaded)
2-element Vector{Int64}:
 3
 3

julia> output
2-element Vector{Int64}:
 3
 3

replicate_particle!(ip,p::SVector{N},box,cl) where N

Internal function or structure - interface may change.

Extended help

Replicates the particle as many times as necessary to fill the computing box.

replicate_system!(
    x::AbstractVector{SVector{N,T}},
    unit_cell_matrix::AbstractMatrix,
    ranges::Tuple
) where {N,T}

Internal function or structure - interface may change.

Extended help

Replicate the system (modifying the original array of coordinates) in all directions defined by the periodic system and by the range of unitary cells of interest. x can be a (N,M) matrix, and the unit cell matrix can be provided instead of the box.

Example

julia> x = rand(SVector{2,Float64},100);

julia> box = Box([1,1],0.1);

julia> CellListMap.replicate_system!(x,box,(0:0,-1:1))
300-element Vector{SVector{2, Float64}}:
 [0.7119987163255118, 0.6788616154460262]
 [0.6188407316804118, 0.8497116428720384]
 [0.21328895963244354, 0.48932085643862977]
 ⋮
 [0.4114499470191678, 1.1034376619603892]
 [0.6094126258851252, 1.2328989485215263]

reset!(cl::CellList{N,T},box,n_real_particles) where{N,T}

Internal function or structure - interface may change.

Extended help

Resets a cell list, by setting everything to zero, but retaining the allocated particles and projected_particles vectors.

set_idxs!(idxs, n_particles, nbatches)

Internal function or structure - interface may change.

Extended help

Sets the indexes of the particles that will be considered for each batch in parallel runs. Modifies the idxs array of ranges, which is usually the aux.idxs array of the the corresponding AuxThreaded structure.

set_number_of_batches!(cl,nbatches::Tuple{Int,Int}=(0,0);parallel=true)

Internal function or structure - interface may change.

Extended help

Functions that set the default number of batches for the construction of the cell lists, and mapping computations. This is of course heuristic, and may not be the best choice for every problem. See the parameter nbatches of the construction of the cell lists for tunning this.

translation_image(x::AbstractVector{<:AbstractVector},unit_cell_matrix,indices)

Translates a complete set of coordinates given a set of indexes of unit-cells. Returns a new set of coordinates.

Example

julia> x = rand(SVector{2,Float64},100);

julia> box = Box([1,1],0.1);

julia> CellListMap.translation_image(x,box.unit_cell.matrix,(1,1))
100-element Vector{SVector{2, Float64}}:
 [1.847791110439223, 1.5989103939725295]
 [1.3493293666090889, 1.4002971843576644]
 [1.4111736701313218, 1.3471780214994182]
 ⋮
 [1.1548437388991908, 1.7034501001177493]
 [1.4066300885242247, 1.2907398318754952]

translation_image(x::SVector{N,T},unit_cell_matrix,indices) where {N,T}

Internal function or structure - interface may change.

Extended help

Translate vector x according to the unit_cell_matrix lattice vectors and the indices provided.

unitcelltype(::Box{T}) where T = T

Returns the type of a unitcell from the Box structure.

Example

julia> box = Box([1,1,1], 0.1)

julia> unitcelltype(box)
OrthorhombicCell

julia> box = Box([1 0 0; 0 1 0; 0 0 1], 0.1)

julia> unitcelltype(box)
TriclinicCell

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)

update_box(
    box::Box{UnitCellType,N,T,TSQ,M};
    unitcell::Union{Nothing,AbstractVector{T},AbstractMatrix{T},Limits,Tuple}=nothing,
    cutoff::Union{Nothing,T}=nothing,
    lcell::Union{Nothing,Int}=nothing
)

Internal function or structure - interface may change.

Function that returns an updated system box in a type-stable manner, given possible variations in the unitcell, cutoff, or lcell parameters.

view_celllist_particles(cl::CellList)

Internal function or structure - interface may change.

Extended help

Auxiliary function to view the particles of a computing box, including images created for computing purposes.

wrap_cell_fraction(x,unit_cell_matrix)

Internal function or structure - interface may change.

Extended help

Obtaint the coordinates of x as a fraction of unit cell vectors, first positive cell. x is a vector of dimension N and cell a matrix of dimension NxN

Example

julia> unit_cell_matrix = [ 10 0
                            0 10 ];

julia> x = [ 15, 13 ];

julia> wrap_cell_fraction(x,unit_cell_matrix)
2-element Vector{Float64}:
 0.5
 0.3

wrap_relative_to(x,xref,sides::AbstractVector)

Internal function or structure - interface may change.

Extended help

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

wrap_relative_to(x, xref, unit_cell_matrix::SMatrix{N,N,T}) where {N,T}

Internal function or structure - interface may change.

Extended help

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

wrap_to_first(x,unit_cell_matrix)

Internal function or structure - interface may change.

Extended help

Wraps the coordinates of point x such that the returning coordinates are in the first unit cell with all-positive coordinates.

Example

julia> unit_cell_matrix = [ 10 0
                            0 10 ];

julia> x = [ 15, 13 ];

julia> wrap_to_first(x,unit_cell_matrix)
2-element Vector{Float64}:
 5.0
 3.0000000000000004

mutable struct PeriodicSystem1{OutputName, V, O, B, C, A} <: CellListMap.PeriodicSystems.AbstractPeriodicSystem{OutputName}

• xpositions::Any

• output::Any

• _box::Any

• _cell_list::Any

• _output_threaded::Vector

• _aux::Any

• parallel::Bool

Structure that carries the information necessary for map_pairwise! computations, for systems with one set of positions (thus, replacing the loops over N(N-1) pairs of particles of the set).

The xpositions, output, and parallel fields are considered part of the API, and you can retrive or mutate xpositions, retrieve the output or its elements, and set the computation to use or not parallelization by directly accessing these elements.

The other fileds of the structure (starting with _) are internal and must not be modified or accessed directly. The construction of the PeriodicSystem1 structure is done through the PeriodicSystem(;xpositions, unitcell, cutoff, output) auxiliary function.

mutable struct PeriodicSystem2{OutputName, V, O, B, C, A} <: CellListMap.PeriodicSystems.AbstractPeriodicSystem{OutputName}

• xpositions::Any

• ypositions::Any

• output::Any

• _box::Any

• _cell_list::Any

• _output_threaded::Vector

• _aux::Any

• parallel::Bool

Structure that carries the information necessary for map_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).

The xpositions, ypositions, output, and parallel fields are considered part of the API, and you can retrive or mutate positions, retrieve the output or its elements, and set the computation to use or not parallelization by directly accessing these elements.

The other fileds of the structure (starting with _) are internal and must not be modified or accessed directly. The construction of the PeriodicSystem1 structure is done through the PeriodicSystem(;xpositions, ypositions, unitcell, cutoff, output) auxiliary function.

PeriodicSystem( 
    xpositions::AbstractVector{<:AbstractVector},
    #or
    xpositions::AbstractVector{<:AbstractVector},
    ypositions::AbstractVector{<:AbstractVector},
    # and
    unitcell::AbstractVecOrMat,
    cutoff::Number,
    output::Any;
    output_name::Symbol,
    parallel::Bool=true,
    nbatches::Tuple{Int,Int}=(0, 0),
    autoswap::Bool = true
)

Function that sets up the PeriodicSystem type given the positions of the particles.

• Positions can be provided as vectors of 2D or 3D vectors (preferentially static vectors from StaticArrays).

• 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), the cutoff used for the construction of the cell lists and the output variable of the calculations.

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. autoswap = false will guarantee that the cell lists will be buitl for the ypositions (by default they are constructed for the smallest set, which is faster).

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.PeriodicSystems, StaticArrays

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

julia> sys = PeriodicSystem(xpositions = positions, unitcell= [1.0,1.0,1.0], cutoff = 0.1, output = 0.0)
PeriodicSystem1 of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 1.0, 0.0, 0.0; 0.0, 1.0, 0.0; 0.0, 0.0, 1.0 ]
      cutoff = 0.1
      number of computing cells on each dimension = [12, 12, 12]
      computing cell sizes = [0.1, 0.1, 0.1] (lcell: 1)
      Total number of cells = 1728
    CellListMap.CellList{3, Float64}
      100 real particles.
      96 cells with real particles.
      164 particles in computing box, including images.
    Parallelization auxiliary data set for: 
      Number of batches for cell list construction: 8
      Number of batches for function mapping: 8
    Type of output variable: Float64

julia> map_pairwise!((x,y,i,j,d2,output) -> output += d2, sys)
0.14556244865996287

Two sets of particles

julia> using CellListMap.PeriodicSystems, StaticArrays

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

julia> sys = PeriodicSystem(
           xpositions = xpositions, 
           ypositions = ypositions,
           unitcell=[1,1,1], 
           cutoff = 0.1, 
           output = 0.0
           )
PeriodicSystem2 of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 1.0, 0.0, 0.0; 0.0, 1.0, 0.0; 0.0, 0.0, 1.0 ]
      cutoff = 0.1
      number of computing cells on each dimension = [12, 12, 12]
      computing cell sizes = [0.1, 0.1, 0.1] (lcell: 1)
      Total number of cells = 1728
    CellListMap.CellListPair{Vector{SVector{3, Float64}}, 3, Float64, CellListMap.Swapped}
       1000 particles in the reference vector.
       97 cells with real particles of target vector.
    Parallelization auxiliary data set for: 
      Number of batches for cell list construction: 8
      Number of batches for function mapping: 8
    Type of output variable: Float64

julia> map_pairwise!((x,y,i,j,d2,output) -> output += d2, sys)
2.3568143238242314

UpdatePeriodicSystem!

Internal function or structure - interface may change.

Updates the cell lists for periodic systems.

_reset_all_output!(output, output_threaded)

Internal function or structure - interface may change.

Function that resets the output variable and the threaded copies of it.

copy_output(x)

Function that 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.PeriodicSystems
# 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`
PeriodicSystems.copy_output(v::Vector{A}) = [ x for x in v ]

# Alternativelly, in this case, one could have defined:
Base.copy(a::A) = a
PeriodicSystems.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 PeriodicSystems.copy_output(v::Vector{T}) where {T<:CustomType} = deepcopy(v).

get_computing_box(sys::AbstractPeriodicSystem)

Internal function or structure - interface may change.

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

map_pairwise!(f::Function, system; show_progress = true, update_lists = 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 must be of the general form:

function f(x,y,i,j,d2,output)
    # operate on particle coordinates, distance and indexes
    # update output
    return output
end

where x and y are the coordinates (adjusted for the minimum image) of the two particles involved, i and j their indices in the original arrays of positions, d2 their squared Euclidean distance, and output the current value of the output variable. The output variable must be updated within this function with the contribution of the two particles involved.

Thread-safety is taken care automatically in parallel executions.

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

If update_lists is false, the cell lists will not be recomputed, this may be useful for computing a different function from the same coordinates.

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 = PeriodicSystem(
           xpositions = rand(SVector{3,Float64},1000), 
           unitcell=[1,1,1], 
           cutoff = 0.1, 
           output = 0.0
           );

julia> map_pairwise((x,y,i,j,d2,output) -> output += 1 / (1 + sqrt(d2)), sys)
1870.0274887950268

reduce_output!(reducer::Function, output, output_threaded)

Internal function or structure - interface may change.

Extended help

Function that defines how to reduce the vector of output variables, after a threaded computation. This function is implemented for output variables that are numbers, and vectors or arrays of number of static arrays, as the sum of the values of the threaded computations, which is the most common application, found in computing forces, energies, etc.

It may be interesting to implement custom PeriodicSystems.reduce_output! function for other types of output variables, considering:

• The arguments of the function must be the return output value and a vector output_threaded of output variables, which is created (automatically) by copying the output the number of times necessary for the multi-threaded computation.

• The function must return the output variable, independently of it being mutable or immutable.

reduce_output is an alias to reduce_output! that can be used for consistency if the output variable is immutable.

Example

In this example we show how to obtain the minimum distance between two sets of particles. This requires a custom reduction function.

using CellListMap.PeriodicSystems, StaticArrays
# Custom output type
struct MinimumDistance
    d::Float64
end
# Custom copy function for `Out`
PeriodicSystems.copy_output(d::MinimumDistance) = MinimumDistance(d.d)
# How to reset an array with elements of type `MinimumDistance`
PeriodicSystems.reset_output!(d::MinimumDistance) = MinimumDistance(+Inf)
# Custom reduction function (keep the minimum distance)
function PeriodicSystems.reduce_output!(
    output::MinimumDistance, 
    output_threaded::Vector{MinimumDistance}
)
    output = reset_output!(output)
    for i in eachindex(output_threaded)
        if output_threaded[i].d < output.d
            output = output_threaded[i]
        end
    end
    return output
end
# Construct the system
sys = PeriodicSystem(;
    xpositions = rand(SVector{3,Float64}, 1000),
    ypositions = rand(SVector{3,Float64}, 1000),
    unitcell = [1,1,1],
    cutoff = 0.1,
    output = MinimumDistance(+Inf),
)

# Obtain the minimum distance between the sets
map_pairwise!((x,y,i,j,d2,output) -> sqrt(d2) < output.d ? MinimumDistance(sqrt(d2)) : output, sys)
# will output something like: MinimumDistance(0.00956913034767034)

reducer!(x,y)

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

The most commont reducer is the sum, and this is how it is implemented for Union{Number, StaticArraysCore.FieldVector, StaticArraysCore.SVector}. For example, when computin 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 carefuly inspected by the user to avoid spurious results.

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), and for AbstractVecOrMat containers of Numbers or SVectors it is implemented as fill!(x, zero(eltype(x)).

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

The reset_output! function must return the output variable, being it mutable or immutable. reset_output is an alias for reset_output! that can be used for consistency if the output variable is immutable.

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).

# Custom data type
struct MinimumDistance d::Float64 end
# How to reset the minimum distance
PeriodicSystems.reset_output!(x::MinimumDistance) = MinimumDistance(+Inf)

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

Resizes the output array and the auxiliary output arrays used for multithreading, if needed because of the system change.

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

julia> using CellListMap.PeriodicSystems, StaticArrays

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

julia> update_cutoff!(sys, 0.2)
PeriodicSystem1 of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 1.0, 0.0, 0.0; 0.0, 1.0, 0.0; 0.0, 0.0, 1.0 ]
      cutoff = 0.2
      number of computing cells on each dimension = [7, 7, 7]
      computing cell sizes = [0.2, 0.2, 0.2] (lcell: 1)
      Total number of cells = 343
    CellListMap.CellList{3, Float64}
      1000 real particles.
      620 cells with real particles.
      1746 particles in computing box, including images.
    Parallelization auxiliary data set for: 
      Number of batches for cell list construction: 8
      Number of batches for function mapping: 12
    Type of output variable: Float64

update_unitcell!(system, unitcell)

Function to update the unit cell of the system. The unicell must be of the same type (OrthorhombicCell or 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 sytem (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.PeriodicSystems, StaticArrays

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

julia> update_unitcell!(sys, [1.2, 1.1, 1.0])
PeriodicSystem1 of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 1.2, 0.0, 0.0; 0.0, 1.1, 0.0; 0.0, 0.0, 1.0 ]
      cutoff = 0.1
      number of computing cells on each dimension = [13, 13, 12]
      computing cell sizes = [0.11, 0.1, 0.1] (lcell: 1)
      Total number of cells = 2028
    CellListMap.CellList{3, Float64}
      1000 real particles.
      633 cells with real particles.
      1703 particles in computing box, including images.
    Parallelization auxiliary data set for: 
      Number of batches for cell list construction: 8
      Number of batches for function mapping: 12
    Type of output variable: Float64

unitcelltype(sys::AbstractPeriodicSystem)

Returns the type of a unitcell from the PeriodicSystem structure.