Single set: Simple outputs

This section shows examples of computing simple outputs (a single scalar or a single array) for a single set of particles.

Potential energy example

For example, here we read the coordinates of Argon atoms from a PDB file. The coordinates are given as vector of SVectors. We then compute an "energy", which in this case is simply the sum of 1/d over all pair of particles, within a cutoff.

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

ParticleSystem1{energy} of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 21.0 0.0 0.0; 0.0 21.0 0.0; 0.0 0.0 21.0 ]
      cutoff = 8.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
    CellListMap.CellList{3, Float64}
      100 real particles.
      8 cells with real particles.
      800 particles in computing box, including images.
    Parallelization auxiliary data set for 8 batch(es).
    Type of output variable (energy): Float64

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

function energy(pair, energy)
    energy += 1 / pair.d
    return energy
end

energy (generic function with 1 method)

pairwise!(energy, system)

207.37593032505288

Because output_name was set to :energy, the system.energy field accesses the resulting value of the computation:

system.energy

207.37593032505288

If the output_name field is not provided, the output value from the system.output field.

The initial value of output

The output = 0.0 in the construction of the example above is used, by default, to set the type of the output variable, here a Float64. This value is stored in the output field, but it is, by default, reset to zero(typeof(output)) when calling the pairwise! function.

To use the initial value provided and accumulate on top of it in the call to pairwise!, the reset=false option must be provided, as:

pairwise!(energy, system; reset=false)

414.7518606501057

Note that the result is twice the previous value, because the initial value of the energy was not reset. If we instead call pairwise! with the default parameters, we recompute the energy from scratch:

pairwise!(energy, system) # reset=true by default

207.37593032505288

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!(pair, forces)
    (; i, j, x, y, d2, d) = pair
    df = (1/d2)*(1/d)*(y - x)
    forces[i] += df
    forces[j] -= df
    return forces
end

update_forces! (generic function with 1 method)

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:

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

ParticleSystem1{forces} of dimension 3, composed of:
    Box{CellListMap.OrthorhombicCell, 3}
      unit cell matrix = [ 21.0 0.0 0.0; 0.0 21.0 0.0; 0.0 0.0 21.0 ]
      cutoff = 8.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
    CellListMap.CellList{3, Float64}
      100 real particles.
      8 cells with real particles.
      800 particles in computing box, including images.
    Parallelization auxiliary data set for 8 batch(es).
    Type of output variable (forces): Vector{StaticArraysCore.SVector{3, Float32}}

A call to pairwise! with the appropriate function definition will update the forces:

pairwise!((pair, forces) -> update_forces!(pair, forces), system)

100-element Vector{StaticArraysCore.SVector{3, Float32}}:
 [0.026493823, 0.18454283, -0.012253985]
 [0.07782604, 0.27910823, 0.21926616]
 [-0.105047956, 0.23932211, -0.27744123]
 ⋮
 [0.11307244, 0.0063534956, -0.05955693]
 [-0.031012032, 0.03543647, 0.031849086]

If we want now to accumulate additional forces, we need to use reset=false in the call to pairwise!.