Marcus N. Bannerman
School of Engineering, University of Aberdeen
m.campbellbannerman@abdn.ac.uk
I distribute an Event-Driven Molecular Dynamics package, called DynamO (dynamomd.org).
EDMD is a general molecular dynamic technique and it is used to study a wide range of fundamental model systems:
A binary hard-sphere fluid, used as an exploratory model for transport in nanofluids (more later).
A coarse-grained simulation of the binding of two heteropolymers. Event-driven simulations are popular in the protein folding/structure literature.
Asymmetric-dumbells or "Snowmen" molecules, used as an model for anisotropic colloids, or Na-Cl crystals, and have a rich phase behavior.
Although it is not as "popular" as time-stepping dynamics, event-driven dynamics is a natural approach to dynamics.
Consider how you would solve this problem:
Event driven dynamics has a similar approach.
Rather than step through time, we calculate the time of "events" exactly (using an approximate model).
Assuming no drag, the exact solution to the projectile motion is a parabola: \begin{align*} r_y(t+\Delta t) = \color{aqua}{r_y(t)} + \color{lime}{v_y(t)}\,\Delta t+g_y\,\frac{\Delta t^2}{2} \end{align*}
\begin{align*} r_y(t+\Delta t) = \color{aqua}{r_y(t)} + \color{lime}{v_y(t)}\,\Delta t+g_y\,\frac{\Delta t^2}{2} \end{align*} When it hits the ground at a time $t+\Delta t$ we have: \begin{align*} r_y(t+\Delta t)=r_{ground} \end{align*}
To find the time of this "event", we are looking for the roots to the equation: \begin{align*} r_y - r_{ground} + v_y\,\Delta t+g_y\,\frac{\Delta t^2}{2} = 0 \end{align*}
In this case, this is a simple quadratic with two roots: \begin{align*} \Delta t = \frac{-v_y \pm \sqrt{v_y^2 - 4\,g_y(r_y-r_{ground})}}{2\,g_y} \end{align*}
Event-detection is driven by root finding.
Consider two dueling catapults:
We calculate $\Delta t$ for all possible events: \begin{align*} \vec{\Delta t}=\left\{\Delta t_{1-ground},\,\Delta t_{2-ground},\,\Delta t_{1-2}\right\} \end{align*} The earliest event is the next event!
Event processing is driven by sorting.
What actually happens at an event?
Consider a two-body event between particle $i$ and particle $j$, which has an energy change of $\Delta E$.
Conservation of momentum requires: \begin{align*} m_i\,\Delta \vec{v}_i = - m_i\,\Delta \vec{v}_i = \Delta \vec{P} \end{align*}
Where the impulse, $\Delta \vec{P}$, is a solution to
energy balance:
\begin{align*}\scriptsize
\frac{1}{2}m_i\, {v}_i^2 + \frac{1}{2}m_j\, v_j^2 + \Delta E = \frac{1}{2}m_i\, \left(\vec{v}_i +\frac{\Delta \vec{P}}{m_i}\right)^2 + \frac{1}{2}m_j\, \left(\vec{v}_j-\frac{\Delta \vec{P}}{m_i}\right)^2
\end{align*}
Event execution is driven by root finding for an impulse.
Note, we do not decide the time of the simulation! We decide how many events to run it for and the time is calculated.
Modern EDMD algorithms are extremely complex and include:
Perhaps the most unusual optimization is the time warp algorithm.
As events only involve small, localized particles, the system can run out of sync. i.e., each particle is left at different points in time:
Initially, the simulation is run in synchronous mode. The time-warp algorithm is then enabled (with the same simulation speed) and only the required updates to the particles are shown.
What is unique to DynamO is its stability.
EDMD is an exact (to machine precision) solution of the dynamics of the model used.
All approximation is explicit in the model used.
(no worries about the time step)
The high-precision of the method conserves detailed
balance.
This is ideal for mixed MC/MD techniques
(e.g., replica exchange, multicanonical).
However, even though the method is analytically exact, there are problems with stability!^{1}
Consider repeated bounces of the projectile: Using a constant coefficient of restitution $e<1$ we have: \begin{align*} \Delta P_{i,y} = -m_i\,(1+e) v_{i,y} \end{align*}
The velocity on impact/event decays to zero, and an infinite number of events happens in a finite time.
This "inelastic collapse" is not an instability, but a phenomena of the model.
Trying to simulate an "inelastic collapse" in a naive way highlights a numerical issue.
Although $r_y$ should equal $r_{ground}$ after an event, limitations of the precision of the computer leave it at $\approx \pm 10^{-14}$.
As $v_y\to0$, numerical precision is lost and the argument of the square-root can turn negative! \begin{align*} \Delta t = \frac{-v_y \pm \sqrt{v_y^2 - 4\,g_y(r_y-r_{ground})}}{2\,g_y} \end{align*}
This gives no real-roots/events and the projectile falls through the ground!
Although this problem is acute in dissipative systems with gravity, it also plagues elastic/molecular systems for certain sequences of events:
Machine precision is not enough to stop particles entering "forbidden" states, which may get rapidly worse.
DynamO has a (small) compile-time computer algebra library and event searching routines:
Variable t;
double deltat = nextroot(rij + t * (vij + aij * t / 2), t);
Additional "stabilizing" events are added if particles are in invalid conditions and this situation is deteriorating.
The library can find all roots of $N$-order polynomials and certain classes of other functions in a stable way, allowing a wide range of models/events to be detected.
Can we simulate "normal" molecular models using EDMD?
Between discontinuities/events/impulses there is no
force:
Exact solution to the motion is possible!
These stepped approximations can closely reproduce the phase behavior of the LJ system:
And the dynamics, for example, the shear viscosity:
But we should only use this approach for gases as for liquids it is not as efficient as standard approaches (as expected).
Transmission electron micrograph of Cu nanoparticles produced by direct evaporation into ethylene glycol, (Eastman et. al. 2001).
Nanofluids are suspensions of nano-sized particles (solid or liquid) in a base fluid.
They are used to attempt to "alloy" the properties of two phases without the macroscopic problems (abrasion, stability).
In particular, metallic nanoparticles can be used to attempt to increase the thermal conductivity of water.
Thermal conductivity enhancement of a ethylene glycol base fluid with $~10$~nm Copper nanoparticles.Taken from Eastman et. al. 2001.
Early attempts to explore the thermal properties of nanofluids were extremely encouraging.
Dispersion of copper nanoparticles in ethylene glycol exhibited enormous thermal conductivity enhancements for low volume fractions of nanoparticle.
These enhancements were well beyond initial predictions, earning the title of an "anomalous" enhancement.
A massive effort, from 30 different organizations, tested identical nanofluid samples with a variety of techniques as part of the International Nanofluid Property Benchmark Exercise.
Finally, in 2010, Eapen et. al. published a paper which highlighted that most results lay within the classical limits of series and parallel conductors (right).
Eapen concluded that there are no "anomalous" results.
This system is an ideal target for EDMD approaches.
The hard-sphere fluid model is sufficiently complex to capture fluid mixing effects (size, mass, heat capacity asymmetry), while being very efficient.
Large asymmetries require large systems, $N>10^5$.
Long simulation times are also required to extract equilibrium transport properties.
Enskog kinetic theory for hard spheres also allows theoretical exploration and validation of the fluid properties.
An equilibrium simulation of a hard-sphere nanofluid at a size ratio of 10:1 and mass ratio of 1000:1.
Thermal conductivity, $L_{\lambda\lambda}$, of a 10:1 nanoparticle fluid as a function of nanoparticle mole fraction, $x_A$, at (a) fixed packing fractions, $\phi$, or (b) fixed pressures, $p$. Symbols are simulation results, lines are predictions from kinetic theory.^{1}
Analyzing the kinetic theory predictions, the dominant contribution has this form: \begin{align*} J_q &\approx \left(C_v^\text{Nano} - C_v^\text{Base}\right)k_B\,T\, L_{A\lambda}\frac{\nabla T}{T} \end{align*}
\begin{align*} J_q &\approx \left(C_v^\text{Nano} - C_v^\text{Base}\right)k_B\,T\, L_{A\lambda}\frac{\nabla T}{T} \end{align*}
A non-equilibrium thermal conductivity simulation of a hard-sphere nanofluid at a size ratio of 10:1 and mass ratio of 1000:1. Thermophoresis causes the larger species to diffuse to the hot plate.
Isotopic 10:1 mass ratio non-equilibrium molecular dynamics measurements of the thermal conductivity. For this system, the equilibrium thermal conductivity is $L_{\lambda\lambda}=2.35$. Packing fraction is 0.1, mole fraction of the heavier component is 0.1.