# ParticleSystem interface

The `ParticleSystem`

interface facilitates the use of `CellListMap`

for the majority of cases.

- This interface requires
`CellListMap.jl`

version`0.8.30`

or greater. - The complete codes of the examples are at the end of this page, with examples of:

The `ParticleSystem`

interface is available since version `0.9.0`

of CellListMap.jl. It replaces the `PeriodicSystems`

interface available in previous versions.

## The mapped function

The purpose of CellListMap is to compute a pairwise-dependent function for all pairs of particles that are closer to each other than a defined cutoff. This pairwise function must be implemented by the user and adhere to the following interface:

```
function f(x, y, i, j, d2, output)
# update output variable
return output
end
```

where `x`

and `y`

are the positions of the particles, already wrapped relative to each other according to the periodic boundary conditions (a minimum-image set of positions), `i`

and `j`

are the indexes of the particles in the arrays of coordinates, `d2`

is the squared distance between the particles, and `output`

is the variable to be computed.

#### Details of the mapped function interface

The input parameters `x`

, `y`

, `i`

, `j`

, and `d2`

must not be modified by the user. They are the the input data that the user may use to update the `output`

variable.

Input Parameter | Type | Meaning |
---|---|---|

`x` | `SVector` | The coordinates of particle `i` of the pair. |

`y` | `SVector` | The coordinates of particle `j` of the pair (minimum-image relative to `x` ). |

`i` | `Int` | Index of first particle in the original array of coordinates. |

`j` | `Int` | Index of second particle in the original array of coordinates. |

`d2` | `<:Real` | Squared distance between the particles. |

`output` | user defined | the value to be updated |

**Notes:** `x`

and `y`

may be 2D or 3D vectors, depending on the dimension of the system. The type of the coordinates of `x`

, `y`

, and of `d2`

are dependent on the input arrays and cutoff, and can be `Float64`

, `Float32`

, unitful quantities, etc.

Return value | Type | Meaning |
---|---|---|

`output` | user defined | the updated value of output. |

The `output`

variable **must** be returned by the function, being it mutable or immutable.

### Basic examples

For example, computing the energy, as the sum of the inverse of the distance between particles, can be done with a function like:

```
function energy(d2,u)
u += 1 / sqrt(d2)
return u
end
```

and the additional parameters required by the interface can be eliminated by the use of an anonymous function, directly on the call to the `map_pairwise`

function:

```
u = map_pairwise(
(x,y,i,j,d2,u) -> energy(d2,u),
system
)
```

(what `system`

is will be explained in the examples below). Note that the `energy`

function does not use the `x`

, `y`

, `i`

, and `j`

input parameters, such that the anonymous function managing the interface could also be written as `(_, _, _, _, d2, u) -> energy(d2, u)`

, making explicit the dummy character of these variables in the example.

Alternatively, the function might require additional parameters, such as the masses of the particles. In this case, we can use a closure to provide such data:

```
function energy(i,j,d2,u,masses)
u += masses[i]*masses[j] / sqrt(d2)
return u
end
const masses = # ... some masses
u = map_pairwise((x,y,i,j,d2,u) -> energy(d2,u,masses), system)
```

Here we reinforce the fact that the `energy`

functions defined above compute the contribution to the energy of the interaction of *a single* pair of particles. This function will be called for every pair of particles within the cutoff, automatically, in the `map_pairwise`

call.

The `output`

of the `CellListMap`

computation may be of any kind. Most commonly, it is an energy, a set of forces, or other data type that can be represented either as a number, an array of numbers, or an array of vectors (`SVectors`

in particular), such as an arrays of forces.

Additionally, the properties are frequently additive (the energy is the sum of the energy of the particles, or the forces are added by summation).

For these types of `output`

data the usage does not require the implementation of any data-type dependent function.

## The ParticleSystem constructor

## Potential energy example

For example, here we read the coordinates of Argon atoms from a PDB file. The coordinates are given as vector of `SVector`

s. We then compute an "energy", which in this case is simply the sum of `1/d`

over all pair of particles, within a cutoff.

The `ParticleSystem`

constructor receives the properties of the system and sets up automatically the most commonly used data structures necessary.

```
julia> using CellListMap, PDBTools
julia> argon_coordinates = coor(readPDB(CellListMap.argon_pdb_file))
julia> system = ParticleSystem(
xpositions=argon_coordinates,
unitcell=[21.0,21.0,21.0],
cutoff = 8.0,
output = 0.0,
output_name = :energy
);
```

- Systems can be 2 or 3-dimensional.
- The
`unitcell`

parameter may be:- a vector, in which case the system periodic boundaries are Orthorhombic, this is faster.
- a matrix, in which case the system periodic boundaries are Triclinic (general).
`nothing`

(by default), in which case no periodic boundary conditions will be used.

`Unitful`

quantities can be provided, given appropriate types for all input parameters.

Now, let us compute the energy of the particles, assuming a simple formula which depends on the inverse of the distance between pairs:

```
julia> function energy(x, y, i, j, d2, energy)
energy += 1 / sqrt(d2)
return energy
end
julia> map_pairwise!(energy, system)
207.37593043370865
```

Note that the first four parameters of `energy`

are not used here but are needed to adhere to the interface. The function input could be written as `(_, _, _, _, d2, energy)`

to make that explicit.

Because `output_name`

was set to `:energy`

, the `system.energy`

field accesses the resulting value of the computation:

```
julia> system.energy
207.37593043370865
```

If the `output_name`

field is not provided, the output value from the `system.output`

field.

## Computing forces

Following the example above, let us compute the forces between the particles. We have to define the function that computes the force between a pair of particles and updates the array of forces:

```
function update_forces!(x,y,i,j,d2,forces)
d = sqrt(d2)
df = (1/d2)*(1/d)*(y - x)
forces[i] += df
forces[j] -= df
return forces
end
```

Importantly, the function *must* return the `forces`

array to follow the API.

Now, let us setup the system with the new type of output variable, which will be now an array of forces with the same type as the positions:

```
julia> using CellListMap, PDBTools
julia> argon_coordinates = coor(readPDB(CellListMap.argon_pdb_file))
julia> system = ParticleSystem(
xpositions=argon_coordinates,
unitcell=[21.0, 21.0, 21.0],
cutoff = 8.0,
output = similar(argon_coordinates),
output_name = :forces
);
```

Let us note that the `forces`

where reset upon the construction of the system:

```
julia> system.forces
1000-element Vector{SVector{3, Float64}}:
[0.0, 0.0, 0.0]
[0.0, 0.0, 0.0]
⋮
[0.0, 0.0, 0.0]
```

A call to `map_pairwise!`

with the appropriate function definition will update the forces:

```
julia> map_pairwise!((x,y,i,j,d2,forces) -> update_forces!(x,y,i,j,d2,forces), system)
100-element Vector{SVector{3, Float64}}:
[0.026493833307357332, 0.18454277989323772, -0.012253902366284965]
[0.07782602581235695, 0.2791082233740261, 0.21926615329195248]
⋮
[0.11307234751448932, 0.006353545239676281, -0.05955687310348302]
[-0.03101200918307673, 0.03543655648545697, 0.031849121630976335]
```

## Computing both energy and forces

In this example we define a general type of `output`

variable, for which custom copy, reset, and reduction functions must be defined. It can be followed for the computation of other general properties from the particle positions.

Interface to be implemented:

Method | Return | What it does |
---|---|---|

`copy_output(x::T)` | new instance of type `T` | Copies an element of the output type `T` . |

`reset_output!(x::T)` | mutated `x` | Resets (usually zero) the value of x to the initial value it must assume before mapping. If `x` is immutable, the function can return a new instance of `T` . |

`reducer(x::T,y::T)` | `mutated x` | Reduces `x` and `y` into `x` (for example `x = x + y` ). If `x` is immutable, returns a new instance of type `T` . |

**Remark:** if the output is an array of an immutable type `T`

, the methods above can be defined for single *instances* of `T`

, which is simpler than for the arrays.

`using CellListMap, StaticArrays, PDBTools`

The computation of energies and forces in a single call is an interesting example for the definition of a custom `output`

type and the required interface functions. Let us first define an output variable containing both quantities:

```
mutable struct EnergyAndForces
energy::Float64
forces::Vector{SVector{3,Float64}}
end
```

Now we need to define what it means to copy, reset, and reduce this new type of output. We overload the default corresponding functions, for our new output type:

The copy method creates a new instance of the `EnergyAndForces`

type, with copied data:

```
function CellListMap.copy_output(x::EnergyAndForces)
return EnergyAndForces(copy(x.energy), copy(x.forces))
end
```

The reset method will zero both the energy and all forces:

```
function CellListMap.reset_output!(output::EnergyAndForces)
output.energy = 0.0
for i in eachindex(output.forces)
output.forces[i] = SVector(0.0, 0.0, 0.0)
end
return output
end
```

The reducer function defines what it means to combine two output variables obtained on independent threads. In this case, we sum the energies and forces. Different reduction functions might be necessary for other custom types (for example if computing minimum distances).

```
function CellListMap.reducer(x::EnergyAndForces, y::EnergyAndForces)
e_tot = x.energy + y.energy
x.forces .+= y.forces
return EnergyAndForces(e_tot, x.forces)
end
```

Note that in the above example, we reuse the `x.forces`

array in the return instance of `EnergyAndForces`

. You must always reduce from right to left, and reuse the possible buffers of the first argument of the reducer (in this case, `x`

).

- All these functions
**must**return the modified`output`

variable, to adhere to the interface. - The proper definition of a reduction function is crucial for correctness. Please verify your results if using the default reducer function, which sums the elements.

Now we can proceed as before, defining a function that updates the output variable appropriately:

```
function energy_and_forces!(x,y,i,j,d2,output::EnergyAndForces)
d = sqrt(d2)
output.energy += 1/d
df = (1/d2)*(1/d)*(y - x)
output.forces[i] += df
output.forces[j] -= df
return output
end
```

To finally define the system and compute the properties:

```
argon_coordinates = coor(readPDB(CellListMap.argon_pdb_file))
system = ParticleSystem(
xpositions = argon_coordinates,
unitcell = [21.0, 21.0, 21.0],
cutoff = 8.0,
output = EnergyAndForces(0.0, similar(argon_coordinates)),
output_name = :energy_and_forces
);
map_pairwise((x,y,i,j,d2,output) -> energy_and_forces!(x,y,i,j,d2,output), system);
```

The output can be seen with the aliases of the `system.output`

variable:

```
julia> system.energy_and_forces.energy
207.37593043370862
julia> system.energy_and_forces.forces
100-element Vector{SVector{3, Float64}}:
[0.02649383330735732, 0.18454277989323772, -0.012253902366284958]
[0.07782602581235692, 0.27910822337402613, 0.21926615329195248]
⋮
[0.11307234751448932, 0.006353545239676281, -0.05955687310348303]
[-0.031012009183076745, 0.03543655648545698, 0.03184912163097636]
```

## Updating coordinates, unit cell, and cutoff

If the `map_pairwise!`

function will compute energy and/or forces in a iterative procedure (a simulation, for instance), we need to update the coordinates, and perhaps the unit cell and the cutoff.

### Updating coordinates

The coordinates can be updated (mutated, or the array of coordinates can change in size by pushing or deleting particles), simply by directly accessing the `xpositions`

field of the system. The `xpositions`

array is a `Vector`

of `SVector`

(from `StaticArrays`

), with coordinates copied from the input array provided. Thus, the coordinates in the `ParticleSystem`

structure must be updated independently of updates in the original array of coordinates.

Let us exemplify the interface with the computation of forces:

```
julia> using CellListMap, StaticArrays
julia> positions = rand(SVector{3,Float64}, 1000);
julia> system = ParticleSystem(
xpositions = positions,
unitcell=[1,1,1],
cutoff = 0.1,
output = similar(positions),
output_name = :forces
);
julia> system.xpositions[1]
3-element SVector{3, Float64} with indices SOneTo(3):
0.6391290709055079
0.43679325975360894
0.8231829019768698
julia> system.xpositions[1] = zeros(SVector{3,Float64})
3-element SVector{3, Float64} with indices SOneTo(3):
0.0
0.0
0.0
julia> push!(system.xpositions, SVector(0.5, 0.5, 0.5))
1001-element Vector{SVector{3, Float64}}:
[0.0, 0.0, 0.0]
[0.5491373098208292, 0.23899915605319244, 0.49058287555218516]
⋮
[0.4700394061063937, 0.5440026379397457, 0.7411235688716618]
[0.5, 0.5, 0.5]
```

The `output`

variable may have to be resized accordingly, depending on the calculation being performed. Use the `resize_output!`

function (do **not** use `Base.resize!`

on your output array directly).

If the `output`

array has to be resized, that has to be done with the `resize_output!`

function, which will keep the consistency of the auxiliary multi-threading buffers. This is, for instance, the case in the example of computation of forces, as the `forces`

array must be of the same length as the array of positions:

```
julia> resize_output!(system, length(system.xpositions));
julia> map_pairwise!((x,y,i,j,d2,forces) -> update_forces!(x,y,i,j,d2,forces), system)
1001-element Vector{SVector{3, Float64}}:
[756.2076075886971, -335.1637545330828, 541.8627090466914]
[-173.02442398784672, -178.782819965489, 4.570607952876692]
⋮
[-722.5400961501635, 182.65287417718935, 380.0394926753039]
[20.27985502389337, -193.77607810950286, -155.28968519541544]
```

In this case, if the `output`

is not resized, a `BoundsError:`

is be obtained, because updates of forces at unavailable positions will be attempted.

### Updating the unit cell

The unit cell can be updated to new dimensions at any moment, with the `update_unitcell!`

function:

```
julia> update_unitcell!(system, SVector(1.2, 1.2, 1.2))
ParticleSystem1 of dimension 3, composed of:
Box{OrthorhombicCell, 3}
unit cell matrix = [ 1.2, 0.0, 0.0; 0.0, 1.2, 0.0; 0.0, 0.0, 1.2 ]
cutoff = 0.1
number of computing cells on each dimension = [13, 13, 13]
computing cell sizes = [0.11, 0.11, 0.11] (lcell: 1)
Total number of cells = 2197
CellListMap.CellList{3, Float64}
1000 real particles.
623 cells with real particles.
1719 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 (forces): Vector{SVector{3, Float64}}
```

The unit cell can be set initially using a vector or a unit cell matrix. If a vector is provided the system is considered Orthorhombic, if a matrix is provided, a Triclinic system is built. Unit cells updates must preserve the system type.

The unit cell of non-periodic systems (initialized with

`nothing`

) cannot be updated manually.It is recommended (but not mandatory) to use static arrays (or Tuples) to update the unitcell, as in this case the update will be non-allocating.

### Updating the cutoff

The cutoff can also be updated, using the `update_cutoff!`

function:

```
julia> update_cutoff!(system, 0.2)
ParticleSystem1 of dimension 3, composed of:
Box{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.
125 cells with real particles.
2792 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 (forces): Vector{SVector{3, Float64}}
julia> map_pairwise!((x,y,i,j,d2,forces) -> update_forces!(x,y,i,j,d2,forces), system)
1000-element Vector{SVector{3, Float64}}:
[306.9612911344924, -618.7375562535321, -607.1449767066479]
[224.0803003775478, -241.05319348787023, 67.53780411933884]
⋮
[2114.4873184508524, -3186.265279868732, -6777.748445712408]
[-25.306486853608945, 119.69319481834582, 104.1501577339471]
```

## Computations for two sets of particles

If the computation involves two sets of particle, a similar interface is available. The only difference is that the coordinates of the two sets must be provided to the `ParticleSystem`

constructor as the `xpositions`

and `ypositions`

arrays.

We will illustrate this interface by computing the minimum distance between two sets of particles, which allows us to showcase further the definition of custom type interfaces:

First, we define a variable type that will carry the indexes and the distance of the closest pair of particles:

```
julia> struct MinimumDistance
i::Int
j::Int
d::Float64
end
```

The function that, given two particles, retains the minimum distance, is:

```
julia> function minimum_distance(i, j, d2, md)
d = sqrt(d2)
if d < md.d
md = MinimumDistance(i, j, d)
end
return md
end
minimum_distance (generic function with 1 method)
```

We overload copy, reset, and reduce functions, accordingly:

```
julia> import CellListMap: copy_output, reset_output!, reducer!
julia> copy_output(md::MinimumDistance) = md
copy_output (generic function with 5 methods)
julia> reset_output!(md::MinimumDistance) = MinimumDistance(0, 0, +Inf)
reset_output! (generic function with 5 methods)
julia> reducer!(md1::MinimumDistance, md2::MinimumDistance) = md1.d < md2.d ? md1 : md2
reducer! (generic function with 2 methods)
```

Note that since `MinimumDistance`

is immutable, copying it is the same as returning the value. Also, resetting the minimum distance consists of setting its `d`

field to `+Inf`

. And, finally, reducing the threaded distances consists of keeping the pair with the shortest distance.

Next, we build the system

```
julia> xpositions = rand(SVector{3,Float64},1000);
julia> ypositions = rand(SVector{3,Float64},1000);
julia> system = ParticleSystem(
xpositions = xpositions,
ypositions = ypositions,
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = MinimumDistance(0,0,+Inf),
output_name = :minimum_distance,
)
```

And finally we can obtain the minimum distance between the sets:

```
julia> map_pairwise((x,y,i,j,d2,md) -> minimum_distance(i,j,d2,md), system)
MinimumDistance(276, 617, 0.006009804808785543)
```

## Additional options

- Turn parallelization on and off
- Displaying a progress bar
- Fine control of the parallelization
- Avoid cell list updating
- Control CellList cell size
- Coordinates as matrices

### Turn parallelization on and off

The use of parallel computations can be tunned on and of by the `system.parallel`

boolean flag. For example, using 6 cores (12 threads) for the calculation of the minimum-distance example:

```
julia> f(system) = map_pairwise((x,y,i,j,d2,md) -> minimum_distance(i,j,d2,md), system)
f (generic function with 1 method)
julia> Threads.nthreads()
8
julia> system.parallel = true
true
julia> @btime f($system)
268.265 μs (144 allocations: 16.91 KiB)
MinimumDistance(783, 497, 0.007213710914619913)
julia> system.parallel = false
false
julia> @btime f($system)
720.304 μs (0 allocations: 0 bytes)
MinimumDistance(783, 497, 0.007213710914619913)
```

### Displaying a progress bar

Displaying a progress bar: for very long runs, the user might want to see the progress of the computation. Use the `show_progress`

keyword parameter of the `map_pairwise!`

function for that.

For example, we execute the computation above, but with much more particles:

```
julia> xpositions = rand(SVector{3,Float64},10^6);
julia> ypositions = rand(SVector{3,Float64},10^6);
julia> system = ParticleSystem(
xpositions = xpositions,
ypositions = ypositions,
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = MinimumDistance(0,0,+Inf),
output_name = :minimum_distance,
);
julia> map_pairwise(
(x,y,i,j,d2,md) -> minimum_distance(i,j,d2,md), system;
show_progress = true
)
Progress: 24%|██████████▏ | ETA: 0:00:29
```

By activating the `show_progress`

flag, a nice progress bar is shown.

### Fine control of the parallelization

The number of batches launched in parallel runs can be tunned by the `nbatches`

keyword parameter of the `ParticleSystem`

constructor. By default, the number of batches is defined as heuristic function dependent on the number of particles, and possibly returns optimal values in most cases. For a detailed discussion about this parameter, see Number of batches.

For example, to set the number of batches for cell list calculation to 4 and the number of batches for mapping to 8, we can do:

```
julia> system = ParticleSystem(
xpositions = rand(SVector{3,Float64},1000),
unitcell=[1,1,1],
cutoff = 0.1,
output = 0.0,
output_name = :energy,
nbatches=(4,8), # use this keyword
);
```

Most times it is expected that the default parameters are optimal. But particularly for inhomogeneous systems increasing the number of batches of the mapping phase (second parameter of the tuple) may improve the performance by reducing the idle time of threads.

### Avoid cell list updating

To compute different properties without recomputing cell lists, use `update_lists=false`

in the call of `map_pairwise`

methods, for example,

```
using CellListMap, StaticArrays
system = ParticleSystem(xpositions=rand(SVector{3,Float64},1000), output=0.0, cutoff=0.1, unitcell=[1,1,1])
# First call, will compute the cell lists
map_pairwise((x,y,i,j,d2,u) -> u += d2, system)
# Second run: do not update the cell lists but compute a different property
map_pairwise((x,y,i,j,d2,u) -> u += sqrt(d2), system; update_lists = false)
```

in which case we are computing the sum of distances from the same cell lists used to compute the energy in the previous example (requires version 0.8.9). Specifically, this will skip the updating of the cell lists, thus be careful to not use this option if the cutoff, unitcell, or any other property of the system changed.

For systems with two sets of particles, the coordinates of the `xpositions`

set can be updated, preserving the cell lists computed for the `ypositions`

, but this requires setting `autoswap=false`

in the construction of the `ParticleSystem`

:

```
using CellListMap, StaticArrays
system = ParticleSystem(
xpositions=rand(SVector{3,Float64},1000),
ypositions=rand(SVector{3,Float64},2000),
output=0.0, cutoff=0.1, unitcell=[1,1,1],
autoswap=false # Cell lists are constructed for ypositions
)
map_pairwise((x,y,i,j,d2,u) -> u += d2, system)
# Second run: preserve the cell lists but compute a different property
map_pairwise((x,y,i,j,d2,u) -> u += sqrt(d2), system; update_lists = false)
```

### Control CellList cell size

The cell sizes of the construction of the cell lists can be controled with the keyword `lcell`

of the `ParticleSystem`

constructor. For example:

```
julia> system = ParticleSystem(
xpositions = rand(SVector{3,Float64},1000),
unitcell=[1,1,1],
cutoff = 0.1,
output = 0.0,
output_name = :energy,
lcell=2,
);
```

Most times using `lcell=1`

(default) or `lcell=2`

will provide the optimal performance. For very dense systems, or systems for which the number of particles within the cutoff is very large, larger values of `lcell`

may improve the performance. To be tested by the user.

The number of cells in which the particles will be classified is, for each dimension `lcell*length/cutoff`

. Thus if the `length`

of the box is too large relative to the `cutoff`

, many cells will be created, and this imposes a perhaps large memory requirement. Usually, it is a good practice to limit the number of cells to be not greater than the number of particles, and for that the cutoff may have to be increased, if there is a memory bottleneck. A reasonable choice is to use `cutoff = max(real_cutoff, length/n^(1/D))`

where `n`

is the number of particles and `D`

is the dimension (2 or 3). With that the number of cells will be close to `n`

in the worst case.

### Coordinates as matrices

Coordinates can also be provided as matrices of size `(D,N)`

where `D`

is the dimension (2 or 3) and `N`

is the number of particles. For example:

```
julia> using CellListMap
julia> system = ParticleSystem(
xpositions=rand(2,100),
ypositions=rand(2,200),
cutoff=0.1,
unitcell=[1,1],
output=0.0,
)
ParticleSystem2{output} of dimension 2, composed of:
Box{OrthorhombicCell, 2}
unit cell matrix = [ 1.0 0.0; 0.0 1.0 ]
cutoff = 0.1
number of computing cells on each dimension = [13, 13]
computing cell sizes = [0.1, 0.1] (lcell: 1)
Total number of cells = 169
CellListMap.CellListPair{Vector{StaticArraysCore.SVector{2, Float64}}, 2, Float64, CellListMap.Swapped}
200 particles in the reference vector.
61 cells with real particles of target vector.
Parallelization auxiliary data set for:
Number of batches for cell list construction: 1
Number of batches for function mapping: 1
Type of output variable (output): Float64
```

This interface less flexible than when the coordinates are input as vectors of vectors, because *the number of particles* cannot be changed, because matrices cannot be resized. Otherwise, matrices can be used as input.

## Complete example codes

### Simple energy computation

In this example, a simple potential energy defined as the sum of the inverse of the distance between the particles is computed.

```
using CellListMap
using StaticArrays
system = ParticleSystem(
xpositions = rand(SVector{3,Float64},1000),
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = 0.0,
output_name = :energy
)
map_pairwise!((x,y,i,j,d2,energy) -> energy += 1 / sqrt(d2), system)
```

### Force computation

Here we compute the force vector associated to the potential energy function of the previous example.

```
using CellListMap
using StaticArrays
positions = rand(SVector{3,Float64},1000)
system = ParticleSystem(
xpositions = positions,
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = similar(positions),
output_name = :forces
)
function update_forces!(x,y,i,j,d2,forces)
d = sqrt(d2)
df = (1/d2)*(1/d)*(y - x)
forces[i] += df
forces[j] -= df
return forces
end
map_pairwise!((x,y,i,j,d2,forces) -> update_forces!(x,y,i,j,d2,forces), system)
```

### Energy and forces

In this example, the potential energy and the forces are computed in a single run, and a custom data structure is defined to store both values.

```
using CellListMap
using StaticArrays
# Define custom type
mutable struct EnergyAndForces
energy::Float64
forces::Vector{SVector{3,Float64}}
end
# Custom copy, reset and reducer functions
import CellListMap: copy_output, reset_output!, reducer
copy_output(x::EnergyAndForces) = EnergyAndForces(copy(x.energy), copy(x.forces))
function reset_output!(output::EnergyAndForces)
output.energy = 0.0
for i in eachindex(output.forces)
output.forces[i] = SVector(0.0, 0.0, 0.0)
end
return output
end
function reducer(x::EnergyAndForces, y::EnergyAndForces)
e_tot = x.energy + y.energy
x.forces .+= y.forces
return EnergyAndForces(e_tot, x.forces)
end
# Function that updates energy and forces for each pair
function energy_and_forces!(x,y,i,j,d2,output::EnergyAndForces)
d = sqrt(d2)
output.energy += 1/d
df = (1/d2)*(1/d)*(y - x)
output.forces[i] += df
output.forces[j] -= df
return output
end
# Initialize system
positions = rand(SVector{3,Float64},1000);
system = ParticleSystem(
xpositions = positions,
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = EnergyAndForces(0.0, similar(positions)),
output_name = :energy_and_forces
)
# Compute energy and forces
map_pairwise((x,y,i,j,d2,output) -> energy_and_forces!(x,y,i,j,d2,output), system)
```

### Two sets of particles

In this example we illustrate the interface for the computation of properties of two sets of particles, by computing the minimum distance between the two sets.

```
using CellListMap
using StaticArrays
# Custom structure to store the minimum distance pair
struct MinimumDistance
i::Int
j::Int
d::Float64
end
# Function that updates the minimum distance found
function minimum_distance(i, j, d2, md)
d = sqrt(d2)
if d < md.d
md = MinimumDistance(i, j, d)
end
return md
end
# Define appropriate methods for copy, reset and reduce
import CellListMap: copy_output, reset_output!, reducer!
copy_output(md::MinimumDistance) = md
reset_output!(md::MinimumDistance) = MinimumDistance(0, 0, +Inf)
reducer!(md1::MinimumDistance, md2::MinimumDistance) = md1.d < md2.d ? md1 : md2
# Build system
xpositions = rand(SVector{3,Float64},1000);
ypositions = rand(SVector{3,Float64},1000);
system = ParticleSystem(
xpositions = xpositions,
ypositions = ypositions,
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = MinimumDistance(0,0,+Inf),
output_name = :minimum_distance,
)
# Compute the minimum distance
map_pairwise((x,y,i,j,d2,md) -> minimum_distance(i,j,d2,md), system)
```

### Particle simulation

In this example, a complete particle simulation is illustrated, with a simple potential. This example can illustrate how particle positions and forces can be updated. Run this simulation with:

```
julia> system = init_system(N=200); # number of particles
julia> trajectory = simulate(system);
julia> animate(trajectory)
```

One important characteristic of this example is that the `system`

is built outside the function that performs the simulation. This is done because the construction of the system is type-unstable (it is dimension, geometry and output-type dependent). Adding a function barrier avoids type-instabilities to propagate to the simulation causing possible performance problems.

```
using StaticArrays
using CellListMap
import CellListMap.wrap_relative_to
# Function that updates the forces, for potential of the form:
# if d < cutoff k*(d^2-cutoff^2)^2 else 0.0 with k = 10^6
function update_forces!(x, y, i, j, d2, forces, cutoff)
r = y - x
dudr = 10^6 * 4 * r * (d2 - cutoff^2)
forces[i] += dudr
forces[j] -= dudr
return forces
end
# Function that initializes the system: it is preferable to initialize
# the system outside the function that performs the simulation, because
# the system (data)type is defined on initialization. Initializing it outside
# the simulation function avoids possible type-instabilities.
function init_system(;N::Int=200)
Vec2D = SVector{2,Float64}
positions = rand(Vec2D, N)
unitcell = [1.0, 1.0]
cutoff = 0.1
system = ParticleSystem(
positions=positions,
cutoff=cutoff,
unitcell=unitcell,
output=similar(positions),
output_name=:forces,
)
return system
end
function simulate(system=init_system(); nsteps::Int=100, isave=1)
# initial velocities
velocities = [ randn(eltype(system.positions)) for _ in 1:length(system.positions) ]
dt = 1e-3
trajectory = typeof(system.positions)[]
for step in 1:nsteps
# compute forces at this step
map_pairwise!(
(x,y,i,j,d2,forces) -> update_forces!(x,y,i,j,d2,forces,system.cutoff),
system
)
# Update positions and velocities
for i in eachindex(system.positions, system.forces)
f = system.forces[i]
x = system.positions[i]
v = velocities[i]
x = x + v * dt + (f / 2) * dt^2
v = v + f * dt
# wrapping to origin for obtaining a pretty animation
x = wrap_relative_to(x, SVector(0.0, 0.0), system.unitcell)
# !!! IMPORTANT: Update arrays of positions and velocities
system.positions[i] = x
velocities[i] = v
end
# Save step for printing
if step % isave == 0
push!(trajectory, copy(system.positions))
end
end
return trajectory
end
using Plots
function animate(trajectory)
anim = @animate for step in trajectory
scatter(
Tuple.(step),
label=nothing,
lims=(-0.5, 0.5),
aspect_ratio=1,
framestyle=:box,
)
end
gif(anim, "simulation.gif", fps=10)
end
```

## Docstrings

`CellListMap.ParticleSystem`

— Method```
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::Any;
output_name::Symbol,
parallel::Bool=true,
nbatches::Tuple{Int,Int}=(0, 0),
autoswap::Bool = true,
validate_coordinates::Union{Nothing,Function}=_validate_coordinates
)
```

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

The `validate_coordinates`

function can be used to validate the coordinates before the construction of the system. If `nothing`

, no validation is performed. By default the validation checks if the coordinates are not missing or NaN.

**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: readPDB, coor
julia> positions = coor(readPDB(CellListMap.argon_pdb_file));
julia> sys = ParticleSystem(
xpositions = positions,
unitcell = [21.0, 21.0, 21.0],
cutoff = 8.0,
output = 0.0,
);
julia> map_pairwise!((x,y,i,j,d2,output) -> output += d2, sys)
43774.54367600001
```

**Two sets of particles**

```
julia> using CellListMap, PDBTools
julia> xpositions = coor(readPDB(CellListMap.argon_pdb_file))[1:50];
julia> ypositions = coor(readPDB(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> map_pairwise!((x,y,i,j,d2,output) -> output += d2, sys)
21886.196785000004
```

`CellListMap.copy_output`

— Method`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 ]
# Alternativelly, 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)`

`CellListMap.map_pairwise!`

— Method```
map_pairwise!(
f::Function, system::AbstractParticleSystem;
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 = ParticleSystem(
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
```

`CellListMap.reducer!`

— Method```
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 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.

**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(readPDB(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:
map_pairwise!((x,y,i,j,d2,output) -> sqrt(d2) < output.d ? MinimumDistance(sqrt(d2)) : output, sys)
MinimumDistance(2.1991993997816563)
```

`CellListMap.reset_output!`

— Method```
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)`

, and for `AbstractVecOrMat`

containers of `Number`

s or `SVector`

s it is implemented as `fill!(x, zero(eltype(x))`

.

Other custom output types must have their `reset_output!`

method implemented.

If the variable is mutable, 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 `reset_output!`

is preferred, by convention for mutating functions.

**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`

.

`CellListMap.resize_output!`

— Method`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.

This function must be implemented by the user if the output variable is a vector whose length is dependent on the number of particles. For example, if the output is a vector of forces acting on each particle, the output vector must be resized if the number of particles changes.

This function *must* be used in that case, to guarantee that the auxiliary arrays used for multi-threading are resized accordingly.

`CellListMap.unitcelltype`

— Method`unitcelltype(sys::AbstractParticleSystem)`

Returns the type of a unitcell from the `ParticleSystem`

structure.

`CellListMap.update_cutoff!`

— Method`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(readPDB(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{output} of dimension 3, composed of:
Box{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
CellList{3, Float64}
100 real particles.
8 cells with real particles.
800 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 (output): Float64
```

`CellListMap.update_unitcell!`

— Method`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 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.

Manual updating of the unit cell of non-periodic systems is not allowed.

**Example**

```
julia> using CellListMap, StaticArrays, PDBTools
julia> xpositions = coor(readPDB(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{output} of dimension 3, composed of:
Box{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
CellList{3, Float64}
100 real particles.
8 cells with real particles.
800 particles in computing box, including images.
Parallelization auxiliary data set for:
Number of batches for cell list construction: 1
Number of batches for function mapping: 1
Type of output variable (output): Float64
```

`CellListMap.ParticleSystem1`

— Type`mutable struct ParticleSystem1{OutputName, V, O, B, C, A, VC} <: CellListMap.AbstractParticleSystem{OutputName}`

`xpositions::Any`

`output::Any`

`_box::Any`

`_cell_list::Any`

`_output_threaded::Vector`

`_aux::Any`

`parallel::Bool`

`validate_coordinates::Any`

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 `ParticleSystem1`

structure is done through the `ParticleSystem(;xpositions, unitcell, cutoff, output)`

auxiliary function.

`CellListMap.ParticleSystem2`

— Type`mutable struct ParticleSystem2{OutputName, V, O, B, C, A, VC} <: CellListMap.AbstractParticleSystem{OutputName}`

`xpositions::Any`

`ypositions::Any`

`output::Any`

`_box::Any`

`_cell_list::Any`

`_output_threaded::Vector`

`_aux::Any`

`parallel::Bool`

`validate_coordinates::Any`

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 `ParticleSystem1`

structure is done through the `ParticleSystem(;xpositions, ypositions, unitcell, cutoff, output)`

auxiliary function.