sciml / nbodysimulator.jl Goto Github PK

@JuliaRegistrator register()

I'm a bit confused - I tried running the example code for gravitational interaction and it didn't generate the expected animation.

This continues a conversation that started here:

SciML/MultiScaleArrays.jl#25 (comment)

HierarchicalMatrices.jl gives a convenient way to implement the fast multipole method, which exploits low rank structure to make matrix multiplication fast. This can be used for efficient simulations: O(n) cost for n bodies.

The current implementation for calculating the potential energy of the system is hardcoded to only recognize a couple of built-in potential functions. NBS supports adding custom potentials for the simulation, but they can't be included in analysis of the energies afterwards. Could this be extended to custom potentials through dynamic dispatch?

Hello and thanks for this project,

is there a reason why simulations are only possible in 3 dimensions?

Such that makedocs warnonly = [:docs_block, :missing_docs] can be removed.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

In a simulation with CubicBoundaryConditions I am seeing bodies go outside of the box after a large number of timesteps. I would expect that get_positions always return coordinates within the range [0, L]

Example script:

using NBodySimulator
using Unitful
using UnitfulAtomic

m = 6.6335209e-26u"kg"
σ = 0.34u"nm"
cutoff = 0.765u"nm"
ϵ = 1.657e-21u"J"

argon_initial_bodies(N::Integer, box_size::Quantity, reference_temp::Quantity) = generate_bodies_in_cell_nodes(N, austrip(m), austrip(√(u"k" * reference_temp / m)), austrip(box_size))
argon_lennard_jones() = Dict(:lennard_jones => LennardJonesParameters(austrip(ϵ), austrip(σ), austrip(cutoff)))
argon_equilibration_thermostat(reference_temp::Quantity, thermostat_prob::AbstractFloat, Δt::Quantity) = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))

N = 8
box_size = 20u"bohr"
reference_temp = 94.4u"K"
thermostat_prob = 0.1
t₀ = 0u"ps"
Δt = 1e-2u"ps"
steps = 20000

initial_bodies = argon_initial_bodies(N, box_size, reference_temp)
potentials = argon_lennard_jones()
system = PotentialNBodySystem(initial_bodies, potentials)
boundary_conditions = CubicPeriodicBoundaryConditions(austrip(box_size))
thermostat = argon_equilibration_thermostat(reference_temp, thermostat_prob, Δt)
simulation = NBodySimulation(system, (austrip(t₀), austrip(t₀ + steps * Δt)), boundary_conditions, thermostat, 1.0)
simulator = VelocityVerlet()

result = run_simulation(simulation, simulator, dt=austrip(Δt))

positions = get_position(result, result.solution.t[end])

if isa(boundary_conditions, CubicPeriodicBoundaryConditions)
    bounded_positions = mod.(positions, boundary_conditions.L)
end

display(positions)
display(bounded_positions

Here we have a box with side length 20, but the resulting positions are the following:

3×8 view(::Matrix{Float64}, :, 1:8) with eltype Float64:
  0.623155  -11.0268   12.3628   -11.9184   25.4363   12.8586  14.1605    2.80354
  2.86595     7.17172   7.84181   -7.04046  12.3267   13.0511  28.6849   23.7666
 17.3475     12.1708   -1.43481   19.108     9.36417  -7.2626   6.38268   6.47773

If I take the coordinates modulo the box size, I get

3×8 Matrix{Float64}:
  0.623155   8.97317  12.3628    8.08165   5.43627  12.8586  14.1605   2.80354
  2.86595    7.17172   7.84181  12.9595   12.3267   13.0511   8.68491  3.76665
 17.3475    12.1708   18.5652   19.108     9.36417  12.7374   6.38268  6.47773

Perhaps this is the right interpretation? Or maybe there is a deeper issue.

An issue for storing graphical stuff.

Warning: Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).

It would be nice, if there were an option to store the forces acting on the bodies separated by interaction type and a function similar to get_positions and get_velocities that allows to analyze the forces when the simulation has finished.

Seems to be a typo in:

i.e. paramteres should be parameters
This doesn't seem to impact the functionality of the code.

@dextorious , have a look at this figures for SI and OpenMM units usage. I tested simulation of 216 particles of liquid argon for 15 ps with timestep at 0.5 fs without thermostating.

I didn't try to set the same initial velocities for both tests, though it is worth to do so for an unbiased experiment/test.

Can you suggest other test systems with already implemented features for comparison calculations in different units? Berendsen thermostat can be used, for example, because it is deterministic.
Well, at least, I can repeat such calculations for water systems.

I get a segfault when getting results out using VelocityVerlet, in a getindex. I'm guessing overagressive inbounds somewhere? I'm too lazy to do an MWE but the below should reproduce the bug. Changing the integrator to something else fixes it. I'm filing against NBS since that's what I'm using, but I can file to DifferentialEquations.jl if that is more appropriate.

using SPICE
using Downloads: download
using NBodySimulator, StaticArrays

# units: kilometers, seconds, kilograms

const LSK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/lsk/naif0012.tls"
const SPK = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp"
# Download kernels
# time kernel
download(LSK, "naif0012.tls")
download(SPK, "de440.bsp")
download("https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/Gravity.tpc", "Gravity.tpc")

furnsh("naif0012.tls")
# ephemeris kernel
furnsh("de440.bsp")
# masses kernel
furnsh("Gravity.tpc")

# Convert the calendar date to ephemeris seconds past J2000
start_date = str2et("Jan 1, 2020")

years_to_observe = 1
dates = collect(range(start_date, start_date + 60*60*24*365*years_to_observe, length=365*years_to_observe))

ref_pos = spkpos.("earth", dates, "J2000", "none", "sun")
ref_pos = [p[1] for p in ref_pos]
ref_pos = reduce(hcat, ref_pos)'

const grav_constant = 6.6743e-11 / 1e9 # m^3 kg^-1 s^-2 -> km^3

planets = ["earth_barycenter", "sun"]
# planets = ["earth_barycenter", "Mercury_barycenter", "Venus_barycenter", "Mars_barycenter", "Jupiter_barycenter", "Saturn_barycenter", "Uranus_barycenter", "Neptune_barycenter", "Sun"]
# planets = ("earth_barycenter", "Mercury_barycenter", "Venus_barycenter", "Mars_barycenter", "Jupiter_barycenter", "Saturn_barycenter", "Uranus_barycenter", "Sun")
planets_to_observe = [1]

masses = bodvrd.(planets, "GM") ./ grav_constant
posvel = [spkezr(p, start_date, "J2000", "none", "sun") for p in planets]

bodies = [MassBody(SVector(posvel[i][1][1:3]...), SVector(posvel[i][1][4:6]...), masses[i][1]) for i = 1:length(planets)]
system = GravitationalSystem(bodies, grav_constant)
tspan = (dates[1], dates[end])
simulation = NBodySimulation(system, tspan)
result = run_simulation(simulation, VelocityVerlet(), dt=60*60*24).solution
pos = [result(t)[:, planets_to_observe] .- result(t)[:, end] for t in dates] # observe wrt sun
pos = cat(pos...; dims=3) # α x p x t

function get_ref_pos(t)
    hcat([spkpos(p, t, "J2000", "none", "sun")[1] for p in planets[planets_to_observe]]...)
end
ref_pos = cat(get_ref_pos.(dates)...; dims=3)

println(maximum(abs, pos-ref_pos))

# using PyPlot
# plot3D(pos[1, 1, :], pos[2, 1, :], pos[3, 1, :])

The current run_simulation interface is used to create and solve the appropriate ODE/SDE problem based on the thermostat algorithm.

The current implementation thus only allows only some some dynamical ODE solvers to be used with run_simulation since any other would use the ODEProblem fallback and fail.

One way to solve this would be via a trait, as indicated in the comment. An other option would be to select only some predefined well performing integrators in the union and indicate to users that they should manually construct and solve the problem for other algorithms.
And a simpler option would be to just use SecondOrderODEProblems for everything and remove the ODEProblem path.

I can start a PR if we decide upon a solution.

When performing a simulation with a large number of timesteps and a non-integral dt, I sometimes get one more tilmestep than expected. The symptoms vary with the exact value of dt, so my first guess is that there is a floating point comparison somewhere that is slightly off. Specifically, if the simulation interval is (0, n * dt) the time series usually has length n + 1 but for some combinations of inputs has length n + 2.

Example script:

using NBodySimulator
using Unitful
using UnitfulAtomic

m = 6.6335209e-26u"kg"
σ = 0.34u"nm"
cutoff = 0.765u"nm"
ϵ = 1.657e-21u"J"

argon_initial_bodies(N::Integer, box_size::Quantity, reference_temp::Quantity) = generate_bodies_in_cell_nodes(N, austrip(m), austrip(√(u"k" * reference_temp / m)), austrip(box_size))
argon_lennard_jones() = Dict(:lennard_jones => LennardJonesParameters(austrip(ϵ), austrip(σ), austrip(cutoff)))
argon_equilibration_thermostat(reference_temp::Quantity, thermostat_prob::AbstractFloat, Δt::AbstractFloat) = AndersenThermostat(austrip(reference_temp), thermostat_prob / Δt)

N = 8
box_size = 20u"bohr"
reference_temp = 94.4u"K"
thermostat_prob = 0.1
t₀ = 0.0
Δt = 413.41373335335163

initial_bodies = argon_initial_bodies(N, box_size, reference_temp)
potentials = argon_lennard_jones()
system = PotentialNBodySystem(initial_bodies, potentials)
boundary_conditions = CubicPeriodicBoundaryConditions(austrip(box_size))
thermostat = argon_equilibration_thermostat(reference_temp, thermostat_prob, Δt)
simulator = VelocityVerlet()

for steps ∈ (10000, 20000, 30000)
    simulation = NBodySimulation(system, (t₀, t₀ + steps * Δt), boundary_conditions, thermostat, 1.0)

    result = run_simulation(simulation, simulator, dt=Δt)

    display(length(result.solution.t))
    # expected output: 10001, 20001, 30001
    # actual output: 10001, 20002, 30001
end

Note that changing Δt to 413.4137 for example causes the output to be flipped (10002, 20001, 30002).

Hi,

I was wondering if it is possible to combine NBodySimulator.jl with DiffEqFlux.jl to learn interaction parameters to reproduce a given RDF. loss is the distance euclidean distance between ground truth and current parameter RDFs to make it differentiable. I have done something similar for stochastic lotka_volterra.

Calculation of RDF is done through a function like this, where x's are radial distances of neighboring atoms.

rkMin = 0 # Ang
rkMax = 10 #Ang
nrk = 80
rk = LinRange(rkMin, rkMax, nrk)

function kernelHist(x, rk=rk)
dr = rk[2] - rk[1]
#x = reshape(x, (1, size(x)[2]size(x)[1]))
dx2 = .-(rk .- x).^2 ./ (2dr^2)
denUn = exp.(dx2 ./ (2dr^2))
denN = denUn ./ sum(denUn .+ 1e-16, dims=1)
totalDenUn = sum(denN, dims = 2)
totalDenUn ./= sum(totalDenUn . dr)
totalDenUn = reshape(totalDenUn, size(totalDenUn)[1])
#totalDenUn
totalDenUn
end

Hi!

(v1.0) pkg> ^C

julia> using NBodySimulator
[ Info: Precompiling NBodySimulator [0e6f8da7-a7fc-5c8b-a220-74e902c310f9]
ERROR: LoadError: LoadError: LoadError: UndefVarError: start not defined
Stacktrace:
 [1] top-level scope at none:0
 [2] include at ./boot.jl:317 [inlined]
 [3] include_relative(::Module, ::String) at ./loading.jl:1041
 [4] include at ./sysimg.jl:29 [inlined]
 [5] include(::String) at /home/username/.julia/packages/NBodySimulator/U07Nf/src/NBodySimulator.jl:3
 [6] top-level scope at none:0
 [7] include at ./boot.jl:317 [inlined]
 [8] include_relative(::Module, ::String) at ./loading.jl:1041
 [9] include at ./sysimg.jl:29 [inlined]
 [10] include(::String) at /home/username/.julia/packages/NBodySimulator/U07Nf/src/NBodySimulator.jl:3
 [11] top-level scope at none:0
 [12] top-level scope at none:2
in expression starting at /home/username/.julia/packages/NBodySimulator/U07Nf/src/boundary_conditions.jl:10
in expression starting at /home/username/.julia/packages/NBodySimulator/U07Nf/src/nbody_simulation.jl:4
in expression starting at /home/username/.julia/packages/NBodySimulator/U07Nf/src/NBodySimulator.jl:10
ERROR: Failed to precompile NBodySimulator [0e6f8da7-a7fc-5c8b-a220-74e902c310f9] to /home/username/.julia/compiled/v1.0/NBodySimulator/TbJpf.ji.

Thanks!

I was trying a simple model of the solar system and I wanted to simulate that for a thousand years. However the NBodySimulator returns the error:

interrupted. Larger maxiters is needed.

It would be nice if there was an option to change the default maximum amount of iterations to a higher value.

Below is the code I used.

using StaticArrays, Plots, NBodySimulator
bodies = [
    MassBody(SVector(5.9133491e10, 0.0 ,0.0), SVector(0.0, 4.74e4 , 0.0), 3.30104e23),
    MassBody(SVector(1.0821141e11, 0.0, 0.0), SVector(0.0, 3.50e4, 0.0), 4.86732e24),
    MassBody(SVector(1.49618773e11, 0.0, 0.0), SVector(0.0, 2.98e4, 0.0), 5.9721986e24),
    MassBody(SVector(2.28931109e11, 0.0, 0.0), SVector(0.0, 2.41e4, 0.0), 6.41693e23),
    MassBody(SVector(7.79323489e11, 0.0, 0.0), SVector(0.0, 1.30e4, 0.0), 1.89813e27),
    MassBody(SVector(1.42881720e12, 0.0, 0.0), SVector(0.0, 9.64e3, 0.0), 5.68319e26),
    MassBody(SVector(2.874165879e12, 0.0, 0.0), SVector(0.0, 6.80e3, 0.0), 8.68103e25),
    MassBody(SVector(4.498418710e12, 0.0, 0.0), SVector(0.0, 5.43e3, 0.0), 1.02410e26),    
    MassBody(SVector(0.0, 0.0, 0.0), SVector(0.0, 0.0, 0.0), 1.989e30)
]
const G = 6.673e-11  # m^3/kg/s^2
const timespan = (0.0, 3.1e10)
 
function nbodysim(nbodies, tspan)
    system = GravitationalSystem(nbodies, G)
    simulation = NBodySimulation(system, tspan)
    run_simulation(simulation)
end
 
simresult = nbodysim(bodies, timespan)