Discovering chemistry with an ab initio nanoreactor.

Lee-Ping Wang, Alexey Titov, Robert McGibbon, Fang Liu, Vijay S. Pande and Todd J. Martínez

Nature Chemistry

6, 1044–1048 (2014)

Contributed by Alán Aspuru-Guzik

When is quantum chemistry going to take the primary discovery role, rather than an explanatory one? This one of the questions asked in the recent Nature Chemistry paper by the groups of Vijay Pande and Todd Martínez at Stanford.

With the advent of fast GPGPU quantum chemistry in codes such as Terachem, it is practical to carry out medium-scale quantum dynamics simulations for a considerable amount of time. In this paper, the authors employ a spherical "piston" to compress and decompress high-temperature reaction mixtures, which in turn result in an accelerated process for triggering chemical reactions. Automated algorithms are employed to recognize the nature of the chemical transformations involved. The authors use this approach to identify novel mechanisms for the Urey-Miller origins of life experiment.

In this modern world of Facebook, Twitter, Youtube, Facebook and blogs, a video is worth more than a thousand words. The video below summarizes the paper very well.

This paper is indeed, one of my top 5 favorite papers in the field of computational chemistry of 2014. 

Todd Martínez explaining the future of the nanoreactor in the QUITEL 2014 conference held in the Galápagos Islands in Ecuador. He suggested, in a joking fashion, that the nanoreactor could eventually evolve a flask with molecules into life ab initio. The virtual end product of the evolution seems to be the contributor of this post.

How Electronic Dynamics with Pauli Exclusion Produces Fermi-Dirac Statistics.

Triet S. Nguyen, Ravindra Nanguneri, and John Parkhill, arxiv:14115324v1 (2014)
Contributed by Alán Aspuru-Guzik

A very interesting preprint by John Parkhill from the University of Notredame and his group provides a solution for a long-standing theoretical problem: Many of the quantum master equations that treat electron-phonon interactions fail to reach a Fermi-Dirac distribution at long times. For example, Ehrenfest dynamics is one of the usual suspects for this behavior. By employing a combination of time-dependent perturbation theory and the extended normal ordering theory of Kutzelnigg and Mukherjee, Parkhill and his co-workers arrive to a novel master equation.

The results from this work show that the Pauli principle "blocks" electron relaxation and slows it down, and in addition it leads to a time-dependent relaxation rate. The plot below shows the correct asymptotic relaxation of a double excitation in a model system and the plot in the right shows a time-dependent relaxation rate.

The paper says that "A paper using these new formulas in a useful all-electron dynamics code is

an immediate follow up.". The contributor is looking forward to it.

Computational Chemistry: 2014 in numbers

In 2014 we learnt that two of the germinal DFT papers (by Becke and Lee, Yang and Parr) are amongst the top ten most cited scientific papers of all time, and of chemistry papers published in the last ten years, the fourth and fifth most cited again relate to computational research (Truhlar and Hess, respectively). In this vein I thought it would be interesting to perform a (pseudo)-scientific analysis of the usage of computation in chemistry research in, and in the years leading up to, 2014 as judged by bibliometric data.

Searching all 2014 chemistry papers in the Web of Science for mention of "computation" or "computational" in either the article title, abstract or keywords suggests that approximately 2.7% of chemistry research involved computation of some variety this year (9,101 of a staggering 331,699 papers). This is most likely an underestimate since searching for more specific phrases such as "DFT" will turn up more hits. The same analysis over previous years reveals a steady increase in the proportion of chemistry research using computation from 0.6% in 1994, to 1.1% in 2004 and 2.2% in 2010.

Around 20% of all the computational chemistry papers published in 2014 emanate from the USA, more than double the closest competitor, China. The top ten nations in terms of publications are USA 19.5%, China 9.3%, Germany 6.1%, India 4.3%, France 4.0%, Italy 3.8%, Spain 3.7%, England 3.6%, Japan 2.7% and Canada 2.4% - making nearly 60% of the total output. A decade ago in 2004 the ten most prolific countries accounted for around 87% of total output, which indicates that recent years have witnessed a greater global involvement  in computational chemistry. Noticeable trends are seen in individual nations share of the computational chemistry pie, with the USA and some European nations effectively halving their fraction of papers between 2004 and 2014, with China's output nearly doubling from 5.4% in 2004 to 9.3% in 2014 and the emergence of India from outside the top ten into fourth place in 2014. It should be borne in mind that the globalization of science will inevitably lead to some over counting of papers here, due to multiple addresses appearing on the same paper.

Two review articles for the general audience

Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

n trying to clean up my in-box of articles for potential posts, I write here about two articles for a more general audience, authored by two of the major leaders in computational organic chemistry.

Ken Houk offers an overview of how computational simulation is a partner with experiment and theory in aiding and guiding our understanding of organic chemistry.¹ The article is written for the non-specialist, really even more for the non-scientist. Ken describes how computations have helped understand relatively simple reactions like pericyclic reactions, that then get more subtle when torquoselection is considered, to metal-catalysis, to designed protein catalysts. If you are ever faced with discussing just what you do as a computational chemist at a cocktail party, this article is a great resource of how to explain our science to the interested lay audience.

Paul Schleyer adds a tutorial on transition state aromaticity.² The authors discusses a variety of aromaticity measures (energetics, geometry, magnetic properties) that can be employed to analyze the nature of transition states, in addition to ground state molecules. This article provides a very clear description of the methods and a few examples. It is written for a more specialized audience than Houk’s article, but is nonetheless completely accessible to any chemist, even those with no computational background.

References

(1) Houk, K. N.; Liu, P. "Using Computational Chemistry to Understand & Discover Chemical Reactions,"Daedalus 2014, 143, 49-66, DOI: 10.1162/DAED_a_00305.

(2) Schleyer, P. v. R.; Wu, J. I.; Cossio, F. P.; Fernandez, I. "Aromaticity in transition structures," Chem. Soc. Rev. 2014, 43, 4909-4921, DOI: 10.1039/C4CS00012A.

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

Benchmarking Ab Initio Binding Energies of Hydrogen-Bonded Molecular Clusters Based on FTIR Spectroscopy

Nicolai Bork, Lin Du, Heidi Reiman, Theo Kurten, and Henrik G. Kjaergaard Journal of Physical Chemistry A 2014, 118, 5316-5322
Contributed by +Jan Jensen

One of the holy grails of computational chemistry is a method that consistently can compute standard binding free energies in bulk solvent to within 1 kcal/mol.  When using electronic structure theory one of the (many) potential problems is the use of the harmonic approximation because low-frequency modes can make a sizable contribution to the vibrational entropy. This paper provides some valuable benchmarking data.

Kjærgaard and co-workers measure the gas phase binding free energy at 295 K of acetonitrile-HCl to be between 1.2 and 1.9 kcal/mol. CCSD(T)/aug-cc-pV(T+d) predicts 2.3 kcal/mol using the harmonic frequencies and 1.9 kcal/mol using the anharmonic scaling factors suggested by Shields and co-workers.  So, at least in this case, the harmonic approximation works quite well, but can be improved quite easily using anharmonic scaling factors.  If this effect is additive and transferable, then making such a correction will be important for systems that bind using several hydrogen bonds of the 1 kcal/mol accuracy-target is to be met.

The authors write "We find that this energy lowering originates almost entirely from the scaling of the two lowest frequency modes (between 28 and 35 cm$^{-1}$, depending on method), emphasizing the importance of calculating the low frequency modes accurately."

Indeed, such low frequency modes are susceptible to numerical noise (e.g. from coarse DFT grids or sloppy optimization convergence criteria) which can even lead to them turning up as imaginary frequencies in the vibrational analysis. Because imaginary frequencies are not included in the vibrational free energy, this can lead to sizable errors in the binding free energy.  For example, a 28 cm$^{-1}$ frequency contributes 1.8 kcal/mol to the free energy at 295 K.

B3LYP-D3, MP2 and PW91/aug-cc-pV(T+d) calculations also yield binding free energies within the experimentally measured range, but in the last two cases this is due for fortunate cancellation of errors in the electronic and free energy components based on comparison to the CCSD(T) calculations.  Having said that, with the exception of B3LYP, all functionals tested in this study are able to predict the binding energy to within 1 kcal/mol using the aug-cc-pV(T+d) basis set.

It is worth noting that the approach by Shields and co-workers find scaling factors by fitting to harmonic vibrational energy and entropy-expressions evaluated using anharmonic fundamental frequencies.  It is not at all obvious (at least to me) that this approach will result in good agreement with experimental binding free energies involving 28 cm$^{-1}$ frequencies, so this (and a previous) study provides very valuable validation of this general approach.



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

Hypercubane: DFT-based prediction of an Oh-symmetric double-shell hydrocarbon

Pichierri, F. Chem. Phys. Lett. 2014, 612, 198-202
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Three-dimensional objects can be projected into four-dimensional objects. So for example a cube can be projected into a hypercube, as in Scheme 1.

Scheme 1.

Pichierri proposes a hydrocarbon analogue of the hypercube. The critical decision is the connecting bridge between the outer (exploded) carbons. This distance is too long to be a single carbon-carbon bond. Pichierri opts to use ethynyl bridges, to give the hypercube 1.¹

Now, unfortunately he does not supply any supporting materials. So I have reoptimized this O_h geometry at B3LYP/6-31G(d), and show this structure in Figure 1. Pichierri does not report much beyond the geometry of 1 and the perfluoronated analogue. One interesting property that might be of interest is the ring strain energy of 1, which I will not take up here.

1

2

But a question I will take up is just what bridges might serve to create the hydrocarbon hypercube. A more fundamental choice might be ethanyl bridges, to create 2. However, the O_h conformer of 2 has 13 imaginary frequencies at B3LYP/6-31G(d). Lowering the symmetry to D₃ give a structure that has only real frequencies, and it’s shown in Figure 1. An interesting exercise is to ponder other choices of bridges, which I will leave for the reader.

1

2

Figure 1. B3LYP/6-31G(d) optimized structures of 1 and 2.
As always, be sure to click on the image to enable Jmol for interactive viewing of these interesting structures!

References

(1) Pichierri, F. "Hypercubane: DFT-based prediction of an O_h-symmetric double-shell hydrocarbon,"Chem. Phys. Lett. 2014, 612, 198-202, DOI: j.cplett.2014.08.032.

InChIs

1: InChI=1S/C40H24/c1-2-26-7-9-29-15-11-27-5-3-25(1)4-6-28-12-16-30(10-8-26)20-23-32(22-19-29)24-21-31(17-13-27,18-14-28)39-35(27)33(25)34(26)37(29,35)40(32,39)38(30,34)36(28,33)39/h1-24H
InChIKey=FFMFUIDOGFAUOP-UHFFFAOYSA-N

2: InChI=1S/C40H48/c1-2-26-7-9-29-15-11-27-5-3-25(1)4-6-28-12-16-30(10-8-26)20-23-32(22-19-29)24-21-31(17-13-27,18-14-28)39-35(27)33(25)34(26)37(29,35)40(32,39)38(30,34)36(28,33)39/h1-24H2
InChIKey=MCSZKKKJCDSRIV-UHFFFAOYSA-N

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

Benchmarks for Characterization of Minima, Transition States, and Pathways in Atomic, Molecular, and Condensed Matter Systems

Samuel T. Chill, Jacob Stevenson, Victor Ruehle, Cheng Shang, Penghao Xiao, James D. Farrell, David J. Wales, and Graeme Henkelman Journal of Chemical Theory and Computation 2014 10, 5476–5482
Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

One of the key tasks that atomistic simulations of all types perform is to characterise a potential energy surface.  This includes diverse problems, such as:

• local optimisation (finding a stable structure given a starting set of coordinates)
• global optimisation (finding the most stable structure for a given set of atoms – a very hard task !)
• transition state searches, both starting from one minimum (single-ended) and interpolating between two minima (double-ended)

These are hugely complex problems, and there are many methods which have been proposed to perform these searches. I have discussed some of these before (here, here,  and here); there is also an overview of searching methods in Chapter 5 of the book.

The paper I’m discussing in this post[1] proposes some benchmarks to test different methods for these optimisation problems.  They also offer a website[2] where the tests can be found, and new results shared.  The problems are all based around simple empirical potentials (which makes the testing relatively light).  They propose the following systems:

• a 38 atom Lennard-Jones cluster
• a 100 atom, two species Lennard-Jones cluster
• water clusters with 10, 15 and 20 atoms modelled with TIP4P
• a seven atom cluster on the (111) surface of a FCC material (Morse potential)
• a bulk sample of the same FCC material

These different problems cover a reasonable set of problems, but there is no reason why more should not be added to the set.  For instance, the present set lacks any examples of covalent bond-breaking systems.

Various methods are compared, including L-BFGS, FIRE and conjugate gradients, nudged elastic band, eigenvector following and basin hopping.  There are so many methods that it would be impossible to cover all, but I’m a little surprised that methods such as AIRSS and metadynamics weren’t tested for the global optimisation problems.

The local optimisation problem offers the simplest problem to characterise, and the conclusions are simplest here: methods that include curvature (L-BFGS and conjugate gradients) are significantly faster.  This is not particularly surprising, as these methods are known to be well optimised, and the test suite avoids any pathological problems.  The conclusions for the other challenges are less clear; the global optimisation problem in particular seems not to distinguish well between the two methods used.

The saddle point searching results surprised me.  I use the nudged-elastic band (NEB) method regularly, and one of the authors is a leading figure in the development of the method, but it did not perform particularly well.  The study did show that the performance is dependent on the initial pathway guessed, and suggested that for double-ended searches an preliminary NEB search should then be refined by the eigenvector following method.  The eigenvector following method performed best on the single-ended searches as well.

What do we learn from the paper ? While the overall conclusions confirm expectations for local optimisation and do not discriminate as much as would be hoped for the other problems, it’s an excellent start.  In particular, in seeking to offer open, standard tests for all the different methods that are available, it offers a way to move forward in this field.  It also suggests an excellent new model of characterisation-based scientific papers: a dynamic database of freely-available results.

[1] J. Chem. Theor. Comput. 10, 5476 (2014) DOI: 10.1021/ct5008718

[2] http://optbench.org/

Gas-Phase Preparation of Carbonic Acid and Its Monomethyl Ester

Reisenauer, H. P.; Wagner, J. P.; Schreiner, P. R.  Angew. Chem. Int. Ed. 2014, 53, 11766-11771
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

I remain amazed at how regularly I read reports of structure determinations of what seem to be simple molecules, yet these structures have eluded determination for decades if not centuries. An example is the recently determined x-ray crystal structure of L-phenylalanine;¹ who knew that growing these crystals would be so difficult?

The paper I want to discuss here is on the gas-phase structure of carbonic acid 1.² Who would have thought that preparing a pure gas-phase sample would be so difficult? Schreiner and co-workers prepared carbonic acid by high-vacuum flash pyrolysis (HVFP) of di-tert-butyl carbonate, as shown in Scheme 1.

Scheme 1

Carbonic acid can appear in three difference conformations, shown in Figure 1. The two lowest energy conformations are separated by a barrier of 9.5 kcal mol^-1 (estimated by focal point energy analysis). These conformations can be interconverted using near IR light. The third conformation is energetically inaccessible.

1cc (0.0)	1ct (1.6)	1tt (10.1)
2cc		2cc

Figure 1. CCSD(T)/cc-pVQZ optimized structures of 1 (and the focal point relative energies in kcal mol^-1) and the CCSD(T)/cc-pVTZ optimized structures of 2.

The structures of these two lowest energy conformations were confirmed by comparing their experimental IR spectra with the computed spectra (CCSD(T)/cc-pVTZ) and their experimental and computed rotational constants.

An interesting added component of this paper is that sublimation of the α- and β-polymorphs of carbonic acid do not produce the same compound. Sublimation of the β-isomorph does produce 1, but sublimation of the α-isomorph produces the methylester of 1, compound 2 (see Figure 1). The structure of 2 is again confirmed by comparison of the experimental and computed IR spectra.

References

(1) Ihlefeldt, F. S.; Pettersen, F. B.; von Bonin, A.; Zawadzka, M.; Görbitz, C. H. "The Polymorphs of L-Phenylalanine," Angew. Chem. Int. Ed. 2014, 53, 13600–13604, DOI: 10.1002/anie.201406886.

(2) Reisenauer, H. P.; Wagner, J. P.; Schreiner, P. R. "Gas-Phase Preparation of Carbonic Acid and Its Monomethyl Ester," Angew. Chem. Int. Ed. 2014, 53, 11766-11771, DOI: 10.1002/anie.201406969.

InChIs

1: InChI=1S/CH2O3/c2-1(3)4/h(H2,2,3,4)
InChIKey=BVKZGUZCCUSVTD-UHFFFAOYSA-N

2: InChI=1S/C2H4O3/c1-5-2(3)4/h1H3,(H,3,4)
InChIKey=CXHHBNMLPJOKQD-UHFFFAOYSA-N

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

Tautomerism in Neutral Histidine

Bermúdez, C.; Mata, S.; Cabezas, C.; Alonso, J. L. Angew. Chem. Int. Ed. 2014, 53, 11015-11018
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

The Alonso group has yet again (see these posts) determined the gas-phase structure of an important, biologically significant molecule using a combination of exquisite microwave spectroscopy and quantum computations. This time they examine the structure of histidine.¹

They optimized four conformations of histidine, as its neutral tautomer, at MP2/6-311++G(d,p). These are schematically drawn in Figure 1. Conformer 1a is the lowest in free energy, likely due to the two internal hydrogen bonds. Its structure is shown in Figure 2.

Figure 1. The four conformers of histidine. The relative free energy (MP2/6-311++G(d,p)) in kcal mol^-1are also indicated.

Figure 2. MP2/6-311++G(d,p) optimized geometry of 1a.

The initial experimental rotation constants were only able to eliminate 1b from consideration. So they then determined the quadrupole coupling constants for the ¹⁴N nuclei. These values strongly implicated 1a as the only structure in the gas phase. The agreement between the experimental values and the computed values at MP2/6-311++G(d,p) was a concern, so they rotated the amine group to try to match the experimental values. This lead to a change in the NHCC dihedral value of -16° to -23° Reoptimization of the structure at MP2/cc-pVTZ led to a dihedral of -21° and overall excellent agreement between the experimental spectral parameters and the computed values.

It is somewhat disappointing the supporting materials does not include the structures of the other three isomers, nor the optimized geometry at MP2/cc-pVTZ.

References

1) Bermúdez, C.; Mata, S.; Cabezas, C.; Alonso, J. L. "Tautomerism in Neutral Histidine," Angew. Chem. Int. Ed. 2014, 53, 11015-11018, DOI: 10.1002/anie.201405347.

InChIs

Histidine: InChI=1S/C6H9N3O2/c7-5(6(10)11)1-4-2-8-3-9-4/h2-3,5H,1,7H2,(H,8,9)(H,10,11)/t5-/m0/s1
InChIKey=HNDVDQJCIGZPNO-YFKPBYRVSA-N

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