Here’s a new paper on ChemRxiv that is very much worth reading if you’re a computational chemist (or work with them). And it makes a larger point that’s applicable to everyone else – not an original point, I fear, but it keeps on coming up.
The computational part first: it has to do with Density Functional Theory (DFT) calculations. (Here’s a site on the topic that does not end up knee-deep in equations quite so swiftly). Note, for example, that “functional” is used here in its mathematical sense, as in “a function of other functions”, note the vernacular sense of “performs some task”. Nonetheless, the task that DFT is trying to perform is to estimate free energies, and the functional is electron density (which is in turn a function of space and time, via the Schrodinger equation, for starters). Walter Kohn (share of the chemistry Nobel in 1998) helped develop the Hohenburg-Kohn theorem, which showed that the total ground state energy of a many-electron system is a functional of that electron density – so if you know that (or can estimate it really well), you can get a value for the energy. You don’t have to deal with every electron individually (a task that ranges from “hideously difficult” to “flat out impossible”); you just have to know how many electrons, on average, are inside a given volume of space: the electron density. This is far more tractable, computationally – you can actually knock the whole thing down to a single Schrodinger-ish equation (the Kohn-Sham equation) and work with that.
We are rapidly approaching the limits of my ability to talk about DFT, sadly, as I’m sure is already apparent to people who really know the stuff. But as you can see from that hand-waving view, a key to doing the calculations is how you divide the space you’re working with as you calculate those electron density numbers: what kind of grid do you use, and how small is each piece of it? You’re going to integrate across all those, and you’d hope that the answer would come out the same as (for example) you rotate the molecule around inside that space. Note: we’re not talking about rotating around any particular bond or anything like that; this is a basic “What if I rotated the whole thing just like I was using the Rotate function in ChemDraw; I should still get the same answer then, right?”
Apparently not. The paper linked above (from Andrea Bootsma and Steven Wheeler at Georgia) shows that some of the common integration grids that people use in DFT calculations are not rotation-invariant, and thus can give very different answers indeed. (An analogy: this is sort of like having a device that can estimate how many potatoes there are in a sack, but then you find that it gives you a significantly different Potato Count if you just flip the sack over). This is a real problem for people who would like to calculate the different free energies of conformers, rotational isomers, transition states of bimolecular reactions, and so on. The paper shows that predicted enantioselectivities and regioselectivities of such reactions can produce almost completely opposite results if you just rotate the transition state structure and run the calculations again with the same integration grid. That ain’t good.
Here, we show that popular grids can lead to large errors in free energies, with harrowing implications for DFT studies of organic and organometallic systems. These errors arise from the sensitivity of the entropic component of the free energy to small variations in the values of low-frequency vibrational modes but are distinct from errors due to the treatment of these modes as harmonic oscillators. Moreover, they occur using different DFT functionals, basis sets, and quantum chemistry codes and require at least a (99,590) grid to be resolved in most cases.
That’s a much finer/denser grid than many people are using (which is apparently often more like 75,302). Those low-frequency vibrational modes contribute to the free energy values in what the authors describe as an “unphysically large” way in the standard calculations, and this doesn’t seem to be generally appreciated. Moving to the 99,590 grid (a recommended minimum, actually) helps a lot, but they recommend even then trying your molecules out in more than one orientation to get a handle on how much variation remains.
And that brings up the other point, the one I promised was unoriginal. Looking under the hood of such systems is essential. Actually, being able to look under the hood is perhaps a better way to put it, and frankly, not all of us are able to do that. I’m not, when it comes to DFT or its cousins, believe me. So I have to work with people who I feel have the competence to do that, and I strongly recommend that everyone else do the same. Every calculation or model of this sort has assumptions built into it, shortcuts that it takes out of necessity, approximations that are set to certain tolerances, and these things can most definitely affect what comes out the chute. If you are expecting Golden Tablets of Truth to come rattling down said chute automatically, you are taking a very substantial risk. Any decent computational chemist will tell you this, and you should make sure as you collaborate with them that they feel able to do that. Open that hood. Kick those tires. It’s the only way.