triatomic molecules [5]. These algebraic solutions for the RSPT equations are found using symbolic mathematics and automatically translated into FortranForm via Mathematica. The subsequent calculation of molecular energies, intensities and properties scales linearly with the order of RSPT [6]. The second is the analytical representation of numerical PESs, where scaling issues seen in first-principles quantum mechanical perturbation theory treatments of intramolecular motion are traced to analytical representations of PESs, accomplished here using singular value decomposition (SVD) [7] For example, using SVD, a 6 th-degree power series expansion of the PES for a 12-atom molecule, having 30 vibrational degrees of freedom, requires the inversion of rank 2-million matrix [2]. However, inversions of high-rank matrices also scale linearly and are accomplished straightforwardly and efficiently using modern parallel supercomputers [8]. In order to describe vibrational and rotational states, the Schrodinger equation for intramolecular nuclear motion is solved in the context of the Born-Oppenheimer approximation. The equation can be cast into a form involving analytic forms for the potential energy operator \(V\) in the Hamiltonian operator. This was accomplished for the ground vibrational state of triatomic molecules by [PERSON] and [PERSON] using variation perturbation theory [9]. It has been solved for any polyatomic molecule using RSPT for the \(I=1\), (3N-6) modes of vibration in terms of normal coordinates \(Q_{i}\), as follows (in atomic units) [2]:  
\[H=T+V \tag{1}\]  
\[H=H_{0}+H_{1}+H_{2} \tag{2}\]  
\[H_{0}=\frac{1}{2}\sum\omega_{i}{p_{i}}^{2}+\frac{1}{2}\sum\omega_{i}{q_{i}}^{2 }=T+V_{0} \tag{3}\]  
\[H_{1}=\sum k_{ijk}q_{i}q_{j}q_{k}=V_{1} \tag{4}\]  
\[H_{2}=\sum\;k_{ijkl}\;q_{i}\;q_{j}\;q_{k}+\frac{1}{8}\;\sum\;\left(I_{\alpha} \right)^{2}\left(4\Pi I_{\alpha}-1\right)=V_{2}+V_{2}t \tag{5}\]  
In Equations (1)-(5) \(T\) is the total kinetic energy operator defined in terms of the harmonic frequencies \(\omega_{i}\), the mass-reduced normal coordinates \(q_{i}=Q_{i}/\surd\omega_{i}\) and the momenta \(p_{l}=-i\;\partial/\partial q_{i}\). The second term in Equation (5) corresponds to the correction due to Coriolis interactions. \(I_{\alpha}\) are the moments of inertia relative to the principal moment of inertia axes \(\alpha=A\), \(B\), and \(C\), and \(\Pi_{\alpha}\) are the vibrational total angular momenta. Normal coordinates can be obtained using a finite differences procedure for calculating the second derivatives of the total energy in terms of mass-weighted cartesian nuclear displacement coordinates, the Hessian matrix. The eigenvectors (L\({}_{\rm{ii}}\)) of the Hessian matrix comprise the linear transformation matrix that defines the normal coordinates in terms of the mass weighted cartesian displacement coordinates. The eigenvalues of the Hessian matrix are the harmonic frequencies [2]. The constants \(k_{ijk}\), and \(k_{ijkl}\) are cubic and quartic force constants that account for anharmonicity. These constants are calculated from an initial power series fit of the PES in terms of internal coordinates by means of SVD, using the normal coordinate transformation matrix L\({}_{\rm{ij}}\)[2]. The total potential energy operator to the second-order of RSPT is then  
\[V=V_{0}+V_{1}+V_{2}+V_{2}^{{}^{\prime}} \tag{6}\]  
Solution of the RSPT equations leads to the following expression for vibrational term energies, the anharmonicity constants being algebraic formulae that only involve the \(\omega_{i}\), \(k_{ijk}\), and \(k_{ijkl}\) constants, obtained as noted above. Any desired line in the vibrational spectrum can be calculated from term energy differences.