Thursday, December 18, 2014

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.

Wednesday, December 17, 2014

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.1

Now, unfortunately he does not supply any supporting materials. So I have reoptimized this Oh 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.


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 Oh conformer of 2 has 13 imaginary frequencies at B3LYP/6-31G(d). Lowering the symmetry to D3 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.


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!


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


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

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

Friday, December 12, 2014

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 (herehere,  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

Wednesday, December 10, 2014

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;1 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.2 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.





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 (see Figure 1). The structure of 2 is again confirmed by comparison of the experimental and computed IR spectra.


(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. 201453, 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. 201453, 11766-11771, DOI: 10.1002/anie.201406969.


1: InChI=1S/CH2O3/c2-1(3)4/h(H2,2,3,4)
2: InChI=1S/C2H4O3/c1-5-2(3)4/h1H3,(H,3,4)

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

Tuesday, December 2, 2014

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.1

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 14N 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.


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


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

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

Sunday, November 23, 2014

Quantum Chemical Approach to Estimating the Thermodynamics of Metabolic Reactions

Contributed by +Jan Jensen 

DOI: 10.1038/srep07022 CC-BY-NC-ND

Computational chemistry is finally in a position to conduct studies on hundreds or even thousands of molecules within a reasonable amount of time.  This has led to benchmark datasets which have been used to develop new semi-empirical methods (e.g. DFT-Dx or HF-3c) as well as high-throughput screening of molecules with desirable properties. 

This paper represents a third kind of study: a concerted effort to fill in missing experimental data needed to quantitatively understand a complex chemical system - in this case cellular metabolism. Standard Gibbs reaction free energies of biochemical reactions are crucial to quantitative modeling of metabolism (metabolic flux analysis) but only a small fraction have been measured experimentally.  

Many of the reactions involve molecule-sizes that can be modeled routinely by electronic structure calculations.  However, there are new challenges compared to, for example, atmospheric chemistry such as the modeling of pH- and solvation-effects and multiply charged ions.  Thus, any conformational search must often be repeated for several different protonation states and multiply charged ions often necessitate the treatment of explicit solvent molecules and, even then, result in relatively solvation energies with relatively large errors.  Also, multiply charged anions may present additional problems for DFT due to incomplete cancellation of the inter-electronic Coulomb repulsion at larger distances. 

This paper represents the first, pioneering, attempt to address these issues with high-throughput in mind and the results are encouraging. For reactions involving non-multiply charged species B3LYP/6-31G* with the COSMO solvation model augmented by explicit water molecules result in MADs less than 5 kcal/mol, while being fast enough to compute more than a hundred reaction energies (involving more than 5000 quantum calculations). 

However, this paper is more a call to action than the definitive answer.  Many new methodological improvements (perhaps such as this one?) are needed in order to increase the accuracy and efficiency needed produce a complete and accurate picture of cellular metabolism. I hope this paper will spur such developments - hopefully in the form of large and open collaborations (perhaps similar to this one?).  

Thanks to Adrian Jinich for alerting me to this paper

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

Wednesday, November 19, 2014

Isotope Effects, Dynamic Matching, and Solvent Dynamics in a Wittig Reaction. Betaines as Bypassed Intermediates

Chen, Z.; Nieves-Quinones, Y.; Waas, J. R.; Singleton, D. A.  J. Am. Chem. Soc. 2014, 136, 13122-13125
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

If you hadn’t noticed, I am a big fan of the work that Dan Singleton is doing concerning the role of dynamics in discerning reaction mechanisms. Dan’s group has reported another outstanding study combining experiments, traditional QM computations, and molecular dynamics – this time on the Wittig reaction.1

The key question concerning the mechanism is whether a betaine intermediate is accessed along the reaction (path A) or whether the reaction proceeds in a concerted manner (path B). Earlier computations had supported the concerted pathway (B).

Experimental determination of the heavy atom kinetic isotope effect was made for Reaction 1.

Reaction 1

Using the 6-31+G(2df,p) basis set, three different density functionals predict three different potential energy surfaces. With M06-2x, the surface indicates path A (stepwise), with the first step rate-limiting. B3P86 also predicts the stepwise reaction, but the second step is rate-limiting. The Lc-wPBE functional predicts a concerted reaction. Using these surfaces, they predicted the carbon isotope effect and compared it to the experimental values. The best agreement is with the M06-2x surface with a weighting of the vibrational energies of the two different TSs. The optimized structures of the two transition states, the betaine intermediate, and the product are shown in Figure 1.




Figure 1. M06-2x/6-31+G(2df,p) optimized geometry of the critical points of Reaction 1.

The agreement of the predicted and experimental KIE is not ideal. So, they performed molecular dynamics computations with the ONIOM approach using M06-2x/6-31G* for Reaction 1 and 53 THF molecules treated at PM3. 360 trajectories were begun in the region of the first transition state (TS1), and they can be organized into 4 groups. The first group (128 trajectories) are reactions that produce product. The second group (76 cases) form the C-C bond but then it ruptures and returns to reactant. The third group (82 cases) have an immediate recrossing back to reactant, and the last group (16 cases) takes product back to the first TS and then returns to product. The predicted KIE using this weighted MD results gives values in outstanding agreement with the experiments.

Of the first group, about 50% pass from TS1 to TS2 in less than 150 fs, or in other words look like a concerted path. But a good number of trajectories reside in the betaine region for 1-2 ps.

In contrast, trajectories initiated from the betaine with equilibrated THF molecules indicate a median of 600 ps to travel from TS1 to TS2 and do not resemble a concerted path.

They argue that this bimodal distribution is in part associated with a solvent effect. When the first TS is crossed the solvent molecules are not equilibrated about the solute, and 10-20% of the trajectories immediately pass through the betaine region due to “dynamic matching” where the entering motion matches with exiting over the second transition state. The longer trajectories result from improper dynamic matching, but faster motion in the solute than motion amongst the solvent needed to stabilize the betaine. So, not only do we need to be concerned about dynamic effects involving the reactants, we need to be concerned about dynamics associated with the solvent too!


(1) Chen, Z.; Nieves-Quinones, Y.; Waas, J. R.; Singleton, D. A. "Isotope Effects, Dynamic Matching, and Solvent Dynamics in a Wittig Reaction. Betaines as Bypassed Intermediates," J. Am. Chem. Soc. 2014136, 13122-13125, DOI: 10.1021/ja506497b.

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