Search Results for "benchmark"

More DFT benchmarks – sugars and “mindless” test sets

Another two benchmarking studies of the performance of DFT have appeared.

The first is an examination by Csonka and French of the ability of DFT to predict the relative energy of carbohydrate conformation energies.1 They examined 15 conformers of α- and β-D-allopyranose, fifteen conformations of 3,6-anydro-4-O-methyl-D-galactitol and four conformers of β-D-glucopyranose. The energies were referenced against those obtained at MP2/a-cc-pVTZ(-f)//B3LYP/6-31+G*. (This unusual basis set lacks the f functions on heavy atoms and d and diffuse functions on H.) Among the many comparisons and conclusions are the following: B3LYP is not the best functional for the sugars, in fact all other tested hybrid functional did better, with MO5-2X giving the best results. They suggest the MO5-2X/6-311+G**//MO5-2x/6-31+G* is the preferred model for sugars, except for evaluating the difference between 1C4 and 4C1 conformers, where they opt for PBE/6-31+G**.

The second, by Korth and Grimme, describes a “mindless” DFT benchmarking study.2
This is really not a “mindless” study (as the term is used by Schaefer and Schleyer3 and discussed in this post, where all searching is done in a totally automated way) but rather Grimme describes a procedure for removing biases in the test set. Selection of “artificial molecules” is made by first deciding how many atoms are to be present and what will be the distribution of elements. In their two samples, they select systems having 8 atoms. The two sets differ by the distribution of the elements. The first set the atoms Na-Cl are one-third as probable as the elements Li-F, which are one-third as probable as H. The second set has the probability distribution similar to those found in naturally occurring organic compounds. The eight atoms, randomly selected by the computer, are placed in the corners of a cube and allowed to optimize (this is reminiscent of the “mindless” procedure of Schaefer and Schleyer3). This process generates a selection of random bonding environments along with open- and closed shell species, and removes (to a large degree) the biases of previous test sets, which are often skewed towards small molecules, ones where accurate experiments are available or geared towards a select group of molecules of interest. Energies are then computed using a variety of functional and compared to the energy at CCSD(T)/CBS. The bottom line is that the functional nicely group along the rungs defined by Perdew:4 LDA is the poorest performer, GGA does much better, the third rung of meta-GGA functionals are slightly better than GGA functionals, hybrids do better still, and the fifth rung functionals (double hybrids) perform quite well. Also of interest is that CCSD(T)/cc-pVDZ gives quite large errors and so Grimme cautions against using this small basis set.


(1) Csonka, G. I.; French, A. D.; Johnson, G. P.; Stortz, C. A., "Evaluation of Density Functionals and Basis Sets for Carbohydrates," J. Chem. Theory Comput. 2009, ASAP, DOI: 10.1021/ct8004479.

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

(3) Bera, P. P.; Sattelmeyer, K. W.; Saunders, M.; Schaefer, H. F.; Schleyer, P. v. R., "Mindless Chemistry," J. Phys. Chem. A 2006, 110, 4287-4290, DOI: 10.1021/jp057107z.

(4) Perdew, J. P.; Ruzsinszky, A.; Tao, J.; Staroverov, V. N.; Scuseria, G. E.; Csonka, G. I., "Prescription for the design and selection of density functional approximations: More constraint satisfaction with fewer fits," J. Chem. Phys. 2005, 123, 062201-9, DOI: 10.1063/1.1904565

DFT &Grimme Steven Bachrach 21 Apr 2009 3 Comments

SN2 and E2 DFT benchmark

Bickelhaupt has reported a broad benchmark study of the prototype SN2 and E2 reactions.1 These are the reactions of ethyl fluoride with fluoride and ethyl chloride with chloride (Scheme A). The critical points were optimized at OLYP/TZ2P and then CCSD(T)/CBS energies are used as benchmark. A variety of different density functionals were then used to obtain single-point energies.

Scheme A

The relative energies of the transition states for the six different reactions are listed for some of the functionals in Table 1. (These are energies relative to separated reactants – and keep in mind that an ion dipole complex is formed between the reactants and the transition states – Bickelhaupt calls this a “reaction complex”.)

Table 1. Relative energies (kcal mol-1) of the transition states for the six reactions shown in Scheme A.






E2 anti

E2 syn


E2 anti

E2 syn


















































There is a lot more data in this paper, along with a summary of the mean absolute errors in the overall and central barriers that mimics the data I show in Table 1. The trends are pretty clear. Generalized gradient approximation (GGA) functions – like BLYP, PW91, and PBE – dramatically underestimate the barriers. The hybrid functionals perform much better. The recently maligned B3LYP functional gets the correct trend and provides reasonable estimates of the barriers. Truhlar’s MO5-2X and MO6-2X functionals do very well in matching up the barrier heights along with getting the correct trends in the relative barriers. Simply looking for the functional with the lowest absolute error is not sufficient; BHandH and MO6-L have small errors but give a wrong trend in barriers, predicting that the SN2 reaction is preferred over the E2 for the fluoride reaction.


(1) Bento, A. P.; Sola, M.; Bickelhaupt, F. M., "E2 and SN2 Reactions of X + CH3CH2X (X = F, Cl); an ab Initio and DFT Benchmark Study," J. Chem. Theory Comput., 2008, 4, 929-940, DOI: 10.1021/ct700318e.

DFT &Substitution Steven Bachrach 22 Jul 2008 No Comments

Using vibrational frequencies to identify stereoisomers

Can vibrational spectroscopy be used to identify stereoisomers? Medel, Stelbrink, and Suhm have examined the vibrational spectra of (+)- and (-)-α-pinene, (±)-1, in the presence of four different chiral terpenes 2-5.1 They recorded gas phase spectra by thermal expansion of a chiral α-pinene with each chiral terpene.

For the complex of 4 with (+)-1 or (-)-1 and 5 with (+)-1 or (-)-1, the OH vibrational frequency is identical for the two different stereoisomers. However, the OH vibrational frequencies differ by 2 cm-1 with 3, and the complex of 3/(+)-1 displays two different OH stretches that differ by 11 cm-1. And in the case of the complex of α-pinene with 2, the OH vibrational frequencies of the two different stereoisomers differ by 11 cm-1!

The B3LYP-D3(BJ)/def2-TZVP optimized geometry of the 2/(+)-1 and 2/(-)-1 complexes are shown in Figure 2, and some subtle differences in sterics and dispersion give rise to the different vibrational frequencies.



Figure 2. B3LYP-D3(BJ)/def2-TZVP optimized geometry of the 2/(+)-1 and 2/(-)-1

Of interest to readers of this blog will be the DFT study of these complexes. The authors used three different well-known methods – B3LYP-D3(BJ)/def2-TZVP, M06-2x/def2-TZVP, and ωB97X-D/def2-TZVP – to compute structures and (most importantly) predict the vibrational frequencies. Interestingly, M06-2x/def2-TZVP and ωB97X-D/ def2-TZVP both failed to predict the vibrational frequency difference between the complexes with the two stereoisomers of α-pinene. However, B3LYP-D3(BJ)/def2-TZVP performed extremely well, with a mean average error (MAE) of only 1.9 cm-1 for the four different terpenes. Using this functional and the larger may-cc-pvtz basis set reduced the MAE to 1.5 cm-1 with the largest error of only 2.5 cm-1.

As the authors note, these complexes provide some fertile ground for further experimental and computational study and benchmarking.


1. Medel, R.; Stelbrink, C.; Suhm, M. A., “Vibrational Signatures of Chirality Recognition Between α-Pinene and Alcohols for Theory Benchmarking.” Angew. Chem. Int. Ed. 2019, 58, 8177-8181, DOI: 10.1002/anie.201901687.


(-)-1, (-)-α-pinene: InChI=1S/C10H16/c1-7-4-5-8-6-9(7)10(8,2)3/h4,8-9H,5-6H2,1-3H3/t8-,9-/m0/s1

(+)-1, (-)-α-pinene: InChI=1S/C10H16/c1-7-4-5-8-6-9(7)10(8,2)3/h4,8-9H,5-6H2,1-3H3/t8-,9-/m1/s1

2, (-)borneol: InChI=1S/C10H18O/c1-9(2)7-4-5-10(9,3)8(11)6-7/h7-8,11H,4-6H2,1-3H3/t7-,8+,10+/m0/s1

3, (+)-fenchol: InChI=1S/C10H18O/c1-9(2)7-4-5-10(3,6-7)8(9)11/h7-8,11H,4-6H2,1-3H3/t7-,8-,10+/m0/s1

4, (-1)-isopinocampheol: InChI=1S/C10H18O/c1-6-8-4-7(5-9(6)11)10(8,2)3/h6-9,11H,4-5H2,1-3H3/t6-,7+,8-,9-/m1/s1

5, (1S)-1-phenylethanol: InChI=1S/C8H10O/c1-7(9)8-5-3-2-4-6-8/h2-7,9H,1H3/t7-/m0/s1

DFT &vibrational frequencies Steven Bachrach 10 Jun 2019 No Comments

C60 Fullerene isomers

The Grimme group has examined all 1812 C60 isomers, in part to benchmark some computational methods.1 They computed all of these structures at PW6B95-D3/def2-QZVP//PBE-D3/def2-TZVP. The lowest energy structure is the expected fullerene 1 and the highest energy structure is the nanorod 2 (see Figure 1).



Figure 1. Optimized structures of the lowest (1) and highest (2) energy C60 isomers.

About 70% of the isomers like in the range of 150-250 kcal mol-1 above the fullerene 1, and the highest energy isomer 2 lies 549.1 kcal mol-1 above 1. To benchmark some computational methods, they selected the five lowest energy isomers and five other isomers with higher energy to serve as a new database (C60ISO), with energies computed at DLPNO-CCSD(T)/CBS*. The mean absolute deviation of the PBE-D3/def2-TZVP relative energies with the DLPNO-CCSD(T)/CBS* energies is relative large 10.7 kcal mol-1. However, the PW6B95-D3/def2-QZVP//PBE-D3/def2-TZVP method is considerably better, with a MAD of only 1.7 kcal mol-1. This is clearly a reasonable compromise method for fullerene-like systems, balancing accuracy with computational time.

They also compared the relative energies of all 1812 isomers computed at PW6B95-D3/def2-QZVP//PBE-D3/def2-TZVP with a number of semi-empirical methods. The best results are with the DFTB-D3 method, with an MAD of 5.3 kcal mol-1.


1) Sure, R.; Hansen, A.; Schwerdtfeger, P.; Grimme, S., "Comprehensive theoretical study of all 1812 C60 isomers." Phys. Chem. Chem. Phys. 2017, 19, 14296-14305, DOI: 10.1039/C7CP00735C.


1: InChI=1S/C60/c1-2-5-6-3(1)8-12-10-4(1)9-11-7(2)17-21-13(5)23-24-14(6)22-18(8)28-20(12)30-26-16(10)15(9)25-29-19(11)27(17)37-41-31(21)33(23)43-44-34(24)32(22)42-38(28)48-40(30)46-36(26)35(25)45-39(29)47(37)55-49(41)51(43)57-52(44)50(42)56(48)59-54(46)53(45)58(55)60(57)59

2: InChI=1S/C60/c1-11-12-2-21(1)31-41-32-22(1)3-13(11)15-5-24(3)34-43(32)53-55-47-36-26-6-16-17-7(26)28-9-19(17)20-10-29-8(18(16)20)27(6)37-46(36)54(51(41)55)52-42(31)33-23(2)4(14(12)15)25(5)35-44(33)58-56(52)48(37)39(29)50-40(30(9)10)49(38(28)47)57(53)59(45(34)35)60(50)58

fullerene &Grimme Steven Bachrach 05 Mar 2018 No Comments

More applications of computed NMR spectra

In this post I cover two papers discussing application of computed NMR chemical shifts to structure identification and (yet) another review of computational techniques towards NMR structure prediction.

Grimblat, Kaufman, and Sarotti1 take up the structure of rubriflordilactone B 1, which was isolated from Schisandra rubriflora. The compound was then synthesized and its x-ray structure reported, however its NMR did not match with the natural extract. It was suggested that there were actually two compounds in the extract, the minor one was less soluble and is the crystallized 1, and a second compound responsible for the NMR signal.

The authors looked at all stereoisomers of this molecule keeping the three left-most rings intact. The low energy rotamers of these 32 stereoisomers were then optimized at B3LYP/6-31G* and the chemical shifts computed at PCM(pyridine)/mPW1PW91/6-31+G**. To benchmark the method, DP4+ was used to identify which stereoisomer best matches with the observed NMR of authentic 1; the top fit (92.6% probability) was the correct structure.

The 32 stereoisomers were then tested against the experimental NMR of the natural extract. DP4+ with just the proton shifts suggested structure 2 (99.8% probability); however, the 13C chemical shifts predicted a different structure. Re-examination of the reported chemical shifts identifies some mis-assigned signals, which led to a higher C-DP4+ prediction. When all 128 stereoisomers were tested, structure 2 had the highest DP4+ prediction (99.5%), but the C-DP4+ prediction remained problematic (10.8%). Analyzing the geometries of all reasonable alternative for agreement with the NOESY spectrum confirmed 2. These results underscore the importance of using all data sources.

Reddy and Kutateladze point out the importance of using coupling constants along with chemical shifts in structure identification.2 They examined cordycepol A 3, obtained from Cordyceps ophioglossoides. They noted that the computed chemical shifts and coupling constants of originally proposed structure 3a differed dramatically from the experimental values.

They first proposed that the compound has structure 3b. The computed coupling constants using their relativistic force field.3 The experimental coupling constants for the proton H1 are 13.4 and 7.1 Hz. The computed values for 3a are 8.9 and 1.6 Hz, and this structure is clearly incorrect. The coupling constants are improved with 3b, but the 13C chemical shifts are in poor agreement with experiment. So, they proposed structure 3c, the epimer at both C1 and C11 of the original structure.

They optimized four conformations of 3c at B3LYP/6-31G(d) and obtained Boltzmann-weighted chemical shifts at mPW1PW91/6-311+G(d,p). The RMS deviation of the computed 13C chemical shifts relative to the experiment is only 1.54 ppm, and more importantly, the computed coupling constants of 13.54 and 6.90 Hz are in excellent agreement with the experiment values.

Lastly, Grimblat and Sarotti present a review of a number of methods for using computed NMR chemical shifts towards structure prediction.4 These methods include CP3, DP4, DP4+ (all of which I have posted on in the past) and an artificial neural network approach of their own design. They discuss a number of interesting cases where each of these methods has been crucial in identifying the correct chemical structure.


1. Grimblat, N.; Kaufman, T. S.; Sarotti, A. M., "Computational Chemistry Driven Solution to Rubriflordilactone B." Org. Letters 2016, 18, 6420-6423, DOI: 10.1021/acs.orglett.6b03318.

2. Reddy, D. S.; Kutateladze, A. G., "Structure Revision of an Acorane Sesquiterpene Cordycepol A." Org. Letters 2016, 18, 4860-4863, DOI: 10.1021/acs.orglett.6b02341.

3. (a) Kutateladze, A. G.; Mukhina, O. A., "Minimalist Relativistic Force Field: Prediction of Proton–Proton Coupling Constants in 1H NMR Spectra Is Perfected with NBO Hybridization Parameters." J. Org. Chem. 2015, 80, 5218-5225, DOI: 10.1021/acs.joc.5b00619; (b) Kutateladze, A. G.; Mukhina, O. A., "Relativistic Force Field: Parametrization of 13C–1H Nuclear Spin–Spin Coupling Constants." J. Org. Chem. 2015, 80, 10838-10848, DOI: 10.1021/acs.joc.5b02001.

4. Grimblat, N.; Sarotti, A. M., "Computational Chemistry to the Rescue: Modern Toolboxes for the Assignment of Complex Molecules by GIAO NMR Calculations." Chem. Eur. J. 2016, 22, 12246-12261, DOI: h10.1002/chem.201601150.


1: InChI=1S/C28H30O6/c1-13-9-20(32-26(13)30)25-14(2)24-17-6-5-15-12-28-21(8-7-16(15)18(17)10-19(24)31-25)27(3,4)33-22(28)11-23(29)34-28/h5-9,14,19-22,24-25H,10-12H2,1-4H3/t14-,19+,20-,21-,22+,24-,25-,28+/m0/s1

2: InChI=1S/C28H30O6/c1-13-9-20(32-26(13)30)25-14(2)24-17-6-5-15-12-28-21(8-7-16(15)18(17)10-19(24)31-25)27(3,4)33-22(28)11-23(29)34-28/h5-9,14,19-22,24-25H,10-12H2,1-4H3/t14-,19-,20-,21-,22+,24+,25-,28+/m0/s1

3c: InChI=1S/C16H28O2/c1-6-11(2)9-14-16(5)12(3)7-8-13(16)15(4,17)10-18-14/h9,12-14,17H,6-8,10H2,1-5H3/b11-9-/t12-,13-,14-,15-,16+/m0/s1

NMR Steven Bachrach 04 Oct 2017 No Comments

A few review articles

A few nice review/opinion pieces have been piling up in my folder of papers of interest for this Blog. So, this post provides a short summary of a number of review articles that computationally-oriented chemists may find of interest.

Holy Grails in computational chemistry

Houk and Liu present a short list of “Holy Grails” in computationally chemistry.1 They begin by pointing out a few technical innovations that must occur for the Grails to be found: development of a universal density functional; an accurate, generic force field; improved sampling for MD; and dealing with the combinatorial explosion with regards to conformations and configurations. Their list of Grails includes predicting crystal structures and structure of amorphous materials, catalyst design, reaction design, and device design. These Grails overlap with the challenges I laid out in my similarly-themed article in 2014.2

Post-transition state bifurcations and dynamics

Hare and Tantillo review the current understanding of post-transition state bifurcations (PTSB).3 This type of potential energy surface has been the subject of much of Chapter 8 of my book and many of my blog posts. What is becoming clear is the possibility of a transition state followed by a valley-ridge inflection leads to reaction dynamics where trajectories cross a single transition state but lead to two different products. This new review updates the state-of-the-art from Houk’s review4 of 2008 (see this post). Mentioned are a number of studies that I have included in this Blog, along with reactions involving metals, and biochemical systems (many of these examples come from the Tantillo lab). They close with the hope that their review might “inspire future studies aimed at controlling selectivity for reactions with PTSBs” (italics theirs). I might offer that controlling selectivity in these types of dynamical systems is another chemical Grail!

The Hase group has a long review of direct dynamics simulations.5 They describe a number of important dynamics studies that provide important new insight to reaction mechanism, such as bimolecular SN2 reactions (including the roundabout mechanism) and unimolecular dissociation. They write a long section on post-transition state bifurcations, and other dynamic effects that cannot be interpreted using transition state theory or RRKM. This section is a nice complement to the Tantillo review.

Benchmarking quantum chemical methods

Mata and Suhm look at our process of benchmarking computational methods.6 They point out the growing use of high-level quantum computations as the reference for benchmarking new methods, often with no mention of any comparison to experiment. In defense of theoreticians, they do note the paucity of useful experimental data that may exist for making suitable comparisons. They detail a long list of better practices that both experimentalists and theoreticians can take to bolster both efforts, leading to stronger computational tools that are more robust at helping to understand and discriminate difficult experimental findings.


1) Houk, K. N.; Liu, F., "Holy Grails for Computational Organic Chemistry and Biochemistry." Acc. Chem. Res. 2017, 50 (3), 539-543, DOI: 10.1021/acs.accounts.6b00532.

2) Bachrach, S. M., "Challenges in computational organic chemistry." WIRES: Comput. Mol. Sci. 2014, 4, 482-487, DOI: 10.1002/wcms.1185.

3) Hare, S. R.; Tantillo, D. J., "Post-transition state bifurcations gain momentum – current state of the field." Pure Appl. Chem. 2017, 89, 679-698, DOI: 0.1515/pac-2017-0104.

4) Ess, D. H.; Wheeler, S. E.; Iafe, R. G.; Xu, L.; Çelebi-Ölçüm, N.; Houk, K. N., "Bifurcations on Potential Energy Surfaces of Organic Reactions." Angew. Chem. Int. Ed. 2008, 47, 7592-7601, DOI: 10.1002/anie.200800918

5) Pratihar, S.; Ma, X.; Homayoon, Z.; Barnes, G. L.; Hase, W. L., "Direct Chemical Dynamics Simulations." J. Am. Chem. Soc. 2017, 139, 3570-3590, DOI: 10.1021/jacs.6b12017.

6) Mata, R. A.; Suhm, M. A., "Benchmarking Quantum Chemical Methods: Are We Heading in the Right Direction?" Angew. Chem. Int. Ed. 2017, ASAP, DOI: 10.1002/anie.201611308.

Dynamics &Houk Steven Bachrach 25 Jul 2017 No Comments

Another procedure for computing NMR chemical shifts

Here’s another take on automating a procedure for using computer 13C chemical shifts to assess chemical structure.1 (Have a look at these previous posts for some alternative methods and applications.) The approach here is to benchmark a few computational methods against a conformationally flexible drug-like molecule, in this case 1. A variety of conformations were optimized using the different computational methods, and 13C chemical shifts evaluated from a Boltzmann-weighted distribution. While the best agreement with the experimental chemical shifts (based on the root-mean-squared deviation) is with ωB97XD/cc-pVDZ, the authors opt for B3LYP/cc-pVDZ for its computational efficiency with only slightly poorer performance. (It should be note that WC04/cc-pVDZ, a functional designed for computing 13 chemical shifts,2 is almost as good as ωB97XD/cc-pVDZ. Also, not mentioned in the article is the dramatically poorer performance of the pcS-2 basis set, despite the fact that it was parametrized3 for NMR computation!)

They apply the procedure to a number of test cases. For example, the HIV-1 reverse transcriptase inhibitor nevirapine hydrolyzes to a compound whose structure has been difficult to identify. The four proposed structures 2a-d were subjected to the computational method, and the 13C chemical shift RMSD for 2d is only 2.3ppm, significantly smaller than for the other 3 structures. Compound 2d was then synthesized and its NMR matches that of the nevirapine hydrolysis product.


1) Xin, D.; Sader, C. A.; Chaudhary, O.; Jones, P.-J.; Wagner, K.; Tautermann, C. S.; Yang, Z.; Busacca, C. A.; Saraceno, R. A.; Fandrick, K. R.; Gonnella, N. C.; Horspool, K.; Hansen, G.; Senanayake, C. H., "Development of a 13C NMR Chemical Shift Prediction Procedure Using B3LYP/cc-pVDZ and Empirically Derived Systematic Error Correction Terms: A Computational Small Molecule Structure Elucidation Method." J. Org. Chem. 2017, ASAP, DOI: 10.1021/acs.joc.7b00321.

2) Wiitala, K. W.; Hoye, T. R.; Cramer, C. J., “Hybrid Density Functional Methods Empirically Optimized for the Computation of 13C and 1H Chemical Shifts in Chloroform Solution,” J. Chem. Theory Comput. 2006, 2, 1085-1092, DOI: 10.1021/ct6001016.

3) Jensen, F., “Basis Set Convergence of Nuclear Magnetic Shielding Constants Calculated by Density Functional Methods,” J. Chem. Theory Comput., 2008, 4, 719-727, DOI: 10.1021/ct800013z.


1: InChI=1S/C24H26F4N2O4S/c1-4-35(33,34)18-7-8-20-15(10-18)9-17(30-20)13-23(32,24(26,27)28)22(2,3)12-14-5-6-16(25)11-19(14)21(29)31/h5-11,30,32H,4,12-13H2,1-3H3,(H2,29,31)/t23-/m0/s1

2d: InChI=1S/C15H16N4O2/c1-9-6-8-17-14(18-10-4-5-10)12(9)19-13-11(15(20)21)3-2-7-16-13/h2-3,6-8,10H,4-5H2,1H3,(H,16,19)(H,17,18)(H,20,21)

NMR Steven Bachrach 10 Jul 2017 No Comments

Progress in DFT development and the density they predict

“Getting the right answer for the right reason” – how important is this principle when it comes to computational chemistry? Medvedev and co-workers argue that when it comes to DFT, trends in functional development have overlooked this maxim in favor of utility.1 Specifically, they note that

There exists an exact functional that yields the exact energy of a system from its exact density.

Over the past two decades a great deal of effort has gone into functional development, mostly in an empirical way done usually to improve energy prediction. This approach has a problem:

[It], however, overlooks the fact that the reproduction of exact energy is not a feature of the exact functional, unless the input electron density is exact as well.

So, these authors have studied functional performance with regards to obtaining proper electron densities. Using CCSD/aug-cc-pwCV5Z as the benchmark, they computed the electron density for a number of neutral and cationic atoms having 2, 4, or 10 electrons. Then, they computed the densities with 128 different functionals of all of the rungs of Jacob’s ladder. They find that accuracy was increasing as new functionals were developed from the 1970s to the early 2000s. Since then, however, newer functionals have tended towards poorer electron densities, even though energy prediction has continued to improve. Medvedev et al argue that the recent trend in DFT development has been towards functionals that are highly parameterized to fit energies with no consideration given to other aspects including the density or constraints of the exact functional.

In the same issue of Science, Hammes-Schiffer comments about this paper.2 She notes some technical issues, most importantly that the benchmark study is for atoms and that molecular densities might be a different issue. But more philosophically (and practically), she points out that for many chemical and biological systems, the energy and structure are of more interest than the density. Depending on where the errors in density occur, these errors may not be of particular relevance in understanding reactivity; i.e., if the errors are largely near the nuclei but the valence region is well described then reactions (transition states) might be treated reasonably well. She proposes that future development of functionals, likely still to be driven by empirical fitting, might include other data to fit to that may better reflect the density, such as dipole moments. This seems like a quite logical and rational step to take next.

A commentary by Korth3 summarizes a number of additional concerns regarding the Medvedev paper. The last concern is the one I find most striking:

Even if there really are (new) problems, it is as unclear as before how they can be overcome…With this in mind, it does not seem unreasonable to compromise on the quality of the atomic densities to improve the description of more relevant properties, such as the energetics of molecules.

Korth concludes with

In the meantime, while theoreticians should not rest until they have the right answer for the right reason, computational chemists and experimentalists will most likely continue to be happy with helpful answers for good reasons.

I do really think this is the correct take-away message: DFT does appear to provide good predictions of a variety of chemical and physical properties, and it will remain a widely utilized tool even if the density that underpins the theory is incorrect. Functional development must continue, and Medvedev et al. remind us of this need.


1) Medvedev, M. G.; Bushmarinov, I. S.; Sun, J.; Perdew, J. P.; Lyssenko, K. A., "Density functional theory is straying from the path toward the exact functional." Science 2017, 355, 49-52, DOI: 10.1126/science.aah5975.

2) Hammes-Schiffer, S., "A conundrum for density functional theory." Science 2017, 355, 28-29, DOI: 10.1126/science.aal3442.

3) Korth, M., "Density Functional Theory: Not Quite the Right Answer for the Right Reason Yet." Angew. Chem. Int. Ed. 2017, 56, 5396-5398, DOI: 10.1002/anie.201701894.

DFT Steven Bachrach 08 May 2017 1 Comment

The strongest base?

The new benchmark has been set for superbases. The previous record holder was LiO, with a computed proton affinity of 424.9 kcal mol-1. A new study by Poad, et al., examines the dianions of the three isomeric phenyldiacetylides: 1o, 1m, and 1p.1 Their computed proton affinities (G4(MP2)-6X) are 440.6, 427.0, and 425.6 kcal mol-1, respectively. The optimized geometries of these dianions are shown in Figure 1.




Figure 1. Optimized geometries of 1o, 1m, and 1p.

The authors also prepared these bases inside a mass spectrometer. All three deprotonate water, but do not deprotonate methane, though that might be a kinetic issue.

The authors speculate that 1o will be hard to beat as a base since loss of an electron is always a concern with small dianions.


1) Poad, B. L. J.; Reed, N. D.; Hansen, C. S.; Trevitt, A. J.; Blanksby, S. J.; Mackay, E. G.; Sherburn, M. S.; Chan, B.; Radom, L., "Preparation of an ion with the highest calculated proton affinity: ortho-diethynylbenzene dianion." Chem. Sci. 2016, 7, 6245-6250, DOI: 10.1039/C6SC01726F.


1o: InChI=1S/C10H4/c1-3-9-7-5-6-8-10(9)4-2/h5-8H/q-2

1m: InChI=1S/C10H4/c1-3-9-6-5-7-10(4-2)8-9/h5-8H/q-2

1p: InChI=1S/C10H4/c1-3-9-5-7-10(4-2)8-6-9/h5-8H/q-2

Acidity Steven Bachrach 21 Feb 2017 1 Comment

Calculating large fullerenes

What is the size of a molecule that will stretch computational resources today? Chan and co-workers have examined some very large fullerenes1 to both answer that question, and also to explore how large a fullerene must be to approach graphene-like properties.

They are interested in predicting the heat of formation of large fullerenes. So, they benchmark the heats of formation of C60 using four different isodesmic reactions (Reaction 1-4), comparing the energies obtained using a variety of different methods and basis sets to those obtained at W1h. The methods include traditional functionals like B3LYP, B3PW91, CAM-B3LYP, PBE1PBE, TPSSh, B98, ωB97X, M06-2X3, and MN12-SX, and supplement them with the D3 dispersion correction. Additionally a number of doubly hybrid methods are tested (again with and without dispersion corrections), such as B2-PLYP, B2GPPLYP, B2K-PLYP, PWP-B95, DSD-PBEPBE, and DSD-B-P86. The cc-pVTZ and cc-pVQZ basis sets were used. Geometries were optimized at B3LYP/6-31G(2df,p).

C60 + 10 benzene → 6 corannulene

Reaction 1

C60 + 10 naphthalene → 8 corannulene

Reaction 2

C60 + 10 phenanthrene → 10 corannulene

Reaction 3

C60 + 10 triphenylene → 12 corannulene

Reaction 4

Excellent results were obtained with DSD-PBEPBE-D3/cc-pVQZ (an error of only 1.8 kJ/mol), though even a method like BMK-D3/cc-pVTZ had an error of only 9.2 kJ/mol. They next set out to examine large fullerenes, including such behemoths as C180, C240, and C320, whose geometries are shown in Figure 1. Heats of formation were obtained using isodesmic reactions that compare back to smaller fullerenes, such as in Reaction 5-8.

C70 + 5 styrene → C60 + 5 naphthalene

Reaction 5

C180 → 3 C60

Reaction 6

C320 + 2/3 C60 → 2 C180

Reaction 7




Figure 1. B3LYP/6-31G(2df,p) optimized geometries of C180, C240, and C320. (Don’t forget that clicking on these images will launch Jmol and allow you to manipulate the molecules in real-time.)

Next, taking the heat of formation per C for these fullerenes, using a power law relationship, they were able to extrapolate out the heat of formation per C for truly huge fullerenes, and find the truly massive fullerenes, like C9680, still have heats of formation per carbon 1 kJ/mol greater than for graphene itself.


(1) Chan, B.; Kawashima, Y.; Katouda, M.; Nakajima, T.; Hirao, K. "From C60 to Infinity: Large-Scale Quantum Chemistry Calculations of the Heats of Formation of Higher Fullerenes," J. Am. Chem. Soc. 2016, 138, 1420-1429, DOI: 10.1021/jacs.5b12518.

fullerene Steven Bachrach 22 Feb 2016 4 Comments

« Previous PageNext Page »