Adventures in Hard Sphere Dynamics

Prolog

After I’ve recently been working on getting the dynamics simulation of hard spheres right, see my previous blog post, I was curious about color coding neighborhood symmetries. To see if my small simulations are already sufficient to see some sort of systematic change. Like can we see “melting” or “crystallization”? 🙂

Turns out yes, for a pretty animation see below. 🙂

You may think this is all settled and researchers have long moved on and everything is understood, spelled out in simple language, organized and filed away for easy grokking. But then you, as I was, would be mistaken. Study of hard spheres seems actually still a matter of ongoing research. There were even experiments done on the space shuttle in 2001, to study the transition between “liquid” and “solid” states (if you are interested in the general concept of “phases” and their transitions, I recommend this wiki page 🙂 ). And more definitive data on liquid-hexatic and hexatic-crystal transition was only collected in 2017. For what hexatic phases are see this wiki page.

Measuring symmetry

So I started researching a bit into how to measure “symmetries”. Various measures could come to mind, but given that hard spheres have been the prime model in statistical mechanics to study many body systems, there ought to be some I can use that have proven useful. If you want to research as well I recommend using the term “bond order parameters”.

I’ve settled on the following

using all 6 nearest neighbors of a hard sphere. So the algorithm is like this:

1. pick a sphere
2. compute the distance vectors to all other spheres
3. compute the lengths of the distance vectors
4. pick the 6 smallest ones
5. for each distance vector take its x and y value to compute the angle theta = arctan(y,x)
6. multiply each angle with 6 and the imaginary unit “i”
7. exponentiate all values
8. compute the mean
9. repeat for the next sphere

and voila you’ve got psi-6 for all your spheres. For me it was most practical to use the absolute value of psi-6.

Visualizing melting

So how could we see “melting” of hard spheres?

By stacking them in a hexagonal grid and then running our molecular dynamics simulation. Easy! This then looks like this

The color scale varies from deep blue = far away from hexagonal order to deep red = perfect hexagonal order.

You may be wondering why the spheres escape upwards and to the right. That’s because I intentionally left some room between the outermost spheres and the simulation box, to encourage things to happen. Because by definition hard spheres don’t have any attractive forces between them they’ll want to fly away.

How to do it yourself

Interested in the code or in how I created the visualization?

I’ve put the code for the above on github for you to check out: https://github.com/eschmidt42/hardspheres-2d/tree/main

You can also install the simulation via uvx and run it from the command line via

uvx --from hardspheres-2d hardspheres2d --help

The visualisation was done using ovito. Scouring the web it just seemed the easiest tool for what I wanted to do. It also has a python package, but I didn’t get done with it what I wanted here. The basic free version of ovito is sufficient for the needs I had as well.

If you just want to dump the molecular dynamics snapshots generated with the above simulation, you can download https://github.com/eschmidt42/hardspheres-2d/blob/d6ea22740fcbf07c2b4d29d88fdc8c2a7f4271ca/hardspheres-hexagon-edmd.xyz and drag it into the ovito GUI.

Once you’ve done that make sure you give the psi-6 column a name that makes sense to you as indicated in the following image:

Then you’ll want to adjust the radius of the spheres / particles a bit like in the following image:

Finally you’ll want to color code each particle using the psi-6 values by adding a “modification” called “Color coding” as indicated in the following image:

That’s it! I hope you enjoyed the read and maybe started tinkering yourself :-).