Calculating accurate energy barriers#

Besides studying the mechanistic details of chemical reactions, ORCA can also be used to calculate accurate energy barriers, and thus predict reaction rates to a certain accuracy.

It must be clear here that it is rather difficult to compare the absolute values of the experimental energy barriers with the calculates ones. First, because these "experimental" values are never measured directly, but obtained from measured data only after a series of assumptions.

Second because, unless your are modeling the solvent explicitly, there will always be a relatively large error associated to solvation. And that is to say the least if one also does not account for conformational entropy, transmission coefficients and so on.

Study case: the Diels-Alder reaction#

Still, it does make sense to compare relative reaction energy barriers and rates, and we will show here how one could do that for a classic Diels-Alder (DA) reaction between cyclopentadiene and some dieneophiles:

That is an important reaction, that conserves all of the atoms from the reactants, and particularly useful to synthesize chiral compounds. Its reaction rate varies drastically with the nature of the reactants, from \(10^{-8}\) to \(10^{2} s^{-1}\) [Guo2012], and its known that electron-withdrawing groups on the dieneophile accelerate the formation of the cyclic products:

Let's try to investigate the relative barrier heights for the reaction between cyclopentadiene and the cyanide compounds above using the NEB-TS method to find the transition states, and investigate the effect of adding proper Correlation energy to it.

First, the transition state#

The transition state (TS) for DA reactions is nothing unusual, except for the fact that three double bonds are being broken and two single bonds are being formed at the same time, which is sometimes hard to track using the usual TS-search algorithms that rely on a single coordinate search.

By using ORCA's NEB-TS method, the TS search is made much simpler, and all that is needed are the structures of the reactant and products. To optimize these using, e.g, B3LYP and a good DEF2-TZVP basis, one can use:

!B3LYP D4 DEF2-SVP OPT FREQ CPCM(TOLUENE)
* XYZFILE 0 1 4CN_reactant.xyz

Here we are using the D4 correction [Grimme2017] to account for dispersion interaction, which is key in this case, and the CPCM solvation model [Truhlar2009] to include solvation effects to some extent.

After checking that both reactant and product are real minima, without any negative frequencies, one can already start the NEB-TS calculation using the optimized structures:

!B3LYP D4 DEF2-SVP NEB-TS FREQ CPCM(TOLUENE)
%NEB
      NEB_END_XYZFILE "4CN_product_optimized.xyz"
END
* XYZFILE 0 1 4CN_reactant_optimized.xyz

The new keywords here are the NEB-TS and FREQ on the main input, and the name of the file containing the product geometry under %NEB. If everything goes well, a TS having a single negative frequency of about \(-368cm^{-1}\) will be automatically found:

Calculating the energy barrier#

With the TS and reactant in hands, one can immediately compute the \(\Delta G^{\ddagger} = G^o_{TS} - G^o_{reactant}\) from the Gibbs free energy difference between them. Here is a table summarizing the results for the three cyanides shown above:

Summary of the energy barriers (\(\Delta G^{\ddagger}\)) calculated to far, in kcal/mol.#
Compound	B3LYP	Exp.
2 CN	9.90	17.7
3 CN	10.69	16.2
4C	10.21	13.6

Well, that is actually not very helpful, since there is not even a clear trend! There is no clear trend resulting from the B3LYP calculation, however, the geometries look fine and the transition states are indeed saddle-points on the potential energy surface. Can we get something better than this?

Using DLPNO-CCSD(T) to correct the electronic energy#

In general DFT is known to reproduce geometries and frequencies with reasonable quality for its low cost, but the energies are definitely a weak point. And they are even worse if the system studied is outside the training set of the functionals, which normally do not included transition states or atypical bonding situations.

There is where higher level ab initio methods should be used to compute the electronic energy of these systems. Correlated methods such as CCSD(T) are unparametrized, based on real physical principles and should truly reproduce the experimental results.

Using those methods together with the The DLPNO scheme developed by the ORCA team allows for fast a accurate calculations of these barriers task. To compute an accurate electronic energy in this case, we could take the DFT geometry and run:

!DLPNO-CCSD(T) DEF2-TZVPP DEF2-TZVPP/C
* XYZFILE 0 1 4CN_reactant_optimized.xyz

for both the reactant and the TS.

Converting \(E_{el}\) to \(G^o\)#

In order to transform the DLPNO-CCSD(T) electronic energy (\(E_{el}\)) into a true \(G^o\), we have to include the vibrational corrections already discussed in the Thermodynamics section. As DFT predicts frequencies relatively well, we can use the correction obtained previously here.

To account for solvation, we can also take the \(\Delta G_{solv}\) presented in the Implicit Solvation Models section, as computed with DFT and add them up to obtained a better final \(\Delta G^{\ddagger}\):

\[ \Delta G^{\ddagger}_{solv} = E_{el}(DLNPO-CC) + \Delta G_{correction} (DFT) + \Delta G^o_{solv} (DFT) \]

Using these results, our updated table is now:

Summary of the energy barriers (\(\Delta G^{\ddagger}\)) calculated, in kcal/mol.#
Compound	B3LYP	CCSD(T)	Exp.
2 CN	9.90	14.71	17.7
3 CN	10.69	12.73	16.2
4 CN	10.31	9.48	13.6

And the patter is clearly there. We can even plot a graphic to show that there is indeed a linear correlation between the calculated and the experimental energy barriers:

That means now we could, in principle, use this model to predict the experimental reaction barriers for unknown compounds with quite good accuracy, by calculating them and adjusting them to the fit!

Important

Again: it is very unlikely that the predicted and measured energy barriers will be equal, unless by coincidence. We not necessarily calculating the same quantity that was measured, but there should be relations like the one demonstrated above.