I have written a notebook on particle simulations using Julia, to be presented at the FortranCon this week:
https://m3g.github.io/2021_FortranCon/
(https://github.com/m3g/2021_FortranCon)
It may be interesting for students learning Julia, or researchers exploring the possibility of writing particle simulation codes.
Any suggestion is of course always appreciated.
Nice, I like this tour of various cool things you can do in Julia.
A few things I noticed reading through:
Point2D for those seems a little strange. Since you already depend on StaticArrays, you might use const Vec2D{T} = SVector{2,T}, or just rename your Point2D struct to Vec2D if you want the x and y fields.random_point would be neater using Distributions.jl and passing the distribution to sample from rather than range. The existing code is good though because it’s self contained; maybe a “further reading” for students which points to Distributions.jl could be helpful.md simulation function is a very simple integrator. Again, this is good because it’s a self-contained tutorial, but much better integrators exist in DifferentialEquations.jl and a note / link to that for further reading could be helpful. Writing integrators is good for learning how things work, but reusing the existing ecosystem of high quality integrators should be encouraged.
Presumably it’s because the md integrator is a low order scheme and not energy conserving. IIRC gravity simulations usually use adaptive high order integrators and the Julia ecosystem has some great ones.We can differentiate everything!
Be careful with differentiating chaotic systems. This would be a good reason to use the DifferentialEquations.jl adjoints. See:
Suggestion taken. I changed to Vec2D. I like to subtype FieldVector because then we can access the fields by keys (i. e. vec.x, vec.y). That makes the plotting codes prettier, in particular.
Indeed, the idea is to not rely to much on packages for simple things, because otherwise the reader may get the impression that Julia is one of many interpreted languages. But, this, in particular, I do not think that Distrubutions help much: the coordinates are just randomly generated in square, thus using something sophisticated to that would seem strange. The velocities could use a Maxwell-Boltzmann distribution, but the last time I checked that was not implemented in Distributions yet (it was a pull request if I remember). Anyway generating that distribution is also simple with randn.
I like the one that is there because it reminds students from the equations they learn in high school for uniformly accelerated motion. They get the concept without distraction. The actual integrator used in most (if not all) molecular dynamics simulation packages is Velocity-Verlet, which in terms of implementation is almost as simple. There has been a lot of research in testing different integration schemes for MD simulations, but AFAIK nothing more sophisticated was worth the cost at the end. What I implement in the true benchmark against NAMD is Velocity-Verlet. It is much more stable than the naive integrator, of course. I think if someone wants to test different integration schemes in that code, that can actually be a research project.
I really didn’t test things to see where the problem is, but very likely it is just that I have taken the average position and velocities from the table, and in reality the distance to the sun and velocity are not constant, so that is not necessarily a realistic starting point. I’ve tested smaller time-steps, but that does not improve the trajectory much.
Yes, of course, both things have to be discussed carefully is a more profound materials on this subject. Note that the simple simulation I have differentiated there (Earth-Sun only) is not chaotic. But I will link that great post Chris mentioned there.
I am adding a final “Remarks” section including general comments concerting all these suggestions.
Thank you both.
Thanks again for the comments. I’ve added some remarks concerning the suggestions, and also fixed some bugs on the code.
Maybe not exactly on-topic for this thread, but: have you tried cell lists with cell size smaller than r_{cut} (e.g. https://arxiv.org/pdf/physics/0311055.pdf)? I had some tests suggesting that there is some improvement. I think LAMMPS uses that approach because of this comment in code. I had some tests for RDF calculation with cell sizes r_{cut}, r_{cut}/2, r_{cut}/3, r_{cut}/4 for LJ with density 0.7 and r_{cut} = 5.0 (IIRC), and r_{cut}/3 gave the best result. But I wasn’t applying other optimization strategies to reduce the number of distance calculations you have.
Btw, reaching NAMD performance is really cool. Even if it is general-purpose etc., the LJ part is probably thoroughly polished to get the best performance in benchmarks.
That can be tuned with a lcell parameter in my implementation. And it does make a difference for very dense systems. For molecular (water-like) systems generally lcell=2 is the best choice. But here that parameter does not make the same difference than in other codes, because I have implemented (sort of) this method, that already reduces significantly the number of unnecessary distance computations.
It looks like you have r^3 in the denominator of your gravitational_force function instead of r^2, could that be why some of the trajectories are off?
No, that is because there is the difference in positions in the numerator, not normalized, so there is one r in plus.
Right, interesting. It makes sense that low order symplectic methods would be the way to go in MD especially if you want to study large scale collective or statistical behavior. The classical equations of motion are only an approximation anyway.
Well planetary orbits should always be closed ellipses in the two-body problem, even though the distance to the sun is not constant over the orbit. It seemed by eye that the orbital radius always increased throughout the simulation which is why I thought it could be an artifact of the integrator. (IIUC Velocity-Verlet is symplectic so should be much better behaved in this respect due to being energy conserving.) Maybe not a big problem for a tutorial but something that immediately caught my eye.
Yes, you are in part correct: most of the error comes from the fact that the period is wrong, but that shift of the initial coordinate was an integration problem. Now I did that more carefully:
The optimization of the initial position effectively fixed the period, but both start and end at the same (for each) positions.
Oh nice, that looks much better!
IIRC, gravity (or in general a \frac{-\gamma}{r} potential ) and harmonic motion are the only two potentials which yield closed orbits regardless of initial conditions. If the orbits aren’t closed there’s something concerning going on.
I wasn’t concerned about the period as that could easily be due to initial conditions and I thought the use of optimization to find something with the correct period was neat.
That’s not quite true BTW. If you benchmark it Verlets are not the most efficient.
https://benchmarks.sciml.ai/html/NBodySimulator/liquid_argon.html
Unless I am mistaken that implementation of Verlets was non-optimal though because it calculated the force twice per time step, i.e. before https://github.com/SciML/OrdinaryDiffEq.jl/pull/1390. The other Verlet variants look like they suffered from this too.
This does need to be re-ran, it will make a dent, but since this is in log scale by eye it doesn’t look it would fully reach McAte5 in some cases for example.
There is something to be said that the definition of order for chaotic N-body problems probably isn’t the best definition though. I think something that converges in terms of ergodic properties faster would be a nice research project in ODE solvers that no one has looked into yet. Numerical ODEs is still such a young field with so much easy wins to explore…
They report the number of force evaluations and, clearly, they are doing 2 evaluations per step, which is not what one would expect. Note that the evaluation of the force is by far the most expensive part of the calculation, so the integrator should really improve on the length of the time step to be able to compensate additional evaluations.
I think exploring these possibilities is very interesting. Certainly people have tried alternatives before and for some reason, which is not my specialty, the current codes all converged to Velocity-Verlet. But maybe this is a consensus that can be revised.
There is something to be said that the definition of order for chaotic N-body problems probably isn’t the best definition though. I think something that converges in terms of ergodic properties faster would be a nice research project in ODE solvers that no one has looked into yet.
I have a preliminary project outline that as I understand is very closely related to the topics discussed here. Despite it being directly not my area of expertise, I have familiarized myself with your lectures @ChrisRackauckas on Parallel Computing and Scientific Machine Learning and I expressed positive words several times in this forum wrt the package developed by @lmiq. I would describe the project as advanced. Should it be potentially to your concern, I would be happy to provide you with a short initial description. Not sure how this forum technically works, if it is possible to send private messages or not, however, should it be to your potential interest and you might have any suggestions please let me know.
I think there are a few reasons:
I agree it’s an interesting area of research. I am surprised there isn’t more overlap between the molecular dynamics and the differential equations communities - maybe Julia can help here.
What’s the tl;Dr from all of those benchmarks?
I would expect about a 25% improvement at most, and what this also shows is that the effect of order is not as significant as the effect of coefficient optimization, which means a different definition of order (one based on ergodic properties) would likely lead to better results.