How Large is the Elephant in the Density Functional Theory Room?

Frank Jensen (2017)

Highlighted by Jan Jensen

Having highlighted this paper it is only right that I highlight Frank Jensen's response. To recap, the previous study used wavelets to compute benchmark energies for PBE and PBE0 functionals and showed that even aug-cc-pV5Z failed to reach chemical accuracy for some atomization energies.

In this paper Jensen shows that this problem goes away when one uses basis sets specifically designed for DFT calculations. At the pentuple-zeta level the maximum errors are reduced by factors of 5 and 10 for segmented contracted and uncontracted basis sets, respectively and the MAE for atomization energies are well below 1 kcal/mol.


This work is licensed under a Creative Commons Attribution 4.0 International License.

Is the Accuracy of Density Functional Theory for Atomization Energies and Densities in Bonding Regions Correlated?

Brorsen, K. R.; Yang, Y.; Pak, M. V.; Hammes-Schiffer, S., J. Phys. Chem. Lett. 2017, 8, 2076-2081
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

I recently blogged about a paper arguing that modern density functional development has strayed from the path of improving density description, in favor of improved energetics. The Medvedev paper¹ was met with a number of criticisms. A potential “out” from the conclusions of the work was that perhaps molecular densities do not fare so poorly with more modern functionals, following the argument that better energies might reflect better densities in bonding regions.

The Hammes-Schiffer group have now examined 14 diatomic molecules with the goal of testing just this hypothesis.² They subjected both homonuclear diatomics, like N₂, Cl₂, and Li₂, and heteronuclear diatomics, like HF, LiF, and SC, to 90 different density functionals using the very large aug-cc-pCVQZ basis set. Using the CCSD density as a reference, they examined the differences in the densities predicted by the various functional both along the internuclear axis and perpendicular to it.

The 20 functionals that do the best job in mimicking the CCSD density are all hybrid GGA functionals, along with the sole double hybrid functional included in the study (B2PLYP). These functionals date from 1993 to 2012. The 20 functionals that do the poorest job include functionals from all rung-types, and date from 1980-2012. A very slight upward trend can be observed in the density error increasing with development year, while the error in the dissociation energy clearly is decreasing over time.

They note that six functionals of the Minnesota-type, those that are highly parameterized and of recent vintage, perform very poorly at predicting atomic densities, but do well with the diatomic densities.

Hammes-Schiffer concludes that their diatomic results support the general trend noted by Medvedev’s atomic results, that density description is lagging in more recently developed functionals. I’d add that this trend is not as dramatic for the diatomics as for atoms.

They pose what is really the key question: “Is the purpose to approximate the exact functional or simply to provide chemists with a useful tool for exploring chemical systems?” Since, as they note, the modern highly parameterized functionals have worked so well for predicting energies and geometries, “the observation that many modern functionals produce incorrect densities could be of no great consequence for many studies”. Nonetheless, “the ultimate goal is still to obtain both accurate densities and accurate energies”.

References

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) Brorsen, K. R.; Yang, Y.; Pak, M. V.; Hammes-Schiffer, S., "Is the Accuracy of Density Functional Theory for Atomization Energies and Densities in Bonding Regions Correlated?" J. Phys. Chem. Lett. 2017, 8, 2076-2081, DOI: 10.1021/acs.jpclett.7b00774.

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Mechanisms and Origins of Periselectivity of the Ambimodal [6 + 4] Cycloadditions of Tropone to Dimethylfulvene

Yu, P.; Chen, T. Q.; Yang, Z.; He, C. Q.; Patel, A.; Lam, Y.-h.; Liu, C.-Y.; Houk, K. N., J. Am. Chem. Soc. 2017, 139 (24), 8251-8258
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Bispericyclic transition states arise when two pericyclic reactions merge to a common transition state. This leads to a potential energy surface with a bifurcation such that reactions that traverse this type of transition state will head towards two different products. The classic example is the dimerization of cyclopentadiene, involving two [4+2] Diels-Alder reactions. Unusual PESs are discussed in my book and in past blog posts.

Houk and coworkers have now identified a bispericyclic transition state involving two [6+4] cycloadditions.¹ Reaching back to work Houk pursued as a graduate student with Woodward for inspiration, these authors examined the reaction of tropone 1 with dimethylfulvene 2. Each moiety can act as the diene or triene component of a [6+4] allowed cycloaddition:

The product with fulvene 2 as the 6 π-e component and tropone as the 4 π-e component [6_F+4_T] is 3, while reversing their participation in the [6_T+4_F] cycloaddition leads to 4. A variety of [4+2] reactions are also possible. All of these reactions were investigated at PCM/M06-2X/6-311+G(d,p)//B3LYP-D3/6-31G(d). The reaction leading to 3 is exothermic by 3.0 kcal mol^-1, while the reaction to 4 endothermic by 1.3 kcal mol^-1.

Interestingly, there is only one transition state that leads to both 3 and 4, the first known bispericyclic transition state for two conjoined [6+4] cycloadditions. The barrier is 27.9 kcal mol^-1. The structures of the two products and the transition state leading to them are shown in Figure 1. 3 and 4 can interconvert through a Cope transition state, also shown in Figure 1, with a barrier of 26.3 kcal mol^-1 (for 4 → 3).

3	4
TS [6+4]	TS Cope

Figure 1. B3LYP-D3/6-31G(d) optimized geometries.

Given that a single transition leads to two products, the product distribution is dependent on the molecular dynamics. A molecular dynamics simulation at B3LYP-D3/6-31G(d) with 117 trajectories indicates that 4 is formed 91% while 3 is formed only 9%. Once again, we are faced with the reality of much more complex reaction mechanisms/processes than simple models would suggest.

References

1) Yu, P.; Chen, T. Q.; Yang, Z.; He, C. Q.; Patel, A.; Lam, Y.-h.; Liu, C.-Y.; Houk, K. N., "Mechanisms and Origins of Periselectivity of the Ambimodal [6 + 4] Cycloadditions of Tropone to Dimethylfulvene." J. Am. Chem. Soc. 2017, 139 (24), 8251-8258, DOI: 10.1021/jacs.7b02966.

InChIs

1: InChI=1S/C7H6O/c8-7-5-3-1-2-4-6-7/h1-6H
InChIKey=QVWDCTQRORVHHT-UHFFFAOYSA-N

2: InChI=1S/C8H10/c1-7(2)8-5-3-4-6-8/h3-6H,1-2H3
InChIKey=WXACXMWYHXOSIX-UHFFFAOYSA-N

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

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

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1−86)

Stefan Grimme, Christoph Bannwarth, and Philip Shushkov (2017)
Highlighted by Jan Jensen

Copyright American Chemical Society (2017)

Inspired by the success of the sTDA-xTB method for predicting electronic spectra, Grimme and co-workers present the tight binding DFT method GFN-xTB for predicting Geometries, vibrational Frequencies and Noncovalent interactions ("x" stands for extensions in the AO basis set and the form of the Hamiltonian).  

The method is implemented in a new program called xtb which includes shared memory parallelization, angeometry optimizer, MD, conformational searches, reaction path optimisation modules, and a GB solvation model supplemented with analytical nuclear gradients.  The software is free to academics and can be obtained by writing to the authors.

Importantly, the method is parameterized for all spd elements, so finally there is an alternative to PM6 for heavier elements. A Fermi smearing technique helps convergence for molecules with small HOMO-LUMO gaps and GFN-xTB gives structures organometallic compounds that are in good agreement with DFT.

As implied by the name, GFN-xTB is not parameterized against reaction energies but the methods gives reasonable results for properties like heats of formation.  However, as the authors write

A key premise of the present method is its special purpose character. In our view, low-cost semiempirical QM methods cannot describe simultaneously very different chemical proper- ties, such as structures and chemical reaction energies, and the GFN-xTB method (as the name conveys) focuses on structural properties. The efficiently computed structures and vibrational frequencies or the conformers obtained from global search procedures can (and should) be used subsequently for more accurate DFT or WFT refinements. We hope that the method can serve as a general tool in quantum chemistry and in particular recommend GFN-xTB optimized structures (and thermostatistical corrections) in a multilevel scheme together with PBEh-3c single-point energies. Large-scale molecular dynamics, screening of huge molecular spaces (libraries), parametrization of force-fields, or providing input for novel machine learning techniques are obvious other fields of application.

I thank Anders Christensen for bringing this paper to my attention


This work is licensed under a Creative Commons Attribution 4.0 International License.

A few review articles

Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

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.¹ 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.²

Post-transition state bifurcations and dynamics

Hare and Tantillo review the current understanding of post-transition state bifurcations (PTSB).³ 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 review⁴ 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.⁵ They describe a number of important dynamics studies that provide important new insight to reaction mechanism, such as bimolecular S_N2 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.⁶ 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.

References

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.

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.

Intrinsic map dynamics exploration for uncharted effective free-energy landscapes

Eliodoro Chiavazzo, Roberto Covino, Ronald R. Coifman, C. William Geard, Anastasia S. Georgiou, Gerhard Hummer, and Ioannis G. Kevrekidis. PNAS June 20, 2017
Contributed by Jesper Madsen

The desire to use enhancing sampling to save computational expense in simulations has been there from the beginning. Broadly speaking, there are two main approaches of enhancing molecular dynamics (MD) simulations in order to determine the free-energy surface (FES) of a computationally expensive Hamiltonian: 1) Trajectory-based enhanced sampling (e.g. temperature replica-exchange) and 2) collective variable (CV)-based methods (e.g. umbrella sampling). It is worth noting that the “zoo” of actual techniques to date is rather large. Either type of method, trajectory-based or CV-based, comes with its own set of advantages and disadvantages and mixing-and-matching is popular. 

Here I highlight a newcomer that is spun off from the rapidly growing field of Machine Learning – a field that most of us is keeping a keen eye on these days. 

The algorithm is called intrinsic map dynamics (iMapD) and it is conceptually simple (See Fig. 1 from the paper). 

1. Run a MD trajectory

2. Figure out, in some abstract sense, what region of configuration space you have sampled.

3. Determine the (non-linear) boundary of the sampled region in the abstract space

4. Initialize new MD trajectories in these boundary areas and explore the uncharted territories of the FES.

It is with the concretization of each step above that machine learning has contributed its ideas and tools. Specifically, the map of the configuration space is data mined and a d-dimensional manifold learning technique called diffusion maps (DMAPs) is applied to find the appropriate manifold and its dimensionality. The (d-1)-dimensional boundary manifold of the explored region is determined by a “wrapping” algorithm (here they use alpha shapes). Outward extrapolation is done by performing local principal component (LPC) analysis in ambient space. New simulations are seeded from these extended initial configurations and the algorithm loops back to the beginning. 

Fig.1: Pictorial illustration of the iMapD exploration procedure with (Left) 1D and (Right) 2D effective FESs. In Left Inset, a good collective coordinate is already available—the collective coordinates in Left and Right are not a priori known. "Copyright (2017) National Academy of Sciences

Since iMapD is a trajectory-based approach, absolutely no prior knowledge about the mechanistics of the process we study is required, which is appealing. One can think of the method as a clever hybrid MD/Monte Carlo algorithm and time will show how it stacks up against other alternative approaches in terms of practical usefulness. 

WInote that the algorithm is pedagogically presented in the paper and if you fancy Indiana Jones, well then there’s a little treat for you in the main text. Enjoy.

The Structure of the Elusive Simplest Dipeptide Gly-Gly

Cabezas, C.; Varela, M.; Alonso, J. L., Angew. Chem. Int. Ed. 2017, 56, 6420-6425
Contributed by Steven Bacharach
Reposted from Computational Organic Chemistry with permission

Continuing their application of laser ablation molecular beam Fourier transform microwave (LA-MB-FTMW) spectroscopy and computational chemistry to biochemical molecules (see these previous posts), the Alonso group reports on the structure of the glycine-glycine dipeptide 1.¹ The microwave spectrum shows three different conformers. MP2/6-311++G(d,p) computations, the same method they have previously utilized for predicting geometries, revealed a number of different conformations. By matching the spectroscopic parameters obtained from the spectrum with those of the computed structures, they proposed the three conformations 1a, 1b, and 1c, shown in Figure 1.

1a

1b

1c

Figure 1. ωb97xd/6-31G(d) optimized structures of the three conformers of 1.
Note that the authors did not report their structures in their supporting materials(!) so I have optimized them.

The structures of conformers 1a and 1b are nearly planar. MP2 predicts a non-planar rotomer of 1a, which brings the carboxyl group out of plane, to be the lowest conformation in terms of electronic energy. With the M06-2x functional, this non-planar rotomer is about isoenergetic with 1a. With all computational levels 1a is the lowest in free energy. The barrier for rotation between the non-planar rotomer and 1a is very small, and this explains why it is not observed in the supersonic expansion.

References

1) Cabezas, C.; Varela, M.; Alonso, J. L., "The Structure of the Elusive Simplest Dipeptide Gly-Gly." Angew. Chem. Int. Ed. 2017, 56, 6420-6425, DOI: 10.1002/anie.201702425.

InChIs

1: InChI=1S/C4H8N2O3/c5-1-3(7)6-2-4(8)9/h1-2,5H2,(H,6,7)(H,8,9)
InChIKey=YMAWOPBAYDPSLA-UHFFFAOYSA-N

'
This work is licensed under a Creative Commons Attribution-NoDerivs 3.0 Unported License.