Highlighted by Jan Jensen

This paper has already been highlighted here and here, so I'll just briefly summarise.

This work is licensed under a Creative Commons Attribution 4.0 International License.

Figure 1B from the paper (CC BY-NC-ND 4.0)

This paper has already been highlighted here and here, so I'll just briefly summarise.

The grid used for the numerical integration in DFT calculations is defined relative to the Cartesian axes, so rotating the molecule will change the integration grid and, hence, the energy. This has been known for some time and, f.eks. Gaussian09 uses a default grid size (Fine, 75,302) where the effect on the electronic energy variation is usually negligible.

Bootsma and Wheeler show that the vibrational entropy and, hence, the free energy is significantly more sensitive to grid size than the electronic energy. Using the Fine grid, the differences in relative free energy changes can be as large as 4 kcal/mol, which could significantly change conclusion regarding mechanisms, etc. The effect comes from the variation in low frequency vibrational modes and the effect can be reduced a little by scaling these frequencies.

However, the errors really only become acceptable when using the UltraFine grid size, which is the default in Gaussian16, especially combined with frequency scaling (which one should do anyway to get consistent results). If you are using Gaussian09 or some other quantum program to compute relative free energies it is definitely a good idea to look at the default grid size and perform some tests.

Note that if you want to perform such tests yourself, you need to re-optimise the molecule after you rotate it because the gradient is also affected by the rotation.

This work is licensed under a Creative Commons Attribution 4.0 International License.

The article in C&EN that picked up on the original preprint provoked me to post a comment, which I reproduce here:

ReplyDeleteI would strongly disagree with the characterisation of the described problem as "flaw inherent to DFT": the problem is entirely a numerical one. As ever in computational science, the translation of a, perhaps even exact, theoretical model into a practical, computationally efficient method involves a plethora of approximations. The accuracy of these approximations is normally controlled by some parameters that need to be carefully chosen by the developer or the user, and the right choice is normally system- and application-dependent. So far so unsurprising.

For the specific case of the reliability/robustness of free energies of "floppy" molecules calculated within the harmonic-oscillator model using DFT, two quite separate issues need to be distinguished:

1) The numerical problems introduced by the integration grid affect mostly low vibrational frequencies. These low frequencies contribute substantially to the vibrational entropy, which is used to compute the free energy. However, the frequencies are obtained within the harmonic-oscillator model, which is invalid for low-frequency vibrational modes. Such modes are significantly anharmonic or even correspond to hindered rotations, rather than vibrations. So notwithstanding their numerical robustness, these frequencies have little physical meaning, but they unduly affect the result.

This is a well-understood and widely known limitation of the standard quantum-chemical approach to calculating free energies. Various solutions have been proposed to eliminate, or at least alleviate, this problem. A simple, effective approach, proposed by Truhlar et al., is to raise all vibrational frequencies below a certain threshold (e.g., 100 cm^–1) up to that threshold. This removes the dependence of the free energy on the exact numerical values of low frequencies, which are physically ill-justified anyway.

Moreover, when things are really "floppy", e.g., two or more initially separate molecules associating into a transition state, the free-energy minimum may differ significantly from the potential-energy minimum. The basic underlying assumption of "static" computational thermochemistry that thermodynamic quantities can be calculated as "single-point corrections" at the potential-energy minimum no longer holds in such cases.

2) Purely numerical accuracy and robustness is quite a separate matter. Irrespective of the physical validity of the approach, the user is entitled to expect that the values obtained (here, the low harmonic frequencies) are numerically accurate and robust. I would therefore agree to this limited extent only that the original article has highlighted an issue that developers should try to address, e.g., by appropriate warnings or automatic tests.