Exact numerical computation of a kinetic energy operator in curvilinear coordinates

The conformation and dynamical behavior of molecular systems is very often advantageously described in terms of physically well-adapted curvilinear coordinates. It is rather easy to show that the numerous analytical expressions of the kinetic energy operator of a molecular system described in terms of n curvilinear coordinates can all be transformed into the following more usable expression: T̂=∑ijf2ij(q)∂2/∂qi∂qj+∑if1i(q)∂/∂qi+ν(q), where f2ij(q), f1i(q), and ν(q) are functions of the curvilinear coordinates q=(…,qi,…). If the advantages of curvilinear coordinates are unquestionable, they do have a major drawback: the sometimes awful complexity of the analytical expression of the kinetic operator T̂ for molecular systems with more than five atoms. Therefore, we develop an algorithm for computing T̂ for a given value of the n curvilinear coordinates q. The calculation of the functions f2ij(q), f1i(q), and ν(q) only requires the knowledge of the Cartesian coordinates and their derivatives in terms of the n curvilinear coordinates. This coordinate transformation (curvilinear→Cartesian) is very easy to perform and is widely used in quantum chemistry codes resorting to a Z-matrix to define the curvilinear coordinates. Thus, the functions f2ij(q), f1i(q), and ν(q) can be evaluated numerically and exactly for a given value of q, which makes it possible to propagate wavepackets or to simulate the spectra of rather complex systems (constrained Hamiltonian). The accuracy of this numerical procedure is tested by comparing two calculations of the bending spectrum of HCN: the first one, performed by using the present numerical kinetic operator procedure, the second one, obtained in previous studies, by using an analytical kinetic operator. Finally, the ab initio computation of the internal rotation spectrum and wave functions of 2-methylpropanal by means of dimensionality reduction, is given as an original application.