Systematic Error Estimation for Chemical Reaction Energies

Gregor N. Simm and Markus Reiher (2016)
Contributed by Jan Jensen

Simm and Reiher present an approach to estimate the error of density functionals parameterized for a specific system.  System-specific parameterization can yield very accurate predictions for the training set but the applicability of the resulting functional to any other systems, including closely related ones, is not guaranteed precisely because the training set is so specific.  However, by Simm and Reiher provide confidence intervals for each result to assess whether the results are reliable. 

Four functional parameters are optimized to best fit 23 CCSD(T) reaction energies of a model system. As one would expect the functional outperforms several other standard functionals for that particular model system. The real question is whether this performance translates to the real system of interest where CCSD(T) benchmarks aren’t available.  To answer this Simm and Reiher generate a set of new values for one of the parameters that lead to a range of predicted reaction energies from which the uncertainty (standard deviations) in the prediction can be computed.  

For the model system the error is within two standard deviations for 21 of the 23 reaction energies. When the functional is then applied to a slightly larger system, for which CCSD(T) reference values also could be computed, the error is within two standard deviations for 17 of the 20 reaction energies. Thus for most of the reaction energies this approach can be used to reliably gauge the accuracy of the results computed using the system dependent density functional.

It should be noted that the current approach only applies to linear parameters, which is the reason the analysis what only performed on one of the four optimized parameters. But it should be possible to extend it to other parameter types, in which case one also needs to address the large number of parameter-combinations that need to be checked.  

This work is licensed under a Creative Commons Attribution 4.0

Beyond DP4: an Improved Probability for the Stereochemical Assignment of Isomeric Compounds using Quantum Chemical Calculations of NMR Shifts

Grimblat, N.; Zanardi, M. M.; Sarotti, A. M. J. Org. Chem. 2015, 80, 12526-12534

Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

Structure determination has been greatly facilitated by the use of computed NMR spectra to compare with experimental spectra. Perhaps the best method for doing this is the DP4 procedure developed by Smith and Goodman.¹ (I have a previous post on their paper.) The basic idea is that if you have an experimental NMR spectrum and a number of potential structures, the computed spectra for each possibility are ranked by a statistical treatment based on the Student t-test.

Grimblat, Zanardi, and Sarotti question a couple of the assumptions embedded within the DP4 method, and offer a revision that they call DP4+.² The two assumptions are (1) that the chemical shifts are computed at B3LYP/6-31G**//MMFF and (2) that the chemical shifts are scaled and then utilized in the analysis.

To test these assumptions, they examine a set of 72 organic compounds comprising 1219 ¹³C shifts and 1123 ¹H shifts. They optimized the structures at B3LYP/6-31G* and computed the chemical shifts of these compounds using the B3LYP and mPW1PW91 functionals with 6 basis sets (6-31G*, 6-31G**, 6-31+G**, 6-311G*, 6-311G**, and 6-311+G**). With all of the combinations, the standard deviation of both the proton and carbon chemical shifts were significantly smaller than with the originally proposed method.

With regards to the second assumption, they define a new probability functions that multiplies the error using scaled chemical shifts with the error using unscaled chemical shifts, and this they call DP4+. Again with all of the computational methods, the DP4+ prediction outperforms the DP4 prediction.

As a test case, they looked at cryptomoscatone D1 and D2 (1), for which the structures were determined with traditional methods. DP4 predicts that both cryptomoscatone D1 and D2 are structure 1d. However, DP4+ correctly predicts that cryptomoscatone D1 is 1b and cryptomoscatone D2 is 1a.

Lin and Tagliatatela-Scafati have reported the use of DP4+ to aid in the structure determination of plakdiepoxide 2.³ ROESY NMR could not provide definitive judgement of the stereochemical relationship about the bond between the two epoxide rings. They computed a number of conformers of the model compounds 2a and 2b at B3LYP/6-31G(d). The computed chemical shifts were then used with the DP4+ procedure to determine that the structure has the stereochemistry of 2b.

References

(1) Smith, S. G.; Goodman, J. M. "Assigning Stereochemistry to Single Diastereoisomers by GIAO NMR Calculation: The DP4 Probability," J. Am. Chem. Soc. 2010, 132, 12946-12959, DOI: 10.1021/ja105035r.

(2) Grimblat, N.; Zanardi, M. M.; Sarotti, A. M. "Beyond DP4: an Improved Probability
for the Stereochemical Assignment of Isomeric Compounds using Quantum Chemical Calculations of NMR Shifts," J. Org. Chem. 2015, 80, 12526-12534, DOI: 10.1021/acs.joc.5b02396.

(3) Chianese, G.; Yu, H.-B.; Yang, F.; Sirignano, C.; Luciano, P.; Han, B.-N.; Khan, S.; Lin, H.-W.; Taglialatela-Scafati, O. "PPAR Modulating Polyketides from a Chinese Plakortis simplex and Clues on the Origin of Their Chemodiversity," J. Org. Chem. 2016, 81 (12), 5135–5143, DOI: 10.1021/acs.joc.6b00695.

InChIs

Cryptomoscatone D1: InChI=1S/C17H20O4/c18-14(10-9-13-5-2-1-3-6-13)11-15(19)12-16-7-4-8-17(20)21-16/h1-6,8-10,14-16,18-19H,7,11-12H2/b10-9+/t14-,15-,16-/m1/s1
InChIKey=GOQOIZFMLWZVMB-OUUZNBFFSA-N

Cryptomoscatone D2: InChI=1S/C17H20O4/c18-14(10-9-13-5-2-1-3-6-13)11-15(19)12-16-7-4-8-17(20)21-16/h1-6,8-10,14-16,18-19H,7,11-12H2/b10-9+/t14-,15+,16+/m0/s1
InChIKey=GOQOIZFMLWZVMB-GLBZDCTLSA-N

plakdiepoxide: InChI=1S/C18H32O4/c1-6-9-10-13(4)12-17(7-2)16(22-17)18(8-3)14(21-18)11-15(19)20-5/h13-14,16H,6-12H2,1-5H3/t13?,14-,16-,17+,18-/m0/s1
InChIKey=YZQCELNTSCPAPG-KSZMJELXSA-N

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Molecular Dynamics of Dimethyldioxirane C–H Oxidation

Yang, Z.; Yu, P.; Houk, K. N.  J. Am. Chem. Soc. 2016, 138, 4237-4242
Contributed by Steven Bachrach
Reposted from Computational Organic Chemistry with permission

Dimethyldioxirane can oxidize alkanes to alcohols. The mechanism for the oxidation has been controversial, ranging from concerted, to radical intermediates to an H-abstraction—O-rebound mechanism. Yang, Yu, and Houk now offer a molecular dynamics examination of the reaction of dimethyldioxirane with isobutane.¹

Gas–phase (U)B3LYP/6-311++G(d,p)//(U)B3LYP/6-31G(d) computations give critical points outlined in the reaction below. The structures of the transition states and the intermediate are shown in Figure 1.

TS1 26.6 (25.9)	INT 8.6 (9.3)
TS2 8.6 (-0.8)

Figure 1. (U)B3LYP/6-31G(d) optimized geometries of TS1, INT, and TS2. Relative free energies (kcal mol^-1) in the gas (top) and solution (bottom) phases

The PES indicates a rebound mechanism, though in acetone solution phase, there was no transition state located for the second step; it appears to be barrierless. It should be noted that the size of the barrier is very small even in the gas phase. The energy given in Figure 1 is for the gas phase structure computed in solution.

Trajectories for both gas and solution phase were computed. For the gas phase, about 90% of the trajectories lead to separated radicals, but in an acetone about 90% of the trajectories lead directly to the alcohol, with only 10% leading to radicals. Even so, the acetone trajectories divide into two types, a dynamically concerted path where the time gap between the formation of the new C-O and O-H bonds is less than 60 fs, and a dynamically stepwise path where the time gap is greater than 60 fs, though for the trajectories that lead to product the gap is typically still less than 150 fs.

References

(1)  Yang, Z.; Yu, P.; Houk, K. N. "Molecular Dynamics of Dimethyldioxirane C–H Oxidation," J. Am. Chem. Soc. 2016, 138, 4237-4242, DOI: 10.1021/jacs.6b01028.

InChIs

Dimethyldioxirane: InChI=1S/C3H6O2/c1-3(2)4-5-3/h1-2H3
InChIKey=FFHWGQQFANVOHV-UHFFFAOYSA-N

Isobutane: InChI=1S/C4H10/c1-4(2)3/h4H,1-3H3
InChIKey=NNPPMTNAJDCUHE-UHFFFAOYSA-N

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Different approaches to creating (meta-GGA) DFT functionals

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Within the DFT community, John Perdew’s idea of the Jacob’s ladder of accuracy[1] starts with LDA, moves to GGAs with the inclusion of the gradient of the electron density, and is further extended (to rungs three and four in the ladder analogy) with meta-GGA, where kinetic energy density is added, and hybrids, where exact exchange plays a role. Although meta-GGAs have been around for ten to fifteen years, they are only starting to become widely used. I will compare two examples which both seem promising, but also encapsulate the two most common approaches to functional creation.

Meta-GGAs are promising because the kinetic energy density allows some discrimination between areas with one or two electrons (in this way, they are similar in some ways to the electron localisation function, or ELF, which can be used to analyse bonding[2]). This gives some hope that they may be able to fit both strong and weak bonding as well as possibly mitigating the self-interaction error that plagues DFT.

The Minnesota meta-GGA MN15-L[3] takes a very complex functional form with 58 parameters and a non-separable form for the exchange-correlation (adding an extra correlation functional to this) and fits it to a portion of a very large database (I estimate that there are at least 900 entries in the database, covering many different chemical properties, specifically including solid state properties, transition barriers and weak interactions). The resulting functional is local (no hybrid or non-local van der Waals terms were included) and produces extremely small errors when compared to those parts of the database which were not used in fitting. Notably, it out-performs functionals with exchange and van der Waals included.

By contrast, the SCAN functional[4] uses only seven parameters (close to, if not at, the minimal number for a meta-GGA) and includes various physically motivated norms and constraints for the electron gas. The early progress in GGAs was made by satisfying important constraints, so this is seen as a good route to reliability. This functional was also tested on various databases, particularly focussing on solid-state and weak interactions. It has excellent agreement with these, though was published before the MN15-L functional so is not compared directly. In a follow-up paper[5] the performance of the SCAN functional for band gaps is found to be good, though it does not calculate Kohn-Sham gaps, but gaps within a generalized Kohn-Sham theory (a distinction which I don’t have time to discuss here; I may write another blog on this, as it is relevant to hybrid functionals among other things.)

What can we learn from this ? Both functionals perform well in the tests which are published. The functional that you choose will depend in part, as always, on which system you wish to study; however, both of these functionals show some promise in being widely applicable. Your choice will also depend on your attitude to fitting[6]: is a reasonable functional form with many parameters something that you trust, or do you prefer to be more prescriptive, and deal with fewer parameters ? Fifty years since its inception, DFT is still developing, communities are still somewhat divided in the approaches that they take to functional development, but there are an increasing number of ways to achieve efficiency and accuracy.

[1] My colleague, Mike Gillan, reckons that we should instead talk about wrestling Jacob when considering how to improve DFT functionals.
[2] J. Chem. Phys. 92, 5397 (1990) DOI: 10.1063/1.458517
[3] J. Chem. Theor. Comput. 12, 1280 (2016) DOI: 10.1021/acs.jctc.5b01082
[4] Phys. Rev. Lett. 115, 036402 (2015) DOI: 10.1103/PhysRevLett.115.036402
[5] Phys. Rev. B 93, 205205 (2016) DOI: 10.1103/PhysRevB.93.205205
[6] The oft-quoted maxim about fitting an elephant with four parameters has been put into practice (see the original paper here and a nice write-up and a python implementation here)