Ecosystem integration

• Agents.jl
• Unitful and units
• Automatic differentiation
• Measurements

Agents.jl

Agents.jl provides a comprehensive framework for simulation, analysis and visualization of agent-based systems. CellListMap can be used to accelerate these simulations, and the integration of the packages is rather simple, particularly using the ParticleSystem interface. A complete integration example can be obtained in the Agents documentation (currently at the development branch).

The example will produce the following animation:

Unitful and units

The functions of CellListMap.jl support the propagation of generic (isbits) types, and thus units and thus automatic differentiation and the use of Unitful. A set of working examples can be found in the generic_types.jl file.

We start illustrating the support for unit propagation. We need to define all involved quantities in the same units:

Using the ParticleSystem interface

The only requirement is to attach proper units to all quantities (positions, cutoff, unitcell, and output variables). Here we compute the square of the distances of the particles within the cutoff, where the particle coordinates are in Angstroms, while the box size and cutoff are defined in nanometers:

using CellListMap, Unitful, PDBTools
positions = coor(read_pdb(CellListMap.argon_pdb_file))u"Å"
system = ParticleSystem(
    positions = positions,
    cutoff = 0.8u"nm",
    unitcell = [2.1,2.1,2.1]u"nm",
    output = 0.0u"nm^2",
    output_name = :sumsqr
)
pairwise!((pair, out) -> out += pair.d2, system)

437.74544424050725 nm^2

Units in neighbor lists

CellListMap.neighborlist also propagates units correctly. Continuing the example above:

cutoff = 0.8u"nm";
neighborlist(positions, cutoff)

857-element Vector{Tuple{Int64, Int64, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}:
 (2, 5, 0.4881879705003039 nm)
 (2, 10, 0.47976206386342274 nm)
 (2, 27, 0.7312039751474433 nm)
 ⋮
 (46, 18, 0.6114386195284445 nm)
 (46, 51, 0.7999472507177153 nm)

Automatic differentiation

Allowing automatic differentiation follows the same principles, meaning that we only need to allow the propagation of dual types through the computation by proper initialization of the input data. However, it is easier to work with the low level interface, which accepts matrices as the input for positions and a more fine control of the types of the variables. Matrices are easier input types for auto diff packages.

The variables are each component of each vector, thus the easiest way to represent the points to interface with differentiation packages is providing the coordinates as a matrix:

julia> x = rand(3,1000)
3×1000 Matrix{Float64}:
 0.186744  0.328719  0.874102  0.503535   …  0.328161  0.0895699  0.917338
 0.176157  0.972954  0.80729   0.624724      0.655268  0.470754   0.327578
 0.648482  0.537362  0.599624  0.0688776     0.92333   0.497984   0.208924

The key here is allow all the types of the parameters to follow the type propagation of the elements of x inside the differentiation routine. The function we define to compute the derivative is, then:

julia> function sumsqr(x, sides, cutoff)
           sys = ParticleSystem(
               positions=x,
               unitcell=eltype(x).(sides),
               cutoff=eltype(x).(cutoff),
               output=zero(eltype(x))
           )
           return  pairwise!((pair, sumsqr) -> sumsqr += pair.d2, sys)
       end

Note that we convert cutoff and sides to the same type of the input x of the function, and set the type of the output variable accordingly. For a simple call to the function this is inconsequential:

julia> cutoff = 0.1; sides = [1,1,1];

julia> sumsqr(x,sides,cutoff)
12.897650398753228

but the conversion is required to allow the differentiation to take place:

julia> ForwardDiff.gradient(x -> sumsqr(x,sides,cutoff),x)
3×1000 Matrix{Float64}:
 -0.132567   0.029865  -0.101301  …   0.249267    0.0486424  -0.0400487
  0.122421   0.207495  -0.184366     -0.201648   -0.105031    0.218342
  0.0856502  0.288924   0.122445     -0.0147022  -0.103314   -0.0862264

Measurements

Propagating uncertainties through the Measurements and other similar packages requires a different strategy, because within CellListMap only isbits types can be used, which is not the case of the type Measurement type.

In cases like this, it is better to bypass all the internals of CellListMap and provide the data to the function that computes pairwise properties directly as a closure. For example:

A vector of particles with uncertainties in their coordinates can be created with:

using CellListMap
using Measurements
using StaticArrays
using LinearAlgebra: norm

x_input = [ SVector{3}(measurement(rand(),0.01*rand()) for i in 1:3) for j in 1:1000 ]

1000-element Vector{StaticArraysCore.SVector{3, Measurements.Measurement{Float64}}}:
 [0.9764 ± 0.0067, 0.6208 ± 0.0029, 0.2386 ± 0.005]
 [0.3477 ± 0.0053, 0.6688 ± 0.0031, 0.3763 ± 0.0066]
 [0.2236 ± 0.0053, 0.6856 ± 0.0015, 0.3659 ± 0.0056]
 ⋮
 [0.65181 ± 0.00011, 0.3951 ± 0.0078, 0.1863 ± 0.0044]
 [0.0335 ± 0.0063, 0.6502 ± 0.0063, 0.4474198 ± 6.7e-6]

The variables within the CellListMap functions will be stripped from the uncertainties, and the ParticleSystem will be built from these coordinates:

x_strip = [ getproperty.(v,:val) for v in x_input ]

1000-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [0.9764189671341497, 0.620836977276611, 0.23864680861354182]
 [0.34774553366697714, 0.6687895016971056, 0.3763242504742369]
 [0.2235575914359884, 0.685594571431364, 0.3659492427810984]
 ⋮
 [0.6518088165679041, 0.3950981395242523, 0.18628429368209154]
 [0.033540411268292325, 0.6501861569080413, 0.4474198219993799]

sys = ParticleSystem(
    positions=x_strip,
    unitcell=[1,1,1],
    cutoff=0.1,
    output=0.0
)

ParticleSystem1{default_output_name} 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 = [13, 13, 13]
      computing cell sizes = [0.1, 0.1, 0.1] (lcell: 1)
      Total number of cells = 2197
    CellListMap.CellList{3, Float64}
      1000 real particles.
      626 cells with real particles.
      1710 particles in computing box, including images.
    Parallelization auxiliary data set for 10 batch(es).
    Type of output variable (default_output_name): Float64

This would be result of the sum of distances:

pairwise!((pair, sumd) -> sumd += pair.d, sys)

162.77337955635903

To compute the same result but propagating the uncertainties, the computation for each pair has to capture the coordinates with uncertainties. To be tread-safe, then, one has to manually guarantee that the sum of squares cannot occur concurrently. Here, then, we introduce a lock:

function sumd_with_uncertainties(sys, x_input)
    lk = ReentrantLock()
    sumd = measurement(0.,0.)
    pairwise!(
        (pair, _) -> begin
            (; i, j) = pair
            x1 = x_input[i]
            x2 = CellListMap.wrap_relative_to(x_input[j], x1, sys.unitcell)
            @lock lk sumd += norm(x2 - x1)
            return 0.0 # not used
        end,
        sys
    )
    return sumd
end
sumd_with_uncertainties(sys, x_input)

\[162.77 \pm 0.37\]

In the function above, the pair.x and pair.y coordinates, which correspond to the coordinates in x_input[pair.i] and x_input[pair.j], but already wrapped relative to each other, are ignored, because they don't carry the uncertainties. We use only the indexes pair.i and pair.j to recompute the relative position of the particles according to the periodic boundary conditions (using the CellListMap.wrap_relative_to function) and their (squared) distance. Since the x_input array carries the uncertainties, the computation of sumsqr will propagate them.