NBodySimulator

Simulating systems of N interacting bodies.

Tutorials and Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

Example

using NBodySimulator
using StaticArrays
using Plots
body1 = MassBody(SVector(0.0, 1.0, 0.0), SVector(5.775e-6, 0.0, 0.0), 2.0)
body2 = MassBody(SVector(0.0, -1.0, 0.0), SVector(-5.775e-6, 0.0, 0.0), 2.0)
G = 6.673e-11
system = GravitationalSystem([body1, body2], G)
tspan = (0.0, 1111150.0)
simulation = NBodySimulation(system, tspan)
sim_result = run_simulation(simulation)
animate(sim_result, "path_to_animated_particles.gif")

<img src="https://user-images.githubusercontent.com/16945627/39958539-d2cf779c-561d-11e8-96a8-ffc3a595be8b.gif" alt="Here should appear a gif of rotating bodies" width="350"/>

NBodySimulation(
    system, tspan, boundary_conditions = InfiniteBox(),
    thermostat = NullThermostat(), kb = 1.38e-23
)

Configuration for an N-body simulation.

Arguments

• system: An NBodySystem describing the bodies and interactions.
• tspan: Tuple with the simulation start and end times.
• boundary_conditions: Boundary-condition model.
• thermostat: Thermostat model.
• kb: Boltzmann constant in the selected unit system.

Fields

• system, tspan: Physical system and integration interval.
• boundary_conditions, thermostat: Simulation constraints and temperature model.
• kb: Boltzmann constant.
• external_electric_field, external_magnetic_field, external_gravitational_field: External-field callables.

Examples

using NBodySimulator, StaticArrays

body = MassBody(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.0)
simulation = NBodySimulation(GravitationalSystem([body], 1.0), (0.0, 1.0))

run_simulation(simulation, algorithm = Tsit5(); kwargs...)

Solve simulation and return an interpolating simulation result.

Arguments

• simulation: Simulation configuration to solve.
• algorithm: Optional compatible OrdinaryDiffEq algorithm.

Keyword Arguments

• kwargs...: Keyword arguments forwarded to solve.

Examples

result = run_simulation(simulation)

Body

Abstract supertype for particle definitions used by N-body systems.

Body implementations provide an initial position, velocity, and mass. The built-in MassBody, ChargedParticle, and MagneticParticle types cover the standard interaction models.

MassBody(r, v, m)

Particle with position r, velocity v, and mass m.

Arguments

• r: Three-dimensional initial position.
• v: Three-dimensional initial velocity.
• m: Particle mass.

Fields

• r: Initial position.
• v: Initial velocity.
• m: Mass.

Examples

using NBodySimulator, StaticArrays

body = MassBody(SVector(0.0, 0.0, 0.0), SVector(0.0, 1.0, 0.0), 1.0)

ChargedParticle(r, v, m, q)

Particle with position r, velocity v, mass m, and charge q.

Arguments

• r: Three-dimensional initial position.
• v: Three-dimensional initial velocity.
• m: Particle mass.
• q: Electric charge.

Fields

• r, v, m: Position, velocity, and mass.
• q: Electric charge.

Examples

using NBodySimulator, StaticArrays

particle = ChargedParticle(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.0, 1.0)

MagneticParticle(r, v, m, mm)

Particle with position r, velocity v, mass m, and magnetic moment mm.

Arguments

• r: Three-dimensional initial position.
• v: Three-dimensional initial velocity.
• m: Particle mass.
• mm: Three-dimensional magnetic moment.

Fields

• r, v, m: Position, velocity, and mass.
• mm: Magnetic moment.

Examples

using NBodySimulator, StaticArrays

particle = MagneticParticle(
    SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.0, SVector(0.0, 0.0, 1.0)
)

generate_bodies_in_cell_nodes(n, m, v_dev, L; rng = MersenneTwister(n))

Generate n equal-mass bodies on a cubic-cell grid with normally distributed velocities.

Arguments

• n: Number of bodies to generate.
• m: Mass assigned to every body.
• v_dev: Standard deviation of each velocity component.
• L: Side length of the cubic cell.

Keyword Arguments

• rng: Random-number generator used for velocities.

Examples

using NBodySimulator

bodies = generate_bodies_in_cell_nodes(8, 1.0, 0.1, 4.0)

PotentialParameters

Abstract interface for potential or force-field parameter types.

PotentialParameters subtypes describe one interaction law in a PotentialNBodySystem. Built-in subtypes include Lennard-Jones, electrostatic, magnetostatic, gravitational, and SPC/Fw water-model parameters.

Interface

• Subtypes should be immutable parameter containers.
• To participate in time stepping, define get_accelerating_function(parameters, simulation).
• The returned acceleration function must support acceleration!(dv, u, v, t, i), mutate only dv, and compute the acceleration contribution for coordinate index i.
• u and v are 3 x n coordinate and velocity arrays. Implementations should not assume Array; use indexing, broadcasting, and eltype generically where possible.

Examples

struct ConstantAcceleration <: PotentialParameters
    a::Float64
end

function NBodySimulator.get_accelerating_function(p::ConstantAcceleration, simulation)
    return (dv, u, v, t, i) -> (dv .= (p.a, 0, 0))
end

get_accelerating_function(parameters::PotentialParameters, simulation::NBodySimulation)

Return the acceleration kernel for a potential parameter object and simulation.

This is the extension interface used by PotentialNBodySystem when converting an NBodySimulation to an ODE problem. Package users can add custom potentials by subtyping PotentialParameters and defining this method for their subtype.

Arguments

• parameters::PotentialParameters: Interaction parameters for one potential.
• simulation::NBodySimulation: Simulation containing bodies, boundary conditions, thermostat, and physical constants.

Returns

• acceleration!: A callable with signature acceleration!(dv, u, v, t, i).

Interface

• acceleration! must mutate dv with the acceleration contribution for coordinate index i.
• acceleration! must treat u and v as read-only.
• u and v are 3 x n coordinate and velocity arrays. Implementations should preserve the element type of dv, u, and v where possible.
• The method may close over precomputed constants from parameters and simulation, but it should not mutate simulation.

Examples

using NBodySimulator, StaticArrays

struct ConstantAcceleration <: PotentialParameters
    a::Float64
end

function NBodySimulator.get_accelerating_function(
        p::ConstantAcceleration,
        simulation::NBodySimulation
    )
    return (dv, u, v, t, i) -> (dv .= SVector(p.a, 0.0, 0.0))
end

GravitationalParameters(G)

Gravitational interaction parameters with gravitational constant G.

Arguments

• G: Gravitational constant.

Fields

• G: Gravitational constant.

Examples

using NBodySimulator

parameters = GravitationalParameters(6.67408e-11)

MagnetostaticParameters(μ_4π)

Magnetostatic interaction parameters storing μ / 4π.

Arguments

• μ_4π: Magnetic interaction constant μ / 4π.

Fields

• μ_4π: Magnetic interaction constant.

Examples

using NBodySimulator

parameters = MagnetostaticParameters(1.0e-7)

ElectrostaticParameters(k, R)

Electrostatic interaction parameters with Coulomb constant k and cutoff R.

Arguments

• k: Coulomb constant.
• R: Interaction cutoff radius. Omit it for no cutoff.

Fields

• k, R: User-supplied electrostatic parameters.
• R2: Cached squared cutoff radius.

Examples

using NBodySimulator

parameters = ElectrostaticParameters(8.9875517923e9, 5.0)

LennardJonesParameters(ϵ, σ, R)

Lennard-Jones potential parameters with energy scale ϵ, length scale σ, and cutoff R.

Arguments

• ϵ: Potential-well depth.
• σ: Particle diameter.
• R: Interaction cutoff radius.

Fields

• ϵ, σ, R: User-supplied potential parameters.
• σ2, R2: Cached squared length scales.

Examples

using NBodySimulator

parameters = LennardJonesParameters(1.0, 1.0, 2.5)

SPCFwParameters(rOH, aHOH, kb, ka)

SPC/Fw water-model parameters for bond length, bond angle, bond stiffness, and angle stiffness.

Arguments

• rOH: Equilibrium oxygen-hydrogen bond length.
• aHOH: Equilibrium hydrogen-oxygen-hydrogen angle in radians.
• kb: Harmonic-bond stiffness.
• ka: Harmonic-angle stiffness.

Fields

• rOH, aHOH: Equilibrium molecular geometry.
• kb, ka: Bond and angle stiffnesses.

Examples

using NBodySimulator

parameters = SPCFwParameters(0.1012, 1.9764, 44315.0, 317.6)

Abstract supertype for thermostat configurations.

Thermostats control a simulation's temperature. Use one of the concrete exported thermostat types when constructing an NBodySimulation.

No thermostat.

AndersenThermostat(T, ν)

Andersen thermostat targeting temperature T with collision frequency ν.

Arguments

• T: Target temperature.
• ν: Collision frequency.

Fields

• T: Target temperature.
• ν: Collision frequency.

Examples

using NBodySimulator

thermostat = AndersenThermostat(300.0, 1.0)

BerendsenThermostat(T, τ)

Berendsen thermostat targeting temperature T with coupling time τ.

Arguments

• T: Target temperature.
• τ: Temperature-coupling time.

Fields

• T: Target temperature.
• τ: Temperature-coupling time.
• γ: Rescaling coefficient derived from τ.

Examples

using NBodySimulator

thermostat = BerendsenThermostat(300.0, 0.1)

NoseHooverThermostat(T, τ)

Nose-Hoover thermostat targeting temperature T with relaxation time τ.

Arguments

• T: Target temperature.
• τ: Relaxation time.

Fields

• T: Target temperature.
• τ: Relaxation time.

Examples

using NBodySimulator

thermostat = NoseHooverThermostat(300.0, 0.1)

LangevinThermostat(T, γ)

Langevin thermostat targeting temperature T with friction coefficient γ.

Arguments

• T: Target temperature.
• γ: Friction coefficient.

Fields

• T: Target temperature.
• γ: Friction coefficient.

Examples

using NBodySimulator

thermostat = LangevinThermostat(300.0, 1.0)

CubicPeriodicBoundaryConditions(L)

Periodic cubic boundary condition with side length L.

Arguments

• L: Cubic-cell side length.

Fields

• L: Cubic-cell side length.

Examples

using NBodySimulator

boundary = CubicPeriodicBoundaryConditions(10.0)

PeriodicBoundaryConditions(boundary)

Periodic boundary conditions over explicit lower and upper bounds for each coordinate.

Arguments

• boundary: Six-element vector (xmin, xmax, ymin, ymax, zmin, zmax).
• L: A scalar side length, expanded to the interval [0, L] on each axis.

Fields

• boundary: Lower and upper bounds for the three Cartesian coordinates.

Examples

using NBodySimulator

boundary = PeriodicBoundaryConditions(10.0)

InfiniteBox()

Boundary condition representing an unbounded three-dimensional domain.

Fields

• boundary: Six infinite lower and upper bounds.

Examples

using NBodySimulator

boundary = InfiniteBox()

PotentialNBodySystem(bodies, potentials)
PotentialNBodySystem(bodies; potentials = Symbol[])

System of bodies governed by a selected set of potential parameters.

Arguments

• bodies: Bodies in the system.
• potentials: A Dict{Symbol, <:PotentialParameters} or a vector of built-in potential names.

Fields

• bodies: Bodies in the system.
• potentials: Potential parameters keyed by interaction name.

Examples

using NBodySimulator, StaticArrays

bodies = [MassBody(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.0)]
system = PotentialNBodySystem(bodies; potentials = [:gravitational])

ChargedParticles(bodies, k)

System of charged particles interacting through an electrostatic constant k.

Arguments

• bodies: Charged particles in the system.
• k: Coulomb constant.

Fields

• bodies: Charged particles.
• k: Coulomb constant.

Examples

using NBodySimulator, StaticArrays

body = ChargedParticle(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.0, 1.0)
system = ChargedParticles([body], 1.0)

GravitationalSystem(bodies, G)

System of massive particles interacting through gravitational constant G.

Arguments

• bodies: Massive particles in the system.
• G: Gravitational constant.

Fields

• bodies: Massive particles.
• G: Gravitational constant.

Examples

using NBodySimulator, StaticArrays

body = MassBody(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.0)
system = GravitationalSystem([body], 1.0)

WaterSPCFw(bodies, mH, mO, qH, qO, lj_parameters, e_parameters, scpfw_parameters)

System representation for SPC/Fw water molecules and their interaction parameters.

Arguments

• bodies: Oxygen sites or explicit water molecules.
• mH, mO: Hydrogen and oxygen masses.
• qH, qO: Hydrogen and oxygen charges.
• lj_parameters, e_parameters, scpfw_parameters: Interaction parameters.

Fields

• bodies: Water molecule definitions.
• mH, mO, qH, qO: Per-site physical constants.
• lj_parameters, e_parameters, scpfw_parameters: Interaction parameters.

Examples

using NBodySimulator, StaticArrays

bodies = [MassBody(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 18.0)]
system = WaterSPCFw(
    bodies, 1.0, 16.0, 0.41, -0.82, LennardJonesParameters(), ElectrostaticParameters(),
    SPCFwParameters(0.1012, 1.9764, 44315.0, 317.6)
)

SimulationResult should provide an interface for working with the properties of a separate particle and with the physical properties of the whole system.

get_position(result::SimulationResult, time, i = 0)

Return all particle positions at time, or the position of particle i when i > 0.

Arguments

• result: Result returned by run_simulation.
• time: Time at which to interpolate the solution.
• i: Particle index; use 0 to return all positions.

Examples

position = get_position(result, 0.5, 1)

get_velocity(result::SimulationResult, time, i = 0)

Return all particle velocities at time, or the velocity of particle i when i > 0.

Arguments

• result: Result returned by run_simulation.
• time: Time at which to interpolate the solution.
• i: Particle index; use 0 to return all velocities.

Examples

velocity = get_velocity(result, 0.5, 1)

get_masses(system::NBodySystem)

Return the particle masses for system in simulation-coordinate order.

Arguments

• system: N-body system whose masses are requested.

Examples

masses = get_masses(simulation.system)

temperature(result::SimulationResult, time)

Compute the system temperature at time.

Arguments

• result: Result returned by run_simulation.
• time: Time at which to interpolate the solution.

Examples

T = temperature(result, 0.5)

rdf(result::SimulationResult)

Compute the radial distribution function from a simulation result.

Arguments

• result: Result for a system with cubic periodic boundary conditions.

Returns

• (r, g): Radial-bin centers and radial-distribution values.

Examples

r, g = rdf(result)

msd(result::SimulationResult)

Compute the mean squared displacement over the saved solution times.

Arguments

• result: Result returned by run_simulation.

Returns

• (t, displacement): Saved times and mean squared displacements.

Examples

t, displacement = msd(result)

initial_energy(simulation::NBodySimulation)

Compute the total energy of the simulation initial condition.

Arguments

• simulation: Simulation whose initial state is evaluated.

Examples

energy = initial_energy(simulation)

kinetic_energy(velocities, masses)
kinetic_energy(result::SimulationResult, time)

Compute kinetic energy from velocity and mass arrays or from a simulation result at time.

Arguments

• velocities: Matrix with one velocity column per particle.
• masses: Particle masses corresponding to the velocity columns.
• result: Result returned by run_simulation.
• time: Time at which to interpolate a result.

Examples

energy = kinetic_energy(result, 0.5)

potential_energy(coordinates, simulation::NBodySimulation)
potential_energy(result::SimulationResult, time)

Compute potential energy from coordinates and simulation parameters or from a result at time.

Arguments

• coordinates: Position matrix with one column per particle.
• simulation: Simulation providing potentials and boundary conditions.
• result: Result returned by run_simulation.
• time: Time at which to interpolate a result.

Examples

energy = potential_energy(result, 0.5)

total_energy(result::SimulationResult, time)

Compute kinetic plus potential energy at time.

Arguments

• result: Result returned by run_simulation.
• time: Time at which to interpolate the solution.

Examples

energy = total_energy(result, 0.5)

Molecules for the SPC/Fw water model can be imported from a PDB file.

Arguments

• path_str: Path to a PDB trajectory written by save_to_pdb.

Examples

molecules = load_water_molecules_from_pdb("water.pdb")

save_to_pdb(result::SimulationResult, path)

Write a simulation result to a Protein Data Bank (PDB) file.

Arguments

• result::SimulationResult: Simulation result to serialize.
• path: File path passed to open(path, "w").

Returns

• Nothing: The file is written for its side effect.