I start every post with a motivation section, and normally have something a bit useful to say. For this post, I have no "useful" motivation: I think event driven dynamics sounds interesting, and I want to give it a go. I think I could shoe-horn it into my research, but not to any great need.

Traditionally, simulation of hard sphere particles is achieved by integrating a set of equations of motion, using some fancy method. This involves a (small) finite timestep. This can result in particles intersecting (slighlty) which is physically unrealisable and makes the simulation useless. To get around this, the "hard sphere" property of the particles is weakened to a "sort of hard sphere" by use of a potential - an equation which gives the force (derivative of the potential energy) experienced by a particle in response to another.

A common potential is the Lennard-Jones 6-12 potential:

\[ V_{LJ} = 4\varepsilon\left[{\left(\frac{\sigma}{r}\right)}^{12} - {\left(\frac{\sigma}{r}\right)}^{6}\right] \]

This potential does a good job: it rapidly increases as \(r_{sep} \to r_{radius}\), giving a continuous function describing the radius of a particle compatible with the finite step. However, I've implemented this in an MD sim (well, my supervisor did) and it still dies when particles intersect. Instead of just being insensible, the force on the intersecting particles becomes huge and the whole sim starts becoming very excited. This is better as it gives an indication of something being wrong, but there is still something wrong.

Enter events. An event driven simulation calculates the time the next event will happen, then jumps to that point. The tricky part is then in sorting the events efficiently. The hard sphere property is maintained in this kind of simualtion, which makes it nice to do simulations with frictional particles as the boundaries of the particle are well defined and it is known when particles are in contact.

Here's where my research could get involved: I do work with suspensions. The prevailing theory on the flow characteristics of dense suspensions of particles involves frictional contact between the particles. I am running experiments to attempt to collect statistics on the frictional interaction. If I could model this, I could get a better idea of what to look for.

In addition, I have a set of 3D printed bearings which are integral to my experimental set up but are a source of mechanical noise in my stress/torque readings: I would like to be able to simulate this (as a matter of curiosity) to see if I can re-create the type of noise I see in my measurements.

Ooft, for a post claiming to not have much motivation, I've really rambled on here!

In this post I'll run through the difficulties and major challenges I faced while building BEARS.

The simulation should:

• [X] be able to model particle interaction in a ballistic (i.e. classical mechanic) manner.
• [X] be able to operate under periodic boundary conditions (where the sim-box is repeated in every dimension to mimic a larger population of particles and to limit edge effects).
• [X] use event driven dynamics to enforce proper hard sphere behaviour.
• [X] model in three dimensions.
• [X] take into account rotational dynamics (even for just spheres: rotational velocity and orientation is important).
• [X] take into account friction.

For this project, I'm defining an event as a collision of two previously separate particles. Every particle has a position and a velocity. These are both described by 3vectors:

\[ \vec{p} = \begin{bmatrix} p^x \\ p^y \\ p^z \end{bmatrix} \] \[ \vec{v} = \begin{bmatrix} v^x \\ v^y \\ v^z \end{bmatrix} \]

Confusingly, they tend to refer to position as \(\vec{r}\) in physics. I will not, position is \(\vec{p}\). It won't collide with momentum, which I will denote with a capital: \(\vec{P}\).

Given these four vectors, \(\vec{p}_1\), \(\vec{p}_2\), \(\vec{r}_1\), and \(\vec{r}_2\), the challenge is to find out if they will collide, the time at which the particles will collide and where they will be.

This is surprisingly simple (for plain ballistic simulation).

Collision at:

\[ \left| \vec{p}_1 - \vec{p}_2 \right| = \left| \Delta\vec{p} \right| = \frac{d_1 + d_2}{2} \]

Where \(d\) is the diameter of a particle. The magnitude of the separation vector is given by:

\[ \left| \Delta\vec{p} \right| = \sqrt{\Delta{}p^x \times \Delta{}p^x + \Delta{}p^x \times \Delta{}p^x + \Delta{}p^x \times \Delta{}p^x} \]

And each position vector after time \(t\) is given by:

\[ \vec{p}_1^t = \vec{p}_1 + \vec{v}_1 \times t \]

Which we can sub into the above and the problem is just finding the roots of a polynomial:

\[ \Delta\vec{v}\cdot\Delta\vec{v}t^2 + 2\Delta\vec{p}\cdot\Delta\vec{v}t + \Delta\vec{p}\cdot\Delta\vec{p} - \left(\frac{d_1 - d_2}{2}\right) = 0 \]

Solving this gives us the time of a collision. Negative discriminant results in complex results: which means that the particles will not collide. Further, if the magnitude of the separation vector is the same as the average diameter, then the particles are already in contact. If the two solutions are positive, then the smallest one is the time to collide (the other being the time to collide and pass through to the other side). If the two solutions are negative, then the collision will not happen (unless time starts to go backwards!)

For a ballistic simulation, we can model the particle interaction by assuming perfect elastic collisions and can use the conservation of both momentum and kinetic energy to calculate the resulting velocity vectors of the two particles.

\[ \vec{P}_1^b + \vec{P}_2^b = \vec{P}_1^a + \vec{P}_2^a \] \[ E_1^b + E_2^b = E_1^a + E_2^a \]

(\(^b\) denotes before a collision and \(^a\) denotes after), we can calculate the velocity of a particle after a collision:

\[ \vec{v}_1^a = \frac{m_1 - m_2}{m_1 + m_2}\vec{v}_1^b + \frac{2m_2}{m_1 + m_2}\vec{v}_2^b \]

(Thanks, wikipedia.)

This meets the requirements of this simple ballistic simulation. What if we go further? What if friction is involved? What if the particles are nolonger "perfectly hard"? What if the particles are not interacting in a vaccuum but are suspended in a fluid? I'll cover these points in future posts as I build up the sim from this starting point.

The structure of the program is fairly simple. We want to have a "containter-ish" sim class which is the box of particles, and manages timestepping and event detection - for now, when we get to parallelisation the particles may be given domain over this matter! Particles will be generated from a configuration file (the particle configuration file, not program settings kind of configuration file), which will be in a simple ascii format (space separated value or tab-separated value).

Particle objects contain their properties:

• position,
• velocity,
• orientation,
• angular velocity,
• mass,
• diameter

And also decide whether they collide with another particle or not.

Particles are built on the heap, pointers to which are owned by the sim (so that pointers can be pased about which makes things nice and fast).

Events are objects too (well, more like structures) and are returned by the particle check_will_collide method. These objects contian information on what particle is involved in the collision and when the collision will happen. (Pointers to) these events are stored in the sim in a container object.

The EventCollection object manages the CollisionEvent objects which are push_back'd on to it. Internally, events are stored as a map, with a std::pair<int,int> of particle IDs as keys. This is a simple way to ensure particles don't interact with each other twice. It allows old events to be cleaned out as they are replaced. However, we need to sort the events by time and choose the smallest one - or, rather, we just need to find the minimum (an \(\mathcak{O}(n)\) operation). This collection becomes the owner of the event pointers pushed on to it: it ~delete~s them on destruction and when replaced.

Trajectory output is important: we want to be able to track sim data (particle positions, kinetic energy, …) over time, so that we may make and test hypotheses using the simulation. Otherwise, what's the point?

During the debug stage, we want to be able to profile the code in some way. The main loop of the run method will need timers to see how long the event detection and output steps take, to identify areas that can be improved and to track how well each improvement performs.

Finally, there will be a "main" file which parses arguments, creates a sim instance, loads the config file, and runs the simulation.

The event finding process is then:

• update event collection
• find next event
• apply event: update all particle positions and transfer momentum of interacting particles
• repeat until no events left, or time has reached limit

The current process works great for small numbers of particles (near instantly for a box of 1000 particles evaporating), but is becomes very computationally expensive as the number increases. The search is \(O(n^2)\). That search order could be better.

We only need to recalculate events for particles which were involved in the previous event, or were going to interact with particles in the previous event.

My supervisor suggested only re-calculating particle positions when they're involved in an event… but we'd need to calculate position anyway when deciding if a collision will take place. Perhaps there's a better way to decide if particles interact?

We can parallelise (not really an efficiency gain, but a performance gain to be sure) by chunking the list of particles and calculating next-events on each thread.

The code for this stage of development is immortalised on my github here. I've made a separate branch for it as I will add features and alter the functionality in the future, but bugfixes may need to be retroactively applied.

I've touched on this earlier, but the goal for this is to have a very performant simulation package for hard spheres. That means I want to be able to simulate upwards of a million particles in a reasonable amount of time on conventional hardware (e.g. my laptop's 8 core i7). I could use 100 or so cores of my department's mad simulation server, but that would be too easy!

Given that end goal, here are my next steps:

• Particles with agency over their interaction. A particle should decide for itself when/if it will interact: this will allow N-threads to pop particles off a queue and calculate its event. This means each particle will
• Simplification of the vector class.
• Simplification of the event calculation.