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.

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.