Rendering for a 3D Movie or Smashing Two Copper Spheres Together at 1 km/s

At one point, I needed to have a general code for finding clusters of atoms. I was going to use a cutoff distance between atoms to determine if there was a bond and thus cluster connectivity. At the time, I didn't realize I was implementing a breadth-first search. Fast forward a few days , I realized that needed I to test it on something. I decided it would be fun to use and exploding/fragmenting structure to get many clusters to test my code on.

For this, I simply smashed two copper spheres into each other at ~ 1 km/s using molecular dynamics. I used a Rydberg potential for copper that I had optimized for mechanical properties. I actually implemented it into LAMMPS and sent the pair potential file and description in their format to the developers so they could add it to pair_style library. That was about three years ago and as of this post, they haven't added it. I guess not enough people use the Rydberg potential for them to take interest.

Anyways, I had these exploding copper spheres via atomic trajectories that had fulfilled their scientific purpose. Just for fun, I decided to render the simulation in POV-Ray and make a movie with FFmpeg. Here is that looked like:

Pretty cool huh? Well that sat on my hard drive for about a year. Meanwhile, I had been trying to convince my adviser, John Kieffer to purchase a 3D TV to replace our office projector for two reasons. The first reason was that 3D TVs could be used to view atomic structures and the second was that the projector exhaust fan would cook whomever was sitting next to it. He finally found the money for this and put me in charge of purchasing a TV. I chose the Sony KDL-70R550A with the advice of the University of Michigan 3D Lab because passive 3D TVs take two side-by-side images and interlace them for the 3D output. It much easier to configure the 3D system when you don't have to worry about the GPU outputing a switching signal for each eye. The disadvantage is that you get half the lateral resolution in 3D mode. So if you have a side-by-side 3D mode on your screen, check out this video:



The funny thing is I don't care for 3D movie showings in theaters, but as a scientific tool they are pretty handy.

Glass Made from Mables

This is an old video, probably my earliest video that I can easily find. I made this video using a collection of scripts and maybe a library that would eventually become Kirke. I had just started learning how to use POV-Ray... so about fall 2010. We didn't have a development server and I didn't have my desktop Mastodon. To get POV-ray, I had ask for it to be installed on Nyx. I didn't even have anything cool to make a video of, I had just joined my research group! So I asked Kathrine Sebeck for some of her data.




This is a soda-lime glass molecular dynamics simulation rendered in POV-Ray with ray tracing. Each atom was had a sphere created for it with some size meant to represent the atom size. I have no idea what size I used for each, probably the atomic radius given on Wikipedia. Anyways, each of the spheres was made to look like a glass marble by changing the index of refraction. After rendering each frame, I combined them all using FFmpeg.

Meshing Densities with Smoothed Binning and Interface Tracking

About two years ago, I was working on growing oxide surfaces on silicon using molecular dynamics. I was basically just creating a bunch of oxygen atoms over the surfaces of a ~ 10 nm thick sheet of silicon. I kept everything at 1200 C and constant pressure so I could get a nice layer thickness versus time curve. The problem is that using binning to compute the layer density is a bad idea because you'll end up with sharp peaks if you reduce to the same size at the lattice spacing or less. There are two solutions to this problem that I know of. One is to use kernel density estimation and the other is to use smoothed bins with weighting functions. At the time, I had the smoothed binning implemented with weighting functions implemented up to 7th order. So I went with that.

I should also point out that kernel density estimation is basically the continuous analytical extension of smoothed bins with the weighting functions. Basically, if the order goes to infinity and the bin size gets infinity small you get the same result. A good reason not to use kernel density estimation is computational complexity. For smoothed binning, it is the number of particles (or measurements positions) Nparticles is multiplied by the order of the weighting function Norder. O(Nparticles*Norder). In higher dimensions, its O(Nparticles*Norder^dimensionality). I think I used 3rd order here.

For kernel density estimation, every output point Noutput has to have a sum over all the position data. So the complexity goes as the number of particle  Nparticle times the number of output points, or mathematically O(Noutput*Nparticles). For most cases, the number of output points is much higher than the order of the binning so its almost always faster to using smoothed binning but not as pretty.

A few words on the weighting functions. Conceivably, you could use whatever function you want to determine what fraction of a particle gets added to its nearest bins. I prefer to use convolutions of the rectangle function. Conceptually, this is like saying that a particle in a 1st order bins has a uniform distribution in that bin (a rectangular function), but giving it a second bin convolutes the two rectangular functions probability distributions giving a triangular function. And indeed, the triangular function gives a fraction of the mass to the two nearest bins linearly weighted for the fractional distance between them. You can go on forever convoluting these rectangular functions getting higher order piecewise polynomials that represent a function that looks increasingly like a Gaussian. There is a second great aspect of the convolutions of the rectangle function: They have simple Fourier transforms. The Fourier transforms of convolutions of the rectangle function are simply the sinc function raised to the power of the order of the convolution. I will say that I am fond of third order weighting functions because the mass is given to a to the cell containing the particle and the two nearest cells. In this scheme, if the particle on on the border of two cells, it is given evenly to those two and no more. If it is in the middle of the cell, it gives ~80% to the middle cell. (I'd have to double check for the exact number) I feel like 3rd order is the most sensible since it shouldn't share with more than the two nearest cells.

But I digress, These first video shows the mass density of the thin sheet of silicon. There are other oxygen atoms but they are removed if they are not bonded to silicon.

The only problem with this method is that people want a layer thickness and my layers are diffuse. I decided that the best course would be to fit a rectangle by using the variance of the height of the oxygen atoms and overlaying it with a rectangle with the same height variance. This rectangle's height would represent my layer thickness. I used total number of oxygen atoms in the layer and the volume of rectangle to compute the oxygen density in the layer. It works very well. Here is a video with the assumed density profiles plotted with binned profiles.

You can see at first the layer grows in density then thickness as the surface reaches full oxygen coverage and the oxide layer grows thicker. Now I can plot my layer thickness and density versus time. I also get two surfaces (more data points) because I use a thin sheet of silicon.

The layer thickness pretty well matches the power law growth predicted by diffusion models. Nice. As for the density, it looks like an exponential decay to a constant value as the surface saturates and then doesn't change much as the layer gets thicker.

Welp, that was cool. I made all the plots with Matplotlib and movies with FFmpeg or whatever they were calling it then. The simulation was done in LAMMPS and post processed with Kirke.

My next post will have less text, more pictures.