Comments on N-body problem

[ChrisRackauckas]

@Mikhail-Vaganov some comments on the proposal someone gave me:

Periodic Boundary Condition (PBC)
Distance

The student says that for the Periodic Boundary condition appropriate particles copies will be generated when it leaves the periodic box. While a nice feature for analysis this does not the defining feature of a simulation with PBC. Rather he has to ensure that any calculated distance follows the minimum image convention (MIC). The MIC says that from any pair i,j in the simulation box and it’s infinite identical copies find to use the shortest distance. See wikipedia. For performance I recommend not using code that uses floor to calculate the distance. Rather use a while loop.

while (x > box_x * .5) x -= box_x;
while (x < -box_x * .5) x += box_x;

The while loop will normally have only one invocation unless the particle jumps several images.

Electro / Magneto Statitics

For PBC simulations electro statics and other long range forces have to be treated specially. As correctly identified by the student one solution is Ewald summation. In simulation without Ewald summation or special screening terms to limit these forces to short-ranges the simulations will give physical wrong results. For PBC conditions if Ewald summation is a stretch goal so should be electro statics.

Thermostats
Nose Hoover Equation
The normal Nose Hoover thermostat does not reproduce the correct equilibrium distribution for a harmonic oscillator ([1] 4.10). One would instead need to use Nose Hoover chains ([1] 4.11). They are more complicated to implement and cannot be converted to a Hamiltonian system

[1] Statistical Mechanic Theory and Simulation
Anderson/Berendson Thermostat

Both only produce approximate results for very large systems and are typically used to equilibrate a system and not for production simulations.
Velocity Rescaling
This algorithm has become a standard thermostat in the MD community. It extends the ideas from Berendson and Anderson by coupling the distribution for selection of random velocities to a stochastic process. This thermostat tends to be relatively easy to implement and gives good results also for small systems. I recommend using it.

[1] Bussi et. al 2007
Langevin Thermostat

This thermostat is interesting for coarse grained simulations where I don’t want to simulate the solvent explicitly. But the random forces in the integrator violate detailed balance, eq 1, leading to different equilibrium distributions.

P(ij)ij=P(ji)ji (1)

A simple example when detailed balance is violated is using a linear increasing force F(j) > F(i) if j>i. Due to the force and the position independent random-force the probability to make a transition from i to j is not equal to the transition from j to i. The only way to minimise this error is to choose a sufficiently small timestep. This can be tested with a harmonic oscillator.

Another problem is the quick divergence of Langevin integrators for stiff potentials. For such a potential the timestep has to be chosen sufficiently small for the simulation not to diverge. Adaptive step size algorithms can help to still achieve good performance for stiff potentials.

An alternative are dynamic Monte-Carlo methods [1]-[3]. A large advantage of MC methods is that you do not need to calculate forces given a potential.

[1] Kotelyanskii et. al 1992
[2] Heyes et. al 1998
[3] Sanz et. al 2010

Particle Collision

Using continuous callbacks to simulate particle collisions sounds similar to a hard sphere simulation. Those can be solved explicitly with a variable timestep [1] Chap 14.2.

For flexibility I would leave this to potential function like the lennard jones potential.

[1] The Art Molecular Dynamics Simulation

Efficient Force/Potential Evaluations

A big problem of MD simulations is that a naiive approach scales as O(N^2) with N the number of particles. This is terrible scaling! We are lucky though that most intersting potential are only short ranged. This allows us to use a domain decomposition algorithm (also known as verlet lists). The idea is that we only need to calculate forces/energies up to a cutoff Rc and therefore only need to calculate distances for particles which are closer than Rc. This can be done by sub-dividing the simulation box into cubes of edge length Rc. Only particles in neighboring cells need to be considered for force/energy calculations. This results in a scaling of O(N). See [1] for an implementation of domain decomposition in python.

[1] CellGrid
Stretch Goals
Other possible interesting stretch goals

Multiple species with different forcefield parameters. (i.e. larger radius)
NPT ensemble simulations

Similar projects

Besides lammps and espresso other big MD engines are GROMACS, OpenMM, Amber, NAMD. The later are specifically optimized for proteins and other biological applications.

MDAnalysis is thinking about taking a student for efficient distance calculations. It might be worth it that both students share experiences for distance calculations in particular.

[Mikhail-Vaganov]

Thanks! Very good comments! I will surely take them into account and change some parts of my project. But before that, I would rethink the comments more thoroughly...

[dextorious]

Just a quick remark on those earlier comments - the Andersen and Berendsen thermostats are in fact fundamentally different in that the Andersen thermostat is local (which means, among other things, that it destroys momentum correlations and is therefore unsuitable for any simulation interested in transport properties - diffusion coefficients, viscosities, etc.), while the Berendsen thermostat (despite having other disadvantages like not quite corresponding to the NVT ensemble) is global and therefore of more general applicability. The main point here is that they both have radically different uses.

More generally, Nose-Hoover variants seem to be the standard at the moment for any simulations involving transport properties and other non-equilibrium phenomena, whereas Langevin type integrators (including various geodesic derivatives) are popular in biological applications where they can increase the timestep by an order of magnitude or so.