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.

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.