Wednesday, July 22, 2015

Simulations of Chemical Reactions with the Frozen Domain Formulation of the Fragment Molecular Orbital Method

Hiroya Nakata, Dmitri G. Fedorov, Takeshi Nagata, Kazuo Kitaura, and Shinichiro Nakamura. Journal of Chemical Theory and Computation 2015, 11, 3053−3064
Contributed by +Jan Jensen

Reprinted (adapted) with permission from J. Chem. Theory Comput. 2015, 11, 3053−3064. Copyright (2015) American Chemical Society

I am highlighting this paper because it is the first paper I have come across that reports the optimized TS, free energy barrier, and IRC for an enzymatic reaction using an all-QM description of the entire enzyme.

The underlying method is the frozen domain with dimers (FDD) formulation of the fragment molecular orbital method (FMO) implemented in the GAMES program.  The FMO is fragmentation method and in the FDD formulation the system is divided into three domain: A) an active domain where the both the geometry and density are allowed to change, B) a polarizable buffer where only the density is allowed to change and F) a frozen domain where neither is allowed to change.  The latter restraints means that the fragments in the frozen domain only have to be calculated once.

This paper presents the implementation of the fully analytic gradient and Hessian at the RHF level and the interface to the IRC methodology in GAMESS.  (The interface to the TS-optimization code was reported previously using a different variant of FMO.)  In this approach only the Hessian of the active domain is computed.

The method is then applied to the keo-enol tautomerization of phosphoglycolohydroxamic acid (PGH) by triosephosphate isomerase (TIM).  The active domain is PGH and the proton donor Glu410 while the polarizable buffer is any residue or water molecule within 5.2 Å of the active domain, which were both treated at the RHF/6-31G(d) + dispersion level of theory.  The frozen domain, i.e. the rest of the TIM monomer, is treated at the RHF/STO-3G level of theory and the entire system is comprised of 9227 atoms.

The workflow was at follows using 128 2.93 GHz Xeon cores:

1. Constrained optimization to obtain TS guess (12 hrs)
2. Hessian calculation (24 hrs)
3. TS optimization (12 hrs)
4. Hessian calculation (24 hrs)
5. IRC calculations (16 points) in each direction (42 hrs)
6. Optimization of reactants and products starting from last IRC points (12 hrs each)
7. Hessian calculation for reactant and products (24 hrs each)

Of course, at this level of theory and with such a small active domain the resulting free energy barrier (21.4 kcal/mol) is considerably different from the experimental value (14 kcal/mol) but this is only a proof-of-concept study.

An important next step for FMO studies of enzyme catalysis will be the implementation of a mixed DFT/RHF-D model where all or part of the active domain can be treated with DFT.  The computational cost-increase associated with using DFT could be off-set by using, for example, HF-3c for the non-reactive part of of the active region as well as the buffer region since this method gives excellent structures as low computational cost.

I thank +Casper Steinmann for alerting me to this article

This work is licensed under a Creative Commons Attribution 4.0

Friday, July 17, 2015

On the tunneling instability of a hypercoordinated carbocation

Kozuch, S. Phys. Chem. Chem. Phys. 2015, 17, 16688-16691
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Last year I wrote a post on the possibility of a stable hypercoordinated carbon in the C(CH3)5+ molecule as proposed by Schleyer and Schaefer.1 Kozuch has re-examined this molecule with an eye towards examining the lifetime of this proposed “fleeting” molecule.2

The computed barriers for either (1) loss of a methane molecule leaving behind the (CH3)2C+CH2CH3cation or (2) loss of an ethane molecule leaving behind the t-butyl cation are small: 1.65 and 1.37 kcal mol-1, respectively. Kozuch employed canonical variational theory with and without small curvature tunneling (SCT). Without the tunneling correction, the pentamethylmethyl cation is predicted to have a long (millennia) lifetime at very low temperatures (<20 K). However, when tunneling is included, the half-life is reduced to 6 and 40 μs for degradation along the two pathways. Clearly, this is not a fleeting molecule – its lifetime is really too short to consider it as anything.

Interestingly, perdeuterating the molecule ((CD3)5C+) substantially increases the half-life to 4 ms, a thousand-fold increase. Tritium substitution would further increase the half-life to 0.2 s – a long enough time to really identify it and perhaps justify the name “molecule”. What is perhaps the most interesting aspect here is that H/D substitution has such a large effect on the tunneling rate even though no C-H bond is broken in the TS! This results from a mass effect (a CH3 vs. a CD3 group is migrating) along with a zero-point vibrational energy effect.


(1) McKee, W. C.; Agarwal, J.; Schaefer, H. F.; Schleyer, P. v. R. "Covalent Hypercoordination: Can Carbon Bind Five Methyl Ligands?," Angew. Chem. Int. Ed. 201453, 7875-7878, DOI: 10.1002/anie.201403314.
(2) Kozuch, S. "On the tunneling instability of a hypercoordinated carbocation," Phys. Chem. Chem. Phys.201517, 16688-16691, DOI: 10.1039/C5CP02080H.


C(CH3)5+: InChI=1S/C6H15/c1-6(2,3,4)5/h1-5H3/q+1

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

Thursday, July 16, 2015

Inverted Carbon Geometries: Challenges to Experiment and Theory

Bremer, M.; Untenecker, H.; Gunchenko, P. A.; Fokin, A. A.; Schreiner, P. R. J. Org. Chem. 2015, 80, 6520–6524
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Inverted carbon atoms, where the bonds from a single carbon atom are made to four other atoms which all on one side of a plane, remain a subject of fascination for organic chemists. We simply like to put carbon into unusual environments! Bremer, Fokin, and Schreiner have examined a selection of molecules possessing inverted carbon atoms and highlights some problems both with experiments and computations.1

The prototype of the inverted carbon is propellane 1. The ­Cinv-Cinv bond distance is 1.594 Å as determined in a gas-phase electron diffraction experiment.2 A selection of bond distance computed with various methods is shown in Figure 1. Note that CASPT2/6-31G(d), CCSD(t)/cc-pVTZ and MP2 does a very fine job in predicting the structure. However, a selection of DFT methods predict a distance that is too short, and these methods include functionals that include dispersion corrections or have been designed to account for medium-range electron correlation.

Figure 1. Optimized Structure of 1 at MP2/cc-pVTZ, along with Cinv-Cinv distances (Å) computed with different methods.

Propellanes without an inverted carbon, like 2, are properly described by these DFT methods; the C-C distance predicted by the DFT methods is close to that predicted by the post-HF methods.

The propellane 3 has been referred to many times for its seemingly very long Cinv-Cinv bond: an x-ray study from 1973 indicates it is 1.643 Å.3 However, this distance is computed at MP2/cc-pVTZ to be considerably shorter: 1.571 Å (Figure 2). Bremer, Fokin, and Schreiner resynthesized 3 and conducted a new x-ray study, and find that the Cinv-Cinv distance is 1.5838 Å, in reasonable agreement with the computation. This is yet another example of where computation has pointed towards experimental errors in chemical structure.

Figure 2. MP2/cc-pVTZ optimized structure of 3.

However, DFT methods fail to properly predict the Cinv-Cinv distance in 3. The functionals B3LYP, B3LYP-D3BJ and M06-2x (with the cc-pVTZ basis set) predict a distance of 1.560, 1.555, and 1.545 Å, respectively. Bremer, Folkin and Schreiner did not consider the ωB97X-D functional, so I optimized the structure of 3 at ωB97X-D/cc-pVTZ and the distance is 1.546 Å.
Inverted carbon atoms appear to be a significant challenge for DFT methods.


(1) Bremer, M.; Untenecker, H.; Gunchenko, P. A.; Fokin, A. A.; Schreiner, P. R. "Inverted Carbon Geometries: Challenges to Experiment and Theory," J. Org. Chem. 201580, 6520–6524, DOI:10.1021/acs.joc.5b00845.
(2) Hedberg, L.; Hedberg, K. "The molecular structure of gaseous [1.1.1]propellane: an electron-diffraction investigation," J. Am. Chem. Soc. 1985107, 7257-7260, DOI: 10.1021/ja00311a004.
(3) Gibbons, C. S.; Trotter, J. "Crystal Structure of 1-Cyanotetracyclo[,7.03,7]decane," Can. J. Chem.197351, 87-91, DOI: 10.1139/v73-012.


1: InChI=1S/C5H6/c1-4-2-5(1,4)3-4/h1-3H2
3: InChI=1S/C11H13N/c12-7-9-1-8-2-10(4-9)6-11(10,3-8)5-9/h8H,1-6H2

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

Wednesday, July 15, 2015

Computational and Experimental Investigations of the Formal Dyotropic Rearrangements of Himbert Arene/Allene Cycloadducts

Pham, H. V.; Karns, A. S.; Vanderwal, C. D.; Houk, K. N.  J. Am. Chem. Soc. 2015,137, 6956-6964
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

Houk and Vanderwal have examined the dyotropic rearrangement of an interesting class of polycyclic compounds using experimental and computational techniques.1 The parent reaction takes the bicyclo[2.2.2]octadiene 1 into the bicyclo[3.2.1]octadiene 3. The M06-2X/6-311+G(d,p)/B3LYP/6-31G(d) (with CPCM simulating xylene) geometries and relative energies are shown in Figure 1. The calculations indicate a stepwise mechanism, with an intervening zwitterion intermediate. The second step is rate determining.





Figure 1. B3LYP/6-31G(d) and relative energies (kcal mol-1) at M06-2X/6-311+G(d,p).

Next they computed the activation barrier for the second TS for a series of substituted analogs of 1, with various electron withdrawing group as R1 and electron donating groups as R2, and compared them with the experimental rates.
Further analysis was done by relating the charge distribution in these TSs with the relative rates, and they find a nice linear relationship between the charge and ln(krel). This led to the prediction that a cyano substituent would significantly activate the reaction, which was then confirmed by experiment. Another prediction of a rate enhancement with Lewis acids was also confirmed by experiment.
A last set of computations addressed the question of whether a ketone or lactone would also undergo this dyotropic rearrangement. The lactam turns out to have the lowest activation barrier by far.


(1) Pham, H. V.; Karns, A. S.; Vanderwal, C. D.; Houk, K. N. "Computational and Experimental Investigations of the Formal Dyotropic Rearrangements of Himbert Arene/Allene Cycloadducts," J. Am. Chem. Soc. 2015,137, 6956-6964, DOI: 10.1021/jacs.5b03718.


1: InChI=1S/C11H11NO/c1-12-10(13)7-9-6-8-2-4-11(9,12)5-3-8/h2-5,7-8H,6H2,1H3
2: InChI=1S/C11H13NO/c1-12-10(13)6-9-3-2-8-4-5-11(9,12)7-8/h4-6,8H,2-3,7H2,1H3

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

Thursday, June 25, 2015

Accurately Modeling Nanosecond Protein Dynamics Requires at least Microseconds of Simulation

Gregory R. Bowman Journal of Computational Chemistry 2015
Contributed by +Jan Jensen

This papers compares computed and experimental order parameters for two proteins, ubiquitin and RNase H, computed using 10, 100, 1000, and 10000 ns (10 $\mu$s) explicit molecular dynamics simulations.

Order parameters quantify the order of particular bonds (typically amide NH and methyl CH) and are typically measured via relaxation-dispersion NMR experiments that are insensitive to dynamics longer than the molecular tumbling time.  For ubiquitin and RNase H the tumbling times are 3.5 and 8.5 ns, respectively so the usual assumption would be that a 100 ns simulation is more than enough for both cases.

Figure 3. Correlation coefficients (r) and RMSDs between experimental backbone order parameters for RNase H and those calculated from simulations with the Amber03 force field. Mean values from 50 bootstrapped samples are shown as a function of the simulation length for the long-time limit approximation (open diamonds) and from the truncated average approximation (closed circles). The error bars represent one standard deviation. (c) 2015 John Wiley and Sons. Reproduced with permission.
Bowman shows that this assumption is not good.  For example, for RNase H (Figure 3 from the paper), 100 ns is barely adequate (especially for r) and a 1 $\mu$s simulation is a minimum requirement to demonstrate convergence. As Bowman points out:
Since order parameters are an ensemble average, they have contributions from all populated states, including those separated from the native state by high-energy barriers, which are unlikely to be accessed during nanoseconds of simulation. 
Using Amber99sb-ILDN and Charmm27 results in similar conclusions. The fact that the agreement with experiment continues to improve as as the length of the simulation increases also suggests that modern force fields predict that the experimentally observed protein structure is in fact a minimum of the free energy surface, which has not always been the case in the past.

This work is licensed under a Creative Commons Attribution 4.0

Tuesday, June 23, 2015

A Practicable Real-Space Measure and Visualization of Static Electron-Correlation Effects

Grimme, S. and Hansen, A., Angew Chem Int Edit, (2015)
Contributed by Tobias Schwabe

The question of how to deal with multireference (MR) cases in DFT has a longstanding history. Of course, the exact functional would also include multireference effects (or non-dynamical/non-local/static electron correlaton, as these effects are also called) and no special care is needed. But when it comes to today's density functional approximations (DFAs) within the Kohn-Sham framework, everything is a little bit more complicated. For example, Baerends and co-workers have shown that is the exchange part in GGA-DFAs that actually accounts for static electron correlation.[1] These studies, among others, led to the conclusion that the (erroneous) electron self-interaction in DFAs accounts for some of the MR character in a system. A good review about how these things are interconnected can be found in Ref. [2].

Instead of searching for better and better DFAs, another approach to the problem is to apply ensemble DFT which introduces the free electron energy and also the concept of entropy into DFT.[3] The key concept here is to allow for fractional occupation numbers in Kohn-Sham orbitals and to look at the system at T > 0 K. In case of systems with MR character which cannot be described with a single Slater determinant fractional occupation will result (for e.g. when computing natural occupation numbers). The interesting thing about ensemble DFT is that it allows to find these numbers directly via a variational approach without computing an MR wavefunction first.

Grimme and Hansen now turned this approach into a tool for a qualitative analysis of molecular systems. They do so by plotting what they call the fractional orbital density (FOD). That is, only those molecular orbitals with non-integer occupation numbers contribute to the density – and only at a finite temperature. This density vanishes completely at T = 0 K. So, the analysis literally shows MR hot spots. Integrating the FOD yields also an absolute scalar which allows to quantify the MR character and to compare different molecules. Due to the authors, this value correlates well with other values which attempts to provide such information.

A great advantage of the approach is that now the MR character can be located (geometrically) within the molecule. The findings presented in the application part of the paper go along well with chemical intuition. The analysis might help to visualize and to interpret MR phenomenon. The tool can provide insight when the nature of the electronic structure is not obvious – for example, when dealing with biradicals in a singlet spin state. It might also be a good starting point to identify relevant regions/orbitals which should be included when one wants to treat a system on a higher level than DFT, for example with WFT-in-DFT based on projector techniques. Last but not least, it can help to identify chemical systems to which standard DFAs should not (or only with great care) be applied.


[1] a) O. V. Grittsenko, P. R. T. Schipper, and E. J. Baerends, J. Chem. Phys. (1997), 107, 5007 b) P. R. T. Schipper, O. V. Grittsenko, and E. J. Baerends, Phys. Rev. A (1998), 57, 1729 c) P. R. T. Schipper, O. V. Grittsenko, and E. J. Baerends, J. Chem. Phys. (1999), 111, 4056

[2] A. J. Cohen, P. Mori-Sánchez, and W. Yang., Chem. Rev. (2012), 112, 289

[3] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press (1989)

Thursday, June 18, 2015

Ring Planarity Problem of 2-Oxazoline Revisited Using Microwave Spectroscopy and Quantum Chemical Calculations

Samdal, S.; Møllendal, H.; Reine, S.; Guillemin, J.-C. J. Phys. Chem. A 2015, 119, 4875–4884
Contributed by Steven Bachrach.
Reposted from Computational Organic Chemistry with permission

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-1with 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 2015119, 4875–4884, DOI: 10.1021/acs.jpca.5b02528.


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

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