# PeriodicSystems interface

The `PeriodicSystems`

interface facilitates the use of `CellListMap`

for the majority of cases. To use it, load the `PeriodicSystems`

module directly, with:

`using CellListMap.PeriodicSystems`

- This interface requires
`CellListMap.jl`

version`0.7.22`

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

## The mapped function

The function to be mapped for every pair of particles within the cutoff follows the same interface as the standard interface. It must be of the form

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

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

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

## Potential energy example

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 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 of `CellListMap.PeriodicSystems`

is the simplest, and does not require the implementation of any data-type dependent function.

For example, let us build a system of random particles in a cubic box, and compute an "energy", which in this case is simply the sum of `1/d`

over all pair of particles, within a cutoff.

The `PeriodicSystem`

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

```
julia> using CellListMap.PeriodicSystems, StaticArrays
julia> system = PeriodicSystem(
xpositions = rand(SVector{3,Float64},1000),
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = 0.0,
output_name = :energy
);
```

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

```
julia> map_pairwise!((x,y,i,j,d2,energy) -> energy += 1 / sqrt(d2), system)
30679.386366872823
```

The `system.energy`

field accesses the resulting value of the computation:

```
julia> system.energy
30679.386366872823
```

because the `output_name`

field was provided. If it is not provided, you can access the output value from the `system.output`

field.

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

parameter may be either a vector, as in the example, or a unit cell matrix, for general boundary conditions. `Unitful`

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

## 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> positions = rand(SVector{3,Float64},1000);
julia> system = PeriodicSystem(
xpositions = positions,
unitcell=[1,1,1],
cutoff = 0.1,
output = similar(positions),
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)
1000-element Vector{SVector{3, Float64}}:
[-151.19529230407284, 159.33819000196905, -261.3055111242796]
[-173.02442398784672, -178.782819965489, 4.570607952876692]
⋮
[-722.5400961501635, 182.65287417718935, 380.0394926753039]
```

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

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:

```
import CellListMap.PeriodicSystems: copy_output
copy_output(x::EnergyAndForces) = EnergyAndForces(copy(x.energy), copy(x.forces))
```

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

```
import CellListMap.PeriodicSystems: reset_output!
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
```

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

```
import CellListMap.PeriodicSystems: reducer
function 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:

```
positions = rand(SVector{3,Float64},1000);
system = PeriodicSystem(
xpositions = positions,
unitcell=[1.0,1.0,1.0],
cutoff = 0.1,
output = EnergyAndForces(0.0, similar(positions)),
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
31696.94766439311
julia> system.energy_and_forces.forces
1000-element Vector{SVector{3, Float64}}:
[-338.1909601911842, 7.7663690656924445, 202.25889647151405]
[33.67299655756128, 282.7581453168999, -79.09639223837306]
⋮
[38.83014327604529, -204.45236278342745, 249.307871211616]
```

## 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 acessing 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 `PeriodicSystem`

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.PeriodicSystems, StaticArrays
julia> positions = rand(SVector{3,Float64}, 1000);
julia> system = PeriodicSystem(
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))
PeriodicSystem1 of dimension 3, composed of:
Box{CellListMap.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.

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

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.PeriodicSystems: 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 = PeriodicSystem(
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 execution options

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

### 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 = PeriodicSystem(
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 paralellization

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

keyword parameter of the `PeriodicSystem`

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 dicussion 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 = PeriodicSystem(
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.PeriodicSystems, StaticArrays
system = PeriodicSystem(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 `PeriodicSystem`

:

```
using CellListMap.PeriodicSystems, StaticArrays
system = PeriodicSystem(
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 constructred 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 `PeriodicSystem`

constructor. For example:

```
julia> system = PeriodicSystem(
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.

## 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.PeriodicSystems
using StaticArrays
system = PeriodicSystem(
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.PeriodicSystems
using StaticArrays
positions = rand(SVector{3,Float64},1000)
system = PeriodicSystem(
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.PeriodicSystems
using StaticArrays
# Define custom type
mutable struct EnergyAndForces
energy::Float64
forces::Vector{SVector{3,Float64}}
end
# Custom copy, reset and reducer functions
import CellListMap.PeriodicSystems: 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 = PeriodicSystem(
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.PeriodicSystems
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.PeriodicSystems: 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 = PeriodicSystem(
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.PeriodicSystems
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 preferrable 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 = PeriodicSystem(
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
```