Here’s one more attempt to discern the failure of DFT to handle simple alkanes (see this earlier post for a previous attempt to answer this question). Tsuneda and co-workers¹ have employed long-range corrected (LC) DFT to the problem of the energy associated with “protobranching”, i.e., from the reaction

CH₃(CH₂)­_nCH₃ + n CH₄ → (n+1) CH₃CH₃

They computed the energy of this reaction for the normal alkanes propane through decane using a variety of functionals, and compared these computed values with experimentally-derived energies. Table 1 gives the mean unsigned error for a few of the functionals. The prefix “LC” indicated inclusion of long-range corrections, “LCgau” indicates the LC scheme with a gaussian attenuation, and “LRD” indicates inclusion of long-range dispersion.

Table 1. Mean unsigned errors of the “protobranching”
reaction energy of various functional compared to experiment.

A number of important conclusions can be drawn. First, with both LC and LRD very nice agreement with experiment can be had. If only LC is included, the error increases on average by over 1 kcal mol^-1. The MO6-2x functional, touted as a fix of the problem, does not provide complete correction, though it is vastly superior to B3LYP and other hybrid functionals. The authors conclude that the need for LC incorporation points out that the exchange functional lacks the ability to account for this effect. Medium-range correlation is not the main source of the problem as large discrepancies in the reaction energy error occur when different functionals are used that are corrected for LC and LRD. Choice of functional still matters, but LC correction appears to be a main culprit and further studies of its addition to standard functionals would be most helpful.

References

(1) Song, J.-W.; Tsuneda, T.; Sato, T.; Hirao, K., "Calculations of Alkane Energies Using Long-Range Corrected DFT Combined with Intramolecular van der Waals Correlation," Org. Lett. 2010, 12, 1440–1443, DOI: 10.1021/ol100082z

Houk has performed a very nice examination of the performance of some density functionals.¹ He takes a quite different approach than what was proposed by Grimme – the “mindless” benchmarking² using random molecules (see this post). Rather, Houk examined a series of simple aldol, Mannich and α-aminoxylation reactions, comparing their reaction energies predicted with DFT against that predicted with CBQ-QB3. The idea here is to benchmark DFT performance for simple reactions of specific interest to organic chemists. These reactions are of notable current interest due their involvement in organocatalytic enantioselective chemistry (see my posts on the aldol, Mannich, and Hajos-Parrish-Eder-Sauer-Wiechert reaction). Examples of the reactions studied (along with their enthalpies at CBS-QB3) are Reaction 1-3.

For the four simple aldol reactions and four simple Mannich reactions, PBE1PBE,
mPW1PW91 and MO6-2X all provided reaction enthalpies with errors of about 2 kcal mol^-1. The much maligned B3LYP functional, along with B3PW91 and B1B95 gave energies with significant larger errors. For the three α-aminoxylation reactions, the errors were better with B3PW91 and B1B95 than with PBE1PBE or MO6-2X. Once again, it appears that one is faced with finding the right functional for the reaction under consideration!

Of particular interest is the decomposition of these reactions into related isogyric, isodesmic
and homdesmic reactions. So for example Reaction 1 can be decomposed into Reactions 4-7 as shown in Scheme 1. (The careful reader might note that these decomposition reactions are isodesmic and homodesmotic and hyperhomodesmotic reactions.) The errors for Reactions 4-7 are typically greater than 4 kcal mol^-1 using B3LYP or B3PW91, and even with MO6-2X the errors are about 2 kcal mol^-1.

Scheme 1.

Houk also points out that Reactions 4, 8 and 9 (Scheme 2) focus on having similar bond changes as in Reactions 1-3. And it’s here that the results are most disappointing. The errors produced by all of the functionals for Reactions 4,8 and 9 are typically greater than 2 kcal mol^-1, and even MO2-6x can be in error by as much as 5 kcal mol^-1. It appears that the reasonable performance of the density functionals for the “real world” aldol and Mannich reactions relies on fortuitous cancellation of errors in the underlying reactions. Houk calls for the development of new functionals designed to deal with fundamental simple bond changing reactions, like the ones in Scheme 2.

Scheme 2

References

(1) Wheeler, S. E.; Moran, A.; Pieniazek, S. N.; Houk, K. N., "Accurate Reaction Enthalpies and Sources of Error in DFT Thermochemistry for Aldol, Mannich, and α-Aminoxylation Reactions," J. Phys. Chem. A 2009, 113, 10376-10384, DOI: 10.1021/jp9058565

(2) Korth, M.; Grimme, S., ""Mindless" DFT Benchmarking," J. Chem. Theory Comput. 2009, 5, 993–1003, DOI: 10.1021/ct800511q

Predicting the pK_a of a compound from first principles remains a challenge, despite all of the many algorithmic and methodological advantages within the sphere of computational chemistry. Predicting the gas-phase deprotonation energy is relatively straightforward as I detail in Section 2.2 of my book. The difficulty is in treating the solvent and the interaction of the acid and its conjugate base in solution. Considering that we are most interested in acidities in water, a very polar solvent, the interactions between water and the conjugate base and the proton are likely to be large and important!

Baker and Pulay report a procedure for determining acidities with the aim of high throughput.¹ Thus, computational efficiency is a primary goal. Their approach is to compute the enthalpy change for deprotonation in solution using a continuum treatment and then employ a linear fit to predict the pK_a with the equation:

pK_a(c) = α_cΔH + β_c

where c designates a class of compound, such as alcohol, carboxylic acid, amine, etc. Fitting constants α_c and β_c need to be found then for each unique class of compound, where the fitting is to experimental pK_as in water. In their test suite, they employed eleven anilines and amines, seven pyridines, nine alcohols and phenols, and seven carboxylic acids.

They test a number of different computational variants: (a) what functional to employ, (b) what basis set to use for optimizing structures, and (c) what basis set to use for the enthalpy computation. They opt to employ COSMO for treating the solvent and quickly reject the use of gas phase structures (and particularly use of geometries obtained with molecular mechanics. Their ultimate model is OLYP/6-311+G**//3-21G(d) with the COSMO solvation model. Mean deviation is less than 0.4 pK units. They do note that use of HF or PW91 provides similar small errors, but ultimately favor OLYP for its computational performance.

While this procedure offers some guidance for future computation of acidity, I find a couple of issues. First, it relies on fitted parameters for every class of compound. If one is interested in a new class, then one must develop the appropriate parameters – and the experimental values may not be available or perhaps an insufficient number of them are experimentally available. Second, the parameters cover-up a great deal of problematic computational sins, like the solvation energy of the proton, small basis sets, missing correlation energies, missing dispersion corrections etc. A purist might hope for a computational algorithm that allows for systematic correction and improvement in the estimation of pK_as. Further work needs to be done to meet this higher goal.

References

(1) Zhang, S.; Baker, J.; Pulay, P., "A Reliable and Efficient First Principles-Based Method for Predicting pKa Values. 1. Methodology," J. Phys. Chem. A 2010, 114 , 425-431, DOI: 10.1021/jp9067069

Another benchmark study of the performance of different functionals – this time looking at the conformations of small alkanes.¹ Martin first establishes high level benchmarks: the difference between the trans and gauche conformers of butane: CCSD(T)/cc-pVQZ, 0.606 kcal mol^-1 and W1h-val, 0.611 kcal mol^-1; and the energy differences of the conformers of pentane, especially the TT and TG gap: 0.586 kcal mol^-1 at CCSD(T)/cc-pVTZ and 0.614 kcal mol^-1 at W1h-val.

They then examine the relative conformational energies of butane, pentane, hexane and a number of branched alkanes with a slew of functionals, covering the second through fifth rung of Perdew’s Jacob’s ladder. The paper has a whole lot of data – and the supporting
materials include Jmol-enhanced visualization of the structures! – but the bottom line is the following. The traditionally used functionals (B3LYP, PBE, etc) overestimate conformer energies while the MO6 family underestimates the interaction energies that occur in GG-type conformers. A dispersion correction tends to overcorrect and leads to wrong energy ordering of conformers. But the new double-hybrid functionals (B2GP-PLYP and B2K-PLYP) with the dispersion correction provide quite nice agreement with the CCSD(T) benchmarks.

Also worrisome is that all the functionals have issues in geometry prediction, particularly in the backbone dihedral angles. So, for example, B3LYP misses the τ₁ dihedral angle in the GG conformer by 5° and even MO6-2x misses the τ₂ angle in the TG conformer by 2.4&deg.

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),"
The Journal of Physical Chemistry A 2009, 113, 11974–11983 , DOI: http://dx.doi.org/10.1021/jp903640h

Here’s another extensive benchmarking study – this time on the use of TD-DFT to predict excitation energies.¹ This study looks at the performance of 28 different functionals, and compares the TD-DFT excitation energies against a data set of (a) computed vertical energies and (b) experimental energies. The performance is generally about the same for both data sets, with many functionals (especially the hybrid functionals) giving errors of about 0.3 eV. Performance can be a bit better when examining subclasses of compounds. For example, PBE0 and mPW1PW91 have a mean unsigned error of only 0.14 eV for a set of organic dyes.

References

(1) Jacquemin, D.; Wathelet, V.; Perpete, E. A.; Adamo, C., "Extensive TD-DFT Benchmark: Singlet-Excited States of Organic Molecules," J. Chem. Theory Comput., 2009, 5, 2420-2435, DOI: 10.1021/ct900298e

Here’s another fine paper from the Alonso group employing laser ablation molecular beam Fourier transform microwave spectroscopy coupled with computation to discern molecular structure. In this work they examine the low-energy tautomers of guanine.1 The four lowest energy guanine tautomers are shown in Figure 1. (Unfortunately, Alonso does not include the optimized coordinates of these structures in the supporting information – we need to more vigorously police this during the review process!) These tautomers are predicted to be very close in energy (MP2/6-311++G(d,p), and so one might expect to see multiple signals in the microwave originating from all four tautomers. In fact, they discern all four, and the agreement between the computed and experimental rotational constants are excellent (Table 1), especially if one applies a scaling factor of 1.004. Once again, this group shows the power of combined experiment and computations!

Figure 1. Four lowest energy (kcal mol^-1, MP2/6-311++G(d,p)) tautomers of guanine.

Table 1. Experimental and computed rotational constants (MHz) of the four guanine tautomers.

References

(1) Alonso, J. L.; Peña, I.; López, J. C.; Vaquero, V., "Rotational Spectral Signatures of Four Tautomers of Guanine," Angew. Chem. Int. Ed. 2009, 48, 6141-6143, DOI: 10.1002/anie.200901462

InChIs

Guanine: InChI=1/C5H5N5O/c6-5-9-3-2(4(11)10-5)7-1-8-3/h1H,(H4,6,7,8,9,10,11)/f/h8,10H,6H2
InChIKey=UYTPUPDQBNUYGX-GSQBSFCVCX

I have written a number of blog posts that deal with the computation of optical activity. Trindle and Altun have now reported TD-DFT computations of circular dichroism of high-symmetry molecules.¹ The employ either B3LYP (with a variety of basis sets, the largest being 6-311++G(2d,2p)) and SOAP/ATZP. For a number of the high symmetry molecules (two examples are shown in Figure 1), the two methods differ a bit in their predictions of the first excited state, with SOAP typically predicting a red shift relative to the B3LYP. However, both methods general give the same sign of the CD signals and their line shapes are similar.

1	2

Figure 1. B3LYP/6-31G(d) optimized structures of 1 and 2 (again due to incomplete supporting materials, I reoptimized these structures)

References

(1) Trindle, C.; Altun, Z., "Circular dichroism of some high-symmetry chiral molecules: B3LYP and SAOP calculations " Theor. Chem. Acc. 2009, 122, 145-155, DOI: 10.1007/s00214-008-0494-8.

InChIs

1: InChI=1/C18H14O2/c19-15-7-11-3-1-4-12-8-16(20)10-14-6-2-5-13(9-15)18(14)17(11)12/h1-6H,7-10H2
InChIKey=DYZSIUYFWKNLHS-UHFFFAOYAB

2: InChI=1/C20H24/c1-13-9-18-7-8-20-12-15(3)19(11-16(20)4)6-5-17(13)10-14(18)2/h9-12H,5-8H2,1-4H3
InChIKey=JTMLLDPOLFRPGJ-UHFFFAOYAC

DFT & Optical Rotation Steven Bachrach 27 Jul 2009 No Comments

I (0.0)	II (4.79)
III (5.81)	IV (5.95)

Figure 1. MP2(FC)/aug-cc-pV(T+d)Z optimized geometries and focal point relative energies (kJ mol^-1) of the four lowest energy conformers of cysteine.¹

The three lowest energy structures found here match up with the lowest two structures found by Alonso and the energy differences are also quite comparable: 4.79 kJ and 5.81 mol^-1 with the focal point method 3.89 and 5.38 kJ mol^-1 with MP4/6-311++G(d,p)// MP2/6-311++G(d,p). So the identification of the cysteine conformers made by Alonso remains on firm ground.

References

(1) Wilke, J. J.; Lind, M. C.; Schaefer, H. F.; Csaszar, A. G.; Allen, W. D., "Conformers of Gaseous Cysteine," J. Chem. Theory Comput. 2009, DOI: 10.1021/ct900005c.

(2) Sanz, M. E.; Blanco, S.; López, J. C.; Alonso, J. L., "Rotational Probes of Six Conformers of Neutral Cysteine," Angew. Chem. Int. Ed. 2008, 4, 6216-6220, DOI: 10.1002/anie.200801337

InChIs

Cysteine:
InChI=1/C3H7NO2S/c4-2(1-7)3(5)6/h2,7H,1,4H2,(H,5,6)/t2-/m0/s1
InChIKey: XUJNEKJLAYXESH-REOHCLBHBU

Back when I was first learning ab initio methods in Cliff Dykstra’s lab, I played a bit with the post-HF method CEPA (couple electron pair approximation). This method fell out of favor over the years with the rise of MP theory and then with DFT. Now, Neese and Grimme and co-workers are resurrecting it.¹ Their Accounts article provides a series of tests of CEPA/1 against benchmark computations (typically CCSD(T)) and lo and behold, CEPA performs remarkably well! It bests B3LYP (no surprise there), B2LYP and MP2 in virtually every category, ranging from reaction energies, hydrogen bond energies, van der Waals interaction energies, and activation barrier heights. As an example, for the isomerization energy of toluene to norbornadiene, CCSD(T) estimates the energy is 42.79 kcal mol^-1. B3LYP does miserably, with an error of nearly 14 kcal mol^-1, but the CEPA/1 estimate is off by only 0.04 kcal mol^-1. Since the computational time of CEPA/1 is competitive with MP2, the authors conclude that CEPA/1 is well-worth reinvestigating as an alternative post-HF methodology.

References

(1) Neese, F.; Hansen, A.; Wennmohs, F.; Grimme, S., "Accurate Theoretical Chemistry with Coupled Pair Models," Acc. Chem. Res. 2009, 42, 641-648 DOI: 10.1021/ar800241t.