Scaling of Rotational Constants

1. Introduction

Rotational spectroscopy, consisting of microwave (MW),[1,2] millimeter wave (MMW),[3] and terahertz (THz)[4,5,6] spectroscopies, is a powerful high-resolution experimental technique that provides unprecedented structural sensitivity for the different structural (including conformational) and constitutional isomerisms of molecules,[7,8,9] as well as for isotopic substitutions and various large-amplitude vibrational motions,[10] such as proton transfer,[11] internal rotation,[12] inversion,[13] and even movements of an entire molecule across another molecule’s surface.[9] This kind of spectroscopic technique can be used in various applications, from monitoring the chemical composition of mixtures[14] and reactions[3,15] to detecting atmospheric[16,17] and interstellar molecular species.[8,18]

In the zeroth approximation, the rotational spectrum of a molecule is given by the rigid rotor approximation, which in the case of non-linear molecules is parametrized using three rotational constants, which are denoted as A, B, and C or, equivalently, as $B_a$, $B_b$, and $B_c$, respectively.[19,20,21] These constants, usually expressed in MHz, are related to the moments of inertia $I_{\alpha }$ along the given $\alpha$-th principal axis of the molecule through the following expression:[19,20]

$$B_{\alpha }={\displaystyle \frac{\hslash }{4\pi I_{\alpha }}},$$

(1)

where $\alpha =a,b,c$ is the given principal axis, $\hslash =1.05\times {10}^{-34}$ J·s is the reduced Planck constant, and the moment of inertia is given as

$$I_{\alpha }=\sum _{n=1}^N{m}_n·r_{\alpha n}^2,$$

(2)

where n enumerates the atoms in the molecule, N is the overall number of atoms in the molecule, $m_n$ is the mass of the n-th atom, and $r_{\alpha n}$ is the distance of the n-th atom to the $\alpha$-th principal axis. In the experiment, the so-called vibrationally averaged rotational constants are obtained, usually in the ground vibrational state, and they are commonly denoted as $A_0$, $B_0$, and $C_0$.[21]

The initial structural assignment of the spectroscopically observed species is usually done based on quantum-chemical (QC) calculations. The candidate structures are optimized at a chosen level of theory, usually with a dispersion-corrected density functional theory (DFT) calculation, and then the theoretical rotational constants are compared with the experimentally determined ones.[22] However, such a comparison does not always yield an unambiguous assignment of the molecular structures, and the reason is two-fold. First of all, the optimized structure corresponds to the so-called equilibrium geometry with corresponding equilibrium rotational constants $A_{\mathrm{e}}$, $B_{\mathrm{e}}$, and $C_{\mathrm{e}}$, in which all vibrational effects are absent.[21,23] Therefore, the experimental $A_0$, $B_0$, $C_0$ values and theoretical $A_{\mathrm{e}}$, $B_{\mathrm{e}}$, $C_{\mathrm{e}}$ are not equal due to the vibrational anharmonic shifts, which usually expands the molecular size, thus increasing the moments of inertia (Equation 2) and consequently decreasing the corresponding rotational constant (Equation 1).[24] In other words, it is generally expected that $B_{0,\alpha }\le B_{\mathrm{e},\alpha }$. The second reason is the quality of the QC approximation, which can distort the equilibrium structure due to complicated underestimation and/or overestimation of various intra- and intermolecular chemical bonds and non-covalent interactions. This systematic error does not have a preferred shift of the rotational constant values with respect to their experimental counterparts and thus can be of any type.[21]

In this work, we propose to systematically improve the inconsistency between experimental ground-state averaged rotational constants ($A_0$, $B_0$, $C_0$) and their theoretical equilibrium counterparts ($A_{\mathrm{e}}$, $B_{\mathrm{e}}$, $C_{\mathrm{e}}$) by applying tabulated scaling factors. Such an approach, where the band shifts due to anharmonic effects and QC approximation failure, was demonstrated to be fruitful in the case of vibrational (e.g., infrared) spectroscopy.[25,26,27,28,29,30,31,32] Therefore, it is interesting to investigate whether a similar systematic improvement for the lower-frequency spectral range can be achieved as well. First, we will introduce the procedure of the scaling, the fitting model, and the training dataset, then we will provide the scaling factors, and in the end, we will give an application example using a few recently studied systems with the usage of the PBE0-D3(BJ)/def2-TZVP level of theory.

2. Scaling Procedure

We propose to perform the scaling of the theoretical equilibrium constants obtained from QC calculations ($A_{\mathrm{e}}^{\mathrm{QC}}$, $B_{\mathrm{e}}^{\mathrm{QC}}$, $C_{\mathrm{e}}^{\mathrm{QC}}$) with a single global scaling factor s for a given QC approximation, such that the adjusted constants $A_{\mathrm{s}}$, $B_{\mathrm{s}}$, and $C_{\mathrm{s}}$ are defined as

$$\left\{\begin{array}{c}{A}_{\mathrm{s}}=s·A_{\mathrm{e}}^{\mathrm{QC}},\hfill \\ B_{\mathrm{s}}=s·B_{\mathrm{e}}^{\mathrm{QC}},\hfill \\ C_{\mathrm{s}}=s·C_{\mathrm{e}}^{\mathrm{QC}}.\hfill \end{array}\right.$$

(3)

Since the rotational constants are inversely proportional to the moments of inertia (see Equations 1 and 2), such scaling procedure is effectively equivalent to a global scaling of the atomic coordinates as

$${\mathbf{r}}_{\mathrm{s},n}={\displaystyle \frac{{\mathbf{r}}_{\mathrm{e},n}}{\sqrt{s}}},$$

(4)

where ${\mathbf{r}}_{\mathrm{e},n}$ and ${\mathbf{r}}_{\mathrm{s},n}$ are the equilibrium and scaled positions of the n-th atom in the molecule (see Equation 2).

To tabulate the scaling factors s for given QC approximations, we need a benchmark dataset. In this case, we chose the set of molecules used for obtaining the scaling factors for harmonic frequencies from Ref. [33]. From the set of 441 neutral singlet molecules, 174 non-linear molecules up to 17 atoms with experimentally available $A_0$, $B_0$, and $C_0$ values were selected. Only the relatively high-quality QC approximations were considered here, namely DF-Dn/def2-mVP ($\mathrm{DF}=\mathrm{B}3\mathrm{LYP},\mathrm{PBE}0$, $n=3\left(\mathrm{BJ}\right),4$, $m=\mathrm{S},\mathrm{TZ}$),[34,35,36,37,38,39,40] PBEh-3c,[41] and ${\mathrm{r}}^2$SCAN-3c[42] levels of theory. This was done since rotational spectroscopy, as a high-resolution technique, generally requires more high-quality data for comparison than rotationally unresolved infrared spectroscopy. Besides, MW spectroscopy, as the lowest frequency gas-phase spectroscopic technique to this date, is practically limited by the size of the systems that can be brought into the gas phase in sufficient amounts and by the spectral resolution, as for large systems the spectra will become dense and uninterpretable, despite of the experimental resolution being of the order of $\Delta \nu /\ni \sim {10}^{-6}-{10}^{-7}$ in the standard arrangement.[43]

As the values of the rotational constants $A\ge B\ge C$ can differ by orders of magnitude, the metrics that use the absolute deviations between the experimental and theoretical rotational constants are essentially useless, as they will mostly fit the A-rotational constants for the small-size molecules. Therefore, an advantageous approach is switching to relative values fitting, which was introduced in Ref. [33]. In the application to rotational constants, the least-squares problem can be written in the following way:

$${\mathrm{rRSMD}}^2\left(s\right)=\sum _{k=1}^{N_{\mathrm{mol}}}\left({\left(s·{\displaystyle \frac{A_{\mathrm{e},k}^{\mathrm{QC}}}{A_{0,k}^{\exp }}}-1\right)}^2+{\left(s·{\displaystyle \frac{B_{\mathrm{e},k}^{\mathrm{QC}}}{B_{0,k}^{\exp }}}-1\right)}^2+{\left(s·{\displaystyle \frac{C_{\mathrm{e},k}^{\mathrm{QC}}}{C_{0,k}^{\exp }}}-1\right)}^2\right)\to min,$$

(5)

where rRMSD denotes the relative root-mean-square deviation (rRMSD) of the rotational constants, k enumerates the molecules in the dataset, and $N_{\mathrm{mol}}$ is the total number of molecules in the dataset. The optimal scaling factors are thus given via equation[33]

$$s_{\mathrm{opt}}=\underset{s}{argmin}\left({\mathrm{rRSMD}}^2\left(s\right)\right)={\displaystyle \frac{{\displaystyle \sum _{k=1}^{N_{\mathrm{mol}}}}\left(\frac{A_{\mathrm{e},k}^{\mathrm{QC}}}{A_{0,k}^{\exp }}+\frac{B_{\mathrm{e},k}^{\mathrm{QC}}}{B_{0,k}^{\exp }}+\frac{C_{\mathrm{e},k}^{\mathrm{QC}}}{C_{0,k}^{\exp }}\right)}{{\displaystyle \sum _{k=1}^{N_{\mathrm{mol}}}}\left({\left(\frac{A_{\mathrm{e},k}^{\mathrm{QC}}}{A_{0,k}^{\exp }}\right)}^2+{\left(\frac{B_{\mathrm{e},k}^{\mathrm{QC}}}{B_{0,k}^{\exp }}\right)}^2+{\left(\frac{C_{\mathrm{e},k}^{\mathrm{QC}}}{C_{0,k}^{\exp }}\right)}^2\right)}},$$

(6)

where argmin denotes the minimal value of s, which minimizes the corresponding function value. The fitting uncertainty of this value is given by the equation

$${\sigma }_{\mathrm{opt}}={\displaystyle \frac{\mathrm{rRSMD}\left(s_{\mathrm{opt}}\right)}{\sqrt{{\displaystyle \sum _{k=1}^{N_{\mathrm{mol}}}}\left({\left(\frac{A_{\mathrm{e},k}^{\mathrm{QC}}}{A_{0,k}^{\exp }}\right)}^2+{\left(\frac{B_{\mathrm{e},k}^{\mathrm{QC}}}{B_{0,k}^{\exp }}\right)}^2+{\left(\frac{C_{\mathrm{e},k}^{\mathrm{QC}}}{C_{0,k}^{\exp }}\right)}^2\right)}}},$$

(7)

where $\mathrm{rRSMD}\left(s_{\mathrm{opt}}\right)$ is the value of rRMSD (Equation 5) with the optimal scaling factor given in Equation 6. Note that we do not use weighting with the standard deviations of the experimental fits here because the theoretical calculations have much larger systematic uncertainties that we cannot account for.

3. Resulting Scaling Factors

The resulting scaling factors for various levels of theory, as well as the rRMSD values (Equation 5) for the unscaled and optimally scaled theoretical rotational constants, are given in Table 1. Several trends can be observed from these results. First, the scaling improves the match between theory and experiment in most cases, except for B3LYP-Dn/def2-TZVP ($n=3$(BJ), 4), which we will discuss later. We also see that the increase of the basis set quality from def2-SVP to def2-TZVP improves the agreement between experiment and theory in both scaled and unscaled cases. Applying either the D3(BJ) or D4 dispersion correction leads to the same scaling factors, which probably points to the equal performance of these corrections. A similar trend was observed for the harmonic frequency scaling factors in Ref. [33]. However, the most unexpected yet predictable result is that the optimal scaling factors for the B3LYP-Dn/def2-TZVP levels of theory are equal to one within the margins of error. This means that the scaling does not significantly improve the predicted rotational constant at this approximation. At the same time, the rRMSD values at B3LYP-Dn/def2-TZVP levels are amongst the best in the dataset. Such behavior matches the popularity of these levels of theory for quantum-chemical computations among the rotational spectroscopy community.

4. Illustrative Cases

The simplest way to illustrate the robustness and generality of the scaling procedure is to apply this procedure to cases outside the training dataset. For this, we demonstrate the applicability of the obtained scaling factors for a popular quantum-chemical approximation, namely PBE0-D3(BJ)/def2-TZVP.[37,38,39] The geometries of the molecules discussed here were optimized at this level of theory using the ORCA 5[44,45] software, and then their rotational constants were taken from the calculations. In the case of isotopic substitutions, the rotational constants of the isotopologues were re-computed from the optimized geometries using the UNEX 1.6 software.[46] To demonstrate the numerical performance of the scaling factor, we compared the rRMSD (Equation 5) and mean absolute deviations (MAD) of the rotational constants from the experimentally determined values for the given molecules in the case of scaled and unscaled theoretical equilibrium rotational constants. The MAD values were calculated according to the expression

$$\mathrm{MAD}\left(s\right)=\sum _{k=1}^{N_{\mathrm{mol}}}\left(\left|s·A_{\mathrm{e},k}^{\mathrm{QC}}-A_{0,k}^{\exp }\right|+\left|s·B_{\mathrm{e},k}^{\mathrm{QC}}-B_{0,k}^{\exp }\right|+\left|s·C_{\mathrm{e},k}^{\mathrm{QC}}-C_{0,k}^{\exp }\right|\right).$$

(8)

The first illustrative set of molecules included 15 linear top molecules from di- to pentatomic molecules. Since linear molecules have only one rotational constant, they were excluded from the training set, and the scaling effect on these systems will be the most clearly visible. The calculation of rRMSD and MAD values (Equations 5 and 8) for linear molecules, thus, included only the B rotational constant. The second illustrative set was the case of isotopologues, which, for simplicity of the analysis, were not included in the training dataset. We chose single-substituted isotopologues of imidazole (${\mathrm{C}}_3{\mathrm{H}}_4{\mathrm{N}}_2$), which had rotational constants of three singly substituted ¹³C and two singly substituted ¹⁵N isotopologues available from the literature.[47] The last example was a set of non-covalently bound molecular systems, namely, water – hydrochloric acid clusters $\mathrm{HCl}{\left({\mathrm{H}}_2\mathrm{O}\right)}_n(n=2-5,7)$, which are examples of hydrogen bond network structures. The rotational constants for these species were taken from Refs. [48] and [49].

We can first take a look at a few exemplary cases of molecular systems from our test dataset, two linear molecules (HCN and HCCCN), one imidazole ¹⁵N-substituted isotopologue, namely ${\mathrm{imidazole-}}^{15}$N(1) (nomenclature adapted from Ref. [47]) and also at the largest of our hydrochloric acid clusters, $\mathrm{HCl}{\left({\mathrm{H}}_2\mathrm{O}\right)}_7$. The structures of these molecules and their rotational constants are given in Figure 1 and Table 2. As one can see, B3LYP-D3(BJ)/def2-TZVP provides a reasonable estimation of the rotational constants, closer to the experimental values than unscaled constants at the PBE0-D3(BJ)/def2-TZVP level of theory. However, PBE0-D3(BJ)/def2-TZVP after scaling becomes as accurate or even more accurate than the B3LYP-D3(BJ)/def2-TZVP-based results. This can be seen by comparing the deviations within the datasets. By looking at the rRMSD and MAD values for the scaled and unscaled rotational constants of these systems at the PBE0-D3(BJ)/def2-TZVP level of theory (Table 3), we observe that the scaling indeed improves the agreement of the theoretical and experimental values.

5. Conclusions

In this work, we introduced the concept of scaling factors for rotational constants. Applying a single tabulated scaling factor for all rotational constants is effectively equivalent to scaling the molecular size to account for systematic errors in the equilibrium structure due to the quantum-chemical approximation and for absent anharmonic effects. The set of scaling factors for ten different DFT approximations, namely DF-Dn/def2-mVP ($\mathrm{DF}=\mathrm{B}3\mathrm{LYP},\mathrm{PBE}0$, $n=3\left(\mathrm{BJ}\right),4$, $m=\mathrm{S},\mathrm{TZ}$) and PBEh-3c and ${\mathrm{r}}^2$SCAN-3c, were produced from the database of 174 non-linear molecules. The applicability of these scaling factors was illustrated for the PBE0-D3(BJ)/def2-TZVP level of theory in the case of linear molecules, isotopologues, and non-covalently bonded systems. Thus, the application of such scaling factors can be recommended for the more accurate identification of species in rotational spectra and to support the assignment of specific molecular species in complicated broadband rotational spectra.