Suppose you are looking at the reaction aA + bB → cC + dD. You can compute each of these molecules at two computational levels; lets call these M1 and M2. Then the reaction energy is

Now, if the two computational methods are reasonably complete, then ΔE^M1 ≈ ΔE^M2. This can also be true if the reaction has been selected such that one might expect very good cancellation of errors. In this case, the overall problems in computing the reactants are similar to the problems computing the products, and so these problems (i.e., errors) will cancel off. So, if we have the latter condition (a reaction constructed to obtain excellent cancellation of errors), then we might be able to exploit this idea in order to obtain energies of large molecules with a large method while avoiding to actually have to do these very large computations!

How does this work? Let’s suppose the largest molecule in the reaction is molecule C. Since

then

and so

So, we can get the energy of the big molecule C at the big method M1 by computing the energy of the big molecule C at the smaller method M2 along with computing all of the other molecules at both levels. If these other molecules are significantly smaller than molecule C, there can be considerable time savings here. This is the idea presented in a recent article by Raghavachari.¹

The key element here is a systematic means for generating appropriate reactions, ones that (a) involve small molecules other than the molecule of interest and (b) get good cancellation of errors. Raghavachari comes up with a systematic way of creating a reaction with ever larger reference molecules. This is analogous with the methodology presented by Wheeler, Schleyer and Allen.² Basically, the method decomposes the molecule of interest into smaller molecules that preserve the immediate chemical environment around each heavy atom, a method they call CBH-2 (connectivity-based hierarchy). The reaction below is an example of the CBH-2 decomposition reaction for methionine. (Note that CBH-2 is essentially a homodesmotic reaction and CBH-3 is essentially the group equivalent reaction I defined years ago.³)

They apply the concept towards computing the energy of larger molecules (having 6-13 heavy atoms) at CCSD(T)/6-31+G(d,p) by only having to compute these large molecules at MP2/6-31+G(d,p). For a set of 30 molecules, the error in the energy of the extrapolated energy vs. the actual CCSD(T) energy is 0.35 kcal mol^-1.

One of the advantages of this approach is that the small molecules are used over and over again, but they need be computed only twice, once at CCSD(T) and once at MP2.

This is certainly an approach that has been implicitly employed by many people for a long time, but here is made explicit and points towards ways to apply it even more widely.

References

(1) Ramabhadran, R. O.; Raghavachari, K. "Extrapolation to the Gold-Standard in Quantum Chemistry: Computationally Efficient and Accurate CCSD(T) Energies for Large Molecules Using an Automated Thermochemical Hierarchy," J. Chem. Theor. Comput. 2013, ASAP DOI: 10.1021/ct400465q.

(2) Wheeler, S. E.; Houk, K. N.; Schleyer, P. v. R.; Allen, W. D. “A Hierarchy of Homodesmotic Reactions for Thermochemistry,” J. Am. Chem. Soc. 2009, 131, 2547-2560, DOI: 10.1021/ja805843n.

(3) Bachrach, S. M. “The Group Equivalent Reaction: An Improved Method for Determining Ring Strain Energy,” J. Chem. Ed. 1990, 67, 907-908, DOI: 10.1021/ed067p90.

The gold standard in quantum chemistry is the method that is considered to be the best, the one that gives accurate reproduction of experimental results. The CCSD(T) method is often referred to as the gold standard, especially when a complete basis set (CBS) extrapolation is utilized. But is this method truly accurate, or simply the highest level method that is within our reach today?

Řezáč and Hobza¹ address the question of the accuracy of CCSD(T)/CBS by examining 24 small systems that exhibit weak interactions, including hydrogen bonding (e.g. in the water dimer and the water^…ammonia complex), dispersion (e.g. in the methane dimer and the methane^…ethane complex) and π-stacking (e.g. as in the stacked ethene and ethyne dimers). Since weak interactions result from quantum mechanical effects, these are a sensitive probe of computational rigor.

A CCSD(T)/CBS computation, a gold standard computation, still entails a number of approximations. These approximations include (a) an incomplete basis set dealt with by an arbitrary extrapolation procedure; (b) neglect of higher order correlations, such as complete inclusion of triples and omission of quadruples, quintuples, etc.; (c) usually the core electrons are frozen and not correlated with each other nor with the valence electrons; and (d) omission of relativistic effects. Do these omissions/approximations matter?

Comparisons with calculations that go beyond CCSD(T)/CBS to test these assumptions were made for the test set. Inclusion of the core electrons within the correlation computation increases the non-covalent bond, but the average omission is about 0.6% of the binding energy. The relativistic effect is even smaller, leaving it off for these systems involving only first and second row elements gives an average error of 0.1%. Comparison of the binding energy at CCSD(T)/CBS with those computed at CCSDT(Q)/6-311G** shows an average error of 0.9% for not including higher order configuration corrections. The largest error is for the formaldehyde dimer (the complex with the largest biding energy of 4.56 kcal mol^-1) is only 0.08 kcal mol^-1. If all three of these corrections are combined, the average error is 1.5%. It is safe to say that the current gold standard appears to be quite acceptable for predicting binding energy in small non-covalent complexes. This certainly gives much support to our notion of CCSD(T)/CBS as the universal gold standard.

An unfortunate note: the authors state that the data associated with these 24 compounds (the so-called A24 dataset) is available on their web site (www.begdb.com), but I could not find it there. Any help?

References

(1) Řezáč, J.; Hobza, P. "Describing Noncovalent Interactions beyond the Common Approximations: How Accurate Is the “Gold Standard,” CCSD(T) at the Complete Basis Set Limit?," J. Chem. Theor. Comput., 2013, 9, 2151–2155, DOI: 10.1021/ct400057w.

Computing the optical rotation of simple organic molecules can be a real challenge. One of the classic problems is methyloxirane. DFT typically gets the wrong sign, let alone the wrong value. Cappelli and Barone¹ have developed a QM/MM procedure where methyloxirane is treated with DFT (B3LYP/aug-cc-pVDZ or CAM-B3LYP/aubg-cc-pVDZ). Then 2000 arrangements of water about methyloxirane were obtained from an MD simulation. For each of these configurations, a supermolecule containing methyloxirane and all water molecules with 16 Å was identified. The waters of the supermolecule were treated as a polarized force field. This supermolecule is embedded into bulk water employing a conductor-polarizable continuum model (C-PCM). Lastly, inclusion of vibrational effects, and averaging over the 2000 configurations, gives a predicted optical rotation at 589 nm that is of the correct sign (which is not accomplished with a gas phase or simple PCM computation) and is within 10% of the correct value. The full experimental ORD spectrum is also quite nicely matched using this theoretical approach.

References

(1) Lipparini, F.; Egidi, F.; Cappelli, C.; Barone, V. "The Optical Rotation of Methyloxirane in Aqueous Solution: A Never Ending Story?," J. Chem. Theor. Comput. 2013, 9, 1880-1884, DOI: 10.1021/ct400061z.

InChIs

(R)-Methyloxirane:
InChI=1S/C3H6O/c1-3-2-4-3/h3H,2H2,1H3/t3-/m1/s1
InChIKey=GOOHAUXETOMSMM-GSVOUGTGSA-N

Alonso and coworkers have again (see this post employed laser-ablation molecular-beam Fourier-transform microwave (LA-MB-MW)spectroscopy to discern the gas phase structure of an important biological compound: cytosine.¹ They identified five tautomers of cytosine 1-5. Comparison between the experimental and computational (MP2/6-311++G(d,p) microwave rotational constants and nitrogen nuclear quadrupole coupling constants led to the complete assignment of the spectra. The experimental and calculated rotational constants are listed in Table 1.

Table 1. Rotational constants (MHz) for 1-5.

The experimental and computed relative free energies are listed in Table 2. There is both not a complete match of the relative energetic ordering of the tautomers, nor is there good agreement in their magnitude. Previous computations² at CCSD(T)/cc-pVQZ//CCSD//cc-pVTZ are in somewhat better agreement with the gas-phase experiments.

Table 2. Relative free energies (kcal mol^-1) of 1-5.

References

(1) Alonso, J. L.; Vaquero, V.; Peña, I.; López, J. C.; Mata, S.; Caminati, W. "All Five Forms of Cytosine Revealed in the Gas Phase," Angew. Chem. Int. Ed. 2013, 52, 2331-2334, DOI: 10.1002/anie.201207744.

(2) Bazso, G.; Tarczay, G.; Fogarasi, G.; Szalay, P. G. "Tautomers of cytosine and their excited electronic states: a matrix isolation spectroscopic and quantum chemical study," Phys. Chem. Chem. Phys., 2011, 13, 6799-6807, DOI:10.1039/C0CP02354J.

InChIs

cytosine: InChI=1S/C4H5N3O/c5-3-1-2-6-4(8)7-3/h1-2H,(H3,5,6,7,8)
InChIKey=OPTASPLRGRRNAP-UHFFFAOYSA-N

Conformational analysis is one of the tasks that computation chemistry is typically quite adept at and computational chemistry is frequently employed for this purpose. Thus, benchmarking methods for their ability to predict accurate conformation energies is quite important. Martin has done this for alkanes¹ (see this post), and now he has looked at a molecule that contains weak intramolecular hydrogen bonds. He examined 52 conformations of melatonin 1.² The structures of the two lowest energy conformations are shown in Figure 1.

1

1a

1b

Figure 1. Structures of the two lowest energy conformers of 1 at SCS-MP2/cc-pVTZ.

The benchmark (i.e. accurate) relative energies of these conformers were obtained at MP2-F12/cc-pVTZ-F12 with a correction for the role of triples: (ECCSD(T)/cc-pVTZ)-E(MP2/cc-pVTZ)). The energies of the conformers were computed with a broad variety of basis sets and quantum methodologies. The root mean square deviation from the benchmark energies is used as a measure of the utility of these alternate methodologies. Of particular note is that HF predicts the wrong ordering of the two lowest energy isomers, as do some DFT methods that use small basis sets and do not incorporate dispersion.

In fact, other than the M06 family or double hybrid functionals, all of the functionals examined here (PBE. BLYP, PBE0, B3LYP, TPSS0 and TPSS) have RMSD values greater than 1 kcal mol^-1. However, inclusion of a dispersion correction, Grimme’s D2 or D3 variety or the Vydrov-van Voorhis (VV10) non-local correction (see this post for a review of dispersion corrections), reduces the error substantially. Among the best performing functionals are B2GP-PLYP-D3, TPSS0-D3, DSD-BLYP and M06-2x. They also find the MP2.5 method to be a practical ab initio alternative. One decidedly unfortunate result is that large basis sets are needed; DZ basis sets are simply unacceptable, and truly accurate performance requires a QZ basis set.

References

(1) Gruzman, D.; Karton, A.; Martin, J. M. L. "Performance of Ab Initio and Density Functional Methods for Conformational Equilibria of CnH2n+2 Alkane Isomers (n = 4-8)," J. Phys. Chem. A 2009, 113, 11974–11983, DOI: 10.1021/jp903640h.

(2) Fogueri, U. R.; Kozuch, S.; Karton, A.; Martin, J. M. L. "The Melatonin Conformer Space: Benchmark and Assessment of Wave Function and DFT Methods for a Paradigmatic Biological and
Pharmacological Molecule," J. Phys. Chem. A 2013, 117, 2269-2277, DOI: 10.1021/jp312644t.

InChIs

1: InChI=1S/C13H16N2O2/c1-9(16)14-6-5-10-8-15-13-4-3-11(17-2)7-12(10)13/h3-4,7-8,15H,5-6H2,1-2H3,(H,14,16)
InChIKey=DRLFMBDRBRZALE-UHFFFAOYSA-N

Truhlar has made a comparison of binding energies and relative energies of five (H₂O)₁₆ clusters.¹ While technically not organic chemistry, this paper is of interest to the readership of this blog as it compares a very large collection of density functionals on a problem that involves extensive hydrogen bonding, a problem of interest to computational organic chemists.

The CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVTZ energies of clusters 1-5 (shown in Figure 1) were obtained by Yoo.² These clusters are notable not just for their size but also that they involve multiple water molecules involved in four hydrogen bonds. Truhlar has used these geometries to compute the energies using 73 different density functionals with the jun-cc-pVTZ basis set (see this post for a definition of the ‘jun’ basis sets). Binding energies (relative to 16 isolated water molecules) were computed along with the 10 relative energies amongst the 5 different clusters. Combining the results of both types of energies, Truhlar finds that the best overall performance relative to CCSD(T) is obtained with ωB97X-D, a hybrid GGA method with a dispersion correction. The next two best performing functionals are LC-ωPBE-D3 and M05-2x. The best non-hybrid performance is with revPBE-D3 and B97-D.

1 (0.0)	2 (0.25)
3 (0.42)	4 (0.51)
5 (0.54)

Figure 1. MP2/aug-cc-pVTZ optimized geometries and relative CCSD(T) energies (kcal mol^-1) of (water)₁₆ clusters 1-5. (Don’t forget to click on any of these molecules above to launch Jmol to interactively view the 3-D structure. This feature is true for all molecular structures displayed in all of my blog posts.)

While this study can help guide selection of a functional, two words of caution. First, Truhlar notes that the best performing methods for the five (H₂O)₁₆ clusters do not do a particularly great job in getting the binding and relative energies of water hexamers, suggesting that no single functional really stands out as best. Second, a better study would also involve geometry optimization using that particular functional. Since this was not done, one can garner little here about what method might be best for use in a typical study where a geometry optimization must also be carried out.

References

(1) Leverentz, H. R.; Qi, H. W.; Truhlar, D. G. "Assessing the Accuracy of Density
Functional and Semiempirical Wave Function Methods for Water Nanoparticles: Comparing Binding and Relative Energies of (H₂O)₁₆ and (H₂O)₁₇ to CCSD(T) Results," J. Chem. Theor. Comput. 2013, ASAP, DOI: 10.1021/ct300848z.

(2) Yoo, S.; Aprà, E.; Zeng, X. C.; Xantheas, S. S. "High-Level Ab Initio Electronic Structure Calculations of Water Clusters (H₂O)₁₆ and (H₂O)₁₇: A New Global Minimum for (H₂O)₁₆," J. Phys. Chem. Lett. 2010, 1, 3122-3127, DOI: 10.1021/jz101245s.

This is a short post mainly to bring to the reader’s attention a couple of recent JCTC papers.

The first is a benchmark study by Hujo and Grimme of the geometries produced by DFT computations that are corrected for dispersion.¹ They use the S22 and S66 test sets that span a range of compounds expressing weak interactions. Of particular note is that the B3LYP-D3 method provided the best geometries, suggesting that this much (and justly) maligned functional can be significantly improved with just the simple D3 fix.

The second paper entails the description of Truhlar and Cramer’s latest iteration on their solvation model, namely SM12.² The main change here is the use of Hirshfeld-based charges, which comprise their Charge Model 5 (CM5). The training set used to obtain the needed parameters is much larger than with previous versions and allows for treating a very broad set of solvents. Performance of the model is excellent.

References

(1) Hujo, W.; Grimme, S. "Performance of Non-Local and Atom-Pairwise Dispersion Corrections to DFT for Structural Parameters of Molecules with Noncovalent Interactions," J. Chem. Theor. Comput. 2013, 9, 308-315, DOI: 10.1021/ct300813c

(2) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. "Generalized Born Solvation Model SM12," J. Chem. Theor. Comput. 2013, 9, 609-620, DOI: 10.1021/ct300900e

Two interesting questions are addressed in a focal-point computational study of t-butyl radical and the t-butyl anion from the Schaefer group.¹ First, is the radical planar? EPR and PES studies from the 1970s indicate a pyramidal structure, with an inversion barrier of only 0.64 kcal mol^-1. The CCSD(T)/cc-pCVTZ optimized structure of t-butyl radical shows it to be pyramidal with the out-of-plane angle formed by one methyl group and the other three carbons of 22.9°, much less than the 54.7° of a perfect tetrahedron. Focal point analysis give the inversion barrier 0.74 kcal mol^-1, in outstanding agreement with experiment.

Second, what is the electron affinity (EA) of the t-butyl radical? Schleyer raised the concern that the alkyl anions may be unbound, and suggested that the electron affinity of t-butyl radical was -9.6 kcal mol^-1; in other words, the anion is thermodynamically unstable. This focal-point study shows just how sensitive the EA is to computational method. The HF/CBS value of the EA is -39.59 kcal mol^-1 (unbound anion), but the MP2/CBS value is +41.57 kcal mol^-1 (bound anion!). The CCSD/aug-cc-pVQZ value is -8.92 while the CCSD(T)/aug-cc-pVQZ value is +4.79 kcal mol^-1. The estimated EA at CCSDT(Q)/CBS is -1.88 kcal mol^-1, and inclusion of correction terms (including ZPE and relativistic effect) gives a final estimate of the EA as -0.48 kcal mol^-1, or a very weakly unbound t-butyl anion. It is somewhat disconcerting that such high-level computations are truly needed for some relatively simple questions about small molecules.

References

(1) Sokolov, A. Y.; Mittal, S.; Simmonett, A. C.; Schaefer, H. F. "Characterization of the t-Butyl Radical and Its Elusive Anion," J. Chem. Theory Comput. 2012, 8, 4323-4329, DOI: 10.1021/ct300753d.

A comprehensive evaluation of how different computational methods perform in predicting the energies of monosaccharides comes to some very interesting conclusions. Sameera and Pantazis¹ have examined the eight different aldohexoses (allose, alltrose, glucose, mannose, gulose, idose, galactose and talose), specifically looking at different rotomers of the hydroxymethyl group, α- vs. β-anomers, pyranose vs. furanose isomers, ring conformations (¹C₄ vs skew boat forms), and ring vs. open chain isomers. In total, 58 different structures were examined. The benchmark computations are CCSD(T)/CBS single point energies using the SCS-MP2/def2-TZVPP optimized geometries. The RMS deviation from these benchmark energies for some of the many different methods examined are listed in Table 1.

Table 1. Average RMS errors (kJ mol^-1) of the 58 different monosaccharide structures for
different computational methods.

Perhaps the most interesting take-home message is that CEPA, MP2, the double hybrid methods and M06-2x all do a very good job at evaluating the energies of the carbohydrates. Given the significant computational advantage of M06-2x over these other methods, this seems to be the functional of choice! The poorer performance of the DFT methods over the ab initio methods is primarily in the relative energies of the open-chain isomers, where errors can be on the order of 10-20 kJ mol^-1 with most of the functionals; even the best overall methods (M06-2x and the double hybrids) have errors in the relative energies of the open-chain isomers of 7 kJ mol^-1. This might be an area of further functional development to probe better treatment of the open-chain aldehydes vs. the ring hemiacetals.

References

(1) Sameera, W. M. C.; Pantazis, D. A. "A Hierarchy of Methods for the Energetically Accurate Modeling of Isomerism in Monosaccharides," J. Chem. Theory Comput. 2012, 8, 2630-2645, DOI:10.1021/ct3002305

I just returned from the Southwest Theoretical Chemistry Conference held at Texas A&M University. My thanks again to Steven Wheeler for the invitation to speak at the meeting and for putting together a very fine program and conference.

Among the many interesting talks was one by Sebastian Kozuch who reported on an interesting double hybrid methodology.^1,2 Working with Jan Martin, they defined a procedure that Kozuch referred to as “putting Stefan Grimme into a blender”. They extend the double hybrid concept first suggested by Grimme that adds on an MP2-like correction functional. Kozuch and Martin substitute a spin-component scaled MP2 (SCS-MP2) model for the original MP2 correction. SCS-MP2 was also proposed by Grimme. Lastly, they add on a dispersion correction, an idea championed by Grimme too. The exchange-correlation term is defined as

E_XC = c_XE_X ^DFT + (1 – c_x)E_x^HF + c_CE_C^DFT + c_OE_O^MP2 + c_SE_S^MP2 + s₆E_D

where c_X is the coefficient for the amount of DFT exchange, c_C the amount of DFT correlation, c_C and c_S the amount of opposite- and same-spin MP2, and s₆ the amount of dispersion. They name this procedure DSD-DFT for Dispersion corrected, Spin-component scaled Double hybrid DFT.

In their second paper on this subject, they propose the use of the PBEP86 functional for the DFT components.² Benchmarking against a variety of standard databases, including kinetic data, thermodynamic data, along with inorganic and weakly interacting systems, this method delivers the lowest mean error among a small set of functionals. Kozuch reported at the conference on a number of other combinations and should have a publication soon suggesting an even better method. Importantly, these DSD-DFT computations can be run with most major quantum codes including Orca, Molpro, Q-Chem and Gaussian (with a series of IOP specifications).

While double hybrid methods don’t have quite the performance capabilities of regular DFT, density fitting procedures offer the possibility of a significant reduction in computational time. These DSD-DFT methods are certainly worthy of fuller explorations.

References

(1) Kozuch, S.; Gruzman, D.; Martin, J. M. L. "DSD-BLYP: A General Purpose Double Hybrid Density Functional Including Spin Component Scaling and Dispersion Correction," J. Phys. Chem. C, 2010, 114, 20801-20808,
DOI: 10.1021/jp1070852

(2) Kozuch, S.; Martin, J. M. L. "DSD-PBEP86: in search of the best double-hybrid DFT with spin-component scaled MP2 and dispersion corrections," Phys. Chem. Chem. Phys. 2011, 13, 20104-20107, DOI: 10.1039/C1CP22592H