Reproducibility in science: calculated kinetic isotope effects for cyclopropyl carbinyl radical.

Previously on the kinetic isotope effects for the Baeyer-Villiger reaction, I was discussing whether a realistic computed model could be constructed for the mechanism. The measured KIE or kinetic isotope effects (along with the approximate rate of the reaction) were to be our reality check. I had used ΔΔG energy differences and then HRR (harmonic rate ratios) to compute[cite]10.5281/zenodo.19272[/cite] the KIE, and Dan Singleton asked if I had included heavy atom tunnelling corrections in the calculation, which I had not. His group has shown these are not negligible for low-barrier reactions such as ring opening of cyclopropyl carbinyl radical.[cite]10.1021/ja1055593[/cite] As a prelude to configuring his suggested programs for computing tunnelling (GAUSSRATE and POLYRATE), it was important I learnt how to reproduce his KIE values.[cite]10.1021/ja1055593[/cite] Hence the title of this post. Now, read on.

cp

I felt I could contribute to the cause by extending the published results in two respects:

  1. The reported[cite]10.1021/ja1055593[/cite] calculations are for the model B3LYP/6-31G(d) but the article does not report the tolerance to e.g. basis set variation (6-31G(d), a modest basis set by 2015 standards),
  2. or to the quantum model used (B3LYP, a veritable DFT method).

These two model chemistries can both be tested by “increasing” their accuracy. The Def2-QZVPP basis set is nearing the CBS, or complete basis set limit. The coupled-cluster CCSD(T) method is regarded as the gold standard for single reference calculations. The CASSCF method tests the response to a multi-reference wave function. Each is applied separately to ensure only one variable is being changed at a time.

Method Expt. KIE[cite]10.1021/ja1055593[/cite] Pred. KIE (my result) Pred. ΔG298 Pred. KIE[cite]10.1021/ja1055593[/cite] KIE + Tunnelling correction[cite]10.1021/ja1055593[/cite]
B3LYP/6-31G(d)[cite]10.14469/ch/191357[/cite],[cite]10.14469/ch/191358[/cite] 1.079295 1.0582 8.0 1.058 1.073
1.163173 1.1067 1.106 1.169
B3LYP/Def2-QZVPP[cite]10.14469/ch/191353[/cite],[cite]10.14469/ch/191352[/cite] 1.079295 1.0563 6.6 1.058 1.073
1.163173 1.1031 1.106 1.169
CASSCF(5,5)/6-31G(d)[cite]10.14469/ch/191361[/cite],[cite]10.14469/ch/191364[/cite] 1.079295 1.0572 8.2 1.058 1.073
1.163173 1.1050 1.106 1.169
CASSCF(5,5)/Def2-TZVPP[cite]10.14469/ch/191363[/cite],[cite]10.14469/ch/191362[/cite] 1.079295 1.0561 7.9 1.058 1.073
1.163173 1.1028 1.106 1.169
CCSD(T)/6-31G(d)[cite]10.14469/ch/191367[/cite],[cite]10.14469/ch/191356[/cite] 1.079295 1.0597 9.7 1.058 1.073
1.163173 1.1099 1.106 1.169

Actually separate ratios of 13C/12C(C-4)/13C/12C(C-3) since C-3 and C-4 are not equivalent in the reactant species because of the methylene group pyramidalisation. The KIE calculation input and outputs are archived.[cite]10.5281/zenodo.19949[/cite]

The first two rows of table are my attempt at an exact replication of the literature. The start point of such a project would be the supporting information or SI[cite]10.1021/ja1055593[/cite] which contains coordinates for the program GAUSSRATE and defines key structures in the form of a double-column, page thrown (broken might be a better word) PDF file. It was going to be a bit of a struggle to reconstitute this format into the structure required for a Gaussian calculation, so I simply constructed the models from scratch and optimised to the ring-opening transition state[cite]10.14469/ch/191358[/cite] and reactant.[cite]10.14469/ch/191357[/cite] I used a more recent version of the Gaussian program (G09/D.01 rather than G03/D.02) to do this, and tightened up some of the criteria to modern cutoff standards. A continuum solvent model could have been specified  (the solvent used in the experiments was 1,2-dichlorobenzene) but since no mention was made of solvent, I assumed a gas phase calculation had originally been done.  The starting geometry of the reactant deliberately had no symmetry, but during optimisation it converged to having a plane of symmetry using the B3LYP/6-31G(d) level of theory (the SI does not note this symmetry, it is implicit). I then used my code[cite]10.5281/zenodo.19272[/cite] to compute the isotope effects. The KIE program used in the original literature calculation was not directly mentioned in the supporting information but is presumed to be Quiver. Dan Singleton has recently sent me these codes, but they still need to be compiled and tested at my end. I ended up with splendid agreement for the KIE as you can see above (top two lines). Its reproducible! Hence the various assumptions I made in achieving this appear justified.

Returning to the geometry of the cyclopropyl carbinyl radical as having a plane of symmetry, two of the other methods, CCSD(T)/6-31G(d) and CASSCF(5,5)/6-31G(d), as well as CASSCF(5,5) at the better Def2-TZVPP basis all predicted that the methylene radical is twisted by about 20° with respect to the Cs plane of the ring.

cp-asymm

It is useful to check whether this twisting has any impact on the predicted KIE. The answer is clear (Table). ALL the methods predict similar KIE to ± 0.003, which is as about as accurate as can be measured experimentally at the 1σ level of confidence. This is a remarkable result; few other computed molecular properties turn out to be so insensitive to the quantum procedure used. The next stage will be to check if the tunnelling corrections required to bring the calculation into congruence with the measured values are similarly insensitive.


The “barrier height” is quoted as 7 kcal/mol[cite]10.1021/ja1055593[/cite]. This is probably NOT the activation free energy.

Tags: , , , , , , ,

2 Responses to “Reproducibility in science: calculated kinetic isotope effects for cyclopropyl carbinyl radical.”

  1. Henry Rzepa says:

    I noted in the post that all the calculations were for a gas phase model. One further useful addition would be to add a continuum solvent model (o-dichlorobenzene at the higher temperatures) to see what effect it has on the B3LYP/6-31G(d) model. At 295K the KIE is 1.0581 (vs 1.0582) and at 173K it is 1.1065 (vs 1.1067). The difference is only in the 4th decimal place. Mind you, the dipole moment of the TS is <1D and solvation only starts to impact at >8D or so. The free energy barrier with solvent is 7.8 kcal/mol, only 0.2 kcal/mol down from the gas phase value.

  2. Dan Singleton says:

    Henry – You have hit upon an aspect of calculated isotope effects that is a bit surprising and yet is general and very important in their utility.

    The Bigeleisen-Mayer theory of isotope effects relates them entirely to the frequencies. Small changes in structure or energetics can have a big effect on the low frequencies, so one would guess that such changes would affect the isotope effects significantly. They do not. When I first starting playing with calculated isotope effects about 1992, I spent a lot of time “messing with” the low-frequencies by putting in constraints on low-energy modes, and found that the effect on the isotope effect was nearly negligible. Later, when we introduced using isotope effects to measure transition state distances, we evaluated the Schramm process of simply fixing distances and calculating the isotope effects from structures with constrained optimizations. From a purist standpoint, there are lots of reasons to think this should fail entirely, but in practice it does not. The Bigeleisen formalism is very forgiving because errors in small frequencies are both minimized in their effect and subject to error canceling.

    Ultimately, isotope effects without tunneling are intimately related to the geometry of the transition state. It is not obvious that this should be the case, but our JACS 2009 paper showed this pretty clearly. The predictions are nearly independent of the computational model, as you have now seen for yourself.

    Outside of high-tunneling regimes, a 1D tunneling analysis based on simply the sharpness of the barrier (as calculated entirely from the imaginary frequency) performs strikingly well in predicting tunneling. The effect of tunneling will only vary significantly when the sharpness of the barrier goes up or down a lot. As long as the overall reaction barrier is approximately correct, the tunneling effect on the KIE will change little. We have reported that barriers that are off by more than 5 kcal/mol can lead to significant errors in the predicted KIES.

Leave a Reply