Neighbor lists

Neighbor lists can be computed, returning all pairs of particles that are found within the cutoff, and the corresponding distances.

Non-periodic systems
Periodic systems
In-place computation of neighbor lists
Options

Note

When computing neighbor lists with cell-lists, it is possible for pairs of particles that are at a distance equal to the cutoff to either be included or excluded due to numerical rounding. As a result, these neighbor lists should only be utilized for calculating properties that vanish at the cutoff distance.

Non-periodic systems

Without periodic boundary conditions, just provide the coordinates and the cutoff:

julia> using CellListMap

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

julia> neighborlist(x,0.05)
376457-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 363, 0.04855594810064624)
 (1, 513, 0.03356381123125866)
 (1, 1209, 0.005159666709130686)
 ⋮
 (6575, 7378, 0.03791567990447959)
 (7378, 3450, 0.01748757015908321)

If the neighbor lists between two sets of points are required, use the following notation, in this case using coordinates as arrays of static arrays:

julia> using StaticArrays

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

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

julia> list = neighborlist(x,y,0.1)
37309-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 971, 0.09867846773727411)
 (1, 567, 0.06630101425431004)
 (1, 3, 0.04103170149300593)
 ⋮
 (10000, 156, 0.08549899843141298)
 (10000, 444, 0.0737386384422871)

where, similarly, the third parameter is the cutoff. The returning array contains tuples with the index of the particle in the first vector, the index of the particle in the second vector, and their distance.

Periodic systems

If periodic boundary conditions are used, the unitcell can be provided explicitly as keyword parameters:

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

julia> neighborlist(x, 0.05; unitcell=[1,1])
392100-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 5, 0.03445098850037766)
 (1, 393, 0.039448810592487206)
 (1, 1632, 0.02276457565643465)
 ⋮
 (9501, 9781, 0.03351665194098955)
 (9501, 5429, 0.04199258248973222)

In the example above, an Orthorhombic cell was assumed, and thus a vector of sides was provided. For general periodic boundary conditions, a unit cell matrix can be provided, for example:

julia> neighborlist(x, 0.05; unitcell=[1.0 0.5; 0.5 1.0])
580693-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 457, 0.03935441952786555)
 (1, 1467, 0.033407692174569875)
 (1, 1767, 0.04490555313598093)
 ⋮
 (3652, 8475, 0.04721628783510375)
 (6260, 8475, 0.04946130971686825)

Note

Positions and unit cells can be 2 or 3-dimensional.

In-place computation of neighbor lists

If neighbor lists are computed within a interactive scenario, it is interesting preallocate all the necessary data and just update the lists at every iteration. This can be achieved by constructing the InPlaceNeighborList object in advance. The performance gain of performing the operations in place might vary and may not be important for single runs, as the allocations do not dominate the computing time.

We will first illustrate the interface for a non-parallel run:

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)
InPlaceNeighborList with types: 
CellList{3, Float64}
Box{OrthorhombicCell, 3, Float64, Float64, 9}
Current list buffer size: 0

Note that the buffer size has size 0. The first time the neighbor lists are computed, the list will be allocated. We will use the neighborlist! (with the bang) function, because it will effectively mutate the system, by allocating all necessary data:

julia> @time list = neighborlist!(system)
  0.017765 seconds (12 allocations: 7.445 MiB)
209190-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 1375, 0.09425551992016712)
 (1, 3076, 0.045320021406080775)
 (1, 3666, 0.07780146666634076)
 ⋮
 (9962, 6983, 0.07355578793348823)
 (9962, 7457, 0.07597724209140656)

Now, if we modify the coordinates, we can update the system and recompute the neighbor lists:

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

julia> @time update!(system, x_new)
  0.003562 seconds
InPlaceNeighborList with types: 
CellList{3, Float64}
Box{OrthorhombicCell, 3, Float64, Float64, 9}
Current list buffer size: 209190

julia> @time list = neighborlist!(system);
  0.012338 seconds

Note

Here we illustrate the behavior of the functions in their second calls, to remove the effects of compilation on the allocation results.
The cutoff and unitcell can be modified by providing additional keyword parameters to the update! function (for example update!(system, x; cutoff=0.1)).
Allocations can occur if the cutoff, unit cell, or number of particles change such that greater buffers are required. The number of allocations tend to diminish as the buffers become large enough to accommodate the possible variations of the computation.

For parallel runs, the allocations are minimal, but some small auxiliary data is required for the launching of multiple threads. We illustrate here the convergence of the allocations to the minimum required for multi-threaded calculations:

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

julia> @time list = neighborlist!(system);
  0.007762 seconds (230 allocations: 18.142 MiB)

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

julia> @time update!(system, x_new)
  0.005283 seconds (20.30 k allocations: 6.200 MiB)
InPlaceNeighborList with types: 
CellList{3, Float64}
Box{OrthorhombicCell, 3, Float64, Float64, 9}
Current list buffer size: 209190

julia> @time neighborlist!(system);
  0.008190 seconds (166 allocations: 6.461 MiB)

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

julia> @time update!(system, x_new);
  0.002723 seconds (221 allocations: 208.922 KiB)

julia> @time neighborlist!(system);
  0.006227 seconds (165 allocations: 2.863 MiB)

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

julia> @time update!(system, x_new);
  0.002396 seconds (275 allocations: 144.078 KiB)

julia> @time neighborlist!(system);
  0.004996 seconds (161 allocations: 15.141 KiB)

Options

Additional optional parameters can be used in a neighborlist call:

Keyword	Values types	Default	About
`parallel`	`Bool`	`true`	turns on and off parallelization
`show_progress`	`Bool`	`false`	turns on and off progress bar
`nbatches`	`Tuple{Int,Int}`	`(0,0)`	Number of batches used in parallelization (see here)
`autoswap`	`Bool`	`true`	(advanced) automatically choose set to construct the cell lists

Docstrings

CellListMap.neighborlist — Method

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.

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(readPDB(CellListMap.argon_pdb_file, "index <= 50"));

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

julia> CellListMap.neighborlist(x, y, 8.0; parallel=false)
439-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 13, 7.0177629626541265)
 (1, 24, 7.976895636774999)
 (1, 29, 3.1770283284856)
 (1, 11, 4.0886518560523095)
 (1, 17, 5.939772807102978)
 (1, 30, 2.457228927063981)
 (1, 44, 5.394713986857875)
 (1, 45, 5.424876588458028)
 (2, 2, 3.9973374888793174)
 (2, 6, 5.355242104704511)
 ⋮
 (50, 27, 6.257296620746054)
 (50, 32, 3.109966559305742)
 (50, 33, 2.9192916949150525)
 (50, 35, 5.043227240567294)
 (50, 10, 3.9736202636890208)
 (50, 20, 6.995405134800989)
 (50, 39, 3.9001540995196584)
 (50, 37, 7.5464903100713)
 (50, 3, 7.232267901564487)

Now, considering the system periodic:

julia> using CellListMap, PDBTools

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

julia> y = coor(readPDB(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, 12, 7.360804915224965)
 (1, 13, 7.017762962654125)
 (1, 24, 7.976895636774997)
 (1, 29, 3.177028328485598)
 (1, 44, 5.394713986857875)
 (1, 45, 5.4248765884580274)
 (1, 11, 4.088651856052312)
 (1, 17, 5.93977280710298)
 (1, 30, 2.457228927063981)
 (1, 28, 6.853834401267658)
 ⋮
 (50, 3, 7.232267901564487)
 (50, 10, 3.9736202636890203)
 (50, 27, 6.257296620746054)
 (50, 32, 3.1099665593057426)
 (50, 33, 2.919291694915052)
 (50, 35, 5.043227240567294)
 (50, 20, 6.995405134800987)
 (50, 37, 7.546490310071297)
 (50, 39, 3.900154099519657)

Note

We set parallel=false in these examples to preserve the order of the pairs in the list. Parallel executions will not guarantee the order of the pairs.

source

CellListMap.neighborlist — Method

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(readPDB(CellListMap.argon_pdb_file));

julia> neighborlist(x, 8.0; parallel=false)
857-element Vector{Tuple{Int64, Int64, Float64}}:
 (1, 20, 3.163779543520692)
 (1, 61, 4.0886518560523095)
 (1, 67, 5.939772807102978)
 (1, 80, 2.457228927063981)
 (1, 94, 5.394713986857875)
 (13, 15, 2.678764267344178)
 (13, 41, 4.408015539900014)
 (13, 44, 6.960112211739117)
 (13, 61, 5.939197673086828)
 (13, 64, 4.560755858407684)
 ⋮
 (46, 18, 6.114385414741209)
 (46, 51, 7.999472795128439)
 (51, 68, 2.200357470957844)
 (51, 22, 6.638020940009152)
 (54, 45, 4.423308377221737)
 (73, 78, 2.853611220891874)
 (73, 88, 6.078711047582372)
 (78, 88, 7.006116541993863)
 (88, 54, 7.933654076149277)

And now, considering the system periodic:

julia> using CellListMap, PDBTools

julia> x = coor(readPDB(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.36387559222989)
 (1, 20, 3.163779543520693)
 (1, 47, 6.243868272153088)
 (1, 63, 7.017762962654125)
 (1, 74, 7.976895636774997)
 (1, 79, 3.177028328485598)
 (1, 94, 5.394713986857875)
 (1, 95, 5.4248765884580274)
 (7, 20, 3.3995637955478935)
 (7, 26, 7.96292025578556)
 ⋮
 (57, 34, 6.536566147450816)
 (57, 84, 7.225401442134547)
 (57, 88, 7.971591246420004)
 (68, 14, 5.2021891545771375)
 (68, 34, 3.955899012866733)
 (68, 84, 5.650943284089833)
 (68, 88, 7.254140403934848)
 (68, 38, 7.4092885623384905)
 (68, 90, 7.875801229081395)

Note

We set parallel=false in these examples to preserve the order of the pairs in the list. Parallel executions will not guarantee the order of the pairs.

source

CellListMap.update! — Method

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)

source

CellListMap.InPlaceNeighborList — Method

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)

source