Getting Isosurfaces from a data grid into POV-Ray

POV-Ray makes for fantastic scientific visualizations but if you want to render isosurfaces from a 3D scalar field, there isn’t a very good option. In the past, I used a meta-ball method to make smooth membranes for isosurfaces. There is a isosurface object that allows you to specify a mathematical function that is evaluated to compute an isosurface, but if your data is just a 3D grid, than this isn’t a very good option. Or is it? My new method makes a Fourier series plane wave reconstruction of the data using an FFT. Since plane waves can be evaluated as function, POV-Ray is effectively doing Fourier interpolation! In this scheme, the discrete Fourier transform value at each k-point is a complex value where the real component is the amplitude of the cosine component and the imaginary component is the amplitude of the sine component. You can implement this easily but immediately there will be a problem for any data grid with more than about 10,000 plane waves. POV-Ray won’t parse it, returning the error: ‘Parse Error: Function too complex!’. Since we the isosurface object is only designed real valued functions, we can use the Hermitian symmetry to effectively increase this upper limit to about 20,000. In other words, because the data is real valued, forwards plane waves and backwards plane waves are the same so the backwards ones can be ignored if the forwards ones have their amplitudes doubled. In my code, any plane wave with a negative x component in k-space is ignored and any with a positive x component is amplitude doubled. We can double the plane wave limit again to 40,000 by merging the sine and cosine components into one cosine wave with a phase shift! Quadrupling our capability is great, but 40,000 plane waves still only represents the FFT of a cube of data 34 cell in edge length. The real trick that makes this method possible is sorting by plane wave amplitude and taking only the strongest amplitudes. For my example system, I am using the electron density of a subphthalocyanine crystal. Below I am showing the plane wave amplitudes sorted in order of decreasing amplitude. I am also showing the cumulative distribution function of these amplitudes, while not directly meaningful since waves have phase too, it helps demonstrate how complete a plane wave reconstruction is for a selection of the strongest waves. What does this look like? For 100 plane waves: For 1000 plane waves For 10000 plane waves:
I think that looks nice enough for publication.

Rendering Voids in a Molecular Dynamics Simulation

It's been a while since I've posted. I still need to make that post about X-ray scattering but in the mean time, I'd like to share a picture that I helped make. Katherine Sebeck wanted to render the voids in her epoxy simulation in smooth organic way. I generalized the way Kirke handles grids and my metaball code to handle density grids. The density grid for her epoxy simulation was generated from my smoothed meshing code in Kirke. I talk a bit more about it here. I'll show a movie if Katie says it's okay, but for now here is a large picture (~41 MiB):

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.

My Problem with Isosurfaces of Molecular Orbitals

When people plot their MO's, they show us an isosurface in electron density. It is common procedure to use a density (or wavefunction) value for this surface and report that (Hopefully). I take issue with this.

By not having a standard, we may or may not be missing features of our MO's. Or worse, undesirable features could be hidden. I propose that the fraction of the MO contained by the isosurface be listed. Because we can create a cumulative distribution function for the MO density, specifying the fraction contained by the surface is a sufficient condition for specifying an isosurface density reproducibly.

To demonstrate this loss of detail, I have two isosurfaces from the same series of calculations. These isosurfaces are of the highest occupied MO (HOMO) for a pair of molecules SubPc and C60. The first animation shows an isosurface containing 80 % of the HOMO as the molecules are separated. This is the percentage that I have used to present at conferences.

80 % of the HOMO looks great! You can clearly see the lobes associated the pi orbitals. Next, I'll render the 99 % isosurface of the HOMO.

You can clearly see that some of the HOMO has tunneled to the C60 and that the amount that has tunneled decreases with separation distance. This has an actual affect on some of my calculations. Since we as field have MO's of varying size, the total fraction contained seems to be a more general metric to allow us to more quickly evaluate the works of others.

Or maybe I am just being pedantic? 

If you are wondering, I detailed how I render MO's in another post. The series of images was animates with ImageMagick. In any case, I am feeling very under the weather so do pardon any grammatical errors in this post.

Total Electron Density in an Organic Crystal

Following on my post about electron density isosurfaces for molecular orbitals (here), I used the same methods to make total electron density isosurfaces from VASP output on a boron subphthalocyanine chloride crystal. Here is the result is very high resolution:

The cool thing is that for some of the atoms you can see where the project augmented waves have left holes in the electron density. This is from the so called muffin-tin approximation, where the nuclear core and inner shell electrons are removed and replaced with a pseudo-potential.

Thermal Ellipsoids

At one point, I was supposed to compute electron/hole mobilities in boron subphthalocyanine chloride using charge transfer integrals. The problem is that these don't give good results unless you have realistic thermal distortions to the molecule. Before getting bogged down in charge transfer physics, I decided to use ab initio molecular dynamics on a single unit cell to get realistic thermal distortions. I used VASP, which has the ability to perform ab initio molecular dynamics at constant temperature via a Nose-Hoover thermostat and at constant volume. I would have preferred to use a constant pressure barostat but there is not option for that in VASP yet. Once they do implement a Parrinello-Rahman barostat, I will be very happy.

Back to the molecular dynamics. Once the system had equilibrated at 300 K, I ran the simulation for about 10 ps. Once I had the trajectory data for each atom, I computed the covariance matrix of the x, y, and z position lists. A singular value decomposition (SVD) can be used to get the principle axes of the multivariate normal distribution. The isosurfaces of probability density are surfaces of constant Mahalanobis distance. We can figure out the Mahalanobis distance that corresponds to any total fraction contained in an iso-probability surface by inverting the radial cumulative distribution function. Which looks like this for 3D multivariate normal:

For an ellipsoid containing the densest 50% of the probability, we can see that we need a Mahalanobis distance of ~ 1.5. So we need scale our principle axes (from the SVD of the covariance matrix) by 1.5 to get our ellipsoid. 

For rendering these ellipsoids in POV-Ray, we can use the matrix transform option on a unit sphere to get our ellipsoid. The principal axis form the column vectors of this transform matrix. Here is what it looks like:

Sorry If changed writing styles much in this post, I've been watching Luther.

How I Render Molecular Orbitals

I want to start by saying that I don't like GaussView. It makes rendering of molecular orbitals well enough to convey information but I want very fine grained control over how I stylistically represent data and it bogs down our old Macs. Also, it is not free, which is odd because it is useless without Gaussian. I know that they are making money off of it, but I feel like the freemium model would work better case.

Anyways, Gaussian is a pretty impressive piece of software with really terrible I/O options. To render my molecular orbitals (MOs), I have to use the formchk utility to make .fchk files for the cubegen utility to make .cube files. I usually make a fine mesh for a nicer rendering.

MOs are usually rendered by their isosurface value of the wavefuction or electron density. I use electron density for my isosurfaces. Once I have my isosurface density chosen, I pick a density above and below that and remove all the cells that are outside that range. This filtering drastically speeds up the next step. I scan through all the cells and compute a linear polynomial fit for each cell and its 8 nearest neighbours. Basically, it uses a cell and neighbours to make linear local density function of the form:

rho = a*x + b*y + c*z + d

Since the desired isosurface forms a plane in this function, we can find the point in the plane that is closest to the center of the cell. If this closet point is inside the cell of interest the point is added to the surface point list. 

For each point in the list, I add a sphere with a size close to that of the cell size. This collection of spheres for representing a surface is blended together in POV-Ray using the blob object. The blob object blends the spheres together using the metaballs method for a smooth look. I color the spheres by the sign of the wavefunction at each cell. Here is an example using the highest occupied MO of subphthalocyanine chloride.

Bonus C60 molecule.