CCP5 Summer School 2015

Research presentation

Marcus N. Bannerman
School of Engineering, University of Aberdeen

Your browser does not support SVG


Marischal College
Aberdeen city center

Aberdeen City

  • Known as the granite city due to its construction.
  • Oil capital of Europe, producing 28% of Scottish GDP with 4% of the population.
  • Unemployment of 0.8% (June 2015), and highest rate of millionaires in the UK (2013).
  • Statistically the coldest city in the UK but still drier than Manchester(!) with 18 hours of daylight in summer (but only 6 in winter).

New Kings
University of Aberdeen

University of Aberdeen

  • Ancient university, established in 1495 (fifth in the UK).
  • 5 Nobel prize winners, including insulin treatment, crystalline electron diffraction, and partition chromatography.
  • Currently has around 16,000 students.
  • A general school of engineering with over 1000 UG students.
  • James Clerk Maxwell was here, but to our infinite shame we fired him in 1860.

Event-Driven Molecular Dynamics

I distribute an Event-Driven Molecular Dynamics package, called DynamO (

Your browser does not support SVG

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.

What is event-driven dynamics?

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:

Q: You are driving from Manchester to Aberdeen, the distance is 360 miles, how long does your journey take?

Event driven dynamics has a similar approach.

Rather than step through time, we calculate the time of "events" exactly (using an approximate model).

Your browser does not support SVG

How are event times calculated?

\definecolor{teal}{RGB}{128,255,128} \definecolor{aqua}{RGB}{128,255,64} \definecolor{lime}{RGB}{0,255,0}

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.

Generalizing for multiple particles

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.

Event execution

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.

The EDMD algorithm

  1. Calculate the times of all possible future events, $\vec{\Delta t}$.
  2. Sort to determine the time of the next event: $\Delta t_{next} = \min\left(\vec{\Delta t}\right)$.
  3. Move the system up to the time of the next event: $t\to t+\Delta t_{next}$
  4. Calculate the impulses, $\Delta \vec{P}$, for all bodies involved and apply.
  5. Repeat.

Note, we do not decide the time of the simulation! We decide how many events to run it for and the time is calculated.

The EDMD algorithm

Modern EDMD algorithms are extremely complex and include:

  • Neighbor lists and particle region tracking.
  • Bounded calendar priority queues for sorting.
  • Lazy deletion methods for optimal event list updates.

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).

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.

Stepped Potentials

Can we simulate "normal" molecular models using EDMD?

Between discontinuities/events/impulses there is no force:
Exact solution to the motion is possible!

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

These stepped approximations can closely reproduce the phase behavior of the LJ system:

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

And the dynamics, for example, the shear viscosity:

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

But we should only use this approach for gases as for liquids it is not as efficient as standard approaches (as expected).

Thomson, C., Lue, L. & Bannerman, MN. . 'Mapping continuous potentials to discrete forms'. J. Chem. Phys, 140, 034105 (2014)

  • EDMD and standard MD are equivalent.
  • EDMD performs exponentially better for low/gas density systems, but is slower for liquids.
  • EDMD is more efficient for:
    • Stiff systems (hard sphere, WCA, spring).
    • Coarse-grained potentials (hard sphere, square-well).
    • Multi-timescale systems (granular systems).

Thermal conductivity of nanofluids

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 of nanofluids

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.

Thermal conductivity of nanofluids

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

Bannerman, M. N. & Lue, L., "Transport properties of highly asymmetric hard-sphere mixtures" J. Chem. Phys., 130, 164507 (2009)

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*}

  • Negative thermophoresis causes the nanoparticles to move towards the hot plate.
  • Conservation of momentum causes an equal mass of base fluid to move towards the cold plate:

\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*}

  • Nanoparticles have $10^3$ times the mass of the base particles, causing a large ``pumping'' effect.
  • As the mass of the nanoparticle is comparatively so high, its mass specific heat capacity is low. \begin{align*} C_v^\text{Nano} &= \frac{5}{2\,m_\text{Nano}} & C_v^\text{Base} &= \frac{5}{2\,m_\text{Base}} \end{align*}
  • Therefore, thermophoresis also causes a ``pumping'' of heat iff nanoparticles and base fluids have different masses/Cv!
  • This effect can oppose or enhance thermal conductivity, depending on the sign of $L_{A\lambda}$.

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.


  • Nanofluids are an exciting development which can revolutionize heat transfer.
  • Anomalous (en/de)hancements of $L_{\lambda\lambda}$ are possible!
  • It arise from a transport of heat due to thermophoresis, provided there is a significant heat capacity difference between the nanoparticle and base fluid.
  • The mechanism is transient which may cause the apparent difficulties in experimental measurements of the thermal conductivity.
  • More data is needed on the thermal diffusivity of nanofluids, to verify this effect.


  • Craig Moir, Aberdeen
  • Chris Thomson, Aberdeen
  • Dr Leo Lue, Strathclyde
  • Severin Strobl, Erlangen-Nurnberg.
  • Dr Dan Serrero, Erlangen-Nurnberg.
  • Prof. Thorsten Poeschel, Erlangen-Nurnberg.