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.

Introduction

Calculation of reaction rates in theoretical chemistry is still nowadays a challenging task. The rate is an intrinsically dynamical quantity, and rigorous methods to compute reaction rates need to rely on dynamics simulations.[1] However, this approach is complicated by the low probability of reactive events occurring in a typical time span of dynamics simulations.[2,3] Moreover, quantum effects, such as tunneling and zero point energy, are recognized to have a significant impact on the rate constant value, especially at low temperatures. Therefore, in the cases where quantum effects are important, the quantum mechanical evolution of the system is mandatory.[4] Obviously, this is a very complicated task,[5−12] and it is not feasible to apply exact quantum methods for practical purposes, for instance, in kinetic modeling applications of complex systems.

Transition state theory (TST) is a clever rate constant approximation that avoids dynamics simulations and delivers rates in terms of static thermodynamics information.[13,14] Since TST is a classical mechanics theory, early theories based on one-dimensional potential approximation were elaborated to account for quantum tunneling, such as Wigner or Eckart corrections.[15,16] Nowadays, more sophisticated approximations have been developed to include, at least to some extent, the effects neglected by 1D approaches to tunneling corrections and limitations of the TST method itself, such as corner cutting, nonseparability of the reaction coordinate, and recrossing.[17−28]

Among these techniques, semiclassical transition state theory (SCTST) initially developed by W. H. Miller in the 1970s and revisited in the 1990s[29−32] has received renewed attention in the past few years.[33−37] What makes SCTST particularly convenient for application is that it requires input quantities that are routinely calculated by quantum chemistry codes. These include the harmonic vibrational frequencies, the height of the reaction barrier, and the anharmonic vibrational coupling constants, which are employed in the context of second order vibrational perturbation theory (VPT2). A user-friendly implementation of SCTST is provided along with the Multiwell program suite.[38−40] Recently, the computational convenience of this program has been enhanced by a parallel implementation of the calculation of the vibrational density of states,[41,42] which are used for the SCTST rate calculations. However, a main computational bottleneck still remains, and it is the calculation of the anharmonic couplings, which we will call below χ. For these reasons approximate reduced dimensionality versions of the SCTST theory have been explored.[43−45]

In this work we propose instead to retain the full dimensional anharmonic couplings matrix to calculate the SCTST rates and develop a convenient parallel implementation to speed up their estimate. The possibility to parallelize this task has already been exploited for spectroscopic applications of VPT2, but it has been applied only to minimum geometries on the potential energy surface (PES).[46] Our implementation extends the method for accelerated anharmonic constant calculations to transition state geometries. In addition, some more specific developments, such as the inclusion of Coriolis ro-vibrational couplings in the calculation of the total anharmonic constants matrix and a detailed treatment of the Fermi resonances along with the deperturbation process, are presented and implemented.

Eventually, our software allows us to compute SCTST rate constants of reactions of medium-high dimensionality at a high level of electronic structure theory with an affordable computational effort. We apply our implementation to the full-dimensional rate calculation of three heavy atom tunneling reactions, respectively composed of 10, 14, and 16 atoms, for which a fitted PES is not available. We carried out our electronic structure calculations at the level of density functional theory (DFT) and second order Møller–Plesset perturbation theory (MP2). For the 10-atom system, we were able to use coupled cluster with a full treatment of singles and doubles and an estimate to the connected triples contribution (CCSD(T)) level of electronic structure theory, while we used CCSD for the 14- and 16-atom systems.

The paper is organized as follows. In Section , we recall the SCTST theory and the expressions of the anharmonic couplings in the case of both a minimum structure and a transition state one. In Section , we describe in details our implementation, and we show the speedup performance of our parallel program. In Section , we describe some applications, and we compare our results with those obtained by other techniques. Eventually, in the Summary and Conclusions section (Section ), we provide our final remarks and anticipate some future perspective.

Semiclassical Transition State Theory

The exact kinetic rate constant can be written as[47−49]where N(E) is the cumulative reaction probability (CRP), h is the Plank constant, Q(T) is the reactants partition function, and β = 1/(kT), where T and k are, respectively, the temperature and the Boltzmann constant. In the SCTST frame the CRP can be written using a multidimensional generalization of the one-dimensional WKB (Wenzel, Kramer, Brillouin) tunneling probability P(E) approximation:[50]where we consider a molecular transition state composed by N atoms with N = 3N – 6 vibrational degrees of freedom. In eq , is a vector of N – 1 quantum integer numbers that defines the vibrational state of the transition state bound coordinates, and ∑ = ∑ ∑... ∑ stands for the sum over all vibrational states. Θ(E, n) is the barrier penetration integral which needs to be determined to calculate N(E) and the rate constant. For nonseparable systems, the total energy E can always be expressed asand in the case of a transition state structure by means of the Bohr Sommerfeld quantization rule,[30] the following identity holds for the Nth reactive modeThis last relation allows to write the total energy as a function of the penetration integral E(n1, ..., n, Θ), and by inversion it is in principle possible to get Θ(E, n).

Miller and Hernandez[31] practically addressed this inversion problem by exploiting the standard perturbative expression[51] for the vibrational energy levels in a minimum of the PES given bywhere V0′ is the potential energy at the stationary point of the PES with the inclusion of a constant term arising from the derivation of this expression in VPT2 context.[35,52,53] ω are the normal-mode frequencies, and χ are the anharmonic constants. Using eq and generalizing eq to the case of a saddle point geometry, we can find an explicit form for the barrier penetration integral as a function of the total energy and the quantum numbers (n1, n2, ..., n) related to the bound degrees of freedom:whereFor the Nth imaginary mode we have that ω = iω̃, χ = −iχ̃, and .

Using eq , we evaluate the sum in eq and get the CRP N(E). Then, by putting the calculated N(E) into eq , we finally get the rate constant.

In this work we compute the rate with the Multiwell program suite by separating the contribution from the different degrees of freedom to the partition functions.[38] The SCTST rate constant is evaluated aswhere Q(T) is the transition state (reactant) translational partition function and Q(T) is the transition state (reactant) rotational one. These are approximated to the corresponding free motion partition functions,where s is the rotational symmetry number, Iα is the moment of inertia along the α = x, y, z axis, and M is the total mass of the reactant or the transition state.

The vibrational partition function of the reactants is fully coupled and anharmonic since it is written in terms of the reactant density of vibrational states (DOS) ρ(E) asFor the numerator in eq , the calculation of the semiclassical N(E) is not trivial. A practical way to address this problem is to divide the energy range of interest into bins of width δE.[33] In this way a certain number # of energy levels, each one identified by a combination of quantum vibrational numbers = , will be found in the jth bin, and the corresponding average reaction probability is defined asTherefore, eq is rewritten aswhere ρ†(E) is the vibrational density of states (DOS) associated with the real-valued frequency vibrations of the TS. Note that, as δE is reduced, the result becomes more accurate. Within this approximation, the sum in eq over the accessible states is replaced by the easier sum over E/δE energy bins with energy lower than or equal to the total energy E. As a result, the rate constant calculation is reduced to the problem of evaluating two vibrational DOSs, ρ(E) for the reactants and ρ†(E) for the TS. To achieve a convenient computational effort to estimate the SCTST rate, some of us have recently implemented a parallel version of the density of vibrational states (DOS).[41,42] However, the main computational bottleneck in the entire approach still remains the evaluation of the anharmonic constants χ.

We now survey how to derive an explicit and convenient form for the anharmonic constants that we will use in our implementation. Equation is a Dunham expansion[54] of the energy for a molecular system where vibrational and rotational motions are coupled. The corresponding rovibrational Hamiltonian, which includes the rotational kinetic energy, has the following form[55]where α and β are the rotational axis indices, μ is the inertia tensor, Π is the vibrational angular momentum operator, k is the normal mode index, and V(Q) is the potential in normal coordinates. Q is obtained as an orthogonal transformation of the mass weighted Cartesian displacement coordinates, and is the conjugated momentum operator. The leading term of this Hamiltonian is the harmonic one and can be written in dimensionless normal coordinates as[52,53]with q = γ1/2Q and , where c is the speed of light and p = (ℏω)−1/2P. Given this harmonic leading term, we are allowed to write the potential as a Taylor expansion in N normal coordinates at the equilibrium point V0. If we truncate this expansion at the 4th order, we get the following quartic force field (QFF) form of the potentialwhere ϕ and ϕ are the force constants. ϕ is related to the third order potential derivatives aswhere f is the third order derivative along the three normal modes:

Analogous formulas hold for the potential fourth derivatives ϕ.

Using the n-Mode coupling Representation (n-MR) notation recalled in Appendix A, we can rewrite eq aswhere V0 is the potential energy at the stationary point of the PES.

The advantage of eq is that it allows us to order the terms of the V() expansion as a function of the number of modes coupled in the potential derivatives. Then, we can truncate the expansion depending on the desired number n of coupled terms in the potential to obtain the n-mode coupling representation of the quartic force field (nMR-QFF). In this way, the potential is better represented for a perturbation treatment, taking the harmonic terms as the zero-order leading ones. The anharmonic terms from the potential will cause the total Hamiltonian matrix to have nonzero off-diagonal terms. Finally, using the van Vleck perturbation theory (VVPT),[53] we can formulate the vibrational perturbation theory of the second order plus resonances (VPT2+K) expression of the anharmonic constants reported in eq at the transition state geometrywith k, k′ = 1, 2, ..., N – 1 and χ = −iχ̃. At the equilibrium geometry, the anharmonic constants are calculated as followswhere now k, k′ = 1, 2, ..., N since all frequencies are real. In eqs and 20, Bα is the rotational constant with respect to the α rotational axis, and ζα is the related Coriolis coupling tensor. To simplify the notation, in eqs and 19 we use the same symbol ω for vibrational frequencies and χ anharmonic couplings for both equilibrium geometries and TS ones. In our code implementation, eqs and 19 are factorized to remove possible resonances, as it will be explained in detail below. The VPT2+K formulation of the anharmonic constants requires a 3MR-QFF, and this is the potential approximation that we will use in the present paper.

Implementation

Our Software Workflow

Our implementation consists in a script that interfaces with any ab initio quantum chemistry software. The script exploits a parallel architecture to compute the VPT2+K anharmonic constants for a transition state or a stable molecular geometry. The program workflow is shown in Figure , where the interface with the ab initio electronic structure code is detailed. Specifically, once the user has provided an initial geometry and the desired level of electronic structure theory (input box), our software automatically composes the appropriate input file for the ab initio code which performs the geometry optimization. The script sets strict optimization thresholds and appropriate grid densities in case of DFT calculations in order to have reliable anharmonic constants, as described in Supporting Information Section 1.[52] Then, our script calculates the inertia moment tensor from the optimized geometry to set the molecule in the Eckart frame.[56] This frame is characterized by two constraints. The first one implies that the origin of the system is placed at the system center of mass. The second condition enforces that the total angular momentum in this frame is zero. Our software first calculates the center of masswhere is the vector of the center of mass coordinates, m is the mass of the ith atom, and = (r, r, r) is the position vector of the ith atom. Then, the coordinates of the molecule are centered in , and the inertia tensor of the molecule is calculated according to Appendix B. The moment of inertia tensor is then diagonalized, and the orthonormal eigenvectors are found. In the last step, the coordinates of the molecule are transformed into the frame defined by these eigenvectors. Within the Eckart frame, the Hessian is calculated either by using the coupled ab initio code and its analytical gradients or by finite differences. The Hessian matrix is then diagonalized. We call the corresponding vibrational eigenvectors matrix l and the eigenvalue diagonal matrix .

Figure 1

General workflow of the program. The two main communicating blocks are the script implemented in this work (green) and the ab initio program (yellow). Plain arrows between different blocks highlight the input–output flow between the script and the ab initio code. Dashed arrows track the information flow to compute Coriolis contributions to the anharmonic constants.

In the next step, the Coriolis couplings are calculated. The Coriolis terms ζα are matrix elements related to the interactions between the rotational and the vibrational motion. These are originated by the contribution arising from the rovibrational coupling as in the kinetic energy of the system[57,58]where φ are the angular velocities of the rotating system of axes andwhere = (Q1, Q, ..., Q) is the transpose vibrational normal coordinate vector and its time derivative vector. ζα are the Coriolis coupling matriceswhere α, σ, and δ are the rotational axes, the indices k and k′ represent the kth and k′th normal mode, and l is the k(k′)m-component of the Hessian eigenvector matrix l defined above. It is not necessary to include in eq the translational and rotational contributions since the system of reference is the Eckart frame.

At this stage, the script is ready for the generation of the displaced geometries to compute by finite differences the third and fourth order derivatives of the potential. This last process is the parallel core of our implementation, and it is detailed in the next paragraph. Once the potential derivatives have been calculated, they are combined with the Coriolis coupling matrices to compute the anharmonic couplings according to eq or 20.

Parallelization of the Derivative Calculations

The specific parallel algorithm used to compute the potential derivatives is shown in Figure . The main idea is to perform parallel single point energy (SPE) calculations after generating a set of displaced geometries according to the finite difference scheme proposed by Yagi et al.[59] for the computation of third and fourth order derivatives of the potential. Specifically, we derived the finite difference formulas from the Fornberg schemes.[60] Mixed derivatives of the potential have been approximated up to the second order of accuracy, while the diagonal terms have been calculated at the fourth order of accuracy. This choice was originally proposed by Allouche et al.[46], and it is convenient because it avoids any additional PES points for the fourth order approximation of the direct (not mixed) derivatives.

Figure 2

Algorithm implemented for the calculation of the third f and fourth f order derivatives of the potential. N is the number of vibrational degrees of freedom for the molecule, D is the displacement matrix, l is the columnwise vibrational eigenvector matrix, δ is a fixed incremental step, M is the diagonal matrix containing the atomic masses, and is the diagonal vibrational eigenvalues matrix. G is the displaced geometry vector, and G0 is the equilibrium geometry vector. After the generation of the geometries, the SPE calculations are launched. Npoints is the total number of energy points required, C is the number of cores available per node, and V is the vector that contains all the energies retrieved from the SPE calculations.

The Fornberg method for the generation of finite difference formulas on spaced grids is based on a simple recursion formula on 1-dimensional grids to determine the weights of the potential points in any derivative order formula and up to any order of accuracy. For our purpose we used an equally spaced grid along each normal mode. For multidimensional derivatives we combined the 1D formulas to get the corresponding multidimensional expressions. In this way, we obtained a finite difference expression for the derivatives f, f, f, f, f, and f reported in Appendix C, which are necessary for anharmonic constants calculation according to eqs and 20.

In a 3MR-QFF formula, the overall number of SPE calculations at different geometries is equal toTo compute the coordinates of each configuration for the derivative calculations, the displacement matrix D is generatedwhere M is the diagonal matrix containing the atomic masses and δ is a fixed displacement. We fixed this displacement to 0.5 to get reliable results according to our tests and the literature.[59] After identifying the equilibrium geometry with a 3N dimensional vector of Cartesian coordinates G0, each displaced geometry is generated aswhere G is the ith geometry vector component corresponding to the application of the ith displacement vector [D] of eq to the equilibrium geometry vector G0. Once all geometries necessary for the finite difference derivative calculation are generated, the SPE calculations are launched in parallel. The total number Npoints of independent SPE ab initio inputs, automatically generated by the program, is divided into launching files, where C is the number of computing cores available per node. The launching scripts are generated such that each node, composed by C physical cores, runs in a parallel fashion C SPE calculations. The energies are then retrieved and saved in a vector from which the third and fourth order derivatives of the potential will be computed using the schemes explained in Appendix C.

Resonance Treatment

The anharmonic constants are calculated using the GVPT2+K theory by implementing eqs and 19. Specifically, in this step it is of fundamental importance to properly treat divergent terms that may arise from the zeroing of the denominators occurring at the resonances. The anharmonic constants can be affected by two different kinds of 1–2 resonances. In Fermi type I resonances between mode i and mode j we havewhile in Fermi type II resonances an additional mode k is involvedThe usual way to deal with these resonances is to set a threshold and see if the difference in wavenumbers between the frequencies is smaller than this threshold. In this case, the resonant terms are set to zero in eqs and 20. Unfortunately, this approach is not very accurate because it completely disregards the entire term, i.e., both numerator and denominator, even if it is only the denominator to be singular. Therefore, we decided to adopt the Martin et al.[61] approach that is based on the evaluation of two parameters. For a Fermi type I resonance the parameterand for a Fermi type II resonanceIf the Δ parameter is greater than a lower threshold (usually set at 1 cm–1) and the frequency difference between the modes involved in the resonance is smaller than a higher threshold (usually set at 200 cm–1), the resonant term is disregarded.

Parallel Implementation Scaling Benchmark

In this paragraph we will show how our program scales with respect to the number of processors. We compare the performance of our script interfaced with Gaussian 16 SPE calculations against Gaussian 16 internal routine for the anharmonic constant calculations.[62] The scaling is rationalized in terms of speedup (S) and efficiency (E) parameters, which are respectively defined aswhere T is the serial execution time and T is the parallel one on P processors. The best speedup that one can get is linear with respect to the number of processors, and the efficiency is 1. However, this is hardly achieved in practice because of the core communication and disk writing/reading latency times present in parallel architectures. In our case the latency time is highly reduced, due to the parallelization strategy in which cores do not communicate with one another. Figure reports our scalability tests for the 10-atom system, using either the MP2 or the DFT method, with aug-cc-pVDZ, jun-cc-pVDZ, and 6-31G* basis sets. Information about the specific architecture used for our benchmarks and further tests on the 14-atom system can be found in Supporting Information Section 2.

Figure 3

(a) Speedup and (b) efficiency plots for the anharmonic constants calculation using the Gaussian code and our program. Calculations are done for the 10 atom cyclobutene molecule at the MP2/aug-cc-pVDZ, MP2/jun-cc-pVDZ, and B3LYP/6-31G* levels of theory.

In Figure , the scaling S is linear for all kinds of computational setups when employing our algorithm. Indeed, the time spent by our script in the nonparallel part is negligible. The same calculation on the same architecture but using the Gaussian 16 software scales almost linearly only for a few cores. For a higher number of cores, the Gaussian 16 software parallel performance deviates form the ideal scaling. The same considerations are valid for the E profile, where the efficiency of the Gaussian 16 software drops with the increasing number of allocated CPUs. These considerations do not depend on the basis set employed. As far as the level of theory is concerned, DFT scales slightly better than MP2 when employing the Gaussian 16 software, especially for high dimensional molecules, but still far from linearity. For example, with the GALILEO IBM NeXtScale architecture by the Italian CINECA HPC center,[63] the best performance for the anharmonicity constant computation of the cis-1,3,5-hexatriene at the MP2/aug-cc-pVDZ level of theory using the Gaussian 16 program was achieved over 12 cores, taking 32 h and 23 min. The same calculation on the same architecture was performed in 19 h with the best possible setup using our program. Furthermore, we point out that the calculation of the cyclobutene anharmonic constants at the CCSD(T)/aug-cc-pVDZ level of theory are not doable with Gaussian 16. Thus, we could not compare our code performance with Gaussian 16 in this case. For this CCSD(T) calculation the maximum number of parallel cores employed is actually limited by the RAM size available on each node. However, we experienced that a 252 core parallelization takes only 20 h to complete the calculation with our code.

Therefore, large molecule anharmonicity calculations are more conveniently addressed with our algorithm rather than with Gaussian 16, when a large amount of cores is available.

Results and Discussion

In this section, we will calculate the rate constants using the anharmonicity matrices computed with our algorithm as input for the Parsctst and Paradensum codes[41,42] of the Multiwell program suite.[38] We decided to apply our algorithm to the study of the kinetic rate constants for the three organic reactions shown in Figure . The reactions involve 10, 14, and 16 atoms. They are the cyclobutene ring opening reaction (R1), the cis-1,3,5-hexatriene electrocyclic ring closing reaction (R2), and the [1,5]-Cope rearrangement of the 1,5-hexadiene molecule (R3). These reactions are known as examples of heavy atom tunneling (HAT) processes. HAT includes all the tunneling contributions to the reaction mechanism involving atoms that do not belong to the first period of the periodic table. These kinds of reactions have attracted a growing interest in recent years, both from the theoretical and the experimental community.[64−67] Specifically, a recent paper by Greer et al.[68] reports theoretical calculations of the tunneling contribution to the rate constant for these three reactions.

Figure 4

Organic reactions simulated in this work: The cyclobutene ring opening reaction (R1), the cis-1,3,5-hexatriene electrocyclic ring closing reaction (R2), and the [1,5]-Cope rearrangement of the 1,5-hexadiene molecule (R3).

We start by validating our method. We compare the anharmonic constants calculated with our algorithm with those obtained using the internal Gaussian 16 subroutine. The aim is to check if possible differences may significantly impact the SCTST rate values. We define the difference and the percentage difference between the anharmonic constants obtained with our algorithm χ and the ones with the Guassian16 software χ aswhere k and k′ are the normal mode indices.

Figure presents our results at the MP2/aug-cc-pVDZ level of theory. The axes in Figure refer to the indices k and k′ of the normal modes, and the color gradient indicates the magnitude according to the label. Other B3LYP/6-31G* and MP2/aug-cc-pVDZ calculations have been carried out for the R1, R2, and R3 reactants and TS geometries, as reported in Supporting Information Section 3, and they show comparable deviations. Few values show more than 100% deviation from the Gaussian 16 ones. These values are all related to small anharmonic couplings. Anyhow, when the rate constant is calculated, we find no significant differences in the values between the two approaches, as shown in Figure where the rates of the R3 reaction are reported. The same comparison of Figure but with reactions R1 and R2 is presented in Supporting Information Section 3, and we find the same degree of numerical agreement between the two approaches, allowing us to conclude that the percentage deviations between the anharmonic coupling values obtained with Gaussian 16 and our implementation are irrelevant for the SCTST reaction rate results.

Figure 5

(a) Difference Δχ between anharmonic constants calculated with Gaussian 16 and our program for the hexatriene molecule. (b) Percentage difference % between the anharmonic constants. In both panels, k and k′ normal mode indices are reported on the axes, and the color gradient indicates the magnitude of the difference. Anharmonic constants are calculated at the MP2/aug-cc-pVDZ level of ab initio theory.

Figure 6

SCTST reaction rate constants k(T) at different temperatures for the [1,5]-Cope rearrangement (R3). The y axis reports the natural logarithm of the rate constants and the abscissa the scaled inverse temperature in Kelvin. The anharmonic constants matrix was calculated using both our program and the internal Gaussian 16 subroutine starting from the same geometry and ab initio calculation setup. DFT calculations have been carried out using the B3LYP/6-31G* setup, while MP2 calculations use the aug-cc-pVDZ basis set.

Reaction Rate Constant Calculations

The reactions of Figure have been previously studied by Greer et al.[68] In their work, they employed as their top-notch method the small curvature tunneling (SCT) approach.[69] In their work they also considered the application of monodimensional tunneling corrections, such as Wigner’s and Bell’s corrections, and applied them to the results obtained with the canonical variational theory (CVT) at the B3LYP/6-31G* level of theory. A good agreement between all the methods was found for temperatures above 250 K, and they suggested that there is a significant tunneling contribution for temperatures up to 420 K. They estimated the HAT contribution to be ≥25%. However, a more precise ab initio level of accuracy is necessary before drawing any conclusion about the presence of HAT at room temperatures.

In the present work, we start from the Greer et al.[68] transition state (TS) geometries, and we optimize them at the level of theory of our anharmonic constants calculation. For each reaction, we double-checked each TS geometry with an intrinsic reaction coordinate (IRC) calculation starting from our TS, obtaining the reactant and product geometries. Geometry details can be found in Supporting Information Section 1.

The TST reaction rate constants are calculated using the following formulawhere ΔV0 is the difference in energy between the TS and the reactant with the inclusion of the harmonic zero point energy (ZPE) and the rotational Q(T) and Q(T) partition functions are calculated using eq , while the vibrational partition functions are calculated with the harmonic approximationSince we deal with unimolecular reactions, N is the vibrational normal mode number of both the reactant and the TS. For the TS, we consider that the Nth mode is the imaginary frequency one. In the first formula of eq , ω is the vibrational frequency of the kth mode of the reactant, while in the second one it is the TS frequency. In eq , we do not indicate the translational partition function ratio because the reactions here considered are unimolecular. To estimate the tunneling contribution we define the following general parameterwhere the transmission coefficient κ is the ratio between the SCTST rate constant computed using eq , which includes tunneling and anharmonic partition function contributions, and the TST one from eq , which does not include tunneling contributions and which is strictly harmonicBefore getting into the details of the rate values, we looked for an experimental validation of our computational setup. Specifically, we compare with the experimental results of the reaction barrier estimate for the R1 reaction obtained at different pressures and in a small temperature range.[70] The comparison is done with the caveat that the experimental values are obtained by enforcing the Arrhenius relation k(T) ∝ exp (−E/kT) to the experimental kinetic constants, where E is the empirical activation energy barrier. Instead, our barrier estimates ΔV0 is the energy difference between the TS and the reactants, both corrected for the harmonic ZPE values.

The values of Table show that the computational results at different levels of theory are not within the experimental error bar interval of confidence. Nevertheless, all experimental and computational values are similar and show an overall good agreement, which validates our ab initio setup.

Table 1

ZPE Corrected Forward Reaction Barriers for the Cyclobutene Ring Opening Reactiona

level of theory/basis set	Ea [kcal/mol]
B3LYP/6-31G*	33.9
MP2/aug-cc-pVDZ	31.2
CCSD(T)/aug-cc-pVDZ	31.7
experiment (8–14 Torr)	32.5 ± 0.5
experiment (100 Torr)	32.5 ± 0.4
experiment (1500 Torr)	32.9 ± 0.7
experiment (5 Torr)	32.7 ± 0.2

The experimental values[70] have been calculated at different pressure conditions and in the temperature range 403–448 K.

The values of the forward reaction barrier used to compute the SCTST rate constants include the anharmonic ZPE correction and the term. This term, although small, can be important in reaction rate calculations, and in our cases it has been determined to be relevant, expecially when evaluating the percentage difference with respect to the TST value. The values of the forward barriers corrected for both anharmonic ZPE and can be found in Supporting Information Section 4.

Before presenting our rate calculation results, we specify how we deal with hindered rotations (HRs). An accurate HR treatment in the partition functions is of fundamental importance for the calculation of reaction rate constants. We identify the HR degrees of freedom using the Ayala schemes[71] implemented in the Gaussian 16 software. After the identification and before the VPT2 calculation, the HR modes are projected out from the Hessian matrix, and their one-dimensional partition functions are calculated. For the 1D HR partition functions, we used the Pitzer and Gewin correction to the quantum harmonic oscillator partition function.[72] With this correction, we can smoothly tune the HR mode treatment from the quantum harmonic oscillator to the free rotor partition function, as the temperature is increased.

The results for the cyclobutene ring opening reaction R1 of Figure are shown in Figure , where panel (a) reports the rate values. An evident bias in the rate estimate is represented by the ab initio theory level, even if post-HF methods give comparable results. These differences should be ascribed mainly to the barrier height, since the exponential function of the rate expression greatly magnifies this dependency. However, when choosing the DFT level of theory, our SCTST results obtained using eq are in close agreement with previous literature results, as highlighted in the inset of panel (a) in Figure . In the same panel, a significant deviation from the Arrhenius behavior is observed below T = 150 K. We observe also that a significant role in the rate estimate is played by the partition function and by the shape of the barrier. Specifically, in Table the MP2 barrier is lower than the CCSD(T) one, while in Figure the CCSD(T) rate is greater than the MP2 one in the tunneling regime. This is due to a thicker MP2 barrier with respect to the CCSD(T) one (see Supporting Information Section 4). The deviation from Arrhenius linearity is more evident from panel (b) of Figure , where the percentage difference with respect to the TST rate is estimated using eq . Even at T = 500 K an amount of 20% rate enhancement is observed, and we think this is mainly due to the presence of tunneling.

Figure 7

Cyclobutene ring opening reaction (R1). (a) TST, SCTST, and SCT(CVT)[68] reaction rate constants. The y axis reports the natural logarithm of the rate constants and the abscissa the scaled inverse temperature in Kelvin. (b) TST, SCTST, and SCT percentage difference between the semiclassical and the TST kinetic rate constant according to eq . For the SCTST rate constant calculations, the anharmonic constant matrix was calculated using our program. DFT calculations have been carried out using B3LYP/6-31G*, while in the MP2 and CCSD(T) calculations we used the aug-cc-pVDZ basis set.

The results for the cis-1,3,5-hexatriene electrocyclic ring closing reaction, R2, are shown in Figure . In both SCTST and TST calculations, we accounted for two HRs in the reactant molecule. The corresponding frequencies and the Pitzer and Gewin parameters employed can be found in Supporting Information Section 5. Figure shows quite different kinetic rate constants depending on the level of ab intitio theory, which goes from the B3LYP/6-31G* to the post-HF MP2/aug-cc-pVDZ and CCSD/jun-cc-pVDZ level of theory. The values of the energy barriers and the TS frequencies can be found in Supporting Information Section 4. All the rate constant calculations show a certain amount of tunneling and a significant deviation from the Arrhenius linear behavior. In particular, the SCTST rates calculated with the DFT potential once again agree quite well with the SCT(CVT) estimates, as shown in the inset of panel (a) in Figure . In panel (b), the percentage difference of each semiclassical rate calculation with respect to the TST one is reported. At low temperatures, we observe for all the semiclassical calculations an almost 100% difference, and we attribute this to tunneling contributions. In the high temperature regime, we think that the anharmonicity of the partition function plays an important role. Specifically, the MP2 barrier is excessively flat and anharmonic if compared to the CCSD and DFT ones, as can be deduced from the higher MP2 rates at any temperature. This flatness generates a significant difference in the TS partition function estimates when using the anharmonic expression of eq versus the harmonic one of eq , and it is responsible for the negative percentage differences observed in panel (b) of Figure . The higher the temperature, the greater the effect of the potential anharmonicity. A higher level of calculation, such as CCSD, confirms that the large negative percentage difference of the MP2 calculations is due to the limitation of the MP2 level of electronic structure theory.

Figure 8

Same as in Figure but for the cis-1,3,5-hexatriene electrocyclic ring closing reaction (R2).

The results for the [1,5]-Cope rearrangement reaction of 1,5-hexadiene, the R3 reaction in Figure , are shown in Figure . Three hindered rotations were considered for the reactant molecule, corresponding to the three normal modes with the lowest frequencies. The results show a strong dependence of the rate constant on the ab initio method used, such as for the R1 and R2 reactions. Specifically, Houk and co-workers[65] showed that extra care should be taken when dealing with pericyclic reaction ab initio calculations, since correlation contributions to the total energy are large and important. To summarize, they concluded after several benchmark calculations of the forward reaction barriers that the MP2 method probably overestimates the correlation energies and underestimates the barrier height. Indeed, we find in our calculations reported in Figure a that the post-HF method CCSD/jun-cc-pVDZ and MP2/aug-cc-pVDZ rates differ by several orders of magnitude. More specifically, the shape of the MP2 potential is so anomalous that the effective imaginary frequecy Ω of eq becomes negative. In this case, the Multiwell program[73] performs the standard TST calculation; thus, we obtain the Arrhenius plot with blue diamonds reported in Figure a. As already observed in the other reactions, the R3 rates shown in the inset of Figure a confirm the previous SCT(CVT) results.[68] Considering the percentage difference between the semiclassical and TST results (see panel (b) of Figure ), we think that this difference is mainly due to tunneling, especially at low temperatures. In this case, the partition function anharmonicity is not as important as in the case of R2, since negative percentage differences are not observed. At T < 150 K, tunneling becomes the only possible way to get to the product side. This is evident also from panel (a) of Figure , where the CCSD rate reaches a plateau. In addition, panel (b) shows how CCSD results predict a larger amount of tunneling than the DFT ones, and that in this case there is a slight difference between the amount of tunneling introduced by the SCT(CVT) approach and our SCTST/DFT one.

Figure 9

[1,5]-Cope rearrangement reaction of the 1,5-hexadiene molecule (R3). (a) TST, SCTST, and SCT(CVT) reaction rate constants. (b) TST, SCTST, and SCT percentage difference with respect to the TST of the kinetic rate constant. For the SCTST rate constant calculations, the anharmonic constant matrix was calculated using our program. DFT calculations have been carried out using B3LYP/6-31G*, MP2 calculations using the aug-cc-pVDZ basis set, and CCSD calculations using the jun-cc-pVDZ basis set.

Summary and Conclusions

In this paper, we calculate at the CCSD(T) and CCSD level of accuracy tunneling rates involving heavy atoms in organic reactions thanks to our fast and efficient algorithm for the computation of anharmonic constants. Our program can be easily interfaced with any ab initio code available, since only SPE calculations are requested. We use the anharmonic constants as input for the semiclassical transition state theory (SCTST) programs Paradensum and Parsctst[41,42] from the Multiwell program suite.[38]

First, we tested the performances of our algorithm, showing that it successfully overcomes the limitations of other software. Then, we have tested our program for the SCTST calculation of reaction rate constants on organic reactions with relevant HAT tunneling contribution even at room temperature, as anticipated by Greer et al.[68] with DFT level of theory estimates. Our code allowed us to extend the HAT phenomena investigation at the post-HF level of theory, up to the CCSD(T) gold standard with consistent basis sets of the type of aug-cc-pVDZ. By varying from the DFT to the MP2 and CCSD(T) or CCSD approaches, we notice conspicuous variations in the reaction rate constants, the tunneling and partition function percentage contributions showing how important the possibility to employ a high level of ab initio theory is. At the B3LYP/6-31G* level, the SCTST approach and the SCT gave similar results for the R1 and R2 reactions, while a slight difference is found for the R3 reaction. We conclude that the use of accurate ab initio post-HF methods is necessary in order to get a reliable PES and anharmonic constants, so that the best possible SCTST reaction rate constant can be calculated.

As far the scaling performances of the code are concerned, we note that the Gaussian 16 software does not actually calculate all the third and fourth order derivatives of the potential but cleverly exploits the symmetry of the molecule to spot zero anharmonicity values and avoid their numerical calculation. We did not implement this symmetry driven strategy in our approach. Nevertheless, the performance of our program overcomes the Gaussian 16 one given the same computational power. A future inclusion of a symmetry tool in our algorithm may be useful to further reduce the computational time, especially for calculations on organic molecules of biological interest.

Gaussian 16[52] employs the Thiel et al.[74] scheme where the gradients are calculated analytically and then the Hessians are differentiated numerically to obtain the needed third and fourth order potential derivatives. This guarantees a higher precision with respect to our approach where the finite differences of the energy are employed. However, as shown above, the scaling and the parallelization efficiency of Gaussian 16 is worse than in our implementation, and the loss of accuracy in estimating the anharmonic couplings with our procedure leads to negligible differences in the rate constant values.

We recall here that couple cluster parallel calculations of anharmonic constants can also be performed with CFOUR ab initio software[75] by setting up an appropriate procedure as explained in the CFOUR user manual.[76] The algorithm presented in our work is a possible alternative to that implementation. Depending on the computational efficiency in parallel architectures of these calculations, one should choose which approach to adopt by estimating the trade-off between computational time and accuracy.

We think the strengths of our code are the high level of automation and that it can be interfaced with any quantum chemistry packages, allowing the user to choose among diverse ab initio methods. In any case, our code is readily modifiable. For instance, it is possible to change the schemes adopted to calculate the anharmonic constants. In the present paper we used a 3MR-QFF for the calculation of the anharmonic constants, but it is possible to further reduce the number of coupled modes in the QFF representation of the potential and find new formulas for the anharmonic constants. For example, the 2MR representation opens the possibility to an alternative approximation of the SCTST rate constants, different from the reduced dimensionality approaches present in the literature.[77] We will further study the implications of the 2MR approximation in the SCTST treatment.

Our algorithm and the related scripts will be made available as part of the Multiwell suite of codes.[38]