Interplay of Rotational and Pseudorotational Motions in Flexible Cyclic Molecules

Solutions to the time-independent nuclear Schrödinger equation associated with the pseudorotational motion of three flexible cyclic molecules are presented and discussed. Structural relaxations related to the pseudorotational motion are described as functions of a pseudorotation angle ϕ which is formulated according to the definition of ring-puckering coordinates originally proposed by Cremer and Pople (J. Am. Chem. Soc.1975, 97 ( (6), ), 1354−1358). In order to take into account the interplay between pseudorotational and rotational motions, the rovibrational Hamiltonian matrices are formulated for the rotational quantum numbers J = 0 and J = 1. The rovibrational Hamiltonian matrices are constructed and diagonalized using a Python program developed by the authors. Suitable algorithms for (i) the construction of one-dimensional cuts of potential energy surfaces along the pseudorotation angle ϕ and (ii) the assignment of the vibrorotational wave functions (which are needed for the automatic calculation of rotational transition energies J = 0 → J = 1) are described and discussed.


INTRODUCTION
The relevance of intramolecular nuclear motion effects on the physical chemical properties of isolated cyclic molecules is well-known. 1−3 A confirmation of the importance of this phenomenon is provided by the analysis of the experimental spectra of cyclic molecular systems. Some of these experimental data were rationalized as consequences of ring puckering (RP) motions (for example, see refs 4−8), which in general are not easily described in terms of Cartesian or primitive internal coordinates (i.e., in terms of bond lengths or dihedral/valence angle variations) of nuclei. Many efforts 9−12 were devoted to the formulation of a curvilinear coordinate system suitable for a general description of RP motions in cyclic molecular systems.
Nowadays, RP motions of N-membered rings (with N equal to or greater than 5) can be elegantly described with Cremer− Pople coordinates (CPCs). 12 With the definitions proposed by Cremer and Pople, the RP of a cyclic moiety comprising N atoms (the substituents are not taken into account in this number, i.e. N = 5 atoms pertains to the cyclic moiety of cyclopentane) can be described with N − 3 CPCs. After the introduction of ring deformation coordinates (RDCs) by Zou, Izotov, and Cremer, 13 the conformation of a cyclic moiety with N atoms can be completely specified employing N − 3 CPCs and 2N − 3 RDCs.
In the special case of a 5-term ring system, RP can be described with 2 CPCs and visualized with a two-dimensional cut of the global potential energy surface (2D-PES). 14 An analysis of the RP motions based on the separation of the RP motions from the other vibrational motions of the molecular system of interest is an approximation. Nevertheless, this approximation is often reliable and useful. In the framework of the Born−Oppenheimer (BO) approximation, the reduction in dimensionality achieved for the PES of a 5-term ring system by employing CPCs can be extended to the solution of the timeindependent nuclear Schrodinger equation (NSE). In other words, the 2D-PES which describes the RP of a 5-term ring system can be employed in the formulation of the NSE associated with the RP of the same system. The picture can be further simplified: as observed in previous works, 15 the minimum-energy path associated with RP motions of 5-term saturated ring systems can be often described with a single CPC, namely, the pseudorotation angle ϕ.
In this work, the solutions of the NSE associated with the Cremer−Pople pseudorotation angle ϕ for the three molecular systems shown in Figure 1 are presented and discussed. First, an algorithm suitable for the construction of a one-dimensional cut of the global PES (1D-PES) along the pseudorotation angle ϕ was devised, implemented, and employed (in order to derive the potential energy term of the nuclear Hamiltonian).
Second, the eigenvalue problem related to the nuclear Hamiltonian was solved with a computational strategy based on the employment of a monodimensional discrete variable representation (DVR). 16 More specifically, the strategy originally proposed by Meyer in refs 17 and 18 was employed.
The solutions of the nuclear problem related to the pseudorotational pathway are provided for two values of the rotational quantum number J; that is, the nuclear Hamiltonian was constructed for J = 0 (pure psuedorotational case) and J = 1. In this way, the interplay between pseudorotational and rotational motions can be numerically evaluated and discussed.
The model employed in this work to reduce the full dimensional NSE to a simpler 1D problem is named f lexible model by Meyer 19 and semirigid model 20 or adiabatic-constraint approach 21,22 by other authors. There are other approaches to dimensionality reduction of the nuclear problem (a list of useful references and a brief historical account concerning development and employment of approximated separation methods for the description of molecular motions can be found in the introduction of chapter 12 of ref 23). In particular, the employment of a BO-like separation between small amplitude motions (SAMs) and large amplitude motions (LAMs) was proposed: 24−26 in this approach, the formulation of an effective nuclear Hamiltonian for the LAM of interest is used to take into account the effects of SAMs (completely neglected in this work) in the solutions of the nuclear problem. Furthermore, higher-dimensional formulation of the nuclear problem can be achieved: in the original article by Meyer, 19 an extension to the 2D problem was proposed, and other extensions suitable for the solution of more complex formulations of the nuclear problems (J > 1 and more than two coupled internal coordinates) are available in literature. 27−29 Nowadays, even the solution of the full dimensional problem in curvilinear coordinates is possible, 30 at least in principle (in practice, to formulate the full dimensional problem the calculation of the global PES must be feasible). An up-to-date picture of the currently available nonperturbative approaches to the formulation of a nuclear Hamiltonian and to the solution of the final eigenvalue problem can be found in ref 31. This article is structured as follows. In section 2, the definition of CPCs is introduced and the method employed for the solution of the NSE associated with the pseudorotation angle ϕ is discussed. In section 3, the computational protocol is outlined. In section 4, computational results are discussed: eigenvalues and eigenvectors of the nuclear Hamiltonians are provided and commented on for each of the compounds shown in Figure 1. In section 5, accuracy, limitations, and possible improvements of the model are discussed. Finally, conclusions are outlined in section 6.

THEORETICAL BACKGROUND
As mentioned in the Introduction, RP motions of flexible cyclic molecular systems can be described with CPCs. The definition of CPCs was proposed for the first time in ref 12. In what follows, the definition of CPCs is provided for the special case of a 5-term ring system. In this special case, only 2 CPCs are required: q (the puckering amplitude) and ϕ (the pseudorotation angle). To define q and ϕ, the Cartesian framework employed to specify the nuclear positions of the molecular system must be translated and rotated. More specifically, the origin of the Cartesian framework must coincide with the geometrical center of the 5 nuclei directly involved in the ring structure: 12 where r j specifies the position of the jth atomic nucleus; then, the z axis must be oriented perpendicularly to a mean plane. The position vectors r′ and r″ can be employed to fix the orientation of the mean plane and the orientation of the zv ersor, which is perpendicular to the mean plane: 12 In a Cartesian framework translated and oriented according to eqs 1−4, the coordinate z j of the jth atomic nucleus is defined as follows: The length q and the pseudorotation angle ϕ are defined in the following manner: 12,32,33   In this work, the conformation of a 5-term ring system is specified through the CPCs q and ϕ. In order to establish a univocal correspondence between the value of the pseudorotation angle ϕ and the conformation of a 5-term ring system, the sequential numbering of the 5 atoms directly involved in the ring structure must be specified (see Figure 2). In other words, each atom is labeled with a specific j value (from 1 to 5), and two consecutive j values must be assigned to adjacent atoms of the ring.
A graphical representation of the correspondence between a couple of values (q, ϕ) and the conformation of a 5-term ring system can be found in Figure 3. The planar conformation (q = 0) corresponds to the center of the circle.
Other coordinate systems were employed by other authors to explore 2D-PES of 5-term ring systems. 34,35 For example, the same conformational space spanned by q and ϕ can be explored defining another two coordinates, q cos ϕ and q sin ϕ (sometimes called ring-twisting and ring-bending coordinates). However, in many cases the minimum-energy path can be described as a one-dimensional pseudorotation pathway: in these cases, a well-defined pseudorotation angle ϕ can be usefully employed to construct a 1D-PES. This simple idea was recognized before the pivotal article of Cremer and Pople 12 (for example, see ref 9), but the general definition of the pseudorotation angle proposed in ref 12 is the starting point for a systematic exploration of pseudorotational pathways in 5term ring systems. Nowadays, the availability of analytic gradients of CPCs 32 allows the construction of a 1D-PES along the pseudorotation angle ϕ defined in eq 7 (for example, see ref 36). In this article, constrained geometry optimizations were performed at different values of ϕ. The constraints were automatically generated with the algorithm proposed in section 1.1 of the Supporting Information: during each optimization, the ϕ value was fixed while all the other nuclear degrees of freedom (including q) were optimized. In this manner, the points of a 1D-PES were calculated.
As mentioned above, dimensionality reduction allows for a formulation of the potential energy term V̂in terms of a limited number of coordinates (typically one or two, i.e., construction of 1D-PES or 2D-PES). This is the starting point for the construction of a nuclear Hamiltonian Ĥemployed to describe nuclear motions in terms of a reduced coordinate space. 37−39 This is an approximated approach: in general, a molecular system is subject to many nuclear motions which cannot be described with a reduced set of one or two curvilinear nuclear coordinates. However, the coupling is often negligible, and a separation between LAMs and SAMs can be introduced. In this work, the pseudorotation of a 5-term ring system is described as a LAM, and its separability from the other internal nuclear degrees of freedom is assumed. In other words: if the number of nuclei in the molecular system is equal to K, the coupling between the pseudorotational motion and the other 3K − 7 internal degrees of freedom is neglected. The problem is therefore reduced to the solution of a onedimensional NSE (1D-NSE) associated with the pseudorotational motion: The nuclear Hamiltonian Ĥcan be partitioned as follows: =̂+Ĥ V T (9) where V̂is the potential energy term and T̂the kinetic energy operator (KEO). The employment of a curvilinear coordinate system is often advantageous for the exploration of the conformational space of a molecular system (i.e., for the construction of V̂). On the other hand, the analytical expression of the KEO in curvilinear coordinates 40−42 can be extremely complicated. In other words, while a simple analytical formulation of the KEO in normal coordinates is well-known, 43−45 attempts to derive an analytical expression for the KEO in curvilinear coordinates led to complex and lengthy formulations. 46−48 Some authors proposed an alternative approach, 49 which was adopted also for this work: the numerical computation of the KEO in curvilinear coordinates. Another possibility is the employment of computer procedures based on symbolic calculations in order to achieve an automatic derivation of an analytical formulation of the KEO. 50,51 The pragmatic approach chosen for this work was proposed and developed by  In particular, the computational strategy proposed in ref 19 was adopted. The construction of the nuclear Hamiltonian matrices is detailed in section 2 of the Supporting Information. In order to take into account the interaction between pseudorotational and rotational motions, two formulations of the nuclear Hamiltonian were employed. The first formulation is suitable for the solution of a pure vibrational problem (i.e., the 1D-NSE related to the pure pseudorotational motion), which is the case of J = 0: where g αβ and g αϕ are elements of the kinetic energy matrix and Pα is a component of the overall angular momentum with respect to molecule-fixed axis.
The 1D-NSE can be viewed as an eigenvalue problem: the eigenvectors are the wave functions (ψ) and the eigenvalues are the energies (E).

COMPUTATIONAL PROTOCOL
Electronic calculations and constrained geometry optimizations were performed using the Gaussian 16 suite of programs. 53 Electronic calculations were carried out with density functional theory (DFT) by employing B3LYP 54−56 as an exchange−correlation functional and maug-cc-pVTZ basis set. 57,58 Each point of the 1D-PES was calculated independently through a constrained optimization. As in ref 33, input geometries were constructed using both Cartesian and primitive internal coordinates (PICs, i.e., bond lengths and dihedral and valence angles) in the same Z-matrix. The positions of the 5 atoms directly involved in the ring structure (e.g., carbon and oxygen atoms in the case of 1,2-dioxolane) were specified with Cartesian coordinates (translated and oriented according to eqs 1−4), while the positions of the other atoms were provided in terms of PICs (e.g., hydrogen atoms in the case of 1,2-dioxolane). In the case of ref 33, both q and ϕ were fixed during the constrained optimization (i.e., 2D-PESs were constructed). In this work, only ϕ is fixed, while q is optimized: therefore, the homemade Python script mentioned in section 3 of ref 33 and employed to automatically write input files corresponding to a specific couple (q, ϕ) was modified accordingly. In practice, the algorithm discussed in section 1.1 of the Supporting Information was devised and implemented. In contrast to the protocol of ref 33 (z j Cartesian coordinates of the 5 atoms directly involved in the ring structure are kept fixed), the following 4 constraints are imposed during the geometry optimization: three ratios z j /z j′ involving 4 of the 5 z j Cartesian coordinates of the atoms directly involved in the ring structure (the z j coordinate with the lower absolute value is excluded) are kept fixed, and the sum of the 5 z j Cartesian coordinates is constrained to 0, according to eq 1. For more details, see section 1.1 of the Supporting Information.
For each of the molecular systems investigated in this work, 360 constrained geometry optimizations were carried out (each one corresponding to a specific value of the pseudorotation angle ϕ, which is kept fixed during the optimization, i.e., for ϕ equal to 1°, 2°, ..., 360°) with extremely tight convergence criteria (RMS values of 1 × 10 −6 hartree/bohr and 4 × 10 −6 bohr and maximum values of 2 × 10 −6 hartree/bohr and 6 × 10 −6 bohr on forces and displacements, respectively). Molecular structures and energies obtained with the constrained geometry optimizations were extracted and employed for the construction of the nuclear Hamiltonian matrices for the cases J = 0 (Ĥv ib 0 ) and J = 1 (Ĥv ibrot 1 ). Data extraction from the Gaussian outputs, construction of the nuclear Hamiltonian matrices (a formulation suitable for an implementation of the matrix elementŝ′ H vib , is provided in eqs SI-42 and SI-79, respectively; see the Supporting Information), solution of the eigenvalue problems, and data analysis were performed with a Python program written by the authors. In order to improve the computational efficiency of the code, Fortran90 subroutines were implemented for the numerical calculation of the nuclear Hamiltonian matrix Ĥv ibrot 1 . F2PY 59 was employed to interface Fortran90 subroutines to the main Python program.
With regard to the solution of the one-dimensional nuclear problem (eqs 8 and 9), the principal axis system (PAS) was adopted as the reference configuration of the Cartesian framework for each optimized molecular structure.

COMPUTATIONAL RESULTS
In this section, solutions of the 1D-NSE related to the pseudorotational angle ϕ are provided and discussed. More specifically, the solutions of the eigenvalue problems associated with the nuclear Hamiltonian matrices Ĥv ib 0 (for the case J = 0) are provided for each of the three molecular systems shown in Figure 1, together with rotational constants and rotational transition energies (for the case J = 0 → J = 1). A previous computational study concerning 1,2-dioxolane, 1,2-oxathiolane, and 3-chloro-1,2-dithiolane (focused on potential energy terms and symmetry considerations) can be found in ref 33 (with regard to 1,2-dioxolane, see also the earlier contributions of Cremer). 14,15 Unless otherwise specified, the mass of the most abundant isotope of the element is assumed for each atom. The computational data provided in this section allow a direct comparison with the experimental observables, especially with regard to experimental transition frequencies measured with far-infrared spectroscopies. An indirect comparison with experimental transition energies measured in the microwave region (i.e., rotational transition energies) is possible as well: in this case, differences of ν-specific rotational constants should be compared with differences of experimental rotational constants associated with different pseudorotational eigenstates.
4.1. 1D-PESs, Eigenvalues, and Eigenvectors of Ĥv ib 0 . Eigenvalues and eigenvectors of Ĥv ib 0 are shown in Figure 4 for each of the molecular systems investigated in this work. Moreover, 1D-PESs are shown: in this manner, the effects of V(ϕ) (i.e., magnitude and shapes of the potential energy barriers) on eigenvalues and eigenfunctions are highlighted.
The 1D-PESs associated with the pseudorotation angle ϕ exhibit a symmetry which is correlated to the symmetry of the planar configuration of a 5-term ring system (this fact was explicitly pointed out by Cremer and co-workers 14,36 and further explored in ref 33). More specifically, in the case of 1,2dioxolane, whose planar configuration pertains to the C 2v symmetry point group, the potential energy term satisfies both V(ϕ) = V(ϕ + π) and V(ϕ) = V(−ϕ); therefore, the 1D-PES can be fitted with a Fourier series containing only the cosine terms of the even harmonics. With regard to 1,2oxathiolane (C s ), only the relation V(ϕ) = V(ϕ + π) is satisfied and the 1D-PES can be fitted with a Fourier series containing both the sine and cosine terms of the even harmonics. 33 Finally, the planar configuration of 3-chloro-1,2-dithiolane pertains to the C 1 symmetry point group; therefore, symmetry cannot be employed to simplify the picture. However, the physical nature of the pseudorotational motion (and the mathematical definition of ϕ) ensures the periodicity of the 1D-PES for each of the molecular systems investigated in this work, i.e., V(ϕ) = V(ϕ + 2π) in all the cases (also in the case of 3-chloro-1,2-dithiolane).
The two minima of the 1D-PES shown in Figure 4a correspond to the same value of the potential energy and to two different structures; these structures form an enantiomeric pair, shown in Figure 5. A similar picture holds in the case of 1,2-oxathiolane: the optimized structures (which correspond to the two minima of the 1D-PES displayed in Figure 4b) are enantiomers and are depicted in Figure 6. The less symmetric 1D-PES of 3-chloro-1,2-dithiolane (see Figure 4c) is characterized by three minima with different energies and structures (reported in Figure 7).  Turning to the solutions of the 1D-NSE, Figure 4 deserves some comments. First, the eigenfunctions plotted in Figure 4 are not normalized, and the zeroes of their amplitudes correspond to the numerical value of the specific eigenvalue associated with the eigenfunction of interest. Second, the sign of an eigenfunction is not relevant, i.e. ψ ν and − ψ ν (where ν labels the eigenstate) denotes exactly the same eigenfunction (this statement is rather obvious when the eigenfunctions are not normalized, i.e., defined up to a numerical constant; however, the statement is true also for normalized eigenfunctions). Finally, it must be pointed out that when two eigenvalues are almost degenerate, their corresponding eigenfunctions can be partially superimposed in Figure 4, and therefore, one of the two can be partially hidden from the other one (e.g., this is the case of the eigenfunctions ψ 4 and ψ 5 of 1,2-dioxolane).
Eigenvalues of Ĥv ib 0 (labeled with E ν ) are reported in Tables  1, 2, and 3 in order to allow quantitative comparisons. The patterns of E ν reported in Tables 1 (1,2-dioxolane) and 2 (1,2oxathiolane) are quite similar: these results can be explained as a consequence of similarities concerning structural features and potential energy profiles. With regard to the structures, the only difference between the two molecules concerns the substitution of one of the two oxygen atoms of 1,2-dioxolane with the sulfur atom found in 1,2-oxathiolane. This substitution increases the mass of one of the atoms directly involved in the central 5-term ring, leading to an increase of the density of states: in the energy interval 0−500 cm −1 , there are 10 eigenstates in the case of 1,2-dioxolane (see Figure 4a and Table 1) and 16 eigenstates in the case of 1,2-oxathiolane (see Figure 4b and Table 2). Moreover, the substitution removes one of the symmetry elements of 1,2-dioxolane, with the consequences already mentioned on the symmetry of the 1D-PESs (symmetric and asymmetric shapes of the potential energy barriers in the cases of 1,2-dioxolane and 1,2oxathiolane, respectively; see Figure 4a,b).
Despite these differences, the undeniable similarities concerning the height of the barrier to pseudorotation (slightly less than 600 and slightly less than 700 cm −1 for 1,2oxathiolane and 1,2-dioxolane, respectively; see Figure 4a,b) and the two minima (which are enantiomers corresponding to the same energy values; see Figures 5 and 6) leads to a pattern which was already observed by other authors in the results obtained for other 5-term ring systems. 60,61 Three regions can be recognized. In the first one (E ν ≲ max[V(ϕ)]), pairs of almost degenerate eigenvalues are obtained (from ν = 0 to ν = 13 in the case of 1,2-dioxolane and from ν = 0 to ν = 17 in the case of 1,2-oxathiolane; see Tables 1 and 2). In the second one (E ν ≈ max[V(ϕ)]), a triplet of eigenvalues is found (ν = 14, 15, 16 in Table 1; ν = 18, 19, 20 in Table 2). Finally, in the third region (E ν > max[V(ϕ)]), the solutions are distributed as pairs of almost degenerate eigenvalues (these solutions are not reported in Tables 1 and 2, but the pairs ν = 17, 18 for 1,2dioxolane and ν = 21, 22 for 1,2-oxathiolane are depicted in panels a and b of Figure 4, respectively). Moreover, in the first region (E ν ≲ max[V(ϕ)]), the energy difference between two almost degenerate eigenvalues increases with the eigenvalue index (e.g., in the case of 1,2-dioxolane the energy difference between E 8 and E 9 is equal to 0.02 cm −1 , while the energy difference between E 10 and E 11 is 0.33 cm −1 ; with regard to the pairs ν = 0, 1, ν = 2, 3, ν = 4, 5, and ν = 6, 7, the values of the energies are equal as a consequence of the limited number of digits reported in the second column of Table 1; that is, the energy differences are smaller than 0.01 cm −1 ).
The eigenvalues of Ĥv ib 0 in the case of 3-chloro-1,2-dithiolane (second column of Table 3) exhibit a pattern which is quite different with respect to the cases discussed above (1,2dioxolane and 1,2-oxathiolane). First, the density of states in the case of 3-chloro-1,2-dithiolane is clearly the highest one among the molecular systems investigated in this work  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article (compare the results shown in Figure 4c with the results depicted in Figure 4a,b); this observation can be explained by the higher molecular mass of 3-chloro-1,2-dithiolane (with respect to 1,2-dioxolane and 1,2-oxathiolane), which is partly due to the presence of a chlorine atom. Second, in the region there is not a pattern of almost degenerate eigenvalues discussed above for the cases of 1,2-dioxolane and 1,2-oxathiolane (compare the values provided in the second column of Table 3 with the values reported in Tables 1 and 2). This is not surprising, because in this case there are three minima with different energies (corresponding to three asymmetric structures; see Figures 4c and 7) and a more complex potential energy profile (see Figure 4). On the other hand, when E ν > max[V(ϕ)] the patterns found for the three molecular systems investigated in this work are similar, with pairs of almost degenerate eigenvalues (only the first pair of eigenstates of this region is depicted in Figure 4; the eigenstates with higher eigenvalues are not reported). From the mathematical point of view, the 1D-NSE associated with the pseudorotation in a 5-term ring system in the J = 0 case is a Sturm−Liouville problem. For completeness, a concise treatment based on this alternative point of view is provided in section 3 of the Supporting Information; more specifically, the existence of exact and approximated analytical solutions to the problem investigated in this subsection is discussed, particularly with regard to the dependence of g ϕϕ , V′, and V (see eq 10) on the pseudorotation angle ϕ. An assessment of the importance of the structural relaxation effect (i.e., the importance of an explicit account of the dependence of g ϕϕ and V′ on ϕ) for the solution of the 1D-NSE in the J = 0 case is also proposed. can be listed in ascending order (on the basis of their numerical value), and the assignment could be carried out automatically. That is, the first element of the list of eigenvalues of Ĥv ib 0 (E 0 ) is associated with the first, the second, and the third elements of the list of eigenvalues of Ĥv ibrot 1 (i.e., these three elements are labeled E 0 [1 01 ], E 0 [1 11 ], and E 0 [1 10 ]); the second element of the list of eigenvalues of Ĥv ib 0 (E 1 ) is associated with the fourth, the fifth, and the sixth elements of the other list, and so on. In practice, this solution works only if the energy differences among the eigenvalues of Ĥv ib 0 are sufficiently large. In this work, this solution works in the case of 3-chloro-1,2-dithiolane but not in the cases of 1,2dioxolane and 1,2-oxathiolane. Therefore, in this work another scheme was employed, which is based on the superposition of the eigenstates of Ĥv ibrot 1 and Ĥv ib 0 (a more detailed explanation can be found in section 2 of the Supporting Information, see eqs SI-70−SI-73): because ⟨ψ vib 0(ν) |ψ vib 0(ν′) ⟩ = 0 when ν ≠ ν′, if the approximation ψ vibrot 1(x) ≈ ψ vib 0(ν) ψ rot 1(νl) holds (i.e., ψ vibrot 1(x) ≡ ψ vibrot 1(νl) ), the superposition between ψ vibrot 1(x) and ψ vib,ext.
The rotational constants A ν , B ν , and C ν are plotted as functions of the index ν in Figures 8−10. Values of the equilibrium rotational constants A e , B e and C e for the optimized geometries (see Figures 5−7) of the three molecular systems investigated in this work are provided in Table 4.
As mentioned in the previous subsection, in the cases of 1,2dioxolane and 1,2-oxathiolane the eigenvalues of the first eigenstates of Ĥv ib 0 are distributed in pairs of almost equal numerical values. Therefore, eigenvalues of the eigenstates of Ĥv ibrot 1 which are assigned to two almost degenerate eigenstates of Ĥv ib 0 are expected to be very similar. This is the case of 1,2oxathiolane: the numerical values of ΔE ν [1 01 ], ΔE ν [1 11 ], and ΔE ν [1 10 ] reported in Table 2 are equal (or almost equal) for each of the following pairs of indices: ν = 0, 1; ν = 2, 3; ν = 4, Table 3. 3-Chloro-1,2-dithiolane: Eigenvalues of Ĥv ib 0 (E ν ), Lowest Rotational Transition Energies (J = 0 → J = 1), and ν-Specific Rotational Constants (A ν , B ν , and C ν )  Table 2). Numerical values of A ν , B ν , and C ν exhibit the same regularities (see Figure 9 Table 2) . Despite the undeniable dependence of A ν , B ν , and C ν on the value of the index ν, in the case of 1,2-oxathiolane limited changes from the original equilibrium values (A e , B e , and C e ) are observed between ν = 0 and ν = 17 (see Tables 2 and 4 With regard to 1,2-dioxolane, rotational constants are plotted in Figure 8. The analysis proposed above for the 1,2oxathiolane molecule is partially valid also for 1,2-dioxolane. However, some differences must be pointed out. First, the values of ΔE ν [1 10 ], A ν and B ν assigned to a pair of almost degenerate eigenstates of Ĥv ib 0 are equal (or almost equal; see Table 1 and Figure 8a,b) 11 ], and C ν assigned to the same pairs of almost degenerate eigenstates of Ĥv ib 0 are different (see Table 1 and  Tables 1 and 4); these intervals are wider than the ones calculated for 1,2-oxathiolane.  The pattern shown in Figure 8c is not the expected one (at least for ν values in the range between 0 and 5) and deserves further comments. In the case of 1,2-dioxolane, the reference system adopted for the Cartesian framework of each optimized molecular structure might be inadequate: as a consequence, the pattern calculated for ΔE ν [1 01 ], ΔE ν [1 11 ], and C ν might be incorrect. As mentioned in Computational Protocol, all the results proposed in this section were calculated employing the PAS reference system. To verify the reliability of the results proposed in Table 1 (and shown in Figure 8) the Hamiltonians Ĥv ib 0 and Ĥv ibrot 1 were computed employing another reference system, named the zeta axis system (ZAS). ZAS is described in section 4 of the Supporting Information, where a comparison of the results obtained with the two reference systems (PAS and ZAS) is proposed. When the ZAS reference system is employed, the expected patterns (first eigenvalues distributed in pairs of equal numerical values) of ΔE ν [1 01 ], ΔE ν [1 11 ], and C ν are obtained (see Table SI 4 and Figure SI 6). Therefore, a possible explanation for the discrepancies between the values of ΔE ν [1 01 ], ΔE ν [1 11 ], and C ν reported in Table 1 and the expected ones is related to the choice of the reference axis system.
The results obtained for 3-chloro-1,2-dithiolane are reported in Table 3 and shown in Figure 10; these results are qualitatively different from the ones depicted in Figures 8  and 9. From lower to higher values of the index ν, three distinct and convergent branches can be identified in each plot (see Figure 10a−c): The first branch consists of the rotational constants corresponding to the following values of the index ν: 0, 1, 2, 3, 5, 7, 8, 10, 12, 13, and 15. The second one comprises the rotational constants corresponding to the values 4, 6, 9, 11, 14, and 16 of the index ν. The third (and the shortest) one involves the rotational constants corresponding to the values 22, 25, and 26 of the index ν. The three branches can be associated with the three minima of the 1D-PES shown in Figure 4c because the values of the ν index, which identify one of the three branches, correspond to eigenstates of Ĥv ib 0 which are prevalently located in the potential well of one of the three minima. This is particularly evident in the cases of the first and the second branches; for the sake of clarity, the relevant portion of the 1D-PES and the eigenstates of Ĥv ib 0 corresponding to a specific branch of Figure 10a−c are shown in Figure 11. With regard to the third branch, it must be underlined that only one eigenstate is prevalently located inside the potential well: this is not surprising and is due to the shallowness of the potential well associated with this branch (see Figure 11c). However, the amplitude of the other two eigenstates corresponding to the third branch is clearly affected by the highest potential well of Figure 11c, and the results presented so far suggest an influence of the highest potential well of Figure 11c also for the values of the rotational constants assigned to the eigenstates 25 and 26. The rotational constants  Table 3) are close to the equilibrium rotational constants reported in Table 4 (associated to the optimized geometries of Figure 7).   The accuracy of the energy values employed to construct the 1D-PESs (see Figure 4) depends on the level of theory chosen for the electronic structure calculations (in this work, B3LYP/ maug-cc-pVTZ). In ref 33, the same basis set was employed in conjunction with a double hybrid exchange−correlation functional 66 and empirical dispersions 67,68 (B2PLYP+D3BJ/ maug-cc-pVTZ); suitable experimental data were available for some of the 5-term ring systems examined in ref 33, allowing an estimation of the accuracy associated with the level of theory employed in that work, which turned out to be satisfactory (particularly for molecules with relevant pseudorotational barriers). The 1D-PESs proposed in this work are not equal, but they are similar to the results reported in ref 33 (see Table 5). Therefore, the 1D-PESs employed in this work should be at least qualitatively correct.
To estimate the accuracy of the approximated 1D flexible model employed in this work to formulate and to solve the nuclear problem, a comparison with the results of DVPT2 calculations 69 (performed with the Gaussian 16 suite of programs) is provided in Tables 6 and 7. A perturbative  Numbers are referred to the eigenvalues ν calculated with the 1D flexible model (see Tables 1−3). b The roman numerals refer to the conformational energies (labeled in ascending order). Calculated employing the ZAS reference system. b Calculated employing the PAS reference system (see Tables 1−3).
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article treatment should be adequate for the characterization of the lowest eigenstates of the molecular systems investigated in this work; each minimum is located at the bottom of a slightly anharmonic potential energy well, and therefore, a VPT2 treatment should provide reliable results (this is no longer true when the splitting due to the tunneling between two adjacent wells is relevant). The data reported in Table 6 show a good agreement between vibrational transition energies calculated with DVPT2 and 1D flexible models and suggest the reliability of the 1D model employed in this work at least for the calculation of the differences between the lowest eigenvalues associated with the pseudorotational motion. 1D models are often employed to reproduce experimental vibrational transition energies associated with a specific LAM (see for example refs 70 and 71): for the molecular systems characterized in this work such a comparison between computational and experimental data would be extremely useful, especially to quantitatively assess the numerical values computed for the eigenvalues of Ĥv ib 0 . With regard to the calculation of rotational constants, it must be pointed out that a 1D model cannot be employed to calculate all the corrections arising from the vibration−rotation interaction, because A ν ≠ A ν , B ν ≠ B ν , and C ν ≠ C ν (where ν ≡ {ν 1 , ν 2 , ..., ν 3K−6 }). In other words, the vibration−rotation interaction term due to a single mode is included in the calculated values of A, B, and C reported in Table 7 (namely, A 0 = A e + ΔA 0 ), while to reproduce the ground-state experimental values A 0 , B 0 , and C 0 the vibration−rotation interaction terms due to all the internal degrees of freedom should be included (namely, A 0 = A e + ∑ i = 1 3K−6 ΔA 0i ). As a consequence, the disagreement between the calculated values of 1,2-dioxolane reported in Table 7 and the experimental  values of the rotational constants taken from ref 65 is expected, because the values correspond to conceptually distinct quantities (i.e., A 0 ≠ A 0 ).
Corrections of the rotational constant values due to the interaction between rotational and pseudorotational motions were calculated with DVPT2 and 1D flexible models (see Table 7). The signs of the calculated corrections are the same for both the models and for each rotational constant, while some differences are observed with regard to the magnitude (especially in the case of 3-chloro-1,2-dithiolane, for which the major discrepancy between DVPT2 and 1D flexible models is found; see Table 7).
ν-specific rotational constants reported in this work could be useful to reproduce the differences between rotational constants ( ); this point should be verified through a comparison with experimental values of the rotational constants associated to excited pseudorotational eigenstates. Unfortunately, the experimental values of A 0 , B 0 , and C 0 are available only for the 1,2-dioxolane molecule, and the rotational constants associated with excited pseudorotational eigenstates are not available in the literature for the three molecular systems investigated in this work.
In the case of 1,2-dioxolane, the calculated values of the equilibrium rotational constants A e , B e , and C e are not far from the experimental values of A 0 , B 0 , and C 0 (see Table 8). With regard to the equilibrium rotational constants, the values calculated at the B2PLYP(D3BJ)/maug-cc-pVTZ level of theory should be more accurate than the values calculated at the B3LYP/maug-cc-pVTZ level of theory (see the values reported in Table 8), as claimed by previous computational studies. 72−75 Moreover, the calculated value of C e obtained at the B3LYP/maug-cc-pVTZ level of theory (see Table 8) can be questioned, because calcd C e < exptl C 0 (generally, the inclusion of all the vibration−rotation interaction terms to perform the calculation of A 0 , B 0 , and C 0 starting from the values of A e , B e , and C e is expected to provide a lower value of the calculated rotational constants, i.e., C 0 < C e ). However, the computational approach suggested in this work should be employed to reproduce the dif ferences between experimental rotational constants, rather than their absolute values. Accurate values of the equilibrium rotational constants are not needed if the efforts are devoted to the calculation of differences between experimental rotational constants; in this case, the accuracy of the numerical values calculated at the B3LYP/maug-cc-pVTZ level of theory should be acceptable.
Finally, a brief summary of the conclusions proposed in this section is provided: Accuracy: For a quantitative assessment of the numerical values reported in this work, an extensive comparison with experimental data is highly desirable. However, for the J = 0 case, the picture provided in this work should be at least qualitatively correct. For the J = 1 case, the protocol employed in this work is not suitable for a direct calculation of the experimentally observed rotational constants A ν , B ν , and C ν because only the interaction between pseudorotational and rotational motions is calculated and taken into account (vibration−rotation interactions involving the other 3K − 5 internal motions cannot be safely neglected). Nevertheless, dif ferences of the ν-specific rotational constant values reported in Tables 1−3 could be a reasonable estimation of the experimental ones, at least for the lowest eigenvalues (in previous contributions, 76,77 the original Meyer's implementation of the 1D flexible model was employed to reproduce the experimental variation of rotational constants associated with different values of the index ν).
Limitations: The dynamical contribution of all the other 3K − 5 internal degrees of freedom is neglected, and the full dimensional NSE is reduced to a 1D-NSE. Furthermore, the rovibrational problem is formulated and solved only in the J = 0 and J = 1 cases, while formulations for the J > 1 cases (available in literature) are not discussed.
Possible Improvements: A proper inclusion of the coupling between the pseudorotational motion and the other internal degrees of freedom would be useful, especially for the computational characterization of pseudorotational pathways with very small potential energy barriers (less than 100 cm −1 ) or for an accurate determination of higher-energy eigenvalues The Journal of Physical Chemistry A pubs.acs.org/JPCA Article of Ĥv ib 0 . Different strategies were proposed and implemented, such as the formulation of an effective 1D nuclear Hamiltonian (BO-like separation between the pseudorotational motion and the other internal degrees of freedom) 24−26 or the formulation of a higher-dimensional nuclear Hamiltonian in curvilinear coordinates (2D, 19 3D, 28 or even full dimensional 30 formulations). With regard to the J = 1 case, a complete calculation of the vibration−rotation interaction (taking into account all the internal motions) would allow the calculation of experimentally observed rotational constants A ν , B ν , and C ν . An extension devoted to the formulation of the nuclear Hamiltonian (and to the solution of the nuclear problem) for higher values of J (J > 1) would be also useful. 27,29 6. CONCLUSIONS In this work, the implementation of an approach suitable for the calculation and the assignment of ν-specific rotational constants A ν , B ν , and C ν is described and the solutions to the 1D-NSEs (for J = 0 and J = 1) associated with the pseudorotational motion of 1,2-dioxolane, 1,2-oxathiolane, and 3-chloro-1,2-dithiolane are reported as test cases. The implementation proposed is suitable for a computational characterization of the pseudorotational motion in 5-term ring systems; in particular, a procedure to perform a relaxed scan along the pseudorotation angle ϕ (defined according to the original definition proposed by Cremer and Pople) without the explicit employment of the analytical gradient was devised, implemented, and applied.
From another point of view, this work is an initial assessment of usefulness and advantages of CPCs for the formulation of the NSE. The usefulness of CPCs and RDCs for the construction of nD-PESs (or to perform conformational studies) is well-known, 33,78−81 and new contributions concerning these coordinate systems and their generalization appeared in the literature in the last two years; 81,82 therefore, in the opinion of the authors the development and the assessment of computational tools aiming at the employment of these coordinates in a quantum mechanical description of the nuclear motions is worthwhile.
Limitations and possible improvements of the 1D model discussed in this contribution are mentioned.