A simplistic computational procedure for tunneling splittings caused by proton transfer

In this manuscript, we present an approach for computing tunneling splittings for large amplitude motions. The core of the approach is a solution of an effective one-dimensional Schrödinger equation with an effective mass and an effective potential energy surface composed of electronic and harmonic zero-point vibrational energies of small amplitude motions in the molecule. The method has been shown to work in cases of three model motions: nitrogen inversion in ammonia, single proton transfer in malonaldehyde, and double proton transfer in the formic acid dimer. In the current work, we also investigate the performance of different DFT and post-Hartree–Fock methods for prediction of the proton transfer tunneling splittings, quality of the effective Schrödinger equation parameters upon the isotopic substitution, and possibility of a complete basis set (CBS) extrapolation for the resulting tunneling splittings.


Introduction
Large amplitude motions (LAM) are ubiquitous in most molecular systems, they are responsible for conformational interconversion and even for chemical reactions, such as in the cases of tautomerization [1][2][3][4][5]. One of the most prominent LAM manifestation is the tunneling splitting of the ground vibrational state, whenever a motion between equivalent minima is involved. Such splittings can be observed, for instance, using microwave (MW) or millimeter wave (MMW) rotational spectroscopy [6][7][8][9][10][11][12][13][14][15]. A knowledge of the LAM and of its characteristics can be helpful for the interpretation of the experimental data [11,14,16,17].
Both these systems have been thoroughly investigated using various theoretical approaches. The low dimensional treatment of proton transfer and some of the coupled motion have been shown to give reasonable splitting estimates for MA [22], while for FAD such applications were less successful [23,24]. Applications of the reaction surface Hamiltonian (RSH) models [25][26][27][28] and of anharmonic normal mode-based treatments [18,[29][30][31][32][33] have yielded in good results for both MA and FAD. As expected, the full-dimensional quantum Monte-Carlo simulations provide results of the experimental quality for MA [29,34,35]. For MA and FAD, instanton approaches have been shown to produce accurate results [36][37][38][39][40].
In this paper, we present a modified version of another computational approach provided in Refs. [17,41,42]. It assumes adiabatic separation of the LAM from the other vibrational modes, that we will further refer to as smallamplitude motions (SAM). The latter are to be treated in the harmonic approximation. Our procedure requires two types of calculations to be done: one-dimensional (1D) relaxed potential energy surface (PES) scans and vibrational frequency calculations in the harmonic approximation. Both computational procedures are known by users of quantumchemical packages, and thus the presented procedure can become an interesting one for members of the highresolution molecular spectroscopy community. We apply the proposed approach to the nitrogen inversion in ammonia to prove the workability of the method, and then we use it to calculate tunneling splittings for the proton transfer motion in MA and FAD at various quantum-chemical approximations employing the def2-TZVPP basis set.

Theoretical framework
The treatment of the LAM that we implemented resembles the procedure introduced in Refs. [41,42]. It is a mixture of the RSH approach [1,28,[67][68][69][70], with the kinetic energy expression taken from the Meyer-Günthard Hamiltonian [71]. The Hamiltonian is expressed as where is the LAM coordinate, p = −i� ∕ = −i� is the momentum conjugated to , V is the electronic potential energy surface (PES) along the LAM coordinate (obtained from a relaxed PES scan), and Ĥ b is the bath of all other 3N − 7 SAM degrees of freedom (N denotes the number of atoms in the molecule). G = −1 is the inverted reduced mass ( ) of the LAM. It is being calculated from the generalized tensor of inertia defined as is the vector of LAM-rotation interaction, and Here, m and denote the atomic masses and positions. 1 G is the last diagonal element of matrix G −1 .
Let = (x 1 , y 1 , z 1 , x 2 …) † be the 3N-sized vector with Cartesian coordinates of atoms, and = (m 1 , m 1 , m 1 , m 2 , …) † 3N-sized vector of atomic masses. Let = √ be massweighted coordinates (MWC) of the system, 2 and F( ) -the Wilson's GF-matrix [72]  where ̂ MWC = −i� . We can remove translations/rotations of the molecule as a whole, as well as the LAM, from F( ) using a projector matrix S( ) , which can be constructed e.g. using the Gram-Schmidt process. As a result, one can get the harmonic bath GF-matrix F 3N−7 ( ) = S † ( )F( )S( ) . We chose one of the structures (at = ref ) as the reference. In our double-well cases, we chose the transition state with ref = 0 . Then we define the vibrational state of the bath using the normal modes of the reference struc- , with 2 ref being a vector of squared angular frequencies) [72,73]. This allows us to write ( − eq ( )) = L ref , where is the position of the nuclei along the normal modes expressed in mass-weighted coordinates. Note that this assumption requires all the structures to be oriented with respect to the reference coordinates. These operations give the harmonic bath Hamiltonian in the form where k enumerates normal modes, P k = −i� Q k , Ω 2 kl are the elements of the matrix 2 ( ) = L † ref F 3N−7 ( )L ref , and k are the zeroth approximation frequencies. We can either set them at ref values, or allow bath vibrations to relax during the LAM by setting k to the respective diagonal values of 2 ( ) . Ĥ 0 denotes the unperturbed vibrational Hamiltonian, and the Ŵ k are the perturbation operators. The mass-weighted P k /Q k (normal modes momentum/ position operators) are related to dimension-less operators p k ,q k corresponding to an oscillator with reference frequency k as We can use them to rewrite the Eq. 3 as (annihilation/ creation operators). 3 If we consider the bath to be in a chosen vibrational , then we can apply perturbation theory (PT) to get the expression of the bath energy. Using the second quantization form of p k and q k , one can compute the bath Hamiltonian energies.
-0 th order.  4 Using the fact that q 2 = (â † ) 2 +â 2 2 +n + 1 2 . For the pure LAM excitations we can assume that the other molecular degrees of freedom are at the ground state ( = 0 ), and therefore, the energy of the bath is the zeropoint vibrational energy (ZPVE) given by an expression If we allow vibrational states relaxation, then k = Ω kk , which gives By combining Eqs. 1 and 6 the final effective Hamiltonian for LAM looks as the following: The parameters for this Hamiltonian (effective masses G( ) and effective PES V eff ( ) ) in this work were computed from relaxed PES scans obtained with Orca 4 using in-house Python scripts [74]. In brief, the script for computing the Hamiltonian parameters works as following.
-The code reads Orca output *.hess files containing vibrational information (molecular geometry, atomic masses, Hessian, etc.). -The molecular geometries and Hessians are oriented with respect to the reference geometry defined as the molecular structure with the lowest value of | | (the closest to the transition state). The orientation corresponds to the minimum of the functional -The effective masses = G −1 are computed as described above. The derivatives quencies are also being set to zero. The resulting Hessians are stored. -The ZPVE correction is computed for each structure by projecting the reference structure normal modes onto the Hessians retrieved in the previous item. -The obtained dependencies for ( ) , V(x) and ZPVE( ) , are extrapolated using cubic splines to a fine uniform grid.
The resulting eigenvalue problem is solved using the discrete variable representation (DVR) method, in particular, sinc-DVR [75].

Proof of concept: ammonia
To test the performance of the implementation, we chose the simplest case of LAM with tunneling: an inversion motion in ammonia (NH 3 ). The LAM coordinate was defined as where NHn denotes a position vector of the hydrogen number n with nitrogen as an origin, r nm = | NHn − NHm | , and s = (r 12 + r 13 + r 23 )∕2 . The NH 3 is basically a height (h) of the tetrahedron constructed on the NHn vectors with an origin at the nitrogen atom, computed through the volume of the tetrahedron V = A ⋅ h∕3 and the area of the base of the tetrahedron (A). The 1D PES scans were performed at B3LYP-D3BJ, CAM-B3LYP, MP2, and B2PLYP levels of theory. The resulting parameters of the LAM Hamiltonian (Eq. 7) are given in Fig. 2. All the methods performed similar to each other. The energy levels of ammonia and some of its isotopologues obtained with CAM-B3LYP are given in Table 1. Concerning simplicity of the calculation, the results are very similar to the experimental values. The values of other methods can be found in the SI. The results for ammonia confirm the workability of the code, which justify testing of the approach on larger molecules.

Tunneling splittings in the parent species
The PESs for proton transfer in MA and FAD were computed at various levels of theory, using both post-Hartree Fock methods (MP2 and SCS-MP2) and different DFT functionals from different steps of Perdew's ladder [76]: GGA (PBE and BLYP), hybrid (B3LYP and PBE0), range-separated hybrid (CAM-B3LYP), and double hybrid (B2PLYP and mPW-2PLYP). FAD is a dimer; therefore, we took into account possible missing dispersion interactions in case of GGA and hybrid functionals via D3BJ correction. The same correction was also tested together with one of the double-hybrid functionals (B2PLYP).
The proton transfer coordinates were as following. In case of MA: and for FAD as The obtained Hamiltonian parameters (PES and effective masses) for these calculations are plotted in Fig. 3, the computed tunneling splittings between ground vibrational states of different parity ( 0 + ∕0 − ) are given in Table 2. As it can be seen from these results, most of the tested methods give reasonable estimates of the tunneling splittings compared to the experimental values determined in Refs. [6,9,15]. The pure GGA functionals with dispersion corrections (PBE-D3BJ and BLYP-D3BJ), as well as the hybrid PBE0-D3BJ, perform not well due to underestimation of the proton transfer barriers. Double hybrids (especially B2PLYP without D3BJ), MP2, and CAM-B3LYP work stable for both MA and FAD. B3LYP-D3BJ performs well for MA, but in case of FAD, it underestimates the barrier height, which might be a result of the intramolecular nature of the bonding forces in this complex. SCS-MP2 showed the worst performance.    The most probable explanation to that is breaking of the adiabaticity assumed in the Hamiltonian given in Eq. 2. Isotopic substitution that leads to an increase of the atomic masses with respect to the parent isotopologue (deuteration, 13 C and 15 N substitution) reduces the vibrational frequencies of the SAM in the bath. This can make these SAM susceptible to vibrational excitation by LAM because of the reduced energy difference between effective potential surfaces. MA has lowenergy vibrational modes, which might lead to this adiabaticity violation. This is supported by the plots of the density of the LAM vibrational manifolds given as V ( ) = V( ) + E ( ) with vibrational energy of the SAM defined by Eq. 5 (see Fig. 4). In case of NH 3 , ND 3 , and the parent MA isotopologue, the effective potential energy V 0 ( ) = V( ) + ZPVE( ) is well separated from the excited vibrational manifolds. However, in the case of MA-D6, the other manifolds get close to the V ( ) at the positions , which are populated in the ground vibrational Table 3 Tunneling splittings for proton transfer in MA and some of its isotopologues. See Fig. 1 and Refs. [7,8]   isotopologue (see Fig. 1 for enumeration). Reaction coordinates for ammonia and MA are given by Eqs. 8 and 9, respectively 1 3 state. This means that multiple vibrational manifolds might be needed for consideration in order to get a good description of the LAM ground vibrational state of the MA-D6.

Complete basis set (CBS) extrapolation of the tunneling splittings
Karlsruhe basis set series [66] can be used for the complete basis set (CBS) extrapolation [77,78], therefore a general possibility of an extrapolation to the complete basis set (CBS) limit for the tunneling splitting values was investigated. For this purpose, we have computed the 0 + ∕0 − tunneling splittings in MA and FAD using def2-nVP and def2-nVPD basis set series (n = S, TZ, QZ, which corresponds to cardinal number = 2, 3, 4 , respectively) with the B3LYP-D3BJ functional. The tunneling splittings were then approximated by an exponential fitting formula ΔE 0 + ∕0 − ( ) = ΔE CBS 0 + ∕0 − + A ⋅ exp(− ⋅ ) , where ΔE 0 + ∕0 − denote the value of the splitting obtained at different basis sets and the CBS-extrapolated value. Due to smaller values of the tunneling splitting of FAD, which give rise to larger numerical instabilities, the parameters (1.5 for def2-nVP and 1.2 for def2-nVPD) for the CBS extrapolation for FAD were taken from the MA fits. The results are shown in Fig. 5. Both def2-nVP and def2-nVPD basis set series converge to similar values of ΔE CBS 0 + ∕0 − . The used basis series consist only of three basis sets, therefore we could not make a proper estimation of the CBS extrapolation errors in each serie, however, the similar trends for the tunneling splittings are clearly visible. The usage of the parameter obtained from MA for FAD tunneling splittings indicates a transferability of the trends from one molecular system to another. Thus, we conclude that the CBS extrapolation is also, in principle, possible for such complicated parameters as tunneling splittings.

Reaction coordinate choice invariance
The proposed approach explicitly relies on the reaction coordinate definition: for calculations of both reduced mass of the motion and of the SAM ZPVE correction. Although, by construction, the Schödinger equation solutions produced by any of the properly chosen coordinates (i.e., by those, that can map a transition from one minimum to another through the vicinity of the transition state) should be the same. However, the numerical implementation of the methods might be rather sensitive to the coordinate choice. In order to check this stability of the algorithm, we have performed the same type of calculation for MA at B2PLYP/def2-TZVPP level of theory with ill-defined reactions coordinates (see Fig. 1 for atomic enumeration): and The value computed with the initial reaction coordinate (Eq. 9) was 20.3 cm−1 (see Table 2). By replacing the to CH and HH the result did not change significantly: 19.9 and 20.1 cm−1 , respectively. Thus, we conclude that the implemented method is rather stable with respect to the reaction coordinate definition.

Calculation using the composite approaches
The additional test of performace was done for the com-  indicates a computationally cheap quantum chemical method, that is being used as a source of relaxed PES scan and harmonic vibrational frequencies, and B is the high accuracy level of theory, that provides accurate single point energies. In order to check the possibility of such schemes, we used the relaxed PES scan done at B3LYP-D3BJ/def2-TZVPP level of theory, and computed single point energies at these geometries using the MP2, B2PLYP, and ae-CCSD(T) levels of theory with the def2-TZVPP basis set. The computed tunneling splittings are given in Table 4. As one can see, all composite schemes gave as reasonable estimates for the splittings, as the single method calculations. However, the scheme employing the most accurate quantum chemical approximation of all tested (ae-CCSD(T)/def2-TZVPP) gave the splitting 15 wavenumbers lower than the expected value. Such decrease of the value might be caused by a different local shape of the PES at this approximation with respect to the B3LYP-D3BJ/def2-TZVPP. This discrepancy requires additional investigation, which is out of the scope of the current manuscript. Nevertheless, we can conclude that the proposed method for computation of the tunneling splittings should also be applicable with the composite quantum chemical calculations.

Limitations of the approach
The introduced approach seems to have three main limitation factors for its performance: -usage of the approximate quantum-chemical methods; -a few steps in the computational procedure, namely, calculation of the effective masses ( ), of the ZPVE, and the interpolation, introduce the numerical noise; -anharmonic effects for the bath vibrational modes.
Despite all three effect, the benchmark calculations given above demonstrate, that with such simple computational procedure we can predict an order of magnitude of the tunneling splitting based only on the quantum-chemical calculations with a relatively low computational cost. Therefore the price-quality ratio of the model given here, seems to be good enough, such as this approach can be used in a routine fashion for augmenting the experimental findings.

Conclusion
We have introduced, implemented, and tested an approach for computing the tunneling splittings for one-dimensional large amplitude motions (LAMs). The computational scheme is based on the adiabatic separation of LAM from the small amplitude motions, which are treated as a vibrational bath. The bath zero-point vibrational energy (ZPVE) is added to the electronic potential energy surface (PES) to yield an effective PES. Combined with an effective mass, we get a one-dimensional Schrödinger equation, and the solutions of the latter allow for determination of the tunneling splitting caused by the LAM between equivalent minima.
The described approach has been shown to work in cases of three model LAMs: nitrogen inversion in ammonia, single proton transfer in MA, and double proton transfer in FAD. Most of the tested methods (B3LYP-D3BJ, CAM-B3LYP, B2PLYP, mPW2PLYP, B2PLYP-D3BJ, and MP2) performed reasonably well and thus can be recommended for general usage. However, the pure GGA functionals tested (PBE-D3BJ and BLYP-D3BJ), as well as the PBE0-D3BJ and SCS-MP2 methods give splitting values that are quite far away from the experimental values and thus should be avoided in the calculations of the proton transfer tunneling splittings. The methodology with certain limitations was found to be applicable for isotopic substitution: an applicability condition is that the adiabatic separation of the LAM from SAM should remain in the isotopologues. We have also demonstrated a possibility of the CBS extrapolation for the tunneling splittings, the of the usage of the composite quantum chemical methods, and the invariance of the calculated spliting with respect to the reaction coordinate choice. In conclusion, the implemented approach can be routinely applied for calculations of the tunneling splittings for the molecules and molecular clusters investigated with MW and MMW spectroscopy.