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.