On the universality of spectroscopic constants of diatomic molecules

Early in the history of molecular spectroscopy, when it became a discipline within chemical physics in the 1920’s  Herzberg (1985), some intriguing empirical relationships between different spectroscopic properties were observed Kratzer (1920); Mecke (1925); Morse (1929). In particular, it was found that the equilibrium distance, $R_{e}$, and the harmonic vibration frequency, $\omega_{e}$, were correlated in diatomic molecules. As the field evolved, the relationship between $R_{e}$ and $\omega_{e}$ became more evident, and more empirical relations between spectroscopic constants were identified Clark (1934); Badger (1933); Clark (1935); Clark and Stoves (1934); Gordy (1946); Guggenheimer (1946); Clark (1941); Clark and Webb (1941). However, these empirical relationships were typically only valid for certain atomic numbers or groups of the constituent atoms. This result motivated the development of realistic diatomic molecular potentials Morse (1929); Linnett (1940); Newing (1940); Linnett (1942); Varshni (1957, 1958) and triggered the physical chemistry community to think about the “periodicity” of diatomic molecules Hefferlin (1989).

The development of quantum chemistry assisted in the understanding of the physics behind the empirical relationships between spectroscopic constants. In particular, thanks to the application of the Hellmann-Feynman theorem, it was possible to connect $\omega_{e}$ directly with the electronic density at $R_{e}$ Salem (1963); Kim and Parr (1964); Borkman and Parr (1968); King (1968). As a result, a first principles-based explanation of the observed empirical relations between spectroscopic constants appeared Borkman et al. (1969); Politzer (1970); Anderson et al. (1969); Anderson and Parr (1970, 1971a); Simons and Parr (1971); And (1972); Gazquez and Parr (1979). Nevertheless, the obtained relations based on the electronic density were only valid for subsets of molecules, and to date it has not been possible to find universal relations between atomic identifiers and for spectroscopic constants of diatomic molecules.

In this letter, we show that the relationship between spectroscopic constants of diatomic molecules is independent of the kind of molecules at hand and hence universal. Our findings rely upon the application of state-of-the-art non-linear machine-learning (ML) models to an orthodox dataset of spectroscopic constants for diatomic molecules. In particular, we apply the Gaussian process (GP) regression model Williams and Rasmussen (2006) to predict $R_{e}$, $\omega_{e}$, and the binding energy, $D_{0}$, as a function of the groups and periods of the constituent atoms. As a result, it is possible to predict those spectroscopic constants with an accuracy of $\lesssim 5\%$. Finally, we show that the spectroscopic constants coming from high-level electronic structure methods (density functional theory (DFT) and ab initio) display the same relationships as the experimental data employed.

The quest for universal relationships between spectroscopic constants is related to the problem of how atomic and molecular properties describe a spectroscopic property of a molecule, $y=f(\mathbf{x})$. Here, $\mathbf{x}=(x_{1},x_{2},...,x_{n})$, consists of different atomic properties of the constituent atoms or molecular properties, as presented in Fig. 1, where $n$ denotes the number of input features relevant for the problem at hand. Unlike traditional (non-)linear regression models, which assume a fixed form of function $f(\mathbf{x})$, GP embraces a Bayesian perspective and presumes a prior distribution over the space of functions $f(\mathbf{x}_{i})\sim\mathcal{GP}(m(\mathbf{x}_{i}),K(\mathbf{x}_{i},\mathbf{x}_{j}))$ with a joint multivariate-Gaussian distribution, centered at $m(\mathbf{x_{i}})$ and characterized by the covariance function $K(\mathbf{x}_{i},\mathbf{x}_{j})$, which specifies the correlation (or “similarity”) between data points Williams and Rasmussen (2006).

In this work, the spectroscopic properties $y$ are modeled as

The employed dataset contains the main spectroscopic constants: $R_{e}$, $\omega_{e}$, and $D_{0}$ for the ground electronic state of heteronuclear diatomic molecules. In particular, it contains the experimental values of $R_{e}$, $\omega_{e}$ for 256 heteronuclear diatomic molecules taken from Refs. Liu et al. ; Huber and Herzberg (1979); Smirnov (2008), whereas the experimentally determined values of $D_{0}$ are only available for 197 of them. Fig. 2 shows the distribution of equilibrium distance of the molecules in the dataset and its ratio to the sum of the atomic radii of the constituent atoms, $R_{1}+R_{2}$. Most of the molecules show an equilibrium distance between 1.4 Å and 3.8 Å, with a most probable value of 1.7 Å. Looking at the values of $R_{e}/(R_{1}+R_{2})$, it is clear that the molecules within the dataset have different bonds: covalent, van der Waals and ionic. In particular, the present dataset contains a majority of covalent and ionic molecules and only a handful of van der Waals molecules, as shown in Fig. 2.

Fueled by the idea of the periodicity of the molecules (see, e.g., Ref. Hefferlin (1989) and references in it), we use groups, $g_{k}$, and periods, $p_{k}$, of the atoms within a molecule, i.e., $k=1,2$, as input features for a GP regression model to predict different combinations of spectroscopic constants, in particular, $R_{e}$, $\omega_{e}$ and $\log\frac{D_{0}}{Re^{3}Z_{1}Z_{2}}$, where $Z_{k}$ stands for the atomic number of the $k$-th atom of the molecule. The last combination of spectroscopic constants was first proposed by Anderson, Parr, and coworkers, by means of the Born-Oppenheimer approximation and assuming that the electron density of an atom within a molecule decays exponentially at the equilibrium distance Anderson and Parr (1971a); Anderson and Parr (1971b); And (1972); Gazquez and Parr (1979). In the GP regression model, as customary in ML, the dataset is divided into training and test sets. However, the present dataset is rather small from a ML perspective. When the dataset is split into training and test sets, the training set may not be representative. This may lead to a bias in the performance of the test set. To solve this problem, we have employed a Monte Carlo (MC) approach, in which the dataset is stratified into $25$ strata based on the level of the true values of the labels ($R_{e}$, $\omega_{e}$, and $\log\frac{D_{0}}{Re^{3}Z_{1}Z_{2}}$ in the present work). In each MC step, the training and test data are randomly selected in each stratum. The stratification helps to minimize the change of the proportions of the dataset compositions upon splitting Raschka (2018). Additionally, in training the GP regression models, the training sets are further split to perform stratified 5-fold cross-validation for hyperparameter optimization to avoid overly optimistic estimates of the model performance.

The performance of a GP regression model may be quantified by the root mean square error (RMSE) as $\text{RMSE}=\sqrt{\frac{1}{N}\sum_{i=1}^{N}{(y_{i}-y^{*}_{i})}^{2}}$, where $y_{i}$ are the true labels (experimental values) and $y^{*}_{i}$ are the predictions, and the normalized error defined as the ratio of the RMSE to the range of the data $r_{E}=\text{RMSE}/(y_{max}-y_{min})$.

The GP regression model predictions of $R_{e}$ as a function of input features $(g_{1},g_{2},p_{1},p_{2})$ in comparison with its true values are displayed in panel (a) of Fig. 3, which shows little dispersion of the predicted values with respect to the true values. To further quantify the GP regression model performance we calculate the average RMSE of the predicted $R_{e}$ on 1000 randomly selected test sets leading to $0.0968~{}\pm~{}0.0070$ Å(Table 1), and $r_{E}$ = $2.80~{}\pm~{}0.20\%$. The model performance increases as the number of molecules in the training set $N$ grows. It is not yet converged for $N=231$, suggesting that the GP regression model can be further improved by learning from more data in the training set SI .

Panel (b) of Fig. 3 shows the comparison between the predicted $\omega_{e}$ and its true value. $({R^{*}_{e}}^{-1},g_{1},g_{2},p_{1},p_{2})$ are used as input features where $R^{*}_{e}$ is the predicted equilibrium distance from $(g_{1},g_{2},p_{1},p_{2})$. The accuracy of GP regression model for $\omega_{e}$ is characterized by RMSE = $227.5~{}\pm~{}4.6$ cm^-1 and $r_{E}=5.56~{}\pm~{}0.11\%$. Thus, the GP regression model accurately predicts $\omega_{e}$, nevertheless for some molecules the prediction is not as good as expected. The overall performance of the GP regression model (Table 1), is improved by utilizing $({R^{*}_{e}}^{-1},g_{1}^{iso},g_{2}^{iso},p_{1},p_{2})$ as input features, where $g_{k}^{iso}$ encodes the information about the hydrogen isotopes of the $k$-th atom in the molecule. Thus, the different isotopologues in the dataset are adequately addressed by input features as explained in the Supplemental Material SI . Further improvement is possible by using the true value of the equilibrium distance, leading to a more precise prediction of $\omega_{e}$, as displayed in Table 1. Despite the improvement on the description of $\omega_{e}$ the outliers shown in panel (b) of Fig. 3 remain, as the ones in panel (c). Indeed, we have identified the majority of these outliers as bi-alkali molecules, for a detailed description see Supplemental Material SI .

The GP regression model prediction of $\log\frac{D_{0}}{R_{e}^{3}Z_{1}Z_{2}}$ vs. its true value is shown in panel (c) of Fig. 3, which shows an outstanding performance. The performance is further supported by an RMSE = $0.451~{}\pm~{}0.007$ and an $r_{E}$ equal to $4.45~{}\pm~{}0.06\%$, as shown in Table 1. In this case, the GP is fed with $({R^{*}_{e}},\bar{g},\bar{p})$ as input features, where $\bar{g}=\frac{g_{1}+g_{2}}{2}$ and $\bar{p}=\frac{p_{1}+p_{2}}{2}$. It is worth mentioning that the accuracy of the GP regression model regression is a factor of 2 better than using standard regression techniques, as shown in the Supplemental Material SI . The learning curve is converged for $n=125$, which suggests that accurate predictions can be made by the models learned from a minimal training set. As a result, with only two atomic properties as input features it is possible to predict the value of $\log\frac{D_{0}}{R_{e}^{3}Z_{1}Z_{2}}$ with an accuracy better than $5\%$. The use of the true value for the equilibrium distance reduced the RMSE, as expected, by $20\%$, as shown in Table 1.

From the GP regression models for the three combinations of spectroscopic constants shown in Fig. 3, $\omega_{e}$ shows the largest error bars regarding the MC sampling technique. This behavior is due to the large variation of the values of $\omega_{e}$ in every iteration. Indeed, by looking at the distribution of $R_{e}^{2}\omega_{e}$ shown in the inset of panel (b) and the distribution of $R_{e}$ in Fig. 2, it is clear that $\omega_{e}$ shows a broad distribution with a multi-modal character. Surprisingly enough, $R_{e}^{2}\omega_{e}=\text{constant}$ seems to be a trend for most of the diatomic molecules, and extensible to any excited state Mecke (1925); Borkman and Parr (1968). However, we do not observe this, as the box plot underneath the inset in panel (b) emphasizes. Indeed, the variation of $R_{e}^{2}\omega_{e}$ may be related to different underlying 2-body potentials for the diatomic molecules considered. Hence, it is another way to show the different bond mechanisms of the molecules within the dataset.

As shown above, the GP regression model for the spectroscopic constants shows a clear universal trend between them. However, the dataset only contains experimentally determined values for spectroscopic constants. Thus, the question arises: do the computationally derived spectroscopic constants lead to the same universal relationships? To answer this question, we have conducted high-level electronic structure calculations using the Molpro package Werner et al. (2019); W. et al. (2012) with the aug-cc-pV5Z basis set Pritchard et al. (2019) to calculate the ground state electronic potential energy curve of $81$ molecules in the dataset. The calculations are performed via coupled-cluster with single, double, and perturbative triple excitations CCSD(T) Raghavachari et al. (1989); Bartlett et al. (1990) and with the B3LYP Stephens et al. (1994) functional as DFT method. From the potential energy curves, we estimate $R_{e}$, and $\omega_{e}$ whose distributions in comparison with the experimental values in the data set are shown in Fig. 4. As a result, the distributions of $R_{e}$ and $\omega_{e}$ from B3LYP and CCSD(T) show the same features as the experimental one. Therefore, the universal relations for the present dataset are equally fulfilled for high-level electronic structure methods, as it is further elaborated in Supplemental Material SI .

In summary, we have shown, using the GP regression model, that the main spectroscopic constants of diatomic molecules are universally related. This result confirms the scenario that Kratzer and Mecke envisioned a century ago Mecke (1925); Kratzer (1920). The relationships are independent of the nature of the chemical bond of the diatomic molecule. In particular, we have demonstrated that merely using the group and period of the atoms within a molecule as input features it is possible to predict particular combinations of spectroscopic constants with an error $r_{E}<5\%$. In other words, the spectroscopic constants of diatomic molecules can be efficiently learned from an appropriate dataset by a GP regression model, and their values accurately predicted without carrying out quantum chemistry calculations. Besides, the high-level electronic structure calculations (ab initio and DFT) for the spectroscopic constants show the same distribution as the experimental ones. We conclude that the computationally-derived spectroscopic constants follow the same universal trends as the experimental ones, which are employed in our GP regression model. From our perspective, having reasonable estimates of the main spectroscopic constants will help to optimize the experimental efforts in performing spectroscopy of diatomic molecules. In the same vein, the present work may motivate data science-driven studies on the field of spectroscopy of diatomic molecules. In particular, it will help to evolve the field of spectroscopy towards the current information era.

Finally, we would like to emphasize that there are around $7000$ heteronuclear molecules, and we only utilize $256$ of these for our GP regression model. The limited availability of spectroscopic data (only around $3\%$ of possible heteronuclear diatomic molecules) shows the vast amount of spectroscopy that can be done within the realm of diatomic molecules. The more data we have, the more accurate will be the GP regression model predictions before reaching convergence of the learning curve, and more knowledgeable the community will be about the fundamental properties of diatomic molecules.

We thank Dr. Matthias Rupp for his comments and suggestions and Drs. Daniel Thomas and Uwe Hergenhahn for carefully reading the manuscript.