Archive for the 'QM Method' Category

Benchmarking π-conjugation

With the proliferation of density functionals, selecting the functional to use in your particular application requires some care. That is why there have been quite a number of benchmark studies (see these posts for some examples). Yu and Karton have now added to our benchmark catalog with a study of π-conjugation.1

They looked at a set of 60 reactions which involve a reactant with π-conjugation and a product which lacks conjugation. A few examples, showing examples involving linear and cyclic systems, are shown in Scheme 1.

Scheme 1.

The reaction energies were evaluated at W2-F12, which should have an error of a fraction of a kcal mol-1. Three of the reactions can be compared with experimental values, and difference in the experimental and computed values are well within the error bars of the experiment. It is too bad that the authors did not also examine 1,3-cyclohexadiene → 1,4-cyclohexadiene, a reaction that is both of broader interest than many of the ones included in the test set and can also be compared with experiment.

These 60 reactions were then evaluated with a slew of functionals from every rung of Jacob’s ladder. The highlights of this benchmark study are that most GGA and meta-GGA and hybrid functionals (like B3LYP) have errors that exceed chemical accuracy (about 1 kcal mol-1). However, the range-separated functionals give very good energies, including ωB97X-D. The best results are provided with double hybrid functionals. Lastly, the D3 dispersion correction does generally improve energies by 10-20%. On the wavefunction side, SCS-MPs gives excellent results, and may be one of the best choices when considering computational resources.


(1) Yu, L.-J.; Karton, A. "Assessment of theoretical procedures for a diverse set of isomerization reactions involving double-bond migration in conjugated dienes," Chem. Phys. 2014, 441, 166-177, DOI: 10.1016/j.chemphys.2014.07.015.

DFT Steven Bachrach 14 Jan 2015 2 Comments

Becke’s Perspective on DFT

The Journal of Chemical Physics has produced a Special Topics issue on Advances in Density Functional Theory. I want to call to your attention the Perspective article by Becke titled “Perspective: Fifty years of density-functional theory in chemical physics”.1 Becke writes a personal account of the history of DFT and makes a number of interesting points and observations. He rightly notes that DFT is exact and we should more properly refer to our actual implementations as Density Functional Approximations (DFA). He also notes that use of the term ab initio as a synonym for wavefunction theory is inappropriate as DFT is just as ab initio as HF and post-HF theories.

A common perception about DFT (well, DFA) is that there is no way to systematically improve functionals. Becke exposes a true underlying logic that has driven much of DFA development.

Lastly, Becke is discouraged by the more recent developments that have included virtual orbitals, such as double hybrid methods. His approach is that true DFT is occupied orbitals only (for which he pointedly does not want to adopt the acronym OOO), and that developments that include the virtual orbitals might toll the “death knell” for DFT.

For those interested in a pretty accessible account of the history of DFT, Becke’s Perspective is an excellent place to get started.


(1) Becke, A. D. "Perspective: Fifty years of density-functional theory in chemical physics," J. Chem. Phys. 2014, 140, 18A301 DOI: 10.1063/1.4869598.

DFT Steven Bachrach 24 Sep 2014 No Comments

Structure of benzene dication

Benzene is certainly one of the most iconic chemical compounds – its planar hexagonal structure is represented often in popular images involving chemists, and its alternating single and double bonds the source of one of chemistry’s most mythic stories: Kekule’s dream of a snake biting its own tail. So while the structure of benzene is well-worn territory, what of the structure of the benzene dication? Jasik, Gerlich and Rithova probe that question using a combined experimental and computational approach.1

The experiment involves generation of the benzene dication at low temperature and complexed
to helium. Then, using infrared predissociation spectroscopy (IRPD), they obtained a spectrum that suggested two different structures.

Next, employing MP2/aug-cc-pVTZ computations, they identified a number of possible geometries, and the two lowest energy singlet dications have the geometries shown in Figure 1. The first structure (1) has a six member ring, but the molecule is no longer planar. Lying a bit lower in energy is 2, having a pentagonal pyramid form. The combination of the computed IR spectra of each of these two structures matches up extremely well with the experimental spectrum.



Figure 1. MP2/aug-cc-pVTZ geometries of benzene dication 1 and 2.


(1) Jašík, J.; Gerlich, D.; Roithová, J. "Probing Isomers of the Benzene Dication in a Low-Temperature Trap," J. Am. Chem. Soc. 2014, 136, 2960-2962, DOI: 10.1021/ja412109h.

Aromaticity &MP Steven Bachrach 08 Apr 2014 2 Comments

Computing OR: norbornenone

Optical activity is a major tool for identifying enantiomers. With recent developments in computational techniques, it is hoped that experiments combined with computations will be a powerful tool for determining absolute configuration. The recent work of Lahiri, et al. demonstrates the scope of theoretical work that is still needed to really make this approach broadly applicable.1

They prepared (1R,4R)-norbornenone 1 and measured its optical rotation in the gas phase and in dilute solutions of acetonitrile and cyclohexane. The specific rotations at three different wavelengths are listed in Table 1. Of first note is that though there is some small differences in solution, as expected, there really is substantial differences between the gas- and solution phases. Thus cautionary point 1: be very careful of comparing solution phase experimental optical activity with computed gas phase predictions.


Table 1. Experimental and computed specific rotation of 1.


355.0 nm

589.3 nm

633.0 nm

Gas phase













Acetonitrile solution









Cyclohexane solution









For the computations, the geometry of 1 was optimized at B3LYP/aug-cc-pVTZ (see Figure 1. The OR was computed at B3LYP with different basis sets, finding that the difference in the predicted specific rotation at 598.3nm differs only quite little (90.6 deg dm-1 (g/mL)-1) between the computations using aug-cc-pVTZ and aug-cc-pVQZ) suggesting that the basis set limit has been reached. The gas –phase computed values at B3LYP and CCSD are also shown in Table 1. Though these computations do get the correct sign of the rotation, they are appreciably off of the experimental values, and the percent error varies with wavelength. Cautionary point 2: it is not obvious what is the appropriate computational method to compute OR, and beware of values that seem reasonable at one wavelength – this does not predict a good agreement at a different wavelength.

Figure 1. Optimized geometry of 1 at B3LYP/aug-cc-pVTZ.

Lastly, computed solution values of the OR were performed with PCM and B3LYP. These values are given in Table 1. Again the agreement is poor. So cautionary point 3: computed (PCM) solution OR
may be in quite poor agreement with experiment.

Often the culprit to poor agreement between computed and experimental OR is attributed to omitted vibrational effects, but in this case, because 1 is so rigid, one might not expect too much error to be caused by the effects of vibrations. So the overall result – we need considerable work towards addressing how to accurately compute optical activity!


(1) Lahiri, P.; Wiberg, K. B.; Vaccaro, P. H.; Caricato, M.; Crawford, T. D. "Large Solvation Effect in the Optical Rotatory Dispersion of Norbornenone," Angew. Chem. Int. Ed. 2014, 53, 1386-1389, DOI: 10.1002/anie.201306339.


1: InChI=1S/C7H8O/c8-7-4-5-1-2-6(7)3-5/h1-2,5-6H,3-4H2/t5-,6+/m1/s1

DFT &Optical Rotation Steven Bachrach 25 Feb 2014 2 Comments

Gas phase structure of 2-deoxyribose

2-deoxyribose 1 is undoubtedly one of the most important sugars as it is incorporated into the backbone of DNA. The conformational landscape of 1 is complicated: it can exist as an open chain, as a five-member ring (furanose), or a six-member ring (pyranose), and intramolecular hydrogen bonding can occur. This internal hydrogen bonding is in competition with hydrogen bonding to water in aqueous solution. Unraveling all this is of great interest in predicting structures of this and a whole host of sugar and sugar containing-molecules.


In order to get a firm starting point, the gas phase structures of the low energy conformers of 1 would constitute a great set of structures to use as a benchmark for gauging force fields and computational methods. Cocinero and Alonso1 have performed a laser ablation molecular beam Fourier transform microwave (LA-MB-FTMW) experiment (see these posts for other studies using this technique) on 1 and identified the experimental conformations by comparison to structures obtained at MP2/6-311++G(d,p). Unfortunately the authors do not include these structures in their supporting materials, so I have optimized the low energy conformers of 1 at ωB97X-D/6-31G(d) and they are shown in Figure 1.

1a (0.0)

1b (4.7)

1c (3.3)

1d (5.6)

1e (8.9)

1f (9.4)

Figure 1. ωB97X-D/6-31G(d) optimized structures of the six lowest energy conformers of 1. Relative free energy in kJ mol-1.

The computed spectroscopic parameters were used to identify the structures responsible for the six different ribose conformers observed in the microwave experiment. To give a sense of the agreement between the computed and experimental parameters, I show these values for the two lowest energy conformers in Table 1.

Table 1. MP2/6-311++G(d,p) computed and observed spectroscopic parameters for the two lowest energy conformers of 1.














B (MHz)





C (MHz)





ΔG (kJ mol-1)





This is yet another excellent example of the symbiotic relationship between experiment and computation in structure identification.


(1) Peña, I.; Cocinero, E. J.; Cabezas, C.; Lesarri, A.; Mata, S.; Écija, P.; Daly, A. M.; Cimas, Á.; Bermúdez, C.; Basterretxea, F. J.; Blanco, S.; Fernández, J. A.; López, J. C.; Castaño, F.; Alonso, J. L. "Six Pyranoside Forms of Free 2-Deoxy-D-ribose," Angew. Chem. Int. Ed. 2013, 52, 11840-11845, DOI: 10.1002/anie.201305589.


1a: InChI=1S/C5H10O4/c6-3-1-5(8)9-2-4(3)7/h3-8H,1-2H2/t3-,4+,5-/m0/s1

MP &sugars Steven Bachrach 16 Dec 2013 No Comments

Extrapolated CCSD(T) Thermochemistry

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

ΔEM1 = cEM1(C) + dEM1(D) – aEM1(A) – bEM1(B)
ΔEM2 = cEM2(C) + dEM2(D) – aEM2(A) – bEM2(A)

Now, if the two computational methods are reasonably complete, then ΔEM1 ≈ ΔEM2. 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



cEM1(C) + dEM1(D) – aEM1(A) – bEM1(B) ≈ cEM2(C) + dEM2(D) – aEM2(A) – bEM2(B)

cEM1(C) ≈ cEM2(C) + d[EM2(D) – EM1(D)] – a[EM2(A) – EM1(A)] – b[EM2(B) – EM1(B)]

and so

cEM1(C) ≈ cEM2(C) + Σ ci(EM2(i) – EM1(i))

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.1

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.2 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.3)

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.


(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.

QM Method Steven Bachrach 02 Dec 2013 2 Comments

Is CCSD(T)/CBS really the gold standard?

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 Hobza1 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 waterammonia complex), dispersion (e.g. in the methane dimer and the methaneethane 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 (, but I could not find it there. Any help?


(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.

QM Method Steven Bachrach 28 May 2013 3 Comments

ORD of methyloxirane

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 Barone1 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.


(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.



DFT &Optical Rotation Steven Bachrach 15 May 2013 2 Comments

Gas-phase structure of cytosine

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.1 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 computations2 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.


























(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.


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

MP &nucleic acids Steven Bachrach 22 Apr 2013 1 Comment

Benchmarking conformations: melatonin

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 alkanes1 (see this post), and now he has looked at a molecule that contains weak intramolecular hydrogen bonds. He examined 52 conformations of melatonin 1.2 The structures of the two lowest energy conformations are shown in Figure 1.




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.


(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.


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)

DFT &MP Steven Bachrach 11 Apr 2013 2 Comments

Next Page »