Archive for the 'QM Method' Category

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 No Comments

Benchmarking Platonic solids and related hydrocarbons

Karton, Schreiner, and Martin have benchmarked the heats of formation of some Platonic Solids and related hydrocarbons.1 The molecules examined are tetrahedrane 1, cubane 2, dodecahedrane 3, trisprismane 4, pentaprismane 5, and octahedrane 6.

The optimized structures (B3LYP-D3BJ/def2-TZVPP) of these compounds are shown in Figure 1.







Figure 1. B3LYP-D3BJ/def2-TZVPP optimized geometries of 1-6.

Using the W1-F12 and W2-F12 composite methods, the estimated the heats of formation of these hydrocarbons are listed in Table 1. Experimental values are available only for 2 and 3; the computed values are off by about 2 kcal mol-1, which the authors argue is just outside the error bars of the computations. They suggest that the experiments might need to be revisited.

Table 1. Heats of formation (kcal mol-1) of 1-6.


ΔHf (comp)

ΔHf (expt)






142.7 ± 1.2



22.4 ± 1










They conclude with a comparison of strain energies computed using isogyric, isodesmic, and homodesmotic reactions with a variety of computational methods. Somewhat disappointingly, most DFT methods have appreciable errors compared with the W1-F12 results, and the errors vary depend on the chemical reaction employed. However, the double hybrid method DSD-PBEP86-D3BJ consistently reproduces the W1-F12 results.


(1)  Karton, A.; Schreiner, P. R.; Martin, J. M. L. "Heats of formation of platonic hydrocarbon cages by means of high-level thermochemical procedures," J. Comput. Chem. 2016, 37, 49-58, DOI: 10.1002/jcc.23963.


1: InChI=1S/C4H4/c1-2-3(1)4(1)2/h1-4H

2: InChI=1S/C8H8/c1-2-5-3(1)7-4(1)6(2)8(5)7/h1-8H

3: InChI=1S/C20H20/c1-2-5-7-3(1)9-10-4(1)8-6(2)12-11(5)17-13(7)15(9)19-16(10)14(8)18(12)20(17)19/h1-20H

4: InChI=1S/C6H6/c1-2-3(1)6-4(1)5(2)6/h1-6H

5: InChI=1S/C10H10/c1-2-5-3(1)7-8-4(1)6(2)10(8)9(5)7/h1-10H

6: InChI=1S/C12H12/c1-2-4-6-5(11-7(1)10(4)11)3(1)9-8(2)12(6)9/h1-12H

QM Method &Schreiner Steven Bachrach 10 May 2016 No Comments

Dispersion in organic chemistry – a review and another example

The role of dispersion in organic chemistry has been slowly recognized as being quite critical in a variety of systems. I have blogged on this subject many times, discussing new methods for properly treating dispersion within quantum computations along with a variety of molecular systems where dispersion plays a critical role. Schreiner1 has recently published a very nice review of molecular systems where dispersion is a key component towards understanding structure and/or properties.

In a similar vein, Wegner and coworkers have examined the Z to E transition of azobenzene systems (1a-g2a-g) using both experiment and computation.2 They excited the azobenzenes to the Z conformation and then monitored the rate for conversion to the E conformation. In addition they optimized the geometries of the two conformers and the transition state for their interconversion at both B3LYP/6-311G(d,p) and B3LYP-D3/6-311G(d,p). The optimized structure of the t-butyl-substituted system is shown in Figure 1.

a: R=H; b: R=tBu; c: R=Me; d: R=iPr; e: R=Cyclohexyl; f: R=Adamantyl; g: R=Ph




Figure 1. B3LYP-D3/6-311G(d,p) optimized geometries of 1a, 2a, and the TS connecting them.

The experiment finds that the largest activation barriers are for the adamantly 1f and t-butyl 1b azobenzenes, while the lowest barriers are for the parent 1a and methylated 1c azobenzenes.

The trends in these barriers are not reproduced at B3LYP but are reproduced at B3LYP-D3. This suggests that dispersion is playing a role. In the Z conformations, the two phenyl groups are close together, and if appropriately substituted with bulky substituents, contrary to what might be traditionally thought, the steric bulk does not destabilize the Z form but actually serves to increase the dispersion stabilization between these groups. This leads to a higher barrier for conversion from the Z conformer to the E conformer with increasing steric bulk.


(1) Wagner, J. P.; Schreiner, P. R. "London Dispersion in Molecular Chemistry—Reconsidering Steric Effects," Angew. Chem. Int. Ed. 2015, 54, 12274-12296, DOI: 10.1002/anie.201503476.

(2) Schweighauser, L.; Strauss, M. A.; Bellotto, S.; Wegner, H. A. "Attraction or Repulsion? London Dispersion Forces Control Azobenzene Switches," Angew. Chem. Int. Ed. 2015, 54, 13436-13439, DOI: 10.1002/anie.201506126.


1b: InChI=1S/C28H42N2/c1-25(2,3)19-13-20(26(4,5)6)16-23(15-19)29-30-24-17-21(27(7,8)9)14-22(18-24)28(10,11)12/h13-18H,1-12H3/b30-29-

2b: InChI=1S/C28H42N2/c1-25(2,3)19-13-20(26(4,5)6)16-23(15-19)29-30-24-17-21(27(7,8)9)14-22(18-24)28(10,11)12/h13-18H,1-12H3/b30-29+

DFT &Schreiner Steven Bachrach 04 Jan 2016 No Comments

Structure of the 2-fluoroethanol trimer

Here is another fine example of the power of combining experiment and computation. Xu and co-worker has applied the FT microwave technique, which has been used in conjunction with computation by the Alonso group (especially) as described in these posts, to the trimer of 2-fluoroethanol.1 They computed a number of trimer structures at MP2/6-311++G(2d,p) in an attempt to match up the computed spectroscopic constants with the experimental constants. The two lowest energy structures are shown in Figure 1. The second lowest energy structure has nice symmetry, but it does not match up well with the experimental spectra. However, the lowest energy structure is in very good agreement with the experiments.



Table 1. MP2/6-311++G(2d,p) optimized structures and relative energies (kJ mol-1) of the two lowest energy structures of the trimer of 2-fluoroethanol. The added orange lines in the lowest energy structure denote the bifurcated hydrogen bonds identified by QTAIM.

Of particular note is that topological electron density analysis (also known as quantum theoretical atoms in a molecule, QTAIM) of the wavefunction of the lowest energy structure of the trimer identifies two hydrogen bond bifurcations. The authors suggest that these additional interactions are responsible, in part, for the stability of this lowest energy structure.


(1) Thomas, J.; Liu, X.; Jäger, W.; Xu, Y. "Unusual H-Bond Topology and Bifurcated H-bonds in the 2-Fluoroethanol Trimer," Angew. Chem. Int. Ed. 2015, 54, 11711-11715, DOI: 10.1002/anie.201505934.


2-fluoroethanol: InChI=1S/C2H5FO/c3-1-2-4/h4H,1-2H2, InChIKey=GGDYAKVUZMZKRV-UHFFFAOYSA-N

Hydrogen bond &MP Steven Bachrach 20 Oct 2015 1 Comment

Keto-enol Benchmark Study

The keto-enol tautomerization is a fundamental concept in organic chemistry, taught in the introductory college course. As such, it provides an excellent test reaction to benchmark the performance computational methods. Acevedo and colleagues have reported just such a benchmark study.1

First, the compare a wide variety of methods, ranging from semi-empirical, to DFT, and to composite procedures, with experimental gas-phase free energy of tautomerization. They use seven such examples, two of which are shown in Scheme 1. The best results from each computation category are AM1, with a mean absolute error (MAE) of 1.73 kcal mol-1, M06/6-31+G(d,p), with a MAE of 0.71 kcal mol-1, and G4, with a MAE of 0.95 kcal mol-1. All of the modern functionals do a fairly good job, with MAEs less than 1.3 kcal mol-1.

Scheme 1

As might be expected, the errors were appreciably larger for predicting the free energy of tautomerization, with a good spread of errors depending on the method for handling solvent (PCM, CPCM, SMD) and the choice of cavity radius. The best results were with the G4/PCM/UA0 procedure, though M06/6-31+G(d,p)/PCM and either UA0 or UFF performed quite well, at considerably less computational expense.


(1) McCann, B. W.; McFarland, S.; Acevedo, O. "Benchmarking Continuum Solvent Models for Keto–Enol Tautomerizations," J. Phys. Chem. A 2015, 119, 8724-8733, DOI: 10.1021/acs.jpca.5b04116.

Keto-enol tautomerization &QM Method Steven Bachrach 12 Oct 2015 No Comments

Domino Tunneling

A 2013 study of oxalic acid 1 failed to uncover any tunneling between its conformations,1 despite observation of tunneling in other carboxylic acids (see this post). This was rationalized by computations which suggested rather high rearrangement barriers. Schreiner, Csaszar, and Allen have now re-examined oxalic acid using both experiments and computations and find what they call domino tunneling.2

First, they determined the structures of the three conformations of 1 along with the two transition states interconnecting them using the focal point method. These geometries and relative energies are shown in Figure 1. The barrier for the two rearrangement steps are smaller than previous computations suggest, which suggests that tunneling may be possible.






Figure 1. Geometries of the conformers of 1 and the TS for rearrangement and relative energies (kcal mol-1)

Placing oxalic acid in a neon matrix at 3 K and then exposing it to IR radiation populates the excited 1tTt conformation. This state then decays to both 1cTt and 1cTc, which can only happen through a tunneling process at this very cold temperature. Kinetic analysis indicates that there are two different rates for decay from both 1tTt and 1cTc, with the two rates associated with different types of sites within the matrix.

The intrinsic reaction paths for the two rearrangements: 1tTt1cTt and → 1cTc were obtained at MP2/aug-cc-pVTZ. Numerical integration and the WKB method yield similar half-lives: about 42 h for the first reaction and 23 h for the second reaction. These match up very well with the experimental half-lives from the fast matrix sites of 43 ± 4 h and 30 ± 20 h, respectively. Thus, the two steps take place sequentially via tunneling, like dominos falling over.


(1) Olbert-Majkut, A.; Ahokas, J.; Pettersson, M.; Lundell, J. "Visible Light-Driven Chemistry of Oxalic Acid in Solid Argon, Probed by Raman Spectroscopy," J. Phys. Chem. A 2013, 117, 1492-1502, DOI: 10.1021/jp311749z.

(2) Schreiner, P. R.; Wagner, J. P.; Reisenauer, H. P.; Gerbig, D.; Ley, D.; Sarka, J.; Császár, A. G.; Vaughn, A.; Allen, W. D. "Domino Tunneling," J. Am. Chem. Soc. 2015, 137, 7828-7834, DOI: 10.1021/jacs.5b03322.


1: InChI=1S/C2H2O4/c3-1(4)2(5)6/h(H,3,4)(H,5,6)

focal point &Schreiner &Tunneling Steven Bachrach 11 Aug 2015 1 Comment

Structure of 2-oxazoline

A recent reinvestigation of the structure of 2-oxazoline demonstrates the difficulties that many computational methods can still have in predicting structure.

Samdal, et al. report the careful examination of the microwave spectrum of 2-oxzoline and find that the molecule is puckered in the ground state.1 It’s not puckered by much, and the barrier for inversion of the pucker, through a planar transition state is only 49 ± 8 J mol-1. The lowest vibrational frequency in the non-planar ground state, which corresponds to the puckering vibration, has a frequency of 92 ± 15 cm-1. This low barrier is a great test case for quantum mechanical methodologies.

And the outcome here is not particularly good. HF/cc-pVQZ, M06-2X/cc-pVQZ, and B3LYP/cc-pVQZ all predict that 2-oxazoline is planar. More concerning is that CCSD and CCSD(T) with either the cc-pVTZ or cc-pVQZ basis sets also predict a planar structure. CCSD(T)-F12 with the cc-pVDZ predicts a non-planar ground state with a barrier of only 8.5 J mol-1, but this barrier shrinks to 5.5 J mol-1 with the larger cc-pVTZ basis set.

The only method that has good agreement with experiment is MP2. This method predicts a non-planar ground state with a pucker barrier of 11 J mol-1 with cc-pVTZ, 39.6 J mol-1 with cc-pVQZ, and 61 J mol-1 with the cc-pV5Z basis set. The non-planar ground state and the planar transition state of 2-oxazoline are shown in Figure 1. The computed puckering vibrational frequency does not reproduce the experiment as well; at MP2/cc-pV5Z the predicted frequency is 61 cm-1 which lies outside of the error range of the experimental value.


Planar TS

Figure 1. MP2/cc-pV5Z optimized geometry of the non-planar ground state and the planar transition
state of 2-oxazoline.


(1) Samdal, S.; Møllendal, H.; Reine, S.; Guillemin, J.-C. "Ring Planarity Problem of 2-Oxazoline Revisited Using Microwave Spectroscopy and Quantum Chemical Calculations," J. Phys. Chem. A 2015, 119, 4875–4884, DOI: 10.1021/acs.jpca.5b02528.


2-oxazoline: InChI=1S/C3H5NO/c1-2-5-3-4-1/h3H,1-2H2

MP &vibrational frequencies Steven Bachrach 15 Jun 2015 1 Comment

Gas phase structure of uridine

To advance our understanding of why ribose takes on the furanose form, rather than the pyranose form, in RNA, Alonso and co-workers have examined the structure of uridine 1 in the gas phase.1


Uridine is sensitive to temperature, and so the laser-ablation method long used by the Alonso group is ideal for examining uridine. The microwave spectrum is quite complicated due to the presence of many photofragments. Careful analysis lead to the identification of a number of lines and hyperfine structure that could be definitively assigned to uridine, leading to experimental values of the rotational constants and the diagonal elements of the 14N nuclear quadrupole coupling tensor for each nitrogen. These values are listed in Table 1.

Table 1. Experimental and calculated rotational constants (MHz), quadrupole coupling constants (MHz) and relative energy (kcal mol-1).
































14N1 χxx







14N1 χyy







14N1 χzz







14N3 χxx







14N3 χyy







14N3 χzz







Rel E







In order to assign a 3-D structure to these experimental values, they examined the PES of uridine with molecular mechanics and semi-empirical methods, before reoptimizing the structure of the lowest 5 energy structures at MP2/6-311++G(d,p). Then, comparison of the resulting rotational constants and 14N nuclear quadrupole coupling constants of these computed structures (see Table 1) led to identification of the lowest energy structure (anti/C2’-endo-g+, see Figure 1) in best agreement with the experiment. Once again, the Alonso group has demonstrated the value of the synergy between experiment and computation in structure identification.

Figure 1. MP2/6-311++G(d,p) optimized structure of 1 (anti/C2’-endo-g+).


(1) Peña, I.; Cabezas, C.; Alonso, J. L. "The Nucleoside Uridine Isolated in the Gas Phase," Angew. Chem. Int. Ed. 2015, 54, 2991-2994, DOI: 10.1002/anie.201412460.


1: Inchi=1S/C9H12N2O6/c12-3-4-6(14)7(15)8(17-4)11-2-1-5(13)10-9(11)16/h1-2,4,6-8,12,14-15H,3H2,(H,10,13,16)/t4-,6-,7-,8-/m1/s1

MP &nucleic acids Steven Bachrach 06 Apr 2015 No Comments

Microsolvated structure of β-propiolactone

The structure of water about a solute remains of critical importance towards understanding aqueous solvation. Microwave spectroscopy and computations are the best tools we have today to gain insight on this problem. This is nicely demonstrated in the Alonso study of the microsolvated structures of β-propiolactone 1.1 They employed chirped-pulse Fourier transform microwave (CP-FTMW) spectroscopy and MP2(fc)/6-311++G(d,p) computations to examine the structure involving 1-5 water molecules.

The computed structures of these microsolvated species are shown in Figure 1. The deviation of the computed and experimental structures (RMS in the atomic positions) is small, though increasing as the size of the cluster increases. The deviation is 0.014 Å for the 1. H2O cluster and 0.244 Å for the 1.(H2O)5 cluster. They identified two clusters with four water molecules; the lower energy structure, labeled as a, is only 0.2 kJ mol-1 more stable than structure b.




1.(H2O)4 a

1.(H2O)4 b


Figure 1. MP2(fc)/6-311++G(d,p) optimized geometries of the hydrates of 1.

Water rings are found in the clusters having four or five water molecules, while chains are identified in the smaller clusters. One might imagine water cages appearing with even more water molecules in the microsolvated structures.


(1) Pérez, C.; Neill, J. L.; Muckle, M. T.; Zaleski, D. P.; Peña, I.; Lopez, J. C.; Alonso, J. L.; Pate, B. H. Angew. Chem. Int. Ed. 2015, 54, 979-982, DOI: 10.1002/anie.201409057.


1: InChI=1S/C3H4O2/c4-3-1-2-5-3/h1-2H2

MP &Solvation Steven Bachrach 24 Feb 2015 1 Comment

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

Next Page »