Event-driven MD and anomalous thermal conductivity

Marcus N. Bannerman
School of Engineering, University of Aberdeen
m.campbellbannerman@abdn.ac.uk

https://DynamOMD.com is an open-source general Event-Driven Particle Dynamics (EDPD) package.
How general can EDPD be?

Hard-spheres are the most well-known EDPD potential and have many direct applications,
E.g. Granular-damper systems in microgravity with a spring forcing the box.
M. N. Bannerman, J. E. Kollmer, A. Sack, M. Heckel, P. Mueller, and T. Poeschel, “Movers and shakers: Granular damping in microgravity,” Phys. Rev. E, 84, 011301 (2011)

Beyond hardness

• Stepped potentials can be used to approximate continuous potentials in EDPD.
• Here, time-steps, $\Delta t$, have been swapped for potential energy steps, $\Delta U$.
• Only faster than traditional methods for low/gas densities (unless accelerated via FPGAs).

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

• What makes DynamO unique is its generality...
  ...but why is that so hard to achieve? Why don't we have barostats, or long-range forces?
• While both ED and TD particle dynamics simulations solve Newton's equation of motion, ED simulations only use models where the integral is piecewise analytic.
• Events are transitions between analytical solutions of the system's motion.

How does EDPD work?

Say we're simulating bouncing hard spheres from a catapult in a gravitational field.

Assuming no aerodynamic 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*}$$

When it hits the ground at a time $t+\Delta t$ we have: $$\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}=r_{ground} \end{align*}$$

The time of this "event" is a root of the overlap equation, $f(t)$: $$\begin{align*} f(t +\Delta t)=r_y - r_{ground} + v_y\,\Delta t+g_y\,\frac{\Delta t^2}{2} = 0 \end{align*}$$ (note $f<0$ is an overlap, $f=0$ is contact, and $f>0$ is not overlapping)

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.

DynamO has root finding algorithms for arbitrary order polynomials, and transcendental functions with bounding functions, but they must be fast AND not miss any roots! This is hard.

Generalising for multiple particles

Consider two duelling 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.

When we update this list, we want events to only update a few particles, to stay fast! This is restrictive (no barostats), and sorting the event list is done in $\mathcal{O}(1)$ time!

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

The impulse, $\Delta \vec{P}$, is a solution to the conservation of 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.

Actually this bit is trivial even for rotating bodies.

The EDPD 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.

Modern EDPD algorithms

Modern EDPD algorithms are quite complex.

• Originally, events were sorted using heaps giving $\mathcal{O}(\log_2(N))$ scaling for updates.
• Modern implementations use bounded-calendar priority queues¹ for sorting, leading to $\mathcal{O}(1)$ costs for updates.
  
  1: Paul, G., 'A Complexity ${O}(1)$ priority queue for event driven molecular dynamics simulations'. J. Comput. Phys. 221, 615-625 (2007)
• Lazy deletion is also used, where invalid events are deleted only when they are "next".
• As events only involve small, localised particles, the system can run out of sync. i.e., each particle is left at different points in time (this time-warp algorithm requires a neighbour list).

Initially, the simulation is run in synchronous mode, then time-warp is enabled (with the same simulation speed). Only the required updates to the particles are shown.

Precision

Much of what is presented so far is well established, what is unique about DynamO?
Stability and generality.

Even though EDPD is analytic, computers are not exact. Round-off error occurs at the limits of the machine precision.

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 inelastic hard-sphere model.

Trying to simulate an "inelastic collapse" highlights a catastrophic numerical issue, the particle falls through the floor.

Although $r_y$ should equal $r_{ground}$ after an event, limitations of the precision of the computer leave it at $r_{ground} \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 (or events) and the projectile falls through the ground!

This is not just an edge case, this problem is systemic.

It even plauges elastic/molecular systems for certain sequences of rare ($P=10^{-10}$) events:

Machine precision is not enough to stop particles entering forbidden/invalid states, which may get rapidly worse.

Previous approaches typically attempt to prevent overlaps from forming using small adjustments, but round-off error is inherent and cannot be easily overcome.

Events, by their nature, move a system approaching an invalid state towards a valid state.

An algorithm that always moves towards a valid state is: Events occur at time, $t$, if particles are in contact, $f(t)=0$, or in an invalid state, $f(t)<0$, and this state is deterioriating $\dot{f}(t)<0$.

This appears to be stable in all cases, and correctly simulates inelastic collapse¹.

The challenge of generality

Applying this stable algorithm requires finding the roots of the overlap function $f(t)$, and its time derivative $\dot{f}(t)$. These functions become increasingly complex. For example,

• Spheres in gravity require quartics to be solved.
• Hard lines and other rotating bodies require transcendental root finding.
• Every variation (compression, shearing, gravity/external-force) has a different overlap function!
• This all must be fast and robust! DynamO uses compile-time computer algebra with homotopy methods to implement this at speed and in general.

We'll now look at an interesting application of DynamO which needs its speed.

Motivation for nanofluids

• Heat transfer is the most important area for research.
• Anything we do must satisfy the second law of thermodynamics: \begin{align*} \oint \frac{{\rm d}Q}{T} \le 0 \end{align*}
• Thus heat transfer ultimately drives (and dissipates from) all processes we carry out.
• It is increasingly the limiting factor in the design of modern devices.
• Although novel techniques are now employed (e.g., heat pipes, right), we are fast reaching the limitations of our heat transfer fluids.

Motivation for nanofluids

Water is an excellent heat-transfer medium, but has a relatively poor thermal conductivity, particularly compared to solids such as copper.

• This came to a head in 2009, where the International Nanofluid Property Benchmark Exercise prepared standard samples and issued them to multiple organizations for testing.

  J. Buongiorno, D. C. Venerus, N. Prabhat, T. McKrell, J. Townsend, R. Christianson, Y. V. Tolmachev, P. Keblinski, Lin-wen Hu, J. L. Alvarado, I. C. Bang, S. W. Bishnoi, M. Bonetti, F. Botz, A. Cecere, Y. Chang, G. Chen, H. Chen, S. J. Chung, M. K. Chyu, S. K. Das, R. Di Paola, Y. Ding, F. Dubois, G. Dzido, J. Eapen, W. Escher, D. Funfschilling, Q. Galand, J. Gao, P. E. Gharagozloo, K. E. Goodson, J. G. Gutierrez, H. Hong, M. Horton, K. S. Hwang, C. S. Iorio, S. P. Jang, A. B. Jarzebski, Y. Jiang, L. Jin, S. Kabelac, A. Kamath, M. A. Kedzierski, L. G. Kieng, C. Kim, J.-H. Kim, S. Kim, S. H. Lee, K. C. Leong, I. Manna, B. Michel, R. Ni, H. E. Patel, J. Philip, D. Poulikakos, C. Reynaud, R. Savino, P. K. Singh, P. Song, T. Sundararajan, E. Timofeeva, T. Tritcak, A. N. Turanov, S. Van Vaerenbergh, D. Wen, S. Witharana, C. Yang, W.-H. Yeh, X.-Z. Zhao, and S.-Q. Zhou
  “A benchmark study on the thermal conductivity of nanofluids”,
  J. Appl. Phys., 106, 094312 (2009)

• All samples were blind, organizations did not know what they were actually testing.
• The size of the author list alone indicates how widespread the issue of repeatability had become.

Conclusions from history:

• Results are surprising, but the majority are within classical/continuum limits (not truly anomalous) and most are within Maxwell's limits (not even anomalous).
• If results can exist beyond classical/continuum limits, they can only be explored and explained by molecular (non-continuum) models.

Questions from history:

• Are anomalous effects possible? If so, why haven't we seen it more conclusively? Why is it so difficult to measure anomalous thermal conductivity?
• Is the scatter in the current anomalous results due to poor experimental procedure, high-difficulty, and/or some phenomena?

The Binary Hard Sphere Model

• For $\sigma_1/\sigma_2=10$ and $m_1/m_2=1000$ (see right), thermal diffusion in response to a heat gradient is evident! First result is that conduction in binary mixtures is both transient and inhomogeneous in concentration.

Challenge: The definition of thermal conductivity

• The standard experimental definition of thermal conductivity is $\vec{Q} = - \frac{A\,k}{L} \Delta T$
• Most experiments use this, does missing thermophoresis explain the scatter in the literature?
• To better understand the effect of thermodiffusion we need hydrodynamics, but multi-component hydrodynamics actually has many definitions of the heat flux, $\vec{q}$: \begin{align*} \vec{q} &= - L_{qq}\,T^{-1}\nabla\,T - \sum_a L_{qa}\, T \nabla\,\frac{\mu_a}{T} & \text{(Mainstream)}\\ &= - L_{qq}'\,T^{-1}\nabla\,T - \sum_a L_{qa}'\, \nabla\,\mu_a & \text{(Prime)}\\ &= - L_{qq}''\,T^{-1}\nabla\,T - \sum_a L_{qa}''\, \left(\nabla\,\mu_a + s_a\,\nabla\,T\right)& \text{(Double prime)}\\ &=\ldots \end{align*} where $L_{qq}$ is the thermal conductivity, $L_{qa}$ is the thermal diffusivity, $\mu_a$ is the chemical potential of species $a$, and the primes denote alternative definitions.
• Temperature terms can be arbitrarily moved in and out of the definitions of $L_{qa}$ and $L_{qq}$, so there is no uniquely defined hydrodynamic thermal conductivity! Which is closest to the experimental value, $k$?

• Considering “steady state” where the mass fluxes in the system are zero, all definitions of the heat flux can be collapsed into to one expression: \begin{align*} \vec{q} &= T^{-1}\left(L_{uu} - \frac{L_{au}^2}{L_{aa}}\right) \nabla\,T \\ &=-\lambda\,\nabla\,T \end{align*} where the “steady-state” thermal conductivity $\lambda$ has been implicitly defined.
• Even though $\lambda$ is unique, it still may depend on temperature and concentration, both of which vary across a conducting system.
• Our first challenge is to use large NEMD simulations to prove $\lambda$ is close to $k$.

• Once system size effects are scaled out, it appears that the steady state¹ thermal conductivity, $\lambda$, corresponds to the observed thermal conductivity, $k$.
  ¹As predicted using Revised Enskog Theory in the third sonine approximation, and verified using equilibrium molecular dynamics simulations. See PhD thesis, Craig Moir, “Thermal Transport in Mixtures,” School of Engineering, University of Aberdeen (2020)
• This is suprising, as it says that strong inhomogeneities in the concentration, temperature, and density all “average” out!
• We can now use kinetic theory to explore the full parameter space looking for truly anomalous results.

• We find that certain mixtures exhibit thermal conductivities not only below the series bound but even below the pure components(!).
• This implies that a mixture can be more insulating than its components.
• Even our simple/fundamental fluid can depart from classical behaviour.
• In fact, the mainstream thermal conductivity $L_{qq}$ also exceeds the parallel bounds (and the components themselves!).
• This could hint that a fluid can switch from anomalously high to anomalously low thermal conductivity over short timescales, but this needs very careful exploration.

• We predict using kinetic theory that certain real gas mixtures will exhibit truly “anomalous” thermal conductivity, below the two components used to make the mixture.
• Helium-Hydrogen appears to be one of the strongest predicted effects, and data is available in the literature.

• Interestingly, the Helium-Hydrogen literature in the 70s also had some controversy over “anomalous” heat transfer with extremely low thermal conductivities, which were later found to be unrepeatable!
  A. G. Shashkov, F.P. Jamchatov, and T. N. Abramenko, “Thermal conductivity of the hydrogen-helium mixture”, J. Eng. Phys., 24, 461-464 (1973)
• “Those who do not learn history are doomed to repeat it.”
  George Santayana

Conclusions

• The steady-state thermal conductivity $\lambda$ is a close approximation of the observed thermal conductivity $k$.
• Binary mixtures of hard-spheres display anomalous thermal conductivities outside classical theory:
  Classical limits can be broken.
• This is possible even in the ideal gas limit:
  Structural arguments (interfacial resistance/layering models) are not relevant for this type of enhancement.
• Helium-Hydrogen systems confirm that anomalous drops in thermal conductivity are possible, but anomalous enhancements have yet to be proven (although they are predicted by Enskog theory).

Appendix

University of Aberdeen

Coupled kinetic theory and hydrodynamic simulation of two heated walls ($k_B\,T=\{1,\,1.5\}$) ten NP diameters apart contacting a size ratio 1:0.5, mass ratio 1:0.5, $\rho=0.2$, $k_B\,T=1$ system.

10:1 size, 1000:1 mass, binary HS at constant volume fraction $\phi$ or constant pressure $p$. Symbols are MD and lines are from kinetic theory.¹

• Truly anomalous results exist in $L_{q,q}$ for hard-spheres, and all nanofluids at least share key features with this system.
• Nanoparticles (in this case) are insulating. If $L_{qq}$ is a useful definition then this is an enhancement like the olive-oil/water system.
• Analyzing the kinetic theory expressions, $L_{qq}$ has a dominant term: \begin{align*} \vec{J}_q \approx\left(C_V^{Nano} - C_V^{Base}\right)k_B\,T\,L_{Nano,q} T^{-1}\nabla T \end{align*}
• Large mass difference (as $C_V=5\,m/2$), coupled with thermophoresis causes heat pumping effect IF $L_{qq}$ is relevant.

Stepped Potentials

Can we simulate "real" contact forces using EDPD?

We recently demonstrated elastic forces can be "stepped".

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

Stepped Potentials

Conversion of a potential first requires the steps to be placed.

Optimal placement was determined to be at fixed changes in the potential energy, $\Delta U$ (analogue of timestep, $\Delta t$!).

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

Stepped Potentials

Optimal energies are given by equating the second virial contribution, which is a volumetric average at $T\to\infty$: $$\begin{align*} \Phi_i^{Volume} &= \frac{3}{r_{i}^3-r_{i+1}^3}\int_{r_{i+1}}^{r_{i}} \Phi \left(r\right) r^2\, {\rm d}r\nonumber \end{align*}$$

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

Stepped Potentials

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

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

Stepped Potentials

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)

Stepped Potentials

Molecular forces are too soft and EDPD is inefficient at high density: harder, granular potentials are more efficient.

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

Event-Driven Spring Dashpot

We have worked on a spring-dashpot force as an example implementation of viscous forces.

Previous work¹ in this direction focussed on single-event collision dynamics and required storing "history". It also demonstrated the importance of considering individual particle trajectories.

A pair of colliding particles under with an elastic spring is simulated to test the trajectory.

^1: Müller, P., Pöschel, T. 'Event-driven molecular dynamics of soft particles'. Phys. Rev. E, 87, 033301 (2013)

For a single collision, the key parameter is where it lands on the final step (inset). $\tau = KE_{final}/(U_{final+1}-U_{final})$

There is an issue in the last step: for vanishing remaining energy the particles become stuck! ("elastic collapse"?).

By enforcing a minimum kinetic energy (immediately bouncing $\tau<\alpha$, where $\alpha=1/9$ in the last step), the average collision time is correct.

Published dissipative forces in EDPD have so far been limited to inelastic coefficients.

To implement arbitrary dissipative forces cheaply, it should be approximated using an impulse/energy loss at each step:

$$\begin{align*} \Delta KE_{diss.} &= \int F_{diss.,ij}\,{\rm d}r_{ij} \\ &\approx F_{diss.,ij}(r(t)) \Delta r \end{align*}$$ Where $\Delta r$ is the width of the step the particle is leaving. No step-change = no energy loss (quasi-static elastic limit/no inelastic collapse).

Even for high inelasticiy ($\gamma$ set such that $e_{exact}=0.4$) this simple approach is accurate with minor time scattering.

Relatively small numbers of steps are acceptable, with collision time showing the most scatter.

Maximum distance is rigidly enforced.

And the model demonstrates the required constant coefficient of restitution for spring-dashpot particles.

DynamO has a compile-time Computer Algebra System (CAS) library. It allows simple definitions of overlap functions:

Variable<'t'> t;
Vector<> rij, vij, aij;
double diam2_ij;		      
/*...*/
delta_t = nextroot(pow<2>(rij+t*(vij+aij*t/2))-diam2_ij, t);
	      

The compiler will transform the overlap function to create derivatives, and select different classes of root finding.

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.

For up to 3rd order polynomials, solution by radicals is used. Solution by radicals is too numerically unstable for 4th order and impossible for 5th order.

For higher order polynomials, roots are bounded using Local Max Quadratic estimates¹, then bisected using Sturm's method to avoid repeated transformations of the polynomial.

For transcendental overlap functions, the roots are bounded using a bounding object with a polynomial overlap function. The bounds are then updated using a "worst case" polynomial approximation constructed using Taylor's theorem. This is accelerated using standard numerical root finding techniques.

Outlook of EDPD

• Complex shapes, forces etc. are easy to implement using the CAS library of DynamO.
• Viscous forces are now available in EDPD (48 years late!).
• Friction forces (Coulomb's law, Cundall and Strack, etc.) should come soon using these ideas.
• With the above improvements in place, new advanced optimisation algorithms for DEM are possible and should be investigated.

"Sleeping particles" algorithm allows free simulation of stationary particles. Possible massive speed-ups in heaps, rolling drums, etc.