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.

Saturday, November 15, 2014

Energy-conserving ab initio molecular dynamics

Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

Energy conservation in molecular dynamics is an important test of the propagators and parameters being used.  In first principles Born-Oppenheimer molecular dynamics (BOMD) this is complicated because of the need for self-consistent iteration to a consistent ground state charge density and potential.  The self-consistency cycle introduces errors since it is not possible to achieve a perfect ground state.  Moreover, to improve efficiency, the starting charge density is often predicted from previous charge densities; while this can results in impressive speed gains, it breaks the time-reversibility of the dynamics and will introduce an unavoidable energy drift.
This blog will discuss a set of papers[1-4] from Anders Niklasson in Los Alamos, which seek to remove this problem, and which have also shown a way to avoid the self-consistency cycle entirely.  The work is perhaps less well-known that it should be, and offers an excellent route to efficiency.  The recent linear scaling DFT molecular dynamics that we have published (going to over 32,000 atoms with DFT)[5] used the time-reversible formulation to conserve energy.
The time reversibility is built on an extended Lagrangian (and is known as XL-BOMD) which introduces a new dynamical variable, or set of variables, for the electronic degrees of freedom (whether charge density or density matrix) alongside the nuclear degrees of freedom, in a trick reminiscent of Car-Parrinello MD[6].  This degree of freedom, written below as P, evolves harmonically in time, depending on a frequency (set by the user) and the difference to the ground state charge density (or whatever is being used).
This degree of freedom, P, can be propagated in a time-reversible manner with an accuracy that only depends on the timestep.  The BOMD can thus be made time-reversible by starting the self-consistency cycle from P – resulting in stability and energy conservation in MD over long timescales.  On its own, this is a very valuable advance.
However, in a set of papers[2-4], Niklasson takes the idea further.  He suggests that, given the right functional for the energy, the self-consistency cycle can be removed.  This involves performing a single diagonalisation/minimisation without self-consistent update, and using P as the input charge density.  The scheme has been developed, with the latest paper[4] generalising the earlier work, and showing that the approach is applicable to systems regardless of gap, provided that an adiabatic separation between electronic and nuclear degrees of freedom can be made.
This work provides techniques to accelerate ab initio MD significantly; we have found that, in combination with linear scaling DFT, we can access picoseconds of MD for tens of thousands of atoms even with a simple initial implementation.  The two approaches together promise to make large-scale, long time accurate MD significantly more affordable.
[1] Phys. Rev. Lett. 100, 123004 (2008) DOI:10.1103/PhysRevLett.100.123004
[2] J. Chem. Phys. 139, 214102 (2013) DOI:10.1063/1.4834015
[3] J. Chem. Phys. 140, 044117 (2014) DOI:10.1063/1.4862907
[4] J. Chem. Phys. 140, 164123 (2014) DOI:10.1063/1.4898803
[5] arxiv 1409.6085 and J. Chem. Theor. Comput. in press (2014) DOI:10.1021/ct500847y

Solution-Phase Dimerization of an Oblong Shape-Persistent Macrocycle

Chu, M.; Scioneaux, A. N.; Hartley, C. S.  J. Org. Chem. 2014, 79, 9009–9017
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

π-π-stacking has been a major theme of my blog, and is discussed in Chapter 3.5.4 in the Second Edition of my book. Most examples involved molecules that are nearly circular (like benzene or triphenylene). Hartley and co-workers discuss the π-π-stacking of the oblong molecule 1, comparing its experimental features with computed features of the model compound 2.1
The key spectroscopic feature associated with assembly of 1 are the changes in the 1H chemical shifts with increasing concentration. For example, the chemical shift of the three protons on the triphenylene unit shift upfield by 0.30 to 0.66 ppm as the concentration increases from 10-5 to 10-2 M.

To see if these NMR shift changes are due to association of 1, they employed a computational approach. First they optimized the structure of model compound 2 at B3LYP/6-31G(d) (Shown in Figure 1a). Then using this fixed geometry, they computed the 1H chemical shifts of the dimer of 2. They explored the stacking distance (ranging from 3.2 to 4.0 Å along with varying the displacement of the two molecules along the major axis from 0.0 to 6.0 Å, finding the best fit to the chemical shifts with a separation of 3.6 Å and a displacement along the major axis of 3.5 Å. Using these two fixed values, they explored displacement of the molecules along the minor axis, along with rotation of the two molecules. The best fit to the experimental chemical shifts was with a displacement of 0.5 Å along the short axis and no rotation. This structure is shown in Figure 1b, with a RMS error of only 0.09 ppm from experiment. Models of the trimer show poorer fit to the experimental data.


Figure 1. B3LYP/6-31G(d) (a) optimized structure of 2 and the (b) structure of the best fit of the dimer of2. (As always, clicking on these images will allow you to manipulate the 3-D structure using JMol – highly recommended for the dimer.)

Using some smaller models and the B97-D functional, they argue that the displacement, which is substantially larger than the displacement found in stacked triphenylene, results from the need to minimize the steric interactions between the alkoxyl chains.


(1) Chu, M.; Scioneaux, A. N.; Hartley, C. S. "Solution-Phase Dimerization of an Oblong Shape-Persistent Macrocycle," J. Org. Chem. 201479, 9009–9017; DOI: 10.1021/jo501260c.


1: InChI=1S/C116H148O10/c1-11-21-31-41-59-117-95-71-87-51-55-91-75-97-99(103-81-111(121-63-45-35-25-15-5)115(125-67-49-39-29-19-9)85-107(103)105-83-113(123-65-47-37-27-17-7)109(79-101(97)105)119-61-43-33-23-13-3)77-93(91)57-53-89-70-90(74-96(73-89)118-60-42-32-22-12-2)54-58-94-78-100-98(76-92(94)56-52-88(69-87)72-95)102-80-110(120-62-44-34-24-14-4)114(124-66-48-38-28-18-8)84-106(102)108-86-116(126-68-50-40-30-20-10)112(82-104(100)108)122-64-46-36-26-16-6/h69-86H,11-50,59-68H2,1-10H3
2: InChI=1S/C76H68O10/c1-11-77-55-31-47-21-25-51-35-57-59(63-41-71(81-15-5)75(85-19-9)45-67(63)65-43-73(83-17-7)69(79-13-3)39-61(57)65)37-53(51)27-23-49-30-50(34-56(33-49)78-12-2)24-28-54-38-60-58(36-52(54)26-22-48(29-47)32-55)62-40-70(80-14-4)74(84-18-8)44-66(62)68-46-76(86-20-10)72(82-16-6)42-64(60)68/h29-46H,11-20H2,1-10H3

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

Thursday, November 6, 2014

Theoretical Foundation for the Presence of Oxacarbenium Ions in Chemical Glycoside Synthesis

Hosoya, T.; Takano, T.; Kosma, P.; Rosenau, T.  J. Org. Chem. 2014, 79, 7889-7894
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Comparing SN2 and SN1 reactions using computational methods is often quite difficult. The problem is that the heterolytic cleavage in the SN1 reaction leads to an ion pair, and in the gas phase this is a highly endothermic process. Optimization of the ion pair in the gas phase invariably leads to recombination. This is disturbingly the result even when one uses PCM to mimic the solvent, which one might have hoped would be sufficient to stabilize the ions.

The computational study of the glycoside cleavage by Hosoya and colleagues offers some guidance towards dealing with this problem.1 They examined the cleavage of triflate from 2,3,4,6-tetra-O-methyl-α-D-glucopyranosyl triflate 1.
Benchmarking the dissociation energy for the cleavage of 1 and considering computational performance, they settled on M06-2X/BS-III//M06-2X/BS-I, where BS-III is aug-cc-pVTZ basis set for O, F, and Cl and cc-pVTZ for H, C, and S and BS-I is 6-31G(d,p) basis sets were employed for H, C, and S, and 6-31+G(d) for O, F, and Cl. Solvent (dichloromethane) was included through PCM.

Attempted optimization of the contact ion pair formed from cleavage of 1 invariably led back to the covalent bound 1. PCM is not capable of properly stabilizing these types of ions in proximity. To solve this problem, they incorporated four explicit dichloromethane molecules. A minor drawback to their approach is that they did not sample much of configuration space and so their best geometries may not be the lowest energy configurations. Nonetheless, with four solvent molecules, they were able to locate contact ion pairs and solvent-separated ion pairs. Representative examples are shown in Figure 1. This method of explicit incorporation of a few solvent molecules seems to be the direction we must take to treat ions or even highly polar molecules in solution.




Figure 1. Representative examples of microsolvated 1, its contact ion pair (CIP) and solvent separated ion pair (SSIP) computed at M06-2X/BS-III//M06-2X/BS-I, and relative energies (kcal mol-1)


(1) Hosoya, T.; Takano, T.; Kosma, P.; Rosenau, T. "Theoretical Foundation for the Presence of Oxacarbenium Ions in Chemical Glycoside Synthesis," J. Org. Chem. 201479, 7889-7894, DOI:10.1021/jo501012s.


1: InChI=1S/C11H19F3O8S/c1-17-5-6-7(18-2)8(19-3)9(20-4)10(21-6)22-23(15,16)11(12,13)14/h6-10H,5H2,1-4H3/t6-,7-,8+,9-,10-/m1/s1

Update: be sure to read the interesting discussion associated with this post on Computational Organic Chemistry.

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

Friday, October 31, 2014

Quantum Monte Carlo Benchmark of Exchange-Correlation Functionals for Bulk Water

Miguel A. Morales, John R. Gergely, Jeremy McMinis, Jeffrey M. McMahon, Jeongnim Kim, and David M. Ceperley Journal of Chemical Theory and Computation 2014,  10, 2355–2362
Contributed by David Bowler
Reposted from Atomistic Computer Simulations with permission

I want to start by noting that the modelling of water is enormously complex and controversial, and that I’m far from an expert; in water, hydrogen bonding, dispersion interactions and quantum nuclear effects are all important.  The paper I want to discuss[1] performs detailed comparisons of the energies of different configurations of water molecules (drawn from a variety of dynamic simulations) calculated with several DFT functionals and quantum Monte Carlo (QMC).  While there have been steps in the direction of doing MD with QMC, it is still too expensive to do this routinely; static QMC calculations have been shown to model ice extremely accurately[2,3].
The configurations that are tested are drawn from classical forcefields and path-integral MD (PIMD) with van der Waals corrected DFT (see here for a recent overview of approaches to vdW in DFT).  The only test that is applied is the difference in total energy between the methods.  This is, probably, the only test that could be applied, as force calculations in QMC are very expensive, but it raises an interesting question: is it the convergence of total energy always the best thing to test ? It’s certainly true in normal DFT calculations that forces and energy differences converge faster than absolute energies, and many calculations use this fact to reduce the computational overhead.
The classical forcefields used rigid water molecules (a standard approach to speed up the simulations) and for these configurations, the DFT errors are around 0.1 mHa (3 meV) per water molecule, with only small variations between functionals.  This goes some way to explaining the success that can be found for some areas of water modelling with DFT (it’s important to note that dispersion effects are important and do reduce the errors).
When the configurations from PIMD are tested, however, these errors become much larger – around 0.3 mHa (10 meV) per water molecule for standard and vdW DFT methods.  Only the hybrid methods maintain a low error (but they are still too expensive to use routinely with MD).  Why does this happen ? Well, it appears that the major sources of error between DFT and QMC comes from the structure of the water molecule itself; when it is allowed to be flexible, configurations are generated that DFT describes poorly, but hybrids model better.
This failure can be explored further, using a correction scheme for modelling water in DFT[4] (where the error between DFT and QMC or coupled cluster calculations is decomposed in a many body expansion).  They find that the major term contributing to the error is the one-body term, and show that the correction scheme gives excellent results.  This confirms the difference between the errors for rigid and flexible water molecules.
The paper gives us an excellent insight into the limitations of DFT, and how it can be improved.  It’s important to note that there are different ways of testing: PIMD simulations with vdW-DF2 have been shown to reproduce the experimental water structure almost perfectly[5].  It’s easy to get trapped into thinking that one particular approach or technique is best, or to become obsessed with methods: but experiment is always what must be tested, though it’s often far from trivial to work out how to compare to experiments.
[1] J. Chem. Theory Comput. 10, 2355 (2014) DOI:10.1021/ct500129p
[2] Phys. Rev. Lett. 107, 185701 (2011) DOI:10.1103/PhysRevLett.107.185701
[3] J. Chem. Phys. 138, 221102 (2013) DOI:10.1063/1.4810882
[4] J. Chem. Phys. 139, 244504 (2013) DOI:10.1063/1.4852182
[5] arXiv:1402.2697

INDO/X: A New Semiempirical Method for Excited States of Organic and Biological Molecules

Alexander A. Voityuk Journal of Chemical Theory and Computation 2014, ASAP
Contributed by +Jan Jensen

This paper presents a significantly more accurate re-parameterization of the popular INDO/S method. The re-parameterization is based on the TBE-2 data set: vertical excitation energies and oscillator strengths computed for the valence excited states of 28 medium-sized organic molecules using multistate-CASPT2 and aug-cc-pVTZ basis set.

The author sums it up succinctly: "The mean absolute deviations of the INDO/X excitation energies relative to the TBE-2 data is 0.26 eV for singlet states and 0.33 eV for triplet states. The corresponding values for INDO/S are 0.56 and 0.64 eV. ...  Unlike INDO/S where two different Hamiltonians are employed for singlet and triplet excitations, both types of excited states are calculated within the same INDO/X scheme."  The MADs are thus similar to the MAD for TD-B3LYP calculations (0.27 eV) and significantly better than TD-BP86 (0.53 eV) and TD-BHLYP (0.48 eV) calculations.

Like INDO/S, INDO/X uses the CIS approach and can therefore be applied to quite large molecules. Furthermore, the author notes that the results are quite insensitive to the chosen active space.

Thanks to @CompChemNews for alerting me to this paper.

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