1 Introduction

The quest for a reliable yet practical modeling of large molecular systems has always played a central role in the field of theoretical and computational chemistry, and this balance between accuracy and feasibility is likely to persist into the future [1]. Combination of different complementary techniques is often the key to determining in an univocal way the molecular structure and bond topology for a wide range of systems, ranging from small molecules to large biomolecules. In particular, high-resolution spectroscopic studies both in the gas-phase (microwave, MW) and in inert matrices (infrared, IR) allow an unbiased disentanglement of intrinsic stereo-electronic effects without any strong perturbing effect from the environment. However, the tuning of experimental spectra by different low-lying structures and the interplay of different contributions can make the interpretation of experimental data troublesome or even impossible without the aid of trustworthy in silico simulations. In fact, ongoing improvements of hardware, software, and, above all, underlying physical–mathematical models are allowing the analysis of experimental data and their interpretation for molecular systems of increasing complexity. Despite the undisputed effectiveness of static structure–property correlations and the fundamental rigid rotor/harmonic oscillator (RRHO) model, results that are directly comparable to experiment can only be achieved through more advanced models: (i) at the electronic level, employing highly-correlated methods; (ii) at the nuclear level, by incorporating anharmonic effects in the description of the nuclear motions. In fact, a strong limitation of the harmonic approximation is that the real vibrations are intrinsically anharmonic, and as a consequence vibrational energies are systematically overestimated. A common empirical scheme for improving the agreement with experimental results rests on the application of scaling factors depending on the employed quantum chemical (QC) method and, possibly, on the range of investigated frequencies. However, those scaling factors, generally obtained from statistical analyses of fundamental transitions, do not improve intensities and are often not suitable for higher-quanta transitions. The introduction of anharmonic effects requires a representation of potential energy surfaces (PESs) around stationary points beyond the quadratic approximation, introducing a series of higher-order terms involving one or more vibrational coordinates. In this context, a crucial point is represented by the inclusion of couplings involving different vibrations, which have to be properly considered in order to obtain accurate estimates of line shapes, frequency shifts, etc. Computational strategies aimed at the study of vibrational properties of molecular systems are usually based on perturbative [2,3,4,5,6,7,8,9] variational [10,11,12] approaches or their combinations [13, 14]. Within variational methods, the energy is minimized starting from an expansion of the wave functions over a basis of known states. This class of methods includes the vibrational self-consistent-field (VSCF) [12, 15,16,17,18,19,20] and vibrational configuration interaction (VCI) [21,22,23]. While the use of accurate variational approaches in conjunction with high-order force fields allows to obtain estimates vibrational energies close to their experimental counterparts, their use is mainly limited to systems containing few atoms due to their unfavorable scaling with the size of the molecule [24]. With the aim of reducing the computational cost, several strategies have been devised, such as the iterative subspace expansion algorithms [25,26,27,28] and protocols for the reduction of the VCI basis [29, 30]. Over the last decades, several variational approaches have been proposed for the resolution of the vibrational problems, not only relying on the Watson’s Hamiltonian and exploiting different sets of nuclear coordinates. In particular, different approaches and sophisticated variational, semi-variational nuclear motion and dynamics theoretical frameworks (and codes) built upon the use of internal coordinates have been discussed by Carrington [31, 32], Csaszar [33, 34], Lauvergnat [35,36,37], Tennyson [38, 39], and Yurchenko [40], among others. Variational approaches exploiting Monte Carlo diffusion methods in internal coordinates have been also investigated [41, 42], as well as the use of optimized and localized coordinates [43, 44], for example, in the vibrational coupled-cluster (VCC) framework [45]. Another option is the parametrization of the vibrational wavefunction in tensor formats such as canonical decomposition [46] and matrix product states. Following that, the wavefunction representations can be optimized through the recent extension of the density matrix renormalization group (DMRG) [47, 48] algorithm to the study of the vibrational problem, the resulting method being referred to as vibrational DMRG (vDMRG) [49, 50].

When comparing the various methods for including anharmonic effects, [2, 3, 16,17,18,19, 21, 22, 51,52,53] those based on perturbation theory applied to the Watson Hamiltonian offer a remarkable balance between accuracy and computational cost provided that resonance effects are taken into proper account. In particular, the vibrational second-order perturbation theory (VPT2), [2, 54] based on a fourth-order polynomial approximation of the potential energy in terms of normal coordinates, allows to study medium-to-large size molecular systems. Furthermore, the formulation of VPT2 based on the Van Vleck contact transformation method [55] enables the extension to a generalized model (GVPT2) [6]. Within this framework, the diagonalization of a limited set of reduced-dimensionality Hamiltonians involving strongly interacting states is carried out, while the other states are corrected by second-order perturbative contributions. At the VPT2 level, also the vibrational dependence of rotational constants on molecular vibrations can be expressed through a perturbative expansion leading to the determination of specific vibro-rotational interaction constants and paving the route for the calculation of accurate molecular structures through the semi-experimental (SE) approach [56].

While VPT2 can be successfully applied in many cases, it is also characterized by some intrinsic drawbacks that prevent its use in some scenarios. As mentioned above, the first problem is related to the presence of resonances (both Fermi and Darling-Dennison), even though several computational strategies as well as robust schemes have been developed over the years [6, 57,58,59,60,61,62] even for the treatment of vibrational intensities. Second, the theoretical framework largely depends on the symmetry of the molecule since the presence of degenerate vibrations implies a distinct derivation of the equations [63, 64]. As a consequence, different derivations and implementations of VPT2 are required for the treatment of systems belonging to different classes of symmetry. Recently, this issue was solved by the formulation of a unified model relying on the proper application of a posteriori transformations to the wavefunctions [65]. Last but not least, large amplitude motions (LAMs) are poorly described at the VPT2 level and more generally with a quartic force field. Unfortunately, variational approaches such as the VCI or more refined perturbative methods (such as VPT4 [66]) become rapidly prohibitive as the size of the molecular systems increases, due to the necessity of more demanding anharmonic force field calculations and the increasing complexity concerning the vibrational calculation. On the other hand, low-dimensionality methods tailored for describing one or a limited number of LAMs [67,68,69,70,71,72] imply their separation from the rest of vibrations, usually referred to as small amplitude motions (SAMs). As shown in a previous study, [73] the definition of a suitable set of internal coordinates able to decouple these two classes of vibration allows to treat each normal mode through the most appropriate method. Starting from the available literature, [74, 75] the VPT2 framework in terms of curvilinear coordinates was developed and successfully applied to both semi-rigid and flexible systems, showing in the latter case that not only LAMs can be effectively decoupled from SAMs, but also that each LAM can be treated independently from the others. As a matter of fact, SAMs can be treated through the internal-based VPT2, while variational approaches can be employed for the calculation of the energy levels of the floppy degrees of freedom. When a single LAM is present, it can be effectively treated through the discrete variable representation (DVR) [76,77,78,79,80,81,82,83]. On the other hand, a higher number of LAMs can be in principle handled through the reaction-surface Hamiltonian (RSH) [84], reaction-volume Hamiltonian (RVH) [85] or VCI-based methodologies.

In this paper, the versatility of internal coordinates and their application in both fields of rotational and vibrational spectroscopy will be described. The paper is organized as follows. First, an overview of the Cartesian-based version of VPT2 will be provided, including the calculation of the vibrational contributions to the rotational constants required in the SE approach. Particular attention will be devoted to the use of different sets of internal coordinates in the nonlinear least-squares fit needed for the calculation of the SE equilibrium molecular geometry. Then, the focus will shift to the extension of VPT2 to the use of internal coordinates, from the calculation of the anharmonic energy levels to the definition of a robust scheme for the automatic identification of Fermi resonances. In this context, the importance of internal coordinates in reducing inter-mode couplings will be underlined through the application to systems of biological and technological interest.

2 Theory

2.1 Overview of the Cartesian-based VPT2

In this section, the unified Cartesian-based, VPT2 framework for asymmetric, symmetric, linear, and spherical tops is described. In this theoretical model, all normal modes are treated independently from one another regardless the presence of vibrational degeneracies. In the latter instance, a series of linear transformations are employed a posteriori to determine the correct energies, wave functions and properties. The discussion is structured to give a broad overview of this methodology up to the derivation of the vibrational energy levels also accounting for the presence of Fermi resonances.

Let us consider a molecule composed of N vibrational modes (N being equal to \(3N_a-6\) for nonlinear molecules and to \(3N_a-5\) for linear ones, where \(N_a\) is the number of atoms). The reference Hamiltonian is obtained by expanding the more general vibro-rotational Hamiltonian proposed by Watson [86] in terms of normal coordinates and then only considering the purely vibrational terms:

$$\begin{aligned} \begin{aligned} \mathcal {H}&= \dfrac{1}{2}\sum _{i=1}^N\omega _i(p_i^2 + q_i^2) + \dfrac{1}{6}\sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N\textsf{f}_{ijk}q_iq_jq_k\\\quad&+\dfrac{1}{24}\sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N\sum _{l=1}^N\textsf{f}_{ijkl}q_iq_jq_kq_l\\&+\sum _{\tau =x,y,z}B^\mathrm {\scriptscriptstyle eq}_\tau \sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N\sum _{l=1}^N\zeta _{ij,\tau }\zeta _{kl,\tau }\sqrt{\dfrac{\omega _j\omega _l}{\omega _i\omega _k}}q_ip_jq_kp_l + {\mathcal{U}} \end{aligned} \end{aligned}$$
(1)

where \(q_i\) is the i-th dimensionless normal coordinate, \(p_i\) is its conjugate momentum, \(\omega _i\) is the corresponding harmonic wavenumber (in cm\(^{-1}\)), \(B^\mathrm {\scriptscriptstyle eq}_\tau\) is the equilibrium rotational constant along the \(\tau\) axis and \(\boldsymbol{\zeta}\) is the anti-symmetric matrix of Coriolis couplings. The terms \(\textsf{f}_{ijk}\) and \(\textsf{f}_{ijkl}\) are, respectively, referred to as cubic and quartic force constants, defined by the following notation,

$$\begin{aligned} \textsf{f}_{ijk\dots } = \left( \dfrac{\partial ^n\mathcal {V}}{\partial q_i\partial q_j\partial q_k\dots }\right) _\mathrm {\scriptscriptstyle eq}\end{aligned}$$
(2)

where \(\mathcal {V}\) represents the potential energy, while \(\mathcal {U}\) is a mass-dependent term which does not contribute to the calculation of the transition energies, and for this reason it will be not considered from now on.

Within the VPT2 framework, the vibrational Hamiltonian is expanded through a perturbative series up to the second order, the resulting Schrödinger equation being typically solved by means of two different approaches, namely the Rayleigh–Schrödinger (RSPT) and Van Vleck (CVPT) [55] perturbation theory. As a result, the VPT2 vibrational energy of a state \(|{\varvec{ v }}_\mathrm{R} \rangle\) is

$$\begin{aligned} \varepsilon _\mathrm{R}= \varepsilon _0 + \sum _{i=1}^N\omega _i v _{{{\scriptscriptstyle R},i}} + \sum _{i=1}^N\sum _{j=i}^N\chi _{ij}\left[ v _{{{\scriptscriptstyle R},i}} v _{{{\scriptscriptstyle R},j}} + \dfrac{1}{2}\left( v _{{{\scriptscriptstyle R},i}}+ v _{{{\scriptscriptstyle R},j}}\right) \right] \end{aligned}$$
(3)

where \(v _{{{\scriptscriptstyle R},i}}\) is the number of quanta associated with the mode i in the state R, and \(\varepsilon _0\) is the anharmonic resonance-free zero-point vibrational energy (ZPVE) [9, 87, 88].

The \({\varvec{{\chi }}}\) matrix contains the anharmonic corrections and is defined by the following expressions,

$$\begin{aligned} 16\chi _{ii}&= \textsf{f}_{iiii} - \frac{5{\textsf{f}_{iii}}^2}{3\omega _i} - \sum _{\begin{array}{c} j=1\\ \left( j\ne i\right) \end{array}}^N \dfrac{\textsf{f}_{iij}^2}{2}\nonumber \\&\quad \left[ \dfrac{4}{\omega _j} + \dfrac{1}{2\omega _i+\omega _j}-\dfrac{1}{2\omega _i-\omega _j}\right] \end{aligned}$$
(4)
$$\begin{aligned}&\begin{aligned} 4\chi _{ij}&= \textsf{f}_{iijj} - \dfrac{{\textsf{f}_{iij}}^2}{2}\left[ \dfrac{1}{2\omega _i+\omega _j}+\dfrac{1}{2\omega _i-\omega _j}\right] \\&\quad - \dfrac{{\textsf{f}_{ijj}}^2}{2}\left[ \dfrac{1}{2\omega _j+\omega _i}+\dfrac{1}{2\omega _j-\omega _i}\right] \\&\quad - \frac{\textsf{f}_{iii}\textsf{f}_{ijj}}{\omega _i} - \frac{\textsf{f}_{jjj}\textsf{f}_{iij}}{\omega _j} \\&\quad + \sum _{\begin{array}{c} k=1\\ \left( k\ne i,j\right) \end{array}}^N \left\{ \dfrac{{\textsf{f}_{ijk}}^2}{2}\left[ \dfrac{1}{\omega _i+\omega _j+\omega _k} + \dfrac{1}{\omega _i-\omega _j+\omega _k}- \dfrac{1}{\omega _i+\omega _j-\omega _k}\right. \right. \\&\left. \left. \quad - \dfrac{1}{\omega _i-\omega _j-\omega _k} \right] - \frac{\textsf{f}_{iik}\textsf{f}_{jjk}}{\omega _k}\right\} \\&\quad + 4\left( \dfrac{\omega _i}{\omega _j} + \dfrac{\omega _j}{\omega _i}\right) \sum _{\tau =x,y,z} B ^{\mathrm {\scriptscriptstyle eq}}_\tau \{\zeta _{ij,\tau }\}^2 \end{aligned} \end{aligned}$$
(5)

Inspection of Eqs. 4 and 5 evidences that the presence of denominators approaching zero would lead to unphysical results. These conditions, collectively known as Fermi resonances (FRs), occur when \(\omega _i\approx 2\omega _j\) or \(\omega _i\approx \omega _j+\omega _k\). The first case is commonly referred to as type I, while the second as type II [89]. Several strategies can be adopted to handle this problem. In the deperturbed VPT2 (DVPT2), a sequential screening of all potentially resonant terms is performed, and an analysis based on one or more criteria allows to state which terms should be removed from the anharmonic calculation. Generally, this procedure is organized into two step. The first one is the evaluation of the energetic proximity of the interacting states at the harmonic level,

$$\begin{aligned} |\omega _i - \omega _j - \omega _k|\le \Delta \omega ^{1-2} \end{aligned}$$
(6)

where i and j can be equal and \(\Delta \omega ^{1-2}\) is a threshold defined a priori. Second, the overall weight of the term is estimated. Several strategies have been proposed for the latter [7, 58, 90]. In this work, the test proposed by Martin and co-workers [90] has been employed, leading to the following condition, [91]

$$\begin{aligned} \dfrac{{\textsf{f}_{ijk}}^4}{64(1+\delta _{jk})^2|\omega _i - \omega _j - \omega _k|^3} \ge K^{1-2} \end{aligned}$$
(7)

where \(K^{1-2}\) is a second threshold required in the procedure. In the DVPT2 scheme, each term fulfilling Eqs. 6 and 7 is labeled as resonant and removed from the anharmonic calculation. If on the one side this method avoids singularities in the calculation of the energy levels, on the other side it can lead to a truncate treatment, since the resonant terms are systematically neglected. In order to prevent this situation, the interaction terms related to FRs can be introduced back in a successive, variational step. For the purpose, a variational matrix \(\tilde{{\varvec{\textbf{H}}}}\) is built starting from the DVPT2 energies and Van Vleck Hamiltonian interaction terms. A new set of energies is obtained by diagonalizing \(\tilde{{\varvec{\textbf{H}}}}\), the full procedure being referred to as GVPT2\(^\text {F}\) (hereafter, simply GVPT2).

2.2 Vibrational corrections

The rotational constants of a target molecule are among the most important parameters obtained from the study of rotational and vibro-rotational spectra, and they are inversely proportional to the principal moments of inertia. In general, the availability of this kind of data for different isotopic species allows for the determination of structural parameters like bond lengths and angles [92]. Nevertheless, since molecules are not rigid rotors and are subject to vibrational motion, the notion of a reference molecular geometry is far from being a simple concept. An important step forward has been done by Pulay and co-workers [56] through the introduction of the SE approach for the determination of the equilibrium structure of a molecule. While more difficult to measure at the experimental level, this kind of structure properly accounts for vibrational effects and it is independent from the isotopic species. Furthermore, it is directly comparable to quantum-mechanical data. Within the SE method, the structural parameters, are obtained through a fit of the principal moments of inertia (or the corresponding rotational constants). Hence, the SE rotational constant \(B^\text {SE}_\tau\) is obtained as follows [93],

$$\begin{aligned} B^\text {SE}_\tau = B^{0,\text {exp}}_\tau - (\Delta B^0_\tau )^\text {QM} \end{aligned}$$
(8)

where \(\tau = a,b,c\) is one of the principal inertial axes, while \(B^{0,\text {exp}}\) and \((\Delta B^0_\tau )^\text {QM}\) are, respectively, the experimental rotational constant and its correction along \(\tau\). In general terms, \((\Delta B^0_\tau )^\text {QM}\) has a double nature. It is composed of a vibrational contribution (\(\Delta B^\text {vib}_\tau\)) stemming from the VPT2 framework, and an electronic contribution (\(\Delta B^\text {ele}_\tau\)) [93,94,95] related to the effect of electron on the moments of inertia. In this context only the vibrational contribution will be analyzed, since the electronic term is often negligible and required only for highly accurate calculations.

The vibrational contribution can be derived from the contact-transformed Hamiltonian, leading to an explicit dependence of rotational constants on vibrational degrees of freedom [94],

$$\begin{aligned} B^{{\varvec{ v }}_\mathrm{R}}_\tau = B^\mathrm {\scriptscriptstyle eq}_\tau - \sum _{i=1}^N\alpha _{\tau ,i}\left( v _{{{\scriptscriptstyle R},i}}+\dfrac{1}{2}\right) \end{aligned}$$
(9)

where \(B^{{\varvec{ v }}_R}_\tau\) and \(B^\mathrm {\scriptscriptstyle eq}_\tau\) are the rotational constants associated with the vibrational state \(|{\varvec{ v }}_\mathrm{R} \rangle\) and the corresponding equilibrium value, respectively, while \(\alpha _{\tau ,i}\) are the vibro-rotational interaction constants:

$$\begin{aligned} \alpha _{i,\tau }= & {} -2(B^\mathrm {\scriptscriptstyle eq}_\tau )^2\left[ \sum _{\eta = x,y,z}\dfrac{3(a_{i,\tau \eta })^2}{4\omega _i I^\mathrm {\scriptscriptstyle eq}_{\eta }} \right. \nonumber \\{} & {} + \left. \sum _{j=1}^N\dfrac{(\zeta _{ij,\tau })^2(3\omega ^2_i + \omega ^2_j)}{\omega _i(\omega ^2_i - \omega ^2_j)} + \pi \sqrt{\dfrac{c}{h}}\sum _{j=1}^N\dfrac{\textsf{f}_{iij}a_{j,\tau \tau }}{\omega ^{3/2}_j}\right] \end{aligned}$$
(10)

where \(\textsf{a}_{i,\tau \eta }\) are the derivatives of the effective moments of inertia tensor evaluated at the equilibrium geometry. When the vibrational ground state is considered, Eq. 9 can be rearranged to define the vibrational contribution \(\Delta B^\text {vib}_\tau = B^0_\tau - B^\mathrm {\scriptscriptstyle eq}_\tau\) as,

$$\begin{aligned} \begin{aligned} \Delta B^\text {vib}_\tau&= B^\mathrm {\scriptscriptstyle eq}_\tau - B^0_\tau \\&= -(B^\mathrm {\scriptscriptstyle eq}_\tau )^2\left[ \sum _{\eta = x,y,z}\sum _{i=1}^N\dfrac{3(a_{i,\tau \eta })^2}{4\omega _i I^\mathrm {\scriptscriptstyle eq}_{\eta }} \right. \\&\quad \left. + \sum _{i=1}^N\sum _{j=1}^N\dfrac{(\zeta _{ij,\tau })^2(3\omega ^2_i + \omega ^2_j)}{\omega _i(\omega ^2_i - \omega ^2_j)}\right. \\&\quad + \left. \pi \sqrt{\dfrac{c}{h}}\sum _{i=1}^N\sum _{j=1}^N\dfrac{\textsf{f}_{iij}a_{j,\tau \tau }}{\omega ^{3/2}_j}\right] \end{aligned} \end{aligned}$$
(11)

As it can be deduced from Eq. 10, the coefficients \(\alpha _{\tau ,i}\) are affected by resonances, occurring when \(\omega _i\approx \omega _j\). However, in the calculation of the vibrational correction this issue is solved by recasting the Coriolis contribution as follows,

$$\begin{aligned} \sum _{i=1}^N\sum _{j=1}^N\dfrac{(\zeta _{ij,\tau })^2(3\omega ^2_i + \omega ^2_j)}{\omega _i(\omega ^2_i - \omega ^2_j)} = - \sum _{i=1}^N\sum _{j=i+1}^N\dfrac{(\zeta _{ij,\tau })^2(\omega _i-\omega _j)^2}{\omega _i\omega _j(\omega _i + \omega _j)} \end{aligned}$$
(12)

Once the set of SE rotational constants is assembled, a nonlinear least-squares fit is carried out in order to determine the molecular geometry, which is described by a proper set of nuclear coordinates. For the purpose, the theoretical framework has been implemented in the MSR (Molecular Structures Refinement) software, [96,97,98] developed in our research group and designed for the determination of equilibrium structures through the SE approach. Over the last years, the MSR code has been employed for the characterization of numerous molecular systems including sulfur-containing [96, 99] and astrochemical systems, [100] biological building blocks [97] and non-covalent complexes [101]. The program has been developed with the target of being as general as possible, and equipped with a series of features at different levels. First, a wide range of optimization algorithms (such as the Gauss–Newton and Levenberg-Marquardt) has been included. Second, the so-called predicate observations [102,103,104] can be added to the set of reference data. Within this approach, the set of SE data is augmented with quantum-mechanical estimates of one or more parameters. This feature can be particularly advantageous when the number of isotopic species is not sufficiently large for the full structural characterization. Third, a detailed error analysis has been implemented, with the possibility to calculate, for example, the standard deviations on the singular parameters, the outliers and the condition number.

2.3 Curvilinear coordinates formalism

As previously anticipated, the choice of a proper set of curvilinear internal coordinates plays a central role in vibro-rotational analyses. For the present discussion, let us consider a set of non-redundant internal coordinates \({\varvec{ s }}\). Since the latter are nonlinear functions of Cartesian coordinates, a strategy commonly adopted is represented by a Taylor-series expansion around the reference structure,

$$\begin{aligned} s_i = s^\mathrm {\scriptscriptstyle eq}_i + \sum _{i=1}^NB_{ij}(x_j - x^\mathrm {\scriptscriptstyle eq}_j) + \dfrac{1}{2}\sum _{j=1}^{3N_a}\sum _{k=1}^{3N_a}B^{'}_{ijk}(x_j - x^\mathrm {\scriptscriptstyle eq}_j)(x_k - x^\mathrm {\scriptscriptstyle eq}_k) + \dots \end{aligned}$$
(13)

where \({\varvec{ x }} = \{x_1,x_2,\dots ,x_{3N_a}\}\) collects the nuclear Cartesian coordinates, \(N_a\) is the number of atoms, and the Wilson \({\varvec{\textbf{B}}}\) matrix and its Cartesian derivative \({\varvec{\textbf{B}}}^{'}\) have been introduced,

$$\begin{aligned} B_{ij} = \left( \dfrac{\partial s_i}{\partial x_j}\right) _\mathrm {\scriptscriptstyle eq}\qquad B^{'}_{ijk} = \left( \dfrac{\partial ^2 s_i}{\partial x_j\partial x_k}\right) _\mathrm {\scriptscriptstyle eq}\end{aligned}$$
(14)

It is worth pointing out that Eq. 13 can be truncated at different levels depending on the area of study. In the present work, a first-order expansion is sufficient to obtain a set of internal coordinates for the structural refinement, while the Wilson \({\varvec{\textbf{B}}}^{'}\) tensor is required in the field of anharmonic vibrational spectroscopy.

Concerning the calculation of accurate molecular structures, Z-matrix internal coordinates (ZICs) represent the most common choice. ZICs are mostly based on chemical intuition; however, they are not completely unambiguous. Because of this, selecting a different set of ZICs could provide a distinct outcome, indicating a substantial user-dependency. Moreover, poorly designed Z-matrices could prevent the optimization process from converging. In order to bypass this problem, a protocol based on molecular symmetry has been adopted. The latter is supposed to be maintained during the entire optimization procedure, and the formulation of suitable geometrical constraints applied to the internal coordinates is typically sufficient when ZICs are the reference set of coordinates. Although this approach has shown to be quite effective, the use several dummy atoms may be necessary when studying complex chemical topologies, which would involve the addition of numerous degrees of freedom. From this perspective, a remarkable improvement can be achieved by incorporating the symmetry in the definition of coordinates itself. More specifically, the totally-symmetric (A\(_1\)) coordinates selected from a non-redundant set of symmetry internal coordinates are used, while all the other ones are kept fixed at their guess values. In the MSR code, the delocalized internal coordinates (DICs) [105,106,107,108,109] are used as a reference set, since they are characterized by simple criteria for the identification of redundancies. In order to build DICs, the starting point is the definition of \(N_r\) redundant internal coordinates, which usually correspond to the so-called primitive internal coordinates (PICs), given by the full list of all bond lengths, valence and dihedral angles. DICs are built as linear combinations of PICs, the transformation matrix \({\varvec{\textbf{U}}}\) ·being generated by the eigenvectors of \({\varvec{\textbf{B}}}{\varvec{\textbf{B}}}^\dagger\) corresponding to non-null eigenvalues,

$$\begin{aligned} {\varvec{B}}{\varvec{B}}^\dagger ({\varvec{U}}{\varvec{R}}) = \begin{pmatrix} \varvec{\Lambda } &{} {\varvec{0}} \\ {\varvec{0}} &{} {\varvec{0}} \end{pmatrix} ({\varvec{U}}{\varvec{R}}) \end{aligned}$$
(15)

while \({\varvec{\textbf{R}}}\) contains the redundancies. The selection of A\(_1\) coordinates is carried out through a multi-step procedure which is shortly outlined here, while a full description is reported in Ref. [97]. First, the guess structure is displaced along each DIC, and the resulting geometry is converted to its Cartesian counterpart. The latter is then subjected to the application of all symmetry operations of the point group, and only if it remains unchanged, the corresponding displacement internal coordinate is marked as totally symmetric. Finally, the matrix product \({\varvec{\textbf{B}}}{\varvec{\textbf{B}}}^\dagger\) expressed in terms of DICs can be cast in a block-diagonal form, each block being related to an irreducible representation of the point group. Once all A\(_1\) coordinates have been detected, their presence in the same block is checked as a further confirmation. The use of A\(_1\) coordinates in the optimization presents two basic advantages. In fact, it is characterized by a black-box definition of a non-redundant set of internal symmetry coordinates. Second, the generation of geometrical constraints is completely automatized in order to handle in a transparent way even large molecular systems. With the aim of highlighting the effectiveness of the mentioned procedure, a complete study at different levels will be presented later. In particular, the stability of the fit with respect to the set of coordinates has been object of interest. The accuracy of the proposed methodology will be shown through the analysis of systems of biological and astrochemical interest. Furthermore, the comparison of our results with advanced, highly accurate protocols will be object of discussion as well.

2.4 Internal-based VPT2 framework

In this section, the main aspects of the internal-based VPT2 framework are described, while interested readers are referred to Ref. [65] for a more detailed discussion.

In order to set up the expressions required for the calculation of VPT2 energies within the internal-based VPT2 framework, the first step is the definition of the vibrational Hamiltonian. At variance with the Cartesian-based formulation, the kinetic energy operator is not diagonal anymore, implying additional terms arising from the perturbative expansion.

Let us consider a set of M internal coordinates \({\varvec{ s }} = \{s_1,s_2,\dots ,s_M\}\) (\(M\ge N\)). The expression of the kinetic energy operator (KEO) in terms of \({\varvec{ s }}\) is well known and can be stated through the introduction of the Wilson \({\varvec{\textbf{G}}}\) matrix,

$$\begin{aligned} {\varvec{\textbf{G}}} = {\varvec{\textbf{B}}}{\varvec{\textbf{M}}}^{-1}{\varvec{\textbf{B}}}^\text {T} \end{aligned}$$
(16)

where \({\varvec{\textbf{M}}}\) is the diagonal matrix of atomic masses and \({\varvec{\textbf{B}}}\) is built over the coordinates \({\varvec{ s }}\). By introducing \(\tilde{\textsf{G}} = \text {det}({{\varvec{\textbf{G}}}})\), the expression of the operator is [35, 110, 111]

$$\begin{aligned} \mathcal {T} = -\dfrac{\hbar ^2}{2}\sum _{i=1}^M\sum _{j=1}^M\dfrac{\partial }{\partial s_i}\textsf{G}_{ij}\dfrac{\partial }{\partial s_j} + \mathcal {V}_\textsf{g} \end{aligned}$$
(17)

where \(\hbar\) is the reduced Planck constant and

$$\begin{aligned} \mathcal {V}_\textsf{g} = \dfrac{\hbar ^2}{32}\sum _{i=1}^M\sum _{j=1}^M\left[ \dfrac{\textsf{G}_{ij}}{\tilde{\textsf{G}}^2}\dfrac{\partial \tilde{\textsf{G}}}{\partial s_i}\dfrac{\partial \tilde{\textsf{G}}}{\partial s_j} - 4\dfrac{\partial }{\partial s_i}\left( \dfrac{\textsf{G}_{ij}}{\tilde{\textsf{G}}}\dfrac{\partial \tilde{\textsf{G}}}{\partial s_j}\right) \right] \end{aligned}$$
(18)

is an inherently quantum-mechanical contribution, commonly known as extra-potential term. Similarly to the Cartesian-based treatment, the unperturbed eigenstates stem from the harmonic theory of vibrations. Within the present framework, the harmonic frequencies and normal modes are obtained through the Wilson GF method: [112]

$$\begin{aligned} {\varvec{{G}}}{\varvec{{F}}}{\varvec{{L}}} = {\varvec{{L}}}\varvec{\Lambda } \end{aligned}$$
(19)

The diagonal elements of \(\varvec{\Lambda }\) are the squared angular frequencies, the matrix \({\varvec{{L}}}\) contains the eigenvectors linking the normal coordinates \({\varvec{ Q }} = \{Q_1,Q_2,\dots ,Q_N\}\) to the vector \({\varvec{ s }}\),

$$\begin{aligned} {\varvec{ s }} = {\varvec{{L}}}{\varvec{ Q }} \end{aligned}$$
(20)

while the Wilson \({\varvec{{F}}}\) matrix is the Hessian of the potential energy in terms of the internal coordinates, and can be directly calculated from its Cartesian counterpart \({\varvec{{H}}}_x\) (see Appendix A). Since the extra-potential term can be generally approximated with its equilibrium value, it can be neglected in the calculation of transition energies. Following the introduction of the customary dimensionless normal coordinates \({\varvec{ q }}\) and their conjugate momenta \({\varvec{ p }}\), the vibrational Hamiltonian \(\mathcal {H}_v\) can be obtained,

$$\begin{aligned} \mathcal {H}_v = \dfrac{1}{2}\sum _{i=1}^N\sum _{j=1}^Np_i\textsf{g}_{ij}p_j + \mathcal {V} \end{aligned}$$
(21)

where \({\varvec{\textsf{g}}}\) is the \({\varvec{{G}}}\) matrix expressed in wavenumbers and \(\mathcal {V}\) is the potential energy operator. The anharmonic contributions due to the kinetic energies can be evaluated through the Taylor-series expansion of the \({\varvec{{g}}}\) matrix,

$$\begin{aligned} \textsf{g}_{ij} = \textsf{g}^\mathrm {\scriptscriptstyle eq}_{ij} + \sum _{k=1}^N\textsf{g}_{ij,k}q_k + \dfrac{1}{2}\sum _{k=1}^N\sum _{l=1}^N\textsf{g}_{ij,kl}q_kq_l \end{aligned}$$
(22)

where \(\textsf{g}^\mathrm {\scriptscriptstyle eq}_{ij} = \omega _i\delta _{ij}\) and \(\textsf{g}_{ij,kl \dots}\) represent the derivatives of the \({\varvec{\textsf{g}}}\) matrix expressed in wavenumbers:

$$\begin{aligned} \textsf{g}_{ij,kl\dots } = \left( \dfrac{\partial ^n\textsf{g}_{ij}}{\partial q_k\partial q_l\dots }\right) _\mathrm {\scriptscriptstyle eq}\end{aligned}$$
(23)

By inserting Eq. 23 into Eq. 21 and considering the expansion of potential energy, the following expression is obtained:

$$\begin{aligned} \mathcal {H}^{(0)}&= \dfrac{1}{2}\sum _{i=1}^N\omega _i(q^2_i + p^2_i) \end{aligned}$$
(24)
$$\begin{aligned} \mathcal {H}^{(1)}&= \dfrac{1}{6}\sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N\left( \textsf{f}_{ijk}q_iq_jq_k + 3\textsf{g}_{ij,k}p_iq_kp_j\right) \end{aligned}$$
(25)
$$\begin{aligned} \mathcal {H}^{(2)}&= \dfrac{1}{24}\sum _{i=1}^N\sum _{j=1}^N\sum _{k=1}^N\sum _{l=1}^N\left( \textsf{f}_{ijkl}q_iq_jq_kq_l + 6\textsf{g}_{ij,kl}p_iq_kq_lp_j\right) \end{aligned}$$
(26)

Analogously to the treatment in terms of Cartesian-based normal coordinates (CNCs), the anharmonic energies can be obtained by either CVPT or RSPT. Furthermore, thanks to the invariance of the harmonic Hamiltonian, the second-quantization formalism can be still applied. The main difference with respect to the Cartesian-based framework is that each vibrational property \(\mathcal {Q}\) is composed of three contributions, namely a purely potential (\(\mathcal {Q}^\mathrm{\mathcal {V}}\)), a purely kinetic (\(\mathcal {Q}^\mathrm{\mathcal {T}}\)) and a cross term (\(\mathcal {Q}^\mathrm{\times }\)). Within the internal-based VPT2, the energy is provided by the same expression as Eq. 3, where the anharmonic \(\varvec{\chi }\) matrix is expressed as the sum of the three contributions mentioned above:

$$\begin{aligned} \varvec{\chi } = \varvec{\chi }^\mathrm{\mathcal {V}}+ \varvec{\chi }^\mathrm{\mathcal {T}}+ \varvec{\chi }^\mathrm{\times }\end{aligned}$$
(27)

As expected, the purely potential term does not show any difference compared with its Cartesian counterpart but the absence of the Coriolis term and the explicit expression of all contributions are reported in the following,

$$\begin{aligned} 16\chi ^\mathrm{\mathcal {V}}_{ii}&= \textsf{f}_{iiii} - \dfrac{5{\textsf{f}_{iii}}^2}{3\omega _i} - \sum _{\begin{array}{c} j=1\\ \left( j\ne i\right) \end{array}}^N\dfrac{{\textsf{f}_{iij}}^2(8\omega ^2_i-3\omega ^2_j)}{\omega _j(4\omega ^2_i - \omega ^2_j)} \end{aligned}$$
(28)
$$\begin{aligned} 4\chi ^\mathrm{\mathcal {V}}_{ij}= \textsf{f}_{iijj} - \dfrac{2{\textsf{f}_{iij}}^2\omega _i}{4\omega ^2_i-\omega ^2_j} - \dfrac{2{\textsf{f}_{ijj}}^2\omega _j}{4\omega ^2_j-\omega ^2_i} - \dfrac{\textsf{f}_{iii}\textsf{f}_{ijj}}{\omega _i} - \dfrac{\textsf{f}_{iij}\textsf{f}_{jjj}}{\omega _j}\nonumber \\ \quad-\sum _{\begin{array}{c} k=1\\ \left( k\ne i,j\right) \end{array}}^N\left[ \dfrac{2\omega _k(\omega ^2_k-\omega ^2_i-\omega ^2_j){\textsf{f}_{ijk}}^2}{(\omega _i+\omega _j+\omega _k)(\omega _i-\omega _j-\omega _k)(\omega _i-\omega _j+\omega _k)(\omega _i+\omega _j-\omega _k)} + \dfrac{\textsf{f}_{iik}\textsf{f}_{jjk}}{\omega _k}\right] \end{aligned}$$
(29)
$$\begin{aligned} 16\chi ^\mathrm{\mathcal {T}}_{ii}&= 2\textsf{g}_{ii,ii} - \dfrac{3{\textsf{g}_{ii,i}^2}}{\omega _i} \nonumber \\&- \sum _{\begin{array}{c} j=1\\ \left( j\ne i\right) \end{array}}^N\dfrac{{\textsf{g}_{ii,j}}^2(8\omega ^2_i-3\omega ^2_j)-4{\textsf{g}_{ij,i}}^2\omega ^2_j + 8\textsf{g}_{ii,j}\textsf{g}_{ij,i}\omega _i\omega _j}{\omega _j(4\omega ^2_i-\omega ^2_j)} \end{aligned}$$
(30)
$$\begin{aligned} 4\chi ^\mathrm{\mathcal {T}}_{ij}&= \textsf{g}_{ii,jj} + \textsf{g}_{jj,ii}\nonumber \\&- \dfrac{2({\textsf{g}_{ii,j}}^2\omega _i + 4{\textsf{g}_{ij,i}}^2\omega _i-2\textsf{g}_{ii,j}\textsf{g}_{ij,i}\omega _j)}{4\omega ^2_i-\omega ^2_j}\nonumber \\&- \dfrac{2({\textsf{g}_{jj,i}}^2\omega _j + 4{\textsf{g}_{ij,j}}^2\omega _j-2\textsf{g}_{jj,i}\textsf{g}_{ij,j}\omega _i)}{4\omega ^2_j-\omega ^2_i}\nonumber \\&- \dfrac{\textsf{g}_{ii,i}\textsf{g}_{jj,i}}{\omega _i}- \dfrac{\textsf{g}_{ii,j}\textsf{g}_{jj,j}}{\omega _j} \end{aligned}$$
(31)
$$\begin{aligned}&-\sum _{\begin{array}{c} k=1\\ \left( k\ne i,j\right) \end{array}}^N\Biggl [\dfrac{2({\textsf{g}_{ij,k}}^2+{\textsf{g}_{ik,j}}^2+{\textsf{g}_{jk,i}}^2)\omega _k(\omega ^2_k-\omega ^2_i-\omega ^2_j)+ 4\textsf{g}_{ij,k}\textsf{g}_{ik,j}\omega _j(\omega ^2_j-\omega ^2_i-\omega ^2_k)}{(\omega _i+\omega _j+\omega _k)(\omega _i-\omega _j-\omega _k)(\omega _i-\omega _j+\omega _k)(\omega _i+\omega _j-\omega _k)}\nonumber \\&+ \dfrac{4\textsf{g}_{ij,k}\textsf{g}_{jk,i}\omega _i(\omega ^2_i-\omega ^2_j-\omega ^2_k)+ 8\omega _i\omega _j\omega _k \textsf{g}_{ik,j}\textsf{g}_{jk,i}}{(\omega _i+\omega _j+\omega _k)(\omega _i-\omega _j-\omega _k)(\omega _i-\omega _j+\omega _k)(\omega _i+\omega _j-\omega _k)}\nonumber \\&+\dfrac{\textsf{g}_{ii,k}\textsf{g}_{jj,k}}{\omega _k}\Biggl ]\nonumber \\ 16\chi ^\mathrm{\times }_{ii}&= -\dfrac{2\textsf{g}_{ii,i}\textsf{f}_{iii}}{\omega _i}\nonumber \\&+\sum _{\begin{array}{c} j=1\\ \left( j\ne i\right) \end{array}}^N\dfrac{8\textsf{g}_{ij,i}\textsf{f}_{iij}\omega _i\omega _j-2\textsf{g}_{ii,j}\textsf{f}_{iij}(8\omega ^2_i-\omega ^2_j)}{\omega _j(4\omega ^2_i-\omega ^2_j)} \end{aligned}$$
(32)
$$\begin{aligned} 4\chi ^\mathrm{\times }_{ij}&= \dfrac{\textsf{g}_{ii,i}\textsf{f}_{ijj}+\textsf{g}_{jj,i}\textsf{f}_{iii}}{\omega _i} + \dfrac{\textsf{f}_{iij}\textsf{g}_{jj,j}+\textsf{g}_{ii,j}\textsf{f}_{jjj}}{\omega _j}\nonumber \\&- \dfrac{4(\textsf{g}_{ij,i}\textsf{f}_{iij}\omega _j-\textsf{g}_{ii,j}\textsf{f}_{iij}\omega _i)}{4\omega ^2_i-\omega ^2_j}\nonumber \\&- \dfrac{4(\textsf{g}_{ij,j}\textsf{f}_{ijj}\omega _i-\textsf{g}_{jj,i}\textsf{f}_{ijj}\omega _j)}{4\omega ^2_j-\omega ^2_i}\nonumber \\&+\sum _{\begin{array}{c} k=1\\ \left( k\ne i,j\right) \end{array}}^N\Biggl [\dfrac{4\textsf{g}_{ik,j}\textsf{f}_{ijk}\omega _{i}(\omega ^2_{i} - \omega ^2_{j} - \omega ^2_{k}) + 4\textsf{g}_{jk,i}\textsf{f}_{ijk}\omega _{j}(\omega ^2_{j} - \omega ^2_{i} - \omega ^2_{k})}{(\omega _i+\omega _j+\omega _k)(\omega _i-\omega _j-\omega _k)(\omega _i-\omega _j+\omega _k)(\omega _i+\omega _j-\omega _k)}\nonumber \\&+ \dfrac{8\omega _i\omega _j\omega _k \textsf{g}_{ij,k}\textsf{f}_{ijk}}{(\omega _i+\omega _j+\omega _k)(\omega _i-\omega _j-\omega _k)(\omega _i-\omega _j+\omega _k)(\omega _i+\omega _j-\omega _k)}\nonumber \\&-\dfrac{\textsf{g}_{ii,k}\textsf{f}_{jjk}+\textsf{f}_{iik}\textsf{g}_{jj,k}}{\omega _k}\Biggl ] \end{aligned}$$
(33)

where few misprints of Eqs. 35 and 37 of Ref. [73] have been corrected. By applying the partial fraction decomposition and introducing the tensors}

$$\begin{aligned} \eta _{ijkl}&= \textsf{f}_{ijkl} + \textsf{g}_{ij,kl} + \textsf{g}_{kl,ij} \end{aligned}$$
(34)
$$\begin{aligned} \sigma _{ijk}&= \textsf{f}_{ijk} - (\textsf{g}_{ij,k} + \textsf{g}_{ik,j} + \textsf{g}_{jk,i}) \end{aligned}$$
(35)
$$\begin{aligned} \rho _{ijk}&= \textsf{f}_{ijk} - (\textsf{g}_{ij,k} - \textsf{g}_{ik,j} - \textsf{g}_{jk,i}) \end{aligned}$$
(36)

the elements of the \({\varvec{{\chi }}}\) matrix can be rewritten in a more compact form, which is also required for the application of the DVPT2 scheme:

$$\begin{aligned} 16\chi _{ii}&= \eta _{iiii} - \dfrac{\sigma _{iii}^2/2 + 9\rho _{iii}^2/2}{3\omega _i} \nonumber \\&\quad - \dfrac{1}{2}\sum _{\begin{array}{c} j=1\\ \left( j\ne i\right) \end{array}}^N\left[ \dfrac{4\rho _{jii}^2}{\omega _j}+\dfrac{\sigma _{iij}^2}{2\omega _i+\omega _j}-\dfrac{\rho _{iij}^2}{2\omega _i-\omega _j}\right] \end{aligned}$$
(37)
$$\begin{aligned} 4\chi _{ij}&= \eta _{iijj} -\dfrac{1}{2}\left[ \dfrac{\sigma _{iij}^2}{2\omega _i+\omega _j}+\dfrac{ \rho _{iij}^2}{2\omega _i-\omega _j}\right] \nonumber \\&\quad -\dfrac{1}{2}\left[ \dfrac{\sigma _{jji}^2}{2\omega _j+\omega _i}+\dfrac{\rho _{jji}^2}{2\omega _j-\omega _i}\right] -\dfrac{\rho _{iii}\rho _{ijj}}{\omega _i}\nonumber \\&\quad -\dfrac{\rho _{jjj}\rho _{jii}}{\omega _j}-\dfrac{1}{2}\sum _{\begin{array}{c} k=1\\ \left( k\ne i,j\right) \end{array}}^N\Biggl [\dfrac{\sigma _{ijk}^2}{\omega _i+\omega _j+\omega _k}\nonumber \\&\quad -\dfrac{ \rho _{ijk}^2}{\omega _i+\omega _j-\omega _k}+\dfrac{\rho _{ikj}^2}{\omega _i-\omega _j+\omega _k}\nonumber \\&\quad -\dfrac{\rho _{jki}^2}{\omega _i-\omega _j-\omega _k}+\dfrac{2\rho _{kii}\rho _{kjj}}{\omega _k}\Biggl ] \end{aligned}$$
(38)

One of the most interesting aspects of Eqs. 37 and 38 is their analogy with the Cartesian-based treatment from the algebraic point of view. In the first place, the internal-based elements of \({\varvec{{\chi }}}\) show the same functional form as those introduced in Eqs. 4 and 5. Consequently, the internal-based \({\varvec{{\chi }}}\) matrix can be interpreted as a full-fledged generalization of the corresponding Cartesian counterpart, implying a remarkable simplification at the implementation level. In fact, the extension of any code implementing the Cartesian-based VPT2 becomes straightforward. Moreover, the diagnostic of Fermi resonances can be carried out through a direct extension of the Martin test. While the first step in the identification of FRs (see Eq. 6) does not require further changes, the second step becomes a straightforward generalization of Eq. 7:

$$\begin{aligned} \dfrac{{\rho _{jki}}^4}{64(1+\delta _{jk})^2|\omega _i - \omega _j - \omega _k|^3} \ge K^{1-2} \end{aligned}$$
(39)

It is evident that Eq. 7 can be generated from Eq. 39 by setting the derivatives of the \({\varvec{{g}}}\) matrix to zero, since in that case \(\rho _{ijk} = \textsf{f}_{ijk}\).

Once the set of FRs has been identified, the corresponding interaction elements of the contact-transformed Hamiltonian can be calculated,

$$\begin{aligned} \langle {\varvec{ v }}_{\mathrm{R}} + 1_i |{\tilde{\mathcal {H}}}|{\varvec{ v }}_{\mathrm{R}} + 1_j + 1_k\rangle = \mathsf {\rho }_{jki}\sqrt{\dfrac{( v _{{{\scriptscriptstyle R},i}}+1)( v _{{{\scriptscriptstyle R},j}}+1)( v _{{{\scriptscriptstyle R},k}}+1+\delta _{jk})}{8(1+3\delta _{jk})}} \end{aligned}$$
(40)

paving the route for the extension of the GVPT2 scheme to the use of internal coordinates. Let us underline that thanks to the analogy between the Cartesian- and internal-based frameworks of VPT2 in terms of the analytical expressions of wave functions, vibrational energies and diagnostic of FRs, the extension of the internal-based VPT2 for asymmetric tops to the treatment of linear, symmetric and spherical tops can be carried out through the procedure described in Ref. [65].

3 Implementation

In this section, the implementation of the algorithms employed in this work is briefly addressed, while a deeper discussion is reported elsewhere [65, 97].

3.1 Internal-based VPT2 framework

A full calculation within the internal-based VPT2 framework is carried out through different steps and exploits an interplay of implementations within a recently devised standalone code and a development version of the quantum-chemistry package Gaussian. In particular, the standalone program performs the harmonic vibrational analysis, the generation of the single-point Gaussian  input files required for the calculation of both potential and kinetic energy derivatives with respect to the normal coordinates by finite differences, and the calculation of such derivatives. On the other side, the tasks carried out through the Gaussian package are the initial optimization, the calculation of Cartesian force constants at both equilibrium and out-of-equilibrium geometries, and the application of the VPT2 framework for the calculation of the anharmonic frequencies. Summarizing, the new program is mainly aimed at generating the anharmonic force field and the Wilson \({\varvec{{g}}}\) derivatives, while the internal-based VPT2 framework has been implemented in the Gaussian package. Even though different quantum-chemistry programs can be employed for the generation of the derivatives required in the anharmonic treatment, in the present work the G16 [113] package has been always employed.

3.2 Calculation of semi-experimental structures

As previously discussed, the main ingredients for the calculation of SE structures are the experimental rotational constants, vibrational corrections, weights and the guess geometry. In general both the guess geometry and vibrational corrections are evaluated through a quantum-chemistry package. Similarly to the application of the VPT2 framework in internal coordinates, the G16 package has been used for this purpose. Once these preliminary operations have been carried out, the whole set of data is used to start the structural refinement by the MSR program, leading to the characterization of the molecular geometry.

Since the new code and MSR share a common set of libraries, they have been recently assembled in a single suit of programs for vibro-rotational analyses.

4 Computational details

Most of the available electronic structure methods implementing analytical second-order derivatives of potential energy have been employed. The hybrid density functional B3PW91 [114] has been used in conjunction with the jul-cc-pVDZ (hereafter julDZ) basis set [115]. Furthermore, tight d functions have been included (julDZd) for the treatment of atoms belonging to the third period. The double-hybrid functional revDSD-PBEP86 [116] and second-order Møller–Plesset perturbation theory (MP2) [117] have been employed in conjunction with the jun-cc-pVTZ (hereafter junTZ) basis set [115]. The latter has been augumented with tight d functions (junTZd) for calculations involving third-period atoms. The empirical dispersion contributions have been systematically accounted for in density functional theory (DFT) computations by means of Grimme’s D3 model with Becke–Johnson damping [118, 119].

5 Results and discussion

Several systems have been analyzed in detail with the objective of validating the theoretical framework developed in the fields of both molecular structure prediction and vibrational spectroscopy. In this section, selected test cases will be discussed to highlight the main aspects characterizing the computational protocols developed so far.

5.1 Determination of SE equilibrium structures

As previously outlined, the choice of the set of coordinates represents one of the main steps needed to obtain reliable SE structures. While the effectiveness of ZICs for this type of calculations is now well recognized, in this context the attention will be focused on symmetry-based optimizations. With the aim of investigating different symmetries, the simplest Criegee intermediate (point group C\(_s\)) and thiophene (C\(_{2v}\)) have been selected as case studies. Furthermore, for both systems detailed analyses are available in the literature, which provided highly accurate SE structures. For this reason, they represent ideal test cases to validate the computational protocol discussed in the present work.

First, let us consider the Criegee intermediate (see Fig. 1), which is a planar system fully described by 7 degrees of freedom.

Fig. 1
figure 1

Molecular structure and atom labeling of the simplest Criegee intermediate

Full size image

A detailed analysis has been performed by McCarthy and co-workers [120], where the experimental rotational constants of nine isotopic species combined with vibrational corrections at the CCSD(T)/ANO1 [121] level of theory have been employed in the prediction of the SE structure. In the present work, the same set of experimental data has been used for the calculation of both the effective (\(r_0\)) structure and, in conjunction with vibrational corrections at the rDSD/junTZ level, its SE counterpart. However, due to the planarity of the system only two of the three moments of inertia are independent. As a matter of fact, 18 rotational constants have been actually included in the optimization. The fit has been performed employing both ZICs (see Section S1 of the Supplementary Information) and A\(_1\)-DICs as nuclear coordinates. Concerning the latter, a set of PICs has been generated from the molecular connectivity, followed by the calculation of DICs and the extraction of the totally-symmetric coordinates. As expected from a C\(_s\) planar systems, the A\(_1\)-DICs correspond to those DICs which do not alter dihedral angles. Since A\(_1\)-DICs are much less chemically intuitive if compared with ZICs, an error propagation has been applied to express both SE geometries in terms of the latter through a double-step procedure. In the first place, the variance–covariance matrix expressed in terms of A\(_1\)-DICs (\({\varvec{{\Sigma }}}^{}_{\text {A}_1}\)) has been converted to the corresponding Cartesian-based counterpart (\({\varvec{{\Sigma }}}^{}_{\text {x}}\)) through the following expression,

$$\begin{aligned} {\varvec{{\Sigma }}}^{}_\text {x} = \lbrace {\varvec{{B}}} _{\text {A}_1} \rbrace^\text {T} {\varvec{{\Sigma }}}^{}_{\text {A}_1}\lbrace {\varvec{{B}}} _{\text {A}_1}\rbrace \end{aligned}$$
(41)

where \({\varvec{{B}}}^\dagger _{\text {A}_1}\) represents the pseudo-inverse of the Wilson \({\varvec{{B}}}\) matrix expressed in terms of A\(_1\)-DICs evaluated at the optimized geometry. Secondly, the \({\varvec{{\Sigma }}}^{}_{\text {x}}\) matrix has been converted to ZICs,

$$\begin{aligned} {\varvec{{\Sigma }}}^{}_\text {Z} = \lbrace {\varvec{{B}}}^ \dagger_{\text {Z}}\rbrace ^\text{T} {\varvec{{\Sigma }}}^{}_{\text {x}}\lbrace {\varvec{{B}}}^\dagger_{\text {Z}}\rbrace \end{aligned}$$
(42)

where \({\varvec{{B}}}_\text {Z}\) is the Wilson \({\varvec{{B}}}\) matrix in terms of ZICs at the optimized geometry. The new standard deviations and confidence intervals have been obtained starting from the square roots of the diagonal elements of \({\varvec{{\Sigma }}}^{}_\text {Z}\). The different sets of structural parameters are reported and compared with those proposed by McCarthy and co-workers in Table 1.

Table 1 Equilibrium molecular structure of the simplest Criegee intermediate (distances in Å, angles in degrees)
Full size table

Inspection of Table 1 shows that that the SE geometries obtained through ZICs and A\(_1\)-DICs exactly coincide both in terms of structural parameters and standard deviation and, what is even more important they present also structural parameters close to their reference counterparts. Consequently, the symmetry-based approach has been applied without the setup of a user-defined Z-matrix, but still reaching the same accuracy.

A further validation of the proposed methodology has been performed by the structure determination of thiophene (see Fig. 2), a C\(_{2v}\) heteroaromatic cycle.

Fig. 2
figure 2

Molecular structure and atom labeling of thiophene

Full size image

A first characterization of this system performed by our research group employed the experimental rotational constants of eight isotopologues combined with vibrational corrections based on the B2PLYP and B3LYP functionals in conjunction with the cc-pVTZ and SNSD basis sets, respectively. [122] More recently, a detailed experimental and theoretical investigation of thiophene has been carried out by Orr and co-workers [123]. In particular, the set of experimental data has been extended to include 26 isotopic species, with the latter being used together with vibrational and electronic corrections at the CCSD(T)/cc-pCPVTZ level of theory for the structural refinement. In agreement with the calculations reported in Ref. [123], 24 isotopic species have been included in the present work. Similarly to the Criegee intermediate, only two rotational constants for each isotopologue have been considered, leading to a set of 48 experimental rotational constants usable in the nonlinear regression. At variance with Ref. [123], only vibrational corrections at the B3P/julDZd level have been used to correct the experimental data, as well as to compute the guess structure for the optimization. A full comparison of the effective and SE structure with those proposed in Ref. [123] is reported in Table 2.

Table 2 Equilibrium molecular structure of thiophene (distances in Å, angles in degrees)
Full size table

In analogy with the Criegee intermediate, both the SE equilibrium structure and the standard deviation are invariant under change of coordinates. The search of totally symmetric coordinates provided the correct number of coordinates (8), and their use reaches the same structure obtained by using the Z-matrix indicated in Ref. [123] and reported in Section S1 of the Supplementary Information. Last but not least, a good agreement between the SE structural parameters retrieved in this work and those proposed by Orr and co-workers has been detected despite the different nature of the corrections to the experimental rotational constants.

5.2 Anharmonic calculations in internal coordinates

With the aim of highlighting the advantages of an internal-based formulation of VPT2, our recently developed engine has been be applied to both semi-rigid and flexible systems. In the former case, the choice of the coordinates set plays a more marginal role even though it leads to a remarkable reduction of inter-mode couplings. This finding paves the way to the systematic use of cost-effective protocols aimed at calculating vibrational properties. In fact, the reduction of couplings between different vibrations results in the calculation of low-dimensionality anharmonic force fields, which represent the bottleneck of an anharmonic calculation. This effect is even stronger when flexible systems are taken into account, especially concerning the couplings involving LAMs.

In the present work, 1,1-difluoroethylene (see Fig. 3) has been selected as test-case for semi-rigid systems.

Fig. 3
figure 3

Molecular structure of 1,1-difluoroethylene

Full size image

The calculation of the fundamental vibrational frequencies at the anharmonic level has been performed by employing all VPT2 schemes discussed in the previous sections within both internal- and Cartesian-based frameworks. A full comparison of the results with the experimental data is reported in Table 3.

Table 3 Comparison of the Cartesian and curvilinear VPT2, DVPT2, and GVPT2 wavenumbers (in cm\(^{-1}\)) of 1,1-difluoroethylene at the MP2/junTZ level with the experimental data
Full size table

As expected, the Cartesian- and internal-based calculations do not show any difference in the resulting transition frequencies. It is also worth specifying that Coriolis contributions have been properly considered in the calculations based on Cartesian coordinates. As a matter of fact, specific second-order terms in the expansion of the KEO yield contributions formally equivalent to the Coriolis terms. Despite the GVPT2 results obtained on the basis of Cartesian and internal coordinates are numerically indistinguishable, the magnitude of the couplings is remarkably affected by the choice of the nuclear coordinates. This pattern has been also observed for terms which do not contribute directly to the calculation of anharmonic transition energies in this context, such as three-mode quartic force constants.

Fig. 4
figure 4

Comparison of the number of cubic (\(\textsf{f}_{ijk}\) (\(i\ne j\ne k\))) and quartic (\(\textsf{f}_{ijkk}\) (\(i\ne j\))) force constants of 1,1-difluoroethylene above a given threshold (in cm\(^{-1}\)) computed at the MP2/junTZ level with Cartesian or curvilinear coordinates

Full size image

Inspection of Fig. 4 shows that the magnitude of inter-mode terms is definitively reduced when switching from the Cartesian-based representation of normal coordinates to that employing curvilinear coordinates. Furthermore, vibrational couplings are not transferred from potential to kinetic energy, since there exist only six first-order derivatives of the Wilson \({\varvec{\textsf{g}}}\) matrix exceeding 100 cm\(^{-1}\), while all second-order derivatives are well below this value.

The AAT conformer of glycolic acid (see Fig. 5) has been selected as a flexible system presenting LAMs.

Fig. 5
figure 5

Molecular structure of the AAT conformer of glycolic acid

Full size image

A recent detailed analysis of the conformers of glycolic acid [126] shows that most of them are characterized by the presence of this kind of vibrations, poorly described at the VPT2 level. In this work, the AAT conformer has been considered, since it is characterized by two distinct LAMs, representing then an ideal system for studying also inter-mode couplings between different LAMs. The simulations have been performed employing a dual-level method based on the substitution approach [73], where the anharmonic calculation at the B3P/julDZ level of theory has been refined by the inclusion of harmonic frequencies evaluated at the rDSD/junTZ level (hereafter rDSD/junTZ//B3P/julDZ). Furthermore, the reduced-dimensionality scheme has been used, so that anharmonic corrections have been applied to all vibrations but LAMs. The set of fundamental wavenumbers obtained within both the Cartesian- and internal-based GVPT2 framework are compared with their experimental counterparts in Table 4.

Table 4 Comparison of the Cartesian and curvilinear GVPT2 fundamental wavenumbers (in cm\(^{-1}\)) of the AAT conformer of glycolic acid with the experimental data
Full size table

The data reported in Table 4 show that both sets of GVPT2 transition frequencies are in remarkable agreement with their experimental counterparts, with a mean absolute error (MAE) always lower than 10 cm\(^{-1}\). Despite the similarity of the results for Cartesian and internal coordinates, the behavior of the inter-mode couplings is very different, especially when interaction terms involving LAMs are examined. In order to highlight this aspect, the contribution of LAMs has been investigated by comparing the quartic force constants involving modes 20 and 21, which are reported in Fig. 6.

Fig. 6
figure 6

Comparison of the Cartesian (top panel) and curvilinear (bottom panel) quartic force constants of the AAT conformer of glycolic acid involving modes 20 and 21 at the rDSD/junTZ//B3P/julDZ level of theory

Full size image

It is apparent that both diagonal and off-diagonal terms are much larger in Cartesian coordinates (exceeding 80,000 cm\(^{-1}\) for \(\textsf{f}_{21,21,21,21}\)), while the values are definitely more reasonable for the internal-based counterparts. Interestingly, this behavior also involves couplings between LAMs and stretching modes. As an example, the values the semi-diagonal quartic force constants \(\textsf{f}_{20,20,1,1}\) and \(\textsf{f}_{21,21,1,1}\) expressed in CNCs are, respectively, −3733 and −9347 cm\(^{-1}\), while the corresponding internal-based counterparts are −27 and −62 cm\(^{-1}\). Hence, the internal-based representation of vibrations allows not only for an improved separation of LAMs from SAMs, but it also lays the foundations for a multi-mode, one-dimensional treatment of LAMs.

6 Conclusion

In this work, I have presented a general engine for vibrational and rotational spectroscopy based on curvilinear internal coordinates.

Concerning the determination of semi-experimental molecular structures, the new methodology provides a transparent way for the a priori definition of the nuclear coordinates required in the optimization process and a unambiguous definition of geometrical constraints without loss of accuracy. The algorithm implemented into the MSR software enables SE nonlinear regression starting from a minimal set of input data without any intermediate, user-dependent step. The applications presented in this context confirmed our initial hypotheses and previous studies, paving the route for the study of larger systems.

The outcomes of a number of test cases demonstrate the reliability of the new GVPT2 engine based on curvilinear internal coordinates and allowing to effectively handle medium-to-large-size molecules. It is worth mentioning that this methodology is completely general for all those electronic structure methods implementing analytical Cartesian force constants. Comparison with the standard Cartesian-based framework shows two main benefits. The first one is the ease of implementation from an existing code based on the Cartesian formulation of VPT2. The second one is a significant improvement in the separation of large amplitude motions and their further treatment by more appropriate methodologies. In fact, a significant reduction in the inherent problems of VPT2 applied to this kind of systems has been verified at different levels. More sophisticated strategies for treating LAMs are currently under study in our research group, with the aim to study molecular systems featured by a growing number of this type of vibrational motions.