Posts Tagged ‘Molecular dynamics’

An Ambimodal Trispericyclic Transition State: the effect of solvation?

Thursday, May 2nd, 2019

Ken Houk’s group has recently published this study of cycloaddition reactions, using a combination of classical transition state location followed by molecular dynamics trajectory calculations,[1] and to which Steve Bachrach’s blog alerted me. The reaction struck me as being quite polar (with cyano groups) and so I took a look at the article to see what both the original[2] experimental conditions were and how the new simulations compared. The reaction itself is shown below.


Turns out that chloroform was used as solvent (also benzene), whilst the transition state calculations and the subsequent molecular dynamics trajectories were modelled for the gas phase. The key observation is that if TS1 is used as the starting point for trajectory calculations, only 87% lead to the product predicted by classical transition state theory (3), as revealed by a classical intrinsic reaction coordinate (IRC) calculation. The remaining 13% lead to 4 and 5, the ambimodal effect. So here, I want to explore what effect including a continuum solvent on the computation of TS1 and its IRC might have on the classical (non-dynamic) model.

Firstly, the model for TS1 as reported[3], ωB97XD/6-31G(d) (FAIR data at DOI: 10.14469/hpc/5590). The basis set is modest by today’s standards, but is largely imposed by the need to use it for very large numbers of trajectory calculations. I was able to copy/paste the coordinates from the reported supporting information and then to replicate the IRC at this level (gas phase), also shown in the SI.

There is a feature in this IRC I want to expand upon (red arrow above). It occurs at an IRC value of ~-0.9 on the axis below.

It can be seen more clearly if the RMS gradient norm is plotted, the value of which drops to almost zero at IRC -0.9. Had it reached exactly 0.0, we would have had an intermediate formed. As it is we have what is called a hidden intermediate. The origins of this intermediate can be more readily inferred from this dipole moment plot along the IRC. At TS1, the dipole moment is > 10D. A rule of thumb I have often used is that if a TS has a DM > 10, then one cannot ignore solvation any more and a gas phase model must be augmented.

Here are the same plots, but now with an added solvent field for water, an extreme polarity.

The IRC now stops at -2, being a high energy true (ionic) intermediate (rather than a hidden one). The IRC stops because gradient norm has now reached 0.0 at IRC -2, again an indication of a true intermediate.

The dipole moment has been increased from 10 in the gas phase to around 15 at TS1 and it continues to increase until the ionic intermediate is reached. These differences from the gas phase plots are induced entirely by applying a continuum solvent model.

Next an intermediate solvent, chloroform, being one of the solvents used for the actual reaction. This time the gradient norm almost reaches a value of 0.0, avoiding it only by a whisker! A barely hidden intermediate.

The dipole moment totters around 13.5D, before finally collapsing as the ionic intermediate itself collapses to a neutral molecule again. Benzene as solvent (not shown here) reaches an intermediate dipole moment of about 12D. It too can stabilize an ionic intermediate noticeably even though it is not ionic itself.

I want to also briefly explore what effect if any the use of a relatively small basis set (6-31G(d)) has on the shape of the IRC. Below is a repeat of the gas phase IRC using the Def2-TZVPP basis, which is about twice the size of the smaller one (and hence is around 16 times slower to compute). The gradient norm shows that the “hidden intermediate” region around IRC -1 is a little more prominent (flatter).

So we see that both a solvent model (as a continuum field) and a larger basis set can increase the degree of “hidden intermediate” character in the classical reaction coordinate for this cycloaddition reaction, to the extent that if water is used as model solvent an actual discrete albeit shallow ionic intermediate forms. As Houk puts it, solvation induces a conversion from an entropic intermediate to an enthalpic one.[4] Molecular dynamics trajectories however have a propensity for not settling into quite shallow intermediates (those with escape barriers of < 3 kcal/mol, as would be the case here).

It will indeed be interesting to see the extent, if any, that either of the augmented models shown above affect the calculated distribution of molecular dynamics trajectories compared to those obtained using a gas phase model.

References

  1. X. Xue, C.S. Jamieson, M. Garcia-Borràs, X. Dong, Z. Yang, and K.N. Houk, "Ambimodal Trispericyclic Transition State and Dynamic Control of Periselectivity", Journal of the American Chemical Society, vol. 141, pp. 1217-1221, 2019. https://doi.org/10.1021/jacs.8b12674
  2. C.Y. Liu, and S.T. Ding, "Cycloadditions of electron-deficient 8,8-disubstituted heptafulvenes to electron-rich 6,6-disubstituted fulvenes", The Journal of Organic Chemistry, vol. 57, pp. 4539-4544, 1992. https://doi.org/10.1021/jo00042a039
  3. https://doi.org/
  4. O.M. Gonzalez-James, E.E. Kwan, and D.A. Singleton, "Entropic Intermediates and Hidden Rate-Limiting Steps in Seemingly Concerted Cycloadditions. Observation, Prediction, and Origin of an Isotope Effect on Recrossing", Journal of the American Chemical Society, vol. 134, pp. 1914-1917, 2012. https://doi.org/10.1021/ja208779k

Tags:, , , , , , , , , , , , ,
Posted in reaction mechanism | 2 Comments »

Reaction coordinates vs Dynamic trajectories as illustrated by an example reaction mechanism.

Monday, March 20th, 2017

The example a few posts back of how methane might invert its configuration by transposing two hydrogen atoms illustrated the reaction mechanism by locating a transition state and following it down in energy using an intrinsic reaction coordinate (IRC). Here I explore an alternative method based instead on computing a molecular dynamics trajectory (MD).

I have used ethane instead of methane, since it is now possible to envisage two outcomes:

An animation of the IRC starting from the located transition state is shown below (DOI: 10.14469/hpc/2331). This is based purely on the computed potential energy surface of the molecule. The IRC is computed from the forces experienced on the atoms as they are displaced from an initial set of coordinates corresponding to the located transition state and then following the direction indicated by the eigenvectors of the negative force constant required of a transition state. Importantly, there is no time component; the path is based entirely on energies and forces.

Next, a molecular dynamics simulation (ωB97XD/6-31G(d,p), DOI: 10.14469/hpc/2330).  This uses the ADMP method, which requests a classical trajectory calculation using the “atom-centered density matrix propagation molecular dynamics model”. This integrates kinetic energy contributions from the molecular vibrations and so explicitly now includes a time component. In this example the evolution of the system from the transition state is charted over a period of 100 femtoseconds (1000 integrated steps). As it happens this is a relatively short period of evolution; sometimes periods of picoseconds may be required to get a realistic model, especially if one is also dealing with explicit solvent molecules (of which perhaps 500 might be required).

Such explicit inclusion of the kinetic energy from molecular vibrations in effect allows the molecule to “jump” over shallow barriers. In this case, the barrier for a [1,2] hydrogen shift from the methyl group to the adjacent carbene (watch atom 8). Simultaneously, the path taken by two hydrogens no longer corresponds to their transposition but to their elimination as a hydrogen molecule. So this quite different outcome from the IRC is very probably also a much more realistic one.

If the MD method is so much more realistic than the IRC, then why is it not always used? The simple answer is computational time! For this very small molecule and using quite a modest basis set (6-31G(d,p)), the relatively short 1000 time steps took about three times as long to compute as the IRC. The factor gets worse as the size of the molecule increases and the number of time steps for a realistic result increases. Perhaps, as technology gets better and new computer architectures emerge, MD simulations of ever increasingly complex reactions will become far more common. In ten years time, I expect most of the examples on this blog will use this method!

Tags:, , , , , , , , , , , , , , , ,
Posted in reaction mechanism | 2 Comments »

Dynamic effects in nucleophilic substitution at trigonal carbon.

Monday, July 16th, 2012

Singleton and co-workers have produced some wonderful work showing how dynamic effects and not just transition states can control the outcome of reactions. Steve Bachrach’s blog contains many examples, including this recent one.

This shows that tolyl thiolate (X=Na) reacts with the dichlorobutenone to give two substitution products in a 81:19 ratio. Singleton and Bogle argue[1] that this arises from a single transition state, and that the two products arise from a statistical distribution of dynamic trajectories bifurcating out of a transition state favouring 2 over 3. Steve puts it very elegantly “I think most organic chemists hold dear to their hearts the notion that selectivity is due to crossing over different transition states“. When I read this, Occam’s razor came to mind: could a simpler (in this case a more conventional) answer in fact be better? 

My thoughts in fact followed a point I have been making here recently, the principle that modelling a complete system is probably better than a partial one. Now, if you look at Figure 1 of the Singleton/Bogle article, captioned “Qualitative energy surface for the reaction of 1 with sodium p-tolyl thiolate” I was struck by something missing; the sodium (X=Na), and possibly also explicit solvent (ethanol). I wanted to see if these missing components may influence the mechanism.

The red arrows follow the proposed mechanism (a), whereas the blue arrows represent a more conventional 1,4-nucleophilic addition to form an intermediate enol anion, this then eliminating to the final product. Singleton & co. explored the potential energy surface using the following computational model: B3LYP/6-31+G(d,p)/PCM(ethanol) for the anionic system (defined by setting an overall charge of -1 during the calculation), finding the potential energy surface corresponded to path (a). They then went on to explore the dynamics of the system emerging out of this single TS, showing that in fact both products would be formed in more or less exactly the ratio observed. 

I thought two things could be considered missing from this model; X+ (the counterion) and explicit solvent (continuum solvent was invoked using the PCM model). On the latter point, I have thought for a little while that there are two types of solvent; those which act via their dielectric field, and those that act via hydrogen bonds. Ethanol does both, and so in this case (I argue) it should be explicitly included (actually, in the first instance it can be approximated using water instead of ethanol). The missing counter-ion is a greater challenge. In what follows I am going to approximate it too, using H+ (Na+ itself I reserve for a future post). The objective is to find out what (if anything) changes when this more complete model is built. It is shown below as the first transition state encountered. Its features include:

First transition state TS1. Click for 3D

  1. Two explicit solvent (water) molecules.
  2. A H+ (I will discuss Na+ in another post), attached to one of the water molecules as a hydronium ion.
  3. The hydronium ion bridges to the carbonyl group (this is the final optimised position; the second water molecule serves only to H-bond to the hydronium ion). 
  4. This overall system is neutral, charge=0 (I like to say it might be found in a bottle or flask; pure anions of course cannot be bottled). 
  5. The model used was B3LYP/6-31+G(d,p)/CPCM(ethanol); I find the CPCM method to be better for calculating intrinsic reaction coordinates (IRC).
  6. Using this transition state to initiate an IRC shows that the presence of this solvent bridge allows X (=H+) to smoothly transfer from sulfur to oxygen as part of a concerted process. This avoids excessive build up of charge separation.
  7. This now forms an intermediate (we are clearly following path (b) and not path (a) now). This is because the enolate anion is stabilised by protonation and a hydrogen bond from the proton to the solvent water, and so this becomes an explicit intermediate in the potential energy surface.

    Intermediate in reaction. Click for 3D

This intermediate now collapses along path (b) to the final product, via the transition state shown below. Again, an IRC shows a solvent bridge allows X to be concertedly transferred, this time from the oxygen to form hydronium chloride and 2 (I have not yet found the equivalent pathway to 3, but given the hydrogen bonds involved it is bound to be different).

Second transition state TS2. Click for 3D 

ΔG† (kcal/mol) along this sequence is 1 (0.0), TS1 (28.2), Int (8.0), TS2 (12.6); the intermediate existing only in a shallow well of 4.6 kcal/mol. The activation barrier is on the high side (the reaction occurs easily at room temperature), and it might be expected that (in part) this might be due to using X=H+ rather than X=Na+ for the model. Watch this space!

What might we conclude from this? That the presence of additional molecules (H3O+ and H2O) can result in structures which can depend on other features of the molecule, in this case the carbonyl group, one that plays little role in mechanism (a). In path (b), the carbonyl group is far from passive, receiving and then releasing X during the course of the reaction. This must mean that the transition state for forming product 2 may indeed be a separate one from the transition state for forming product 3, since the relationship of these two to the carbonyl is different. To re-quote Steve again “I think most organic chemists hold dear to their hearts the notion that selectivity is due to crossing over different transition states”.

Perhaps the explanation might indeed be due to different transition states rather than different dynamics? Clearly, more research needs to be done; I for one do not regard the case as closed on this example just yet.

References

  1. X.S. Bogle, and D.A. Singleton, "Dynamic Origin of the Stereoselectivity of a Nucleophilic Substitution Reaction", Organic Letters, vol. 14, pp. 2528-2531, 2012. https://doi.org/10.1021/ol300817a

Tags:, , , , , , , ,
Posted in Interesting chemistry | 3 Comments »