Heavy Atom Tunneling in Organic Reactions at Coupled Cluster Potential Accuracy with a Parallel Implementation of Anharmonic Constant Calculations and Semiclassical Transition State Theory

We describe and test on some organic reactions a parallel implementation strategy to compute anharmonic constants, which are employed in semiclassical transition state theory reaction rate calculations. Our software can interface with any quantum chemistry code capable of a single point energy estimate, and it is suitable for both minimum and transition state geometry calculations. After testing the accuracy and comparing the efficiency of our implementation against other software, we use it to estimate the semiclassical transition state theory (SCTST) rate constant of three reactions of increasing dimensionality, known as examples of heavy atom tunneling. We show how our method is improved in efficiency with respect to other existing implementations. In conclusion, our approach allows SCTST rates and heavy atom tunneling at a high level of electronic structure theory (up to CCSD(T)) to be evaluated. This work shows how crucial the possibility to perform high level ab initio rate evaluations can be.


Electronic Structure Calculations and Optimized Geometries
In our paper, we performed Single Point electronic structure Energy (SPE) calculations and geometry optimizations using the Gaussian16 software. 1 The geometry optimizations are carried out requiring the lowest possible residual gradient in the optimization process. The following input keywords are set as default in our program for the geometry optimization step: • SCF=(XQC,Tight) and Int=UltraFine (for DFT calculations) • Opt=(VeryTight,CalcAll,MaxCycles=200) (add TS option for the transition state opt) For the CCSD(T) optimization we started from the CCSD optimized geometry, for which we have analytic gradients available. Then we switched to numerical gradients for CCSD (T) calculations with the eigenvalue-following algorithm for the optimization process. A first calculation is done using the standard setup, with a maximum step size of 0.1 Bohr. Once the Tight convergence is reached, new optimization parameters have been used for the VeryTight convergence criteria. The following overlay 1 options have been used for the R1 reactant

Accuracy Tests
We report here the results from the direct comparison of our anharmonic constants with the ones obtained using Gaussian16. The definitions of the difference ∆χ and the percentage difference % dif f are given in the article in eq. 30.
(a) DFT difference ∆χ (a) DFT difference ∆χ (a) DFT difference ∆χ (a) DFT difference ∆χ (a) DFT difference ∆χ (a) DFT difference ∆χ harmonic constants calculated with our program and the Gaussian software is shown in the following results for the three reactions studied in our paper. Figure 9: SCTST Reaction Rate Constants at different temperatures for the R1 reaction. The anharmonic constants matrix was calculated using our program and the Gaussian16 software starting from the same geometry and ab initio calculation options. DFT calculations have been carried out using the B3LYP functional with the 6-31G* basis set while in the MP2 calculations we used the aug-cc-pVDZ one. Figure 10: SCTST Reaction Rate Constants at different temperatures for the R2 reaction. The anharmonic constants matrix was calculated using our program and the Gaussian16 software starting from the same geometry and ab initio calculation options. DFT calculations have been carried out using the B3LYP functional with the 6-31G* basis set, while in the MP2 calculations we used the aug-cc-pVDZ one. Figure 11: SCTST Reaction Rate Constants at different temperatures for the R3 reaction. The anharmonic constants matrix was calculated using our program and the Gaussian16 software starting from the same geometry and ab initio calculation options. DFT calculations have been carried out using the B3LYP functional with the 6-31G* basis set, while in the MP2 calculations we used the aug-cc-pVDZ basis set.
The barriers here reported are the classical barriers corrected for the harmonic Zero Point k=1 ω k and corrected for the anharmonic ZPE: where N is the number of vibrational degrees of freedom, ω k are the harmonic vibrational frequencies, G R 0 and G T S 0 are the constant terms coming from the VPT2 perturbation expansion and R and TS refer to the formula for the reactant and the transition state.   For the hindered rotational modes (HRM) treatment, involved in the low frequency modes of the R2 and R3 reactants, we used the Pitzer and Gewin (PG) coefficients for the correction of the harmonic oscillator partition function.
The overall correction values for the total harmonic partition functions associated with the HRM are listed in Tables (9) and (10).  The absolute values of the rate constants for the three reactions at different temperatures and ab initio levels of theory are here shown in the reported tables 1 . Table 11: TST and SCTST reaction rate constants at the MP2/aug-cc-pVDZ, CCSD(T)/aug-cc-pVDZ and B3LYP/6-31G* level of theory for the R1 reaction.  Table 12: TST and SCTST reaction rate constants at the MP2/aug-cc-pVDZ, B3LYP/6-31G* and CCSD/jun-cc-pVDZ level of theory for the R2 reaction.  The values of the % difference between the SCTST and the TST rate constants are reported in the following tables.