A general formulation of the quasiclassical trajectory method for reduced-dimensionality reaction dynamics calculations

Abstract

Dimension reduction by freezing the unimportant coordinates is widely used in intramolecular and reaction dynamics calculations when the solution of the accurate full-dimensional nuclear Schrödinger equation is not feasible. In this paper we report on a novel form of the exact classical internal-coordinate Hamiltonian for full and reduced-dimensional vibrational motion of polyatomic molecules with the purpose of using it in quasiclassical trajectory (QCT) calculations. The derivation is based on the internal to body-fixed frame transformation, as in the t-vector formalism, however it does not require the introduction of rotational variables to allow cancellation of non-physical rotations within the body-fixed frame. The formulas needed for QCT calculations: normal mode analysis and state sampling as well as for following the dynamics and normal-mode quantum number assignment at instantaneous states are presented. The procedure is demonstrated on the CH₄, CD₄, CH₃D and CHD₃ isotopologs of methane using three reduced-dimensional models, which were previously used in quantum reactive scattering studies of the CH₄ + X → CH₃ + HX type reactions. The reduced-dimensional QCT methodology formulated this way combined with full-dimensional QCT calculations makes possible the classical validation of reduced-dimensional models that are used in the quantum mechanical description of the nuclear dynamics in reactive systems [A. Vikár et al., J. Phys. Chem. A, 2016, 120, 5083–5093].

---

1. Introduction

Dimension reduction is often used in modeling phenomena in chemical physics to reduce the complexity of the model. By selecting the degrees of freedom that are relevant to the investigated properties of the system, one can concentrate the effort on a model whose smaller size allows one to perform a simulation at a higher level of sophistication. In molecular physics, reduced-dimensional (RD) models have been used in the description of nuclear motion both in rovibrational spectroscopy^1,2 as well as in molecular^3–5 and reaction dynamics simulations.^6–17 In both cases, the possibility of simplification is offered by the separation of time scales and the weakness of the coupling between the various modes of nuclear motion.

In vibrational spectroscopy, the semi-rigid modes are well described by the normal mode approximation^18–20 in which the potential and kinetic energies are considered as quadratic functions of Cartesian or internal coordinates and of the conjugate momenta, resp. The description of large amplitude motion (LAM), in floppy molecules, such as hindered rotor type modes, however, require a more sophisticated treatment, because the quadratic approximation does not work, often within the space swept by the zero-point motion. The frequencies of such modes are generally much smaller than those of stiffer vibrations, and the potential coupling between the fast and slow degrees of freedom is often also limited. A way to achieve a numerically feasible description of vibration of molecules with LAMs is that one reduces the dimensionality of the problem to those of the strongly anharmonic low-frequency modes by freezing the fast vibrations.^1,2

In molecular dynamics simulation of, for example, biomolecules in water, when a large number of solvent molecules are present, the high-frequency OH stretching vibrations require the integration time step to be small. However, their instantaneous phase has no effect on the much slower conformational motion, thus one can freeze them,^4,5 which allows a significantly faster, yet realistic simulation of the system.

In reaction dynamics, the reactive event often concerns a small number of internal coordinates involving only a few atoms, called active modes, and the rest are considered as “spectators”. Dimension reduction is possible because the potential- and often the kinematic coupling between the active and spectator modes are small. In the corresponding scattering calculation the active degrees of freedom are treated explicitly, while the coordinates of spectator modes are frozen at their values at the saddle point of the potential energy surface (PES) or at the equilibrium geometry of reactant molecules. In quantum and state-to-state quasiclassical dynamical calculations it is necessary to know the quantum states of the reactants and products, for which one needs to properly characterize the vibrational motion involving only the active degrees of freedom of the reactant and product molecules, respectively. This is equivalent to the reduced-dimensional vibrational spectroscopic problem. RD models are mostly used in quantum mechanical simulations of reactions, because of the exponential growth of the computational effort with the number of degrees of freedom,^1,2,6–16 but only rarely in quasiclassical trajectory (QCT) simulations, where the growth is linear.

In the QCT method,^20,21 the motion of the atoms is described classically and the only nuclear quantum effect considered is that the rovibrational energies of the reactant molecules are discrete. Accordingly, each rovibrational quantum state of the reactant molecules is simulated by an ensemble of semiclassically quantized classical states (i.e. coordinates and momenta). In RD reaction dynamical models the spectator degrees of freedom are frozen during the reaction. This means that the number of vibrational degrees of freedom of the reactant and the product molecules is also reduced. Consequently, the problem of semiclassical quantization also arises when state-to-state reactivity parameters or simply state distributions of product molecules are to be determined. In what follows, the generation of classical states corresponding to a quantum state will be called the direct problem, and the determination of the quantum state corresponding to a given classical state will be referred to as the inverse problem.

Application of dimension reduction in QCT simulations is rather scarce. One of the reasons for this is that there is no general theory for the vibrational analysis, initial state preparation and final state analysis is available. Among the few reduced-dimensional trajectory calculations performed so far, in a set only trivial reductions were applied: some of the Cartesian coordinates were simply disregarded. Lu and Hase²² applied RD models of benzene, obtained by constraining and truncating the molecule to planar C₃H₃ and C₃H moieties to prevent zero-point energy leakage from neighboring high-frequency modes during the simulation of intramolecular vibrational energy redistribution (IVR). Klossika and Schinke²³ investigated the photodissociation of HNCO induced by NH vibrational excitation, by constraining the atoms into a plane using a reduced-dimensional analytic potential energy surface calculated only for planar arrangement of atoms.

In another set of dynamical studies the vibrations of molecules or fragments were completely frozen and only their relative motion was simulated. Raff and coworkers^24,25 applied rigid-body dynamics to investigate rotational energy transfer between CO₂ and He as well as H₂. Rotational dynamics in collisions of H₂O and H₂ with frozen stretch vibrations were studied by Faure et al.²⁶ In some reactive scattering calculations bond lengths and angles were frozen, focusing again on rotational dynamics in the capture step of some bimolecular reactions (Maergoiz et al.,²⁷ Faure et al.,²⁸ Harding et al.²⁹). Harding et al.³⁰ investigated the roaming dynamics of the photochemical decomposition of CH₃CHO, where they froze the vibrations of the fragments to avoid zero-point energy leakage and the need for constructing a high-dimensional analytic PES.

More complex constraints were considered by the authors of the present paper in a recent study comparing the results of reduced- and full-dimensional (FD) quasiclassical trajectory calculations. The purpose of that work was to assess the accuracy of the Palma–Clary RD quantum dynamical model of the CH₄ + H → CH₃ + H₂ reaction,³¹ whose FD counterpart is computationally too expensive to solve. In the present paper the general theory of reduced-dimensional QCT calculations used in that study is described. We demonstrate how the QCT method, including initial condition generation and final state can be consistently applied to RD models involving arbitrary constraints.

In the following, first we describe the three fundamental coordinate systems used in this work (Section 2.1); then in Section 2.2 we derive the vibrational Lagrangian in body-fixed Cartesian coordinates and then the vibrational Lagrangian and Hamiltonian in internal coordinates (Section 2.3). We discuss the connection of our formulation to the t-vector and s-vector formalisms in both full and reduced-dimensionality (Section 2.4). As applications, normal mode analysis (Section 2.5) and normal mode sampling (NMS, Section 2.6) in internal coordinates and the subsequent transformation of states to laboratory frame are described. In Section 2.7, two methods of RD trajectory integration are presented and compared. Section 2.8 is devoted to the inverse problem, where first the classical state given in laboratory frame is transformed to the internal coordinate system, and the normal coordinate displacements and momenta are calculated, from which the normal mode quantum numbers are determined. The equations presented in Section 2 (apart from Section 2.4) hold not only for reduced-dimensional models but are also applicable in full dimensionality and can also serve as a basis in the derivation of RD quantum Hamiltonians. In Section 3, as a proof of principle, the method is applied to a hierarchy of three RD models of the methane molecule CZ₃Y, in each of which the CZ₃ group is constrained to maintain C_3v symmetry. A complete reaction dynamics study based on this theory has been presented in ref. 31.

2. Theory

2.1 Frames and coordinate systems

Derivation of the classical vibrational Hamiltonian in internal coordinates starts from the full Lagrangian expressed in Cartesian coordinates in a space-fixed (a.k.a. laboratory) frame:where for an N-atomic system X = (X_1x, X_1y, X_1z,…,X_Nz)^T and Ẋ are the 3N-component coordinate and corresponding velocity column vectors, respectively, which are composed of the corresponding atomic coordinate R_i = (X_ix, X_iy, X_iz)^T and velocity Ṙ_i vectors. M = diag(m₁,…,m_N) is the diagonal 3N × 3N mass matrix, containing the atomic masses. The classical mechanical state of the system can be given either by the coordinates and the velocities, (X,Ẋ) or by the coordinates and the conjugate momenta, (X,P_X), the Cartesian momenta being defined as P_X = MẊ.

For the description of molecular vibrations, internal coordinates such as valence coordinates (bond lengths, bond and torsion angles etc.) are much more meaningful than Cartesians as the forces acting between atoms are inherently intramolecular, i.e., they do not depend on the position and orientation of the molecule. In addition, the force constants defined in terms of valence coordinates can be rationalized using chemical intuition (for example, they are roughly transferable between molecules).³² Furthermore, the use of internal coordinates is advantageous also when approximations to the Hamiltonian (e.g. quadratic- or quartic-order) are applied as they can describe large amplitude curvilinear motion more effectively than Cartesians.

In general, an internal coordinate is such a function of Cartesian coordinates, whose value is invariant under displacement and rotation of the molecule, thus it is necessarily formulated by using scalar (a·b = a^Tb), and vector (a × b) products of atomic Cartesian coordinate vector differences R_i − R_j. Consequently, functions of internal coordinates are also internal coordinates. It is worth noting that those internal coordinates that are defined in terms of a cross product (e.g., torsion angles), change sign under mirroring the molecule through a plane (they are pseudoscalars).

For the description of the vibrating molecule one needs to define n independent, otherwise arbitrary internal coordinates y(X) = (y₁(X),…,y_n(X))^T in terms of Cartesian coordinates. n is less than or may be equal to f = 3N − 6, the number of internal degrees of freedom. If n equals f then the set of internal coordinates is complete and the model is called full-dimensional. In reduced-dimensional models the set of irredundant internal coordinates is incomplete (n < f) and the remaining f − n internal degrees of freedom are constrained by fixing f − n functions y_n+1(X),…,y_f(X) at values y_n+1,0,…,y_f,0. These constraint functions are generally expressed in terms of the usual valence coordinates. For example, such a function can measure the deviation from some desired symmetry, e.g., it may be the difference of two bond lengths, which is constrained to zero. Variables y_n+1(X),…,y_f(X) expressing the constraints are formally internal coordinates, because their values should also be independent of the position and the orientation of the molecule. Note that these constrained variables are not included in vector y.

In full-dimensional models of vibrating molecules, the kinetic energy in internal coordinates is given with the help of the Wilson B matrix evaluated at the instantaneous geometry:¹⁹

The row vectors of the (3N − 6) × 3N dimensional B matrix are called vibrational s-vectors^1,33 (s_i = dy_i/dX). The nonredundant internal coordinates y are defined for all values of laboratory coordinates X, so that the inverse mass matrix

G_y,vib = BM⁻¹B^T (3)

properly assigns masses to the internal coordinates and the vibrational kinetic energy in internal coordinates written as,will be exact. However, in reduced-dimensional models (n < f), the X coordinates are interrelated by the constraints and the B matrix defined in (2) with reduced number of rows lacks the information on the corresponding constrained internal coordinates, which is required to disentangle them from the free ones. Consequently, such a reduced-dimensional B matrix cannot be used for the construction of the exact reduced-dimensional kinetic energy expression, unless the frozen internal coordinates are orthogonal in some sense to the free ones (see ref. 1 and Section 2.4).

In order to circumvent this problem, we use the inverse transformation, which converts f internal coordinates y to 3N lab Cartesian coordinates X. The inverse function X(y) and its partial derivatives by definition take into account the constraints because they are evaluated under the condition y_j(X) = y_j,0 for j = n + 1,…,f. However, the internal coordinates do not determine the position and orientation of the molecule in the Cartesian system. To locate and orient the molecule one can utilize an intermediate body-fixed frame and an attached Cartesian coordinate system. In this auxiliary body-fixed frame the Cartesian coordinate 3N-vector and the coordinate 3-vector of atom i will be denoted by x = (x_1x, x_1y, x_1z,…,x_Nz)^T and r_i = (x_ix, x_iy, x_iz,)^T, respectively. The body-fixed coordinates X are connected to the space-fixed frame by an instantaneous translation and rotation, summarized in the function x(X). The definition of the intermediate frame allows us to derive the Lagrangian in internal coordinates via converting the kinetic and potential energy expressions (1) from space-fixed to body-fixed frame and (2) from body-fixed to internal frame using the inverse transformations X(x) and x(y), respectively.

In what follows we proceed on this route in two steps. First we describe the body-fixed frame and its connection to the internal coordinates, x(y) and then its relationship to the space-fixed Cartesian frame, X(x) and show how a classical state given in internal coordinates can be transformed into the space-fixed frame.

There are many possible ways of defining a body-fixed Cartesian frame, and it depends on the system which one is the most favorable. One of the simplest possibilities is that one places the origin either at the center of mass of the molecule or at one of the atoms and selects four non-coplanar atoms within the molecule and orthonormalizes three vectors pointing from one of them to the other three to obtain the basis vectors e (α = x, y, z).

The definition of the body-fixed Cartesian coordinates as a function of internal ones is given by a vector–vector function

x = x(y), (5)

whose time derivative connects the internal coordinate and the Cartesian velocities:

ẋ = C(y)ẏ. (6)

The columns of the 3N × n matrix C(y) = dx/dy are known as vibrational t-vectors^1,33,34 (t_i = ∂x/∂y_i).

When one knows the x coordinates, the Cartesian coordinates of the atoms in the space-fixed frame can be obtained by considering that the body-fixed and space-fixed Cartesian frames can be brought into overlap by a translation and a rotation, i.e., the r_i coordinate vectors need to be rotated by matrix O_frame and shifted by vector R_frame to get the space-fixed coordinates R_i:

R_i = R_frame + O_framer_i. (7)

The atomic velocity Ṙ_i is obtained from that in the body-fixed frame (ẋ) by differentiating eqn (7) with respect to time:

Ṙ_i = Ṙ_frame + Ȯ_framer_i + O_frameṙ_i, (8)

where Ṙ_frame is the velocity of the body-fixed frame with respect to the space-fixed one. To get Ȯ_frame, one needs to take into account that the body-fixed frame may rotate around its origin with angular velocity ω_frame. The total time derivative of the rotation matrix O_frame is then Ȯ_frame = ω_frame × O_frame, where the cross product of ω_frame and matrix O_frame needs to be evaluated column by column. With this, the atomic velocity R_i in the space-fixed Cartesian frame is:

Ṙ_i = Ṙ_frame + (ω_frame × O_frame)r_i + O_frameṙ_i. (9)

Unless special care is taken, the body-fixed frame does not move together with the molecule: its instantaneous linear (Ṙ_frame) and angular velocities (ω_frame) differ from those of the molecule in the space-fixed frame. For example, when the body-fixed frame for a vibrating water molecule (H₂O) is selected to be centered at the O atom with the x-axis parallel to H¹H², then the linear velocity of the origin of the frame is not identical to that of the center of mass, and the antisymmetric OH stretch vibration generates angular motion of the molecule (see Fig. 1).

Fig. 1 Schematic drawing of the motion of a water molecule (H¹OH²) in a specific body-fixed frame, which is centered at the O atom with the x-axis parallel to the H¹H² line. (a) The equilibrium geometry with space-fixed displacements (blue arrows) due to antisymmetric stretch vibration. (b) The distorted molecule and the body-fixed frame aligned according to the new H–H axis. (c) The distorted molecule when the body-fixed frame is aligned as in (a), showing the corresponding atomic displacements in the body-fixed frame, which result in a clockwise rotated structure whose center of mass (CM) is displaced upward and to the left in the body-fixed frame.

Obviously, the physically correct description of motion requires both the displacements of the atoms of the molecule in the body-fixed frame and translation + rotation of the body-fixed frame in the lab frame. Consequently, the transformation of internal coordinates and momenta to body-fixed frames according to eqn (5) and (6) usually generates unphysical (often referred to as spurious³⁵) translation and rotation.

When the focus is on the vibrational motion of a molecule, however, the general procedure is that the motion of the body-fixed frame is disregarded and extra steps are made to eliminate the spurious rotation and translation which would falsify the effective masses assigned to internal coordinates. In the t-vector formalism the cancellation of spurious rotation is achieved by the introduction of rotational coordinates when vibrational energy levels are calculated.^1,2

In classical mechanics, when a vibrational state is generated by normal mode sampling, the molecule can artificially be cleared of these unwanted velocities.³⁶ In general, it is desirable to avoid the appearance of the unphysical translational and rotational terms, a novel way of which is proposed in the following sections.

2.2 Vibrational Lagrangian in body-fixed Cartesian coordinates

To derive the pure vibrational Lagrangian L_x,vib(x,ẋ) in body-fixed Cartesian coordinates, first the Lagrangian L_x(x,ẋ) for the non-translating and non-rotating body-fixed frame is obtained. To achieve this, one substitutes eqn (7) and (9) into eqn (1) and eliminates the rotational and translational motion of the body-fixed frame by setting Ṙ_frame = 0 and ω_frame = 0. Exploiting that matrix M is diagonal and matrix O_frame is unitary, the kinetic energy function in the new coordinates can be transformed into the same form in body-fixed Cartesian coordinates as it was in the lab Cartesian frame, in accordance with the expectations. The form of the potential energy function will also remain the same because it is a function of internal coordinates only, which are defined with dot and cross products (brief common notation: ), that are also left unchanged by O_frame.The curly bracket refers to the sets of products involving the atomic indices (j, k, l, m) that are used in the definition of internal coordinates.

We would like to obtain the pure vibrational Lagrangian. For this we need to decompose the kinetic energy into separate vibrational as well as translational and rotational parts. If one introduces mass-scaled coordinates x̃ = M^1/2x, the quadratic form of the kinetic energy expressed in the body-fixed Cartesian frame (T_x) reduces to the square of the mass-scaled velocity vector . This form is advantageous, because it allows one to decompose the instantaneous mass-scaled velocity vector into orthogonal translational, rotational and vibrational parts (we show later how). Once these components, , and are available, the kinetic energy can also be broken down into the corresponding T_trans, T_rot and T_vib terms:

Consequently, the vibrational part can be obtained by eliminating the components of the mass-scaled velocity that correspond to translation and rotation, i.e., projecting the velocities onto the complementary and orthogonal, instantaneous (geometry-dependent) vibrational subspace. If the 3N × 3N matrix that performs the desired projection in the space of mass-scaled velocities is denoted by P_vib(x), the pure vibrational Lagrangian can be formally written as:

Utilizing the fact that the orthogonal projectors onto the translational and rotational subspaces, P_trans and P_rot(x), can be easily found (see below) and their sum is complementary to P_vib(x), we define P_vib(x) in terms of them:

P_vib(x) = E − P_trans − P_rot(x), (13)

where E is the 3N × 3N unit matrix. In the following we define the basis vectors of the translational and rotational subspaces and using them, we construct the corresponding projectors.

The translational subspace of mass-scaled displacements and velocities is spanned by the three 3N-component translational basis vectors u^trans,α (α = x, y, z):

u^trans,α_iγ denotes the γ (γ = x, y, z) component of the mass-scaled displacement of atom i during translation of the whole molecule along axis α (see Fig. 2a–c). δ_αγ is the Kronecker symbol, which is evaluated using x = 1, y = 2, z = 3 assignments to the possible values of its indices α and γ. The translational basis vectors are related to translational t-vectors (t^trans,α) by mass-scaling. The translational basis vectors only depend on the masses and are independent of the geometry and also of the choice of the origin of the body-fixed frame. Consequently, the translational subspace is the same at all geometries and thus it includes all finite mass-scaled translations of the molecule.

Fig. 2 The three translational (a–c), two rotational (d and e) and one vibrational (f) normalized basis vectors of a homonuclear diatomic molecule (H₂, O₂, etc.) at a given orientation. The Cartesian coordinate system and the indices of atoms are shown in the leftmost panel. In panels (d and e), the center of mass (CM) is shown in blue, and the unit vectors along the principal axes are shown in red.

The rotational subspace is spanned by the three 3N-component rotational basis vectors u^rot,α(x), which describe the relative magnitude of the mass-scaled displacements of atoms due to an infinitesimal rotation of the molecule around the α principal axis (PA) of the instantaneous moment of inertia tensor (α = 1, 2, 3, but not x, y, z). Rotational basis vectors can be calculated from the orthonormal unit vectors of principal axes e^PA_α (given in the body-fixed frame) and the instantaneous position vectors of atoms ρ_i := r_i − r_CM (i = 1,…,N) relative to the center of mass of the molecule, r_CM:

u^rot,α_iγ characterizes the relative magnitude of the γ component (γ = x, y, z) of the mass-scaled displacement of atom i when the molecule is rotated infinitesimally around the α (= 1, 2, 3) principal axis (see Fig. 2d and e). The angle of the rotation around e^PA_α is denoted by φ^PA_α. ε_γστ is the Levi-Civita tensor which is evaluated using assignments x = 1, y = 2, z = 3 regarding the possible values of its indices γ, σ and τ. The rotational basis vectors are related to rotational t-vectors (t^rot,α) corresponding to rotations around the principal axes by mass-scaling. The rotational basis vectors depend on the geometry. Thus they can be used to describe only infinitesimal mass-scaled displacements of atoms during the rotation of the molecule. The translational and the rotational basis vectors are orthogonal to each other, and by normalizing them an orthonormal set of translational and rotational basis vectors, u^trans,α₀ and u^rot,α₀(x) can be obtained. The proof of orthogonality and the derivation of the squared norms of the basis vectors are presented in the Appendix.

The instantaneous vibrational subspace, which is orthogonal to the translational and rotational subspaces, also depends on the geometry, thus it will span only infinitesimal mass-scaled displacements during the vibration of the molecule. When the shape and the orientation of the molecule do not change during its motion, for example, during the symmetric stretch vibration of CH₄ or during the translation of any molecule, the rotational and vibrational subspaces will not change either. A possible orthonormal basis (u^vib,i₀(x_e), i = 1,…,f) of the vibrational subspace at equilibrium geometry (x_e) is formed by the vibrational normal mode eigenvectors (see Fig. 2f), which are obtained from harmonic vibrational analysis done in Cartesian frame at x_e and are mass-scaled by definition. The orthonormal basis vectors of the translational, rotational and vibrational subspaces of a homonuclear diatomic molecule are summarized in Fig. 2. The orthogonal projection matrices P_trans and P_rot(x) can be defined as dyadic products of the relevant basis vectors:

These are the projectors to be used in eqn (13) to generate the orthogonal projection matrix P_vib(x) onto the complementary, vibrational subspace. Matrix P_vib(x), being an orthogonal projector, is idempotent (P_vib²(x) = P_vib(x)) and symmetric while the mass matrix M is diagonal. Consequently, by introducing an effective vibrational mass matrix:

M_vib(x) = M^1/2P_vib(x)M^1/2, (18)

which is geometry dependent and dense as opposed to matrix M, the Lagrangian in eqn (12) can be rewritten in the form:Matrix M_vib(x) is singular, since it assigns non-zero masses only to the motion within the vibrational subspace, which has fewer dimensions than 3N. Momenta p_x,vib canonically conjugate to coordinates x are obtained by differentiating the Lagrangian:

While the velocity vector ẋ may describe translation and rotation in addition to vibration of the molecule in the body-fixed frame, the momentum vector p_x,vib describes only vibrations, because it is obtained by the singular mass matrix M_vib(x). Thus, the Euler–Lagrange equations of motion cannot describe translation and rotation within the body-fixed frame.

We note that the present derivations are similar to the projection method proposed by Miller et al.³⁷ and Szalay used projections for the approximate decomposition of the kinetic energy in the Eckart frame.³⁸ The difference is that we apply the exact decomposition of the instantaneous kinetic energy.

2.3 Vibrational Hamiltonian in internal coordinates

The vibrational Lagrangian in internal coordinates y and velocities ẏ can be obtained from eqn (19) using the transformations in eqn (5) and (6):Here we introduced the potential energy V_y(y) = V(x(y)) as well as the n × n vibrational mass matrix, M_y,vib as a function of internal coordinates y:

M_y,vib(y) = C^T(y)M^1/2P_vib(x(y))M^1/2C(y). (22)

At this point, it becomes obvious that by projecting onto the vibrational subspace, we cancel the spurious translation and rotation in the body-fixed frame and the masses corresponding to them will not contaminate the mass matrix assigned to the internal coordinates, which would unphysically increase the matrix elements. Without the projection, for example, harmonic vibrational analysis would give incorrect, reduced frequencies for some of the normal modes.

Eqn (22) applies as it is regardless whether n = f or n < f. Its form implies that in reduced-dimensional models the reduced-dimensional M_y,vib(y) matrix can also be obtained from the full-dimensional M_y,vib by simply deleting the rows and columns corresponding to the constrained internal coordinates.

Momenta canonically conjugate to the internal coordinates are obtained as:

and the velocity ẏ as a function of y and p_y will be

ẏ = M_y,vib⁻¹(y)p_y = G_y,vib(y)p_y. (24)

Here, the n × n dimensional G_y,vib(y) matrix is the inverse of the non-singular M_y,vib(y) matrix. Applying Legendre transformation to the Lagrangian, one can obtain the vibrational Hamiltonian in internal coordinates:This form is correct in any nonredundant set of internal coordinates, be it reduced (n < f) or complete (n = f).

2.4 Connection with the t-vector and the s-vector formalisms

In the previous two sections we derived the vibrational kinetic energy in internal coordinates by introducing the body-fixed frame x and using the inverses of the arbitrarily defined x(X) and y(x) transformations. In this section we show the relationship of our method to those that are generally used for this purpose in the literature, one those based on the X → y transformation^1,33,39 (s-vector formalism) as well as on the y → x transformation^1,40 (t-vector formalism).

The s-vector formalism provides the exact vibrational Hamiltonian for the full-dimensional vibrational problem. To derive the reduced-dimensional kinetic energy and the mass matrices using the X → y transformation, it is necessary to extend the set of variables to a full set of internal variables and first construct the full-dimensional G_y,vib (see eqn (3)). Then one can either calculate the full-dimensional M_y,vib by inversion, then delete rows and columns and invert again, or directly correct the block of G_y,vib corresponding to the kept variables (ref. 19, Appendix IX) to account for the effect of the constrained ones.

The t-vector formalism requires the definition of a body-fixed frame (e.g., Eckart frame⁴¹), which is rarely the absolute co-rotating frame. To compensate for the arising spurious rotation during the change of internal coordinates of the molecule, rotational coordinates need to be introduced and a rovibrational Hamiltonian has to be constructed to make possible the exact description of the vibrational problem. The reduced-dimensional Hamiltonian is directly obtained by using only those vibrational t-vectors which belong to non-constrained internal variables.

The method presented in this work is analogous to the t-vector formalism also in the sense that both are based on the dx/dy derivatives and allow the direct construction of a reduced-dimensional Hamiltonian. However, instead of introducing rotational variables, our method exploits the orthogonality of the mass-scaled translational and (instantaneous) rotational basis vectors (3N-component) to vibrational ones so that rotation and translation are exactly removed. This approach provides a pure vibrational kinetic energy expression equivalent with the one provided with the s-vector formalism.

It can be shown that the inverse mass matrix G_y,vib in the s-vector formalism (eqn (3)) is the inverse of the M_y,vib mass matrix in eqn (22) and that this does not hold for the reduced-dimensional case as the pure-vibrational infinitesimal mass-scaled Cartesian displacement vectors corresponding to the various internal coordinates are not orthogonal to each other in general. This also implies that the s-vector formalism in full internal dimensionality is in fact a reduced-dimensional approximation, because it considers only vibrations and translational and rotational coordinates are simply omitted. This raises the question how it can be exact for vibrations. The reason for this is the inherent orthogonality of vibrations to rotations and translations in mass-scaled infinitesimal displacement space.

2.5 Normal mode analysis in internal coordinates

During normal mode analysis the vibration of the molecule is approximately decomposed into independent harmonic oscillators using the harmonic approximation to the kinetic and potential energy expressions. The harmonic approximation to the Lagrangian in internal coordinates in the neighborhood of a stationary point y₀ can be obtained by replacing the M_y,vib(y) matrix function with its value at y₀, M_y,vib,0 = M_y,vib(y₀), approximating the potential energy to second order around y₀, and setting its zero level to that at geometry y₀:Here F_y,0 = V_y′′(y₀) is the force constant matrix. This generally proves to be a good approximation to the Lagrangian at low energies in semi-rigid molecules, where no internal rotations or other large amplitude motion can take place.

From here on one can follow the standard procedure of normal mode analysis. After introducing mass-scaled vectors of deformation ỹ = M_y,vib,0^1/2(y − y₀) and velocity and the F̃_y,0 = G_y,vib,0^1/2F_y,0G_y,vib,0^1/2 mass-scaled force constant matrix (where G_y,vib,0 = M_y,vib,0⁻¹) one solves its F̃_y,0U = UΛ eigenproblem. Matrix F̃_y,0 is symmetric, thus all n eigenvalues λ_i (in matrix Λ = diag(λ₁,…,λ_n)) are real and the eigenvectors (i.e. columns of U) can be chosen to be orthonormal. If y₀ is a potential minimum, all eigenvalues are positive, whereas for a k^th-order saddle point k of them are negative. The vector of normal-mode deformation coordinates Q = (Q₁,…,Q_n)^T and the canonically conjugate momentum vector P = (P₁,…,P_n)^T (in fact, the normal-mode velocity vector Q̇) are defined as:

Q = U^Tỹ = U^TM_y,vib,0^1/2(y − y₀), (27)

P = Q̇ = U^Tỹ = U^TM_y,vib,0^1/2ẏ. (28)

In normal coordinates both the Lagrangian and the Hamiltonian (i.e. energy) decompose into sums of n harmonic oscillator (HO) terms:The normal mode frequencies can be obtained as ω_i = λ_i^1/2.

2.6 Normal mode sampling

The purpose of normal mode sampling is the generation of a set of classical states corresponding to a preselected vibrational state of a reactant molecule that will serve as initial conditions for collision or intramolecular trajectories. In QCT calculations, the initial states of molecules are usually generated by assuming that vibration and rotation are separable and the vibrations are well described by the normal mode approximation. The procedure for normal mode sampling is well known for full-dimensional models, but without the proper RD normal mode analysis, it cannot be used in reduced dimensionality. In normal mode sampling it is assumed that the vibrating molecule can be well approximated as a set of independent normal oscillators with n distinct quantum numbers v₁,…,v_n (n constants of motion). In quasiclassical quantization, normal coordinate and momentum amplitudes (Q_i,max, P_i,max = Q̇_i,max) of each normal mode oscillator are set so that the energy of the oscillator (E^1D_harm,i) matches that of the corresponding quantum harmonic oscillator in a given v_i quantum state:For the generation of classical states (Q,P), a random initial phase φ_i,0 is selected from a uniform distribution in the [0,2π) interval for each vibrational mode, and then the normal mode coordinate and velocity of normal oscillator i can be calculated according to:

Q_i = Q_i,max cos φ_i,0, (31)

P_i = Q̇_i = −Q̇_i,max sin φ_i,0. (32)

The corresponding classical state (y,ẏ) in internal coordinates can be obtained by inverting eqn (27) and (28):

y = y₀ + G_y,vib,0^1/2UQ, (33)

ẏ = G_y,vib,0^1/2UP = G_y,vib,0^1/2UQ̇. (34)

The resulting classical state (y,ẏ) is the appropriate initial state if one intends to follow the time evolution of the system using the Euler–Lagrange or Newton equations of motion. If one would like to describe the same dynamics in the more convenient Hamiltonian formalism (see Section 2.7.1), then the initial state should be expressed in internal coordinates and the conjugate momenta (y,p_y), where the momenta are obtained from eqn (23) with the mass matrix M_y,vib(y) calculated at the instantaneous geometry y.

When the trajectory integration is to be performed in Cartesian coordinates, the state (y,ẏ) needs to be transformed to the body-fixed Cartesian frame and one should ensure that the unphysical translation and rotation generated during transformations in eqn (5) and (6) are removed. To obtain the pure vibrational classical state (x,p_x,vib), the body-fixed Cartesian coordinates x are calculated from y using eqn (5); the corresponding Cartesian momenta describing vibration only are found by transforming ẏ to ẋ with matrix C(y) and finally to p_x,vib with matrix M_vib(x) using eqn (6) and (20), respectively:

p_x,vib = M_vib(x)ẋ = M_vib(x(y))C(y)ẏ. (35)

The point here is that the geometry-dependent mass matrix M_vib(x), which assigns mass only to vibrations, guarantees that the resulting classical states (x,p_x,vib) will have zero center-of-mass momentum and angular momentum. The coordinate and velocity conversion equations also ensure that Cartesian coordinates x, velocities ẋ and momenta p_x,vib fulfill the equations of geometrical constraints (y_n+1(x) = y_n+1,0,…,y_f(x) = y_f,0) and their time derivatives (see Section 2.7.2).

The ensemble of classical states (x,p_x,vib) in Cartesian coordinates, which corresponds to the pre-selected vibrational quantum state of the reduced-dimensional model of the reactant molecule is generated by carrying out the sampling procedure many times with different random phases in eqn (31) and (32). It should be noted that the ensemble obtained this way is not monoenergetic. In QCT calculations the nonrotating molecules are randomly oriented before collision. If a rovibrational state is to be generated, then the molecule's angular momentum is also set to that of the desired quantum state. The obtained classical states are in laboratory Cartesian frame and thereby they can be used in the definition of initial conditions of the molecule (X,P_X), for QCT calculations.

The standard method to prepare monoenergetic ro-vibrating ensembles by iterative rescaling of deformations and momenta and angular momentum vector adjustment after normal mode sampling has been proposed by Hase and coworkers.⁴² That method can be generalized to reduced-dimensional models. What remains the same is that the scaling factor is determined in the lab frame. The important difference is that in the full-dimensional model the scaling factor can be applied to scale the lab-frame deformations and momenta. In RD models the rescaling step needs to be performed in internal coordinates, because otherwise it would violate the constraints. It should be noted that the ensembles generated this way are not stationary when allowed to evolve on the real, anharmonic PES.²¹ Monoenergetic, stationary and accurately quantized ensembles of classical states representing a rovibrational quantum state can be generated by applying the adiabatic principle of classical mechanics. A generalized version of the adiabatic switching method, which accurately includes anharmonicity and coupling of vibrations, has recently been shown to perform well for polyatomic molecules.⁴³ The method has also been extended to generate ensembles corresponding to rovibrational quantum states.

2.7 Dynamics in reduced dimensionality

In reduced-dimensional classical trajectory calculations the equations of motion for the intramolecular motion of molecules have to guarantee the fulfillment of the constraints prescribed by the model. In the following, we discuss two choices: the application of equations of motion in internal coordinates and the integration of the equations of motion in lab Cartesian frame, supplemented by constraint forces. Integration in internal coordinates is more appropriate for the description of a vibrating molecule, whereas the equations presented for integration in 3N Cartesians are equally applicable both to pure bound motion and to scattering problems.

2.7.1 Equations of motion in internal coordinates. Hamiltonian equations of the pure vibrational motion in internal coordinates can be derived from the Hamiltonian in eqn (25).Rank-3 tensor G_y,vib′(y) and vector V_y′(y) are the coordinate derivatives of the inverse mass matrix G_y,vib(y) and the potential energy V_y(y). The initial conditions of the motion can be obtained viaeqn (33) and (34) and formula p_y = M_y,vib(y)ẏ. When vibrational dynamics are simulated, internal to body-fixed frame transformation and projection onto the vibrational subspace have to be performed several times at every integration step for the evaluation of G_y,vib(y), G_y,vib′(y) and V_y′(y) (see eqn (5), (6), (22) and (24)). The generally involved calculation of the rank-3 G_y,vib′(y) derivative tensor can be sped up by using some analytical transformations (see ESI†). A more serious drawback of using internal coordinates is that they may become indeterminate and in that neighborhood they change very steeply. At these points the coordinate transformation in eqn (5) or its inverse has singularities, which can cause significant numerical errors during the solution of the equations of motion. For example, such a singularity arises for 3D polar coordinates (r, θ, φ) at small (θ = 0°) and large (θ = 180°) polar angles where φ may change very quickly during dynamics. With careful selection of the body-fixed frame and the set of internal coordinates one can achieve that the singularities of the transformation equations (eqn (5) and (6)) will be at highly deformed geometries which are not sampled by normal mode state preparation and not visited during the vibration of the molecule. It is important to emphasize that eqn (36) and (37) describe the vibration of non-rotating molecules only.

2.7.2 Equations of motion with constraints in space-fixed Cartesian coordinates. During the internal motion of floppy molecules, as well as in loose clusters and chemical reactions, where bonds get broken and formed, highly-deformed geometries are visited and the moieties of reactant molecules can have arbitrary orientation with respect to each other. For such systems, it is preferable to simulate the reduced-dimensional dynamics by integrating the equations of motion in the full set of 3N lab Cartesian coordinates (X) and enforcing the restrictions in the internal degrees of freedom by constraint forces. To derive constraint forces, we introduce functions g_i(X) (i = 1,…,f − n) as the deviations of functions y_n+i(X) from their corresponding frozen values y_n+i,0. The equations of constraints are obtained by equating the g_i(X) functions to zero.

g_i(X) := y_n+i(X) − y_n+i,0 = 0 i = 1,…,f − n. (38)

These constraints are holonomic as they depend only on the position coordinates (but not on their time derivatives) and are scleronomic as they do not depend on time explicitly. They are expected to be fulfilled all along a trajectory, implying that their time derivatives must also be zero:

ġ_i(X,Ẋ) = ∇^Tg_iẊ = ∇^Tg_iM⁻¹P_X = 0 i = 1,…,f − n, (39)

which serve as constraint equations for the velocities Ẋ and momenta P_X. Consequently, a mechanical state, which is fully characterized either by (X,Ẋ) or (X,P_X), should fulfill eqn (38) and (39) simultaneously. The 3N-component constraint force arising from constraint i in eqn (39) is necessarily parallel to gradien t∇g_i, because it confines the allowed motion (velocity Ẋ) to a (3N − 1)-dimensional surface orthogonal to ∇g_i. Consequently, each constraint force F^constr_i is expressed as being proportional to ∇g_i which, multiplied by Lagrange multipliers λ_i (i = 1,…,f − n) are added to Hamilton's equations for momenta, supplementing there the potential forces:The Hamiltonian equations of motion for the coordinates remain the same (Ẋ = M⁻¹P_X) since the constraints in eqn (38) are holonomic. The Lagrange multipliers need to be calculated at every time step of trajectory integration via a set of linear equations. An excellent description of how the Lagrange multipliers are determined in practice can be found in ref. 3, 4 and 20. For completeness, the equations with the present notations are summarized in the ESI.†

2.7.3 Comparison of the computational aspects of the two descriptions. Usually the most expensive part of trajectory simulations is the evaluation of potential gradients. If no analytical gradients of the potential energy are available, then they need to be evaluated numerically with finite difference formulae, whose computational cost scales linearly with the number of gradient components. When nonredundant internal coordinates are used, some computational savings come from the reduced number (f < 3N − 6) of gradient components in the V_y′(y) potential gradient compared to the full-dimensional Cartesian problem (3N). On the other hand, simulation of the dynamics in irredundant internal coordinates requires repeated evaluation of rank-3 tensor G_y,vib′(y), which can be expensive unless analytical first and second derivatives (matrices C(y) and C′(y)) of the coordinate transformation x(y) are available. Analytical derivation of matrix C can be complicated and requires non-negligible human effort even when computer algebra packages are employed, and even when the derivatives are available, their evaluation may require numerous algebraic operations.

In the alternative method, integration of the equations of motion in Cartesian coordinates under the control of constraints, at least 6 more potential gradient components need to be evaluated; in addition, the constraint forces need to be determined. Yet, the application of Cartesian coordinates together with constraints can be overall cheaper than using internal coordinates, especially when analytic first derivatives of the PES, and first and second analytic derivatives of the constraints are available.

An additional aspect is the accuracy to which the constraints prescribed by the reduced-dimensional model are fulfilled. Integrating in internal coordinates automatically satisfies these constraints. On the other hand, when full-dimensional Cartesian coordinates are used, the constraints are enforced numerically, so that their fulfillment (eqn (38) and (39)) depends on the accuracy of the numerical procedure. Accordingly, in long simulations it is desirable to check regularly how well the constraints are met. If any of them is violated significantly, fulfillment of them can be reestablished by minimizing to zero the sum of the properly scaled squares of the left hand sides of eqn (38) and (39).

2.8 Assignment of normal mode quantum numbers to reduced-dimensional classical states

Solution of the inverse problem, determination of normal mode quantum numbers corresponding to a reduced-dimensional classical state of the molecule is needed when product-state resolved properties are calculated for bimolecular reactions or when the classical evolution of vibrationally excited states of molecules is of interest. Depending on which of the two sets of Hamiltonian equations of motion (Section 2.7) were used for trajectory integration, the classical state of the system is provided either in nonredundant internal coordinates and momenta (y,p_y) or in space-fixed Cartesian coordinates and momenta (X,P_X) fulfilling constraints in eqn (38) and (39).

In the former case, since the normal mode analysis of the reduced-dimensional model is done in internal coordinates (see Section 2.5), the state expressed as (y,ẏ = G_y,vib(y)p_y) can be directly used for quantum number assignment as described below. However, when the state is given in Cartesians, X and P_X need to be converted to internal coordinates and velocities. Coordinates y can be computed directly from the X coordinates, which fulfill the f − n constraints given in eqn (38), using the definition of the internal coordinates y(X). Velocities of internal coordinates can be obtained by either numerical derivation along the trajectory or analytical differentiation utilizing the reduced-dimensional variant of Wilson's B matrix (see eqn (2)) and the Cartesian velocities Ẋ = M⁻¹P_X:

It is important to point out that, while the reduced-dimensional B matrix cannot be used for the construction of G_y,vib due to the non-orthogonality of internal coordinates (see Section 2.4), it can be applied to transform Cartesian velocities to internal coordinate velocities if and only if they fulfill the constraints.

Once the state is given as (y,ẏ), the normal mode coordinates and momenta can be calculated by inverting eqn (33) and (34):

Q = U^TM_y,vib,0^1/2(y − y₀), (42)

P = Q̇ = U^TM_y,vib,0^1/2ẏ. (43)

The E^1D_harm,i mode energies can be calculated from the instantaneous normal coordinates and momenta viaeqn (29). The quantum numbers v_i are obtained from the energy corresponding to each vibrational mode using the inverse of the harmonic oscillator quantization rule (eqn (30)).

Note that the vibrational quantum numbers ν_i obtained this way are in general not integer numbers and because of this, the assignment of a quantum state to the given classical one is not unequivocal, and various tricks are generally used to do it (see, e.g.ref. 44 and 45).

3. Results

3.1 Three reduced-dimensional models of methane

In this section we apply the methods described above to normal mode analysis and sampling in various reduced-dimensional models of methane used in ref. 31 to study RD QCT reaction dynamics calculations of the reaction, which is the most complicated type of reaction for which exact RD quantum mechanical scattering calculations are available.^13,46,47 The RD model with the fewest restrictions is the one proposed by Palma and Clary (in the following, PC):^9,10 the only constraint is that the CZ₃ group keeps C_3v symmetry. The model and a series of its versions with additional constraints have been used in numerous quantum scattering simulations. In order to assess the performance of the RD models, Vikár et al.³¹ compared the results of (quasi)classical trajectory simulations obtained with the Palma–Clary and the full-dimensional models of the CH₄ + H → CH₃ + H₂ reaction with various CH₄ isotopologs. We just mention in passing that the PC model was found to give good results in many cases, but its performance proved to depend on the quality of the potential energy surface and the mass combination.³¹

Starting from the PC model, one can design a hierarchy of reduced-dimensional models for methane, in each of which the geometry of the CZ₃ moiety is constrained to maintain C_3v point group symmetry and by freezing additional types of motion, its internal degrees of freedom are reduced from f = 9 to 5, 4 or 3. Since the three Z atoms are treated as equivalent, such models are suitable for the description of methane isotopologs CH_nD_4−n with n = 0, 1, 3, 4 but not for CH₂D₂. In Fig. 3 the models with their reduced sets of internal variables and the attached body-fixed Cartesian frame are shown.

Fig. 3 The reduced-dimensional (a) Palma–Clary 5D (PC-5D) and (b) 4D (PC-4D) as well as (c) the rotating bond umbrella (RBU-3D) models of CZ₃Y molecule (composed of C, Y and identical Z¹, Z², Z³ atoms) with their respective rectilinear internal variables printed in blue. Point Z₃ denotes the geometrical center of the three Z atoms. In all the three models the CZ₃ moiety is constrained to C_3v symmetry. Internal variables in orange show the additional constraints in the PC-4D and RBU-3D models compared to the PC-5D model.

The latter was defined for all models in the same way: the origin was placed at the carbon atom and the axes were determined by the constrained CZ₃ moiety: z axis is along the C_3v axis, x axis within the YCZ¹ plane, forming an acute angle with the CZ¹ vector, and the y axis is chosen to ensure that the coordinate system is right-handed.

In the “full” Palma–Clary model⁹ (Fig. 3a) the only constraint is that the CZ₃ group is restricted to C_3v symmetry, resulting in 5 degrees of freedom for CZ₃Y; we abbreviate this model as PC-5D. In the 4D model (called PC-4D, Fig. 3b) considered here, the additional constraint is that the C–Z bond length is frozen at the equilibrium value (r_eq) in CZ₃Y as in a frequently applied restricted variant of the Palma–Clary model.^9,11,16 In the 3D model (Fig. 3c) we study, the motion of the Y atom is also constrained to one of the symmetry planes (for example in this case to the y plane) of the CZ₃ group. This is identical to the three-dimensional methane model in the rotating bond umbrella (RBU) method^14,48,49 and will be referred to as RBU-3D. The coordinates as well as the constraint equations used in the internal coordinate and Cartesian representations of the models are shown in Table 1. The common in the applied constraints is that all are defined by sums of dot products, which are very practical because their first and second derivatives have simple analytical forms as shown in the ESI.†

Table 1 Properties of three reduced-dimensional models (PC-5D, PC-4D, RBU-3D) of CZ₃Y and the lists of internal coordinates and Cartesian constraints used for simulating vibrational dynamics. For PC-4D and RBU-3D models simpler alternatives of the 3rd and 4th constraints are shown in the brackets. Vector R_AB denotes difference vector R_B − R_A, which points from atom A to B

---

In the following we report on normal mode analysis and the time evolution of ensembles of classical states of methane isotopologs CH_nD_4−n (n = 0, 1, 3, 4) generated by RD NMS using the three models, performed on the ZBB3 PES.⁵⁰ At the T_d equilibrium geometry the C–H bond length is 2.0579 a₀.

3.2 Reduced-dimensional normal mode analysis for methane isotopologs

The frequencies and the nature of normal modes obtained from harmonic vibrational analysis using the procedure described in Section 2.5 are shown in Table 2; the corresponding displacement vectors can be seen in Fig. 4. Vibrations of full-dimensional CH₄ and CD₄ are characterized by four frequencies, corresponding to one nondegenerate fully symmetric stretch mode (A₁ irreducible representation), a pair of doubly degenerate deformation (E) modes, a triply degenerate stretch (T₂) and a triply degenerate deformation (T₂). Symmetry reduction from T_d to C_3v by distinguishing one of the ligand atoms in CZ₄ from the other three, but not changing its mass and the forces it feels, decomposes each of the T₂ modes into a totally symmetric A₁ and a doubly degenerate E symmetry mode (according to T₂ = A₁ ⊕ E) while leaves their frequencies unchanged. The arising A₁ stretching mode will be an antisymmetric combination of the C–Y and the symmetric C–Z₃ stretching modes, while the new A₁ deformation will be the umbrella mode, which is the symmetric deformational mode of the CZ₃ moiety. The original totally symmetric A₁ stretching mode will preserve its character and frequency.

Table 2 The irreducible representation (irrep), degeneracy, type and frequency of vibrational normal modes of four methane isotopologs determined for the full-dimensional (FD) and the PC-5D, PC-4D and RBU-3D reduced-dimensional (RD) models on the ZBB3 PES

Isotopolog, sym. group	Irrep : degeneracy				Type				Harmonic frequency/cm^−1
	FD	5D	4D	3D	FD	5D	4D	3D	FD	5D	4D	3D
H3CH	T2 : 3	A1 : 1	1	1	Deg. def.	Umbrella	Umbrella	Umbrella	1335.6	1335.6	1336.7	1336.7
FD : Td	E : 2	E : 2	2	1	Deg. def.	Rocking	Rocking	Rocking	1547.8	1436.4	1436.3	1436.3
RD : C3v	A1 : 1	A1 : 1	—	—	Sym. stre.	Sym. stre.	—	—	3026.1	3026.2	—	—
	T2 : 3	A1 : 1	1	1	Deg. stre.	Antisym.stre.	CH stre.	CH stre.	3166.7	3166.8	3129.4	3129.4
D3CH	A1 : 1	A1 : 1	1	1	CD3 sym. def.	Umbrella	Umbrella	Umbrella	1018.1	1018.1	1018.1	1018.1
FD : C3v	E : 2	—	—	—	CD3 deg. def.	—	—	—	1045.7	—	—	—
RD : C3v	E : 2	E : 2	2	1	CD3 rock.	Rocking	Rocking	Rocking	1305.7	1298.4	1298.3	1298.3
	A1 : 1	A1 : 1	—	—	CD3 sym. stre.	CD3 stre.	—	—	2186.5	2186.5	—	—
	E : 2	—	—	—	CD3 deg. stre.	—	—	—	2343.2	—	—	—
	A1 : 1	A1 : 1	1	1	CH stre.	CH stre.	CH stre.	CH stre.	3136.1	3136.2	3125.9	3125.9
H3CD	E : 2	E : 2	2	1	CH3 rock.	Rocking	Rocking	Rocking	1178.3	1218.1	1218.0	1218.0
FD : C3v	A1 : 1	A1 : 1	1	1	CH3 sym. def.	Umbrella	Umbrella	Umbrella	1330.5	1330.5	1331.8	1331.8
RD : C3v	E : 2	—	—	—	CH3 deg. def.	—	—	—	1487.9	—	—	—
	A1 : 1	A1 : 1	1	1	CD stre.	CD stre.	CD stre.	CD stre.	2286.7	2286.7	2290.5	2290.5
	A1 : 1	A1 : 1	—	—	CH3 sym. stre.	CH3 stre.	—	—	3067.6	3067.7	—	—
	E : 2	—	—	—	CH3 deg. stre.	—	—	—	3166.4	—	—	—
D3CD	T2 : 3	A1 : 1	1	1	Deg. def.	Umbrella	Umbrella	Umbrella	1010.3	1010.3	1010.4	1010.4
FD : Td	E : 2	E : 2	2	1	Deg. def.	Rocking	Rocking	Rocking	1094.9	1051.8	1051.8	1051.8
RD : C3v	A1 : 1	A1 : 1	—	—	Sym. stre.	Sym. stre.	—	—	2140.6	2140.6	—	—
	T2 : 3	A1 : 1	1	1	Deg. stre.	Antisym.stre.	CD stre.	CD stre.	2343.1	2343.1	2286.3	2286.3

---

Fig. 4 Atomic displacements in vibrational normal modes of the CZ₃Y molecule with Z = Y = H, determined using the PC-5D, PC-4D and RBU-3D reduced-dimensional models. The Y atom points upward. The types and the harmonic frequencies of the vibrational modes are also shown.

When symmetry reduction in CZ₄ methane isotopologs (i.e. Z = Y, top and bottom rows of Table 2) is caused by constraining the symmetry of the CZ₃ group to C_3v (leading to the PC-5D model), the fully symmetric stretch mode (involving all four C–Z bonds) remains again intact, whereas the A₁ modes resulting from symmetry reduction (the antisymmetric stretch and the umbrella bend) together with their frequencies remain essentially unchanged while the asymmetric CZ₃ stretch modes of E symmetry disappear. The original and the resulting doubly-degenerate E-symmetry deformational modes will mix to conform to the constrained C_3v symmetry and result in a doubly-degenerate E-symmetry rocking mode with an intermediate frequency. The two combinations violating the C_3v symmetry disappear. In CZ₃Y methane isotopologs (i.e. Z ≠ Y), whose symmetry is reduced to C_3v due to the difference of isotopes Z and Y, the decomposition of the triply-degenerate modes is accompanied by frequency changes and mode mixing. The resulting antisymmetric combination of the local C–Y and symmetric C–Z₃ stretching modes of A₁ symmetry will mix with the totally symmetric stretching mode (present in CZ₄) and gives essentially local A₁-symmetry C–Z₃ and C–Y stretching modes with significantly different frequencies. Furthermore, a pair of doubly-degenerate rocking modes appears due to mode mixing. When the symmetry of the CZ₃ group is constrained to C_3v in mixed methane isotopologs according to the PC-5D model, the doubly-degenerate bending and stretching modes which violate the C_3v symmetry disappear. The A₁-symmetry modes (umbrella and two local stretches) remain essentially intact, while the frequency of the doubly-degenerate rocking mode changes slightly.

Further reduction of dimensionality to four by freezing the C–Z bond lengths of the CZ₃ moiety (in the PC-4D model) in CZ₄ isotopologs causes the mixing of the symmetric and antisymmetric stretches and results in a local C–Z stretch, whereas in CZ₃Y isotopologs removes the C–Z₃ stretch, and the frequency of the remaining modes changes slightly. Confining the motion of the Y atom into one of the symmetry planes of CZ₃ (in the RBU-3D model) freezes one of the rocking vibrations, without causing any additional frequency change.

3.3 Evolution of reduced dimensional ensembles of classical states of methane isotopologs

The reduced-dimensional vibrational Hamiltonian in internal coordinates and the equations of motion derived in Section 2.7 allow us to monitor the time evolution of two parameters that proved to be important in the QCT simulations of the abstraction reaction between methane isotopologs (CZ₃Y) and H atoms,^21,31 namely, the C–Z and C–Y bond lengths.

A single classical state corresponding to the ground state of the PC-5D model of the CH₄ molecule (Z, Y = H) has been prepared by normal mode sampling (Section 2.6) and propagated in time. Fig. 5 shows the two bond lengths as a function of time up to 10⁴τ₀ (1τ₀ = atomic time unit ≈0.0242 fs). The equations of motion for the PC-5D model have been integrated in two sets of coordinates, internal (continuous lines) and space-fixed Cartesian (dotted lines, with the constraints being taken into account with the Lagrange multipliers method), yielding identical trajectories (see Fig. 5). One can see that the evolution of the two bond lengths is coupled. At the beginning of this specific trajectory, the amplitude of the bond length oscillation decreases for the C–Y bond and increases for the C–Z bond, indicating energy transfer between them. Later, the coherence of the energy exchange between the two oscillations is reduced, which indicates that other modes are also coupled to these key reaction dynamical parameters.

Fig. 5 Evolution of the C–Z₃ (upper panel) and C–Y (lower panel) bond lengths for a reduced-dimensional model of a single ground-state CH₄ molecule obtained by normal-mode sampling. The equations of motion corresponding to the PC-5D model were integrated using (i) the reduced set of internal coordinates (continuous grey lines) and (ii) Cartesian coordinates with constraints (dotted lines).

More important is the behavior of an ensemble of classical states that represents a quantum state. Ensembles of 10⁴ classical states corresponding to the ground state of CH₄ molecules were prepared by normal mode sampling as described in Section 2.6 for all four methane isotopologs and all three models. No deformation and momentum scaling was applied to make the ensembles monoenergetic. The evolution of the ensemble was simulated using the space-fixed Cartesian frame description with constraints, which were always fulfilled precisely, and the obtained ensemble average C–Y and C–Z bond lengths and normal mode quantum numbers are plotted in Fig. 6 and 7. The oscillations of the average bond length plots indicate the same kind of coherence in the evolution of the members of the ensemble as was reported in ref. 21 for the FD evolution of methane isotopologs. The phenomenon is caused in both cases by the failure of the normal mode approximation. In QCT simulations of the CH₄ + H → CH₃ + H₂ reaction and its isotopologs,²¹ the temporal oscillation of the average C–Y bond length of the ensemble of methane molecules prepared with normal mode sampling was shown to give rise to “spatial” oscillation of the calculated reaction probabilities and cross sections as a function of initial separation of the reactants. The conversion of the temporal oscillation to spatial can be understood by considering that the ensemble evolves during the initial “free” flight of the reactants as it adapts to the anharmonic PES.

Fig. 6 Evolution of the ensemble average C–Y (black line) and C–Z₃ (red line) bond lengths (equilibrium value: 2.0579a₀) for ground-state methane isotopologs CH_nD_4−n (n = 0, 1, 3, 4) within the PC-5D, PC-4D and RBU-3D models. The initial states of the 10⁴-member ensembles were generated by normal mode sampling.

Fig. 7 Evolution of the ensemble average normal mode quantum numbers for four ground-state methane isotopologs within the PC-5D, PC-4D and RBU-3D reduced-dimensional models. The initial states of the 10⁴-member ensembles were generated by normal mode sampling.

The relaxation dynamics of the methane ensembles generated by normal-mode sampling on the anharmonic PES provides information on the couplings and energy flow between vibrational modes. The coarse-grained pattern of the evolution of the two mean bond lengths in the CZ₄ isotopologs indicate initially periodic energy exchange between the two stretching modes of the PC-5D model. The other two models also show beating in their ensemble average bond length oscillation, which suggest significant energy exchange with the bending modes. Eventually, the amplitude of the oscillations decreases as the ensemble adapts to the anharmonic PES and the phases of stretching vibration of the individual trajectories decohere.

The lower the dimensionality of the model, the slower is the decay of coherence in the stretching oscillations, which is consistent with the lower-dimensional, simpler structure of the phase space and suggests that the chaotic character of the vibration is reduced. This is also reflected in the coarse-grained structure of oscillations, which shows regularity and is symmetric with respect to a mean value in the case of all investigated models and isotopologs, except for the evolution of CH₄ and CD₄ in the PC-5D model.

More informative on the intramolecular energy transfer is the evolution of the ensemble average normal mode quantum numbers (see Fig. 7), which also represents a stringent test of the applicability of our method to assigning a vibrational quantum state to the molecule in motion. For the doubly degenerate rocking mode present in the PC-5D and PC-4D models, the sum of the two ensemble averaged quantum numbers was taken, as the degenerate modes are kinematically strongly coupled, thus their individual values are not meaningful. Similarly to the full-dimensional case,²¹ long-term oscillations can be observed and the patterns hint at significant mode-to-mode energy transfer. The initial value of each quantum number is zero at the NMS preparation, which immediately changes when the atoms are allowed to move.

The pattern of oscillation is in agreement with that of the C–Z and C–Y bond lengths. In the PC-5D model, the CZ₄ isotopologs show strong coupling between the symmetric stretch and the rocking modes. In the CZ₃Y isotopologs stronger coupling emerges between the “internal” modes of the CZ₃ group, i.e. the bond stretching and the umbrella mode. It is obvious, however, that in all cases every mode participates in the energy exchange. In the PC-4D and RBU-3D models where the stretching of C–Z bonds is frozen, intense energy exchange can be observed between the rocking mode(s) and the C–Y stretch modes for CZ₄ isotopologs, which is a consequence of Coriolis coupling. Furthermore, these two models in the case of isotopolog CH₃D show a very long-period (80–160 × 10³) beating-like energy exchange between the C–Y stretch and the umbrella mode. The quantum numbers obtained with our formalism show the expected pattern, confirming the applicability of the formulas.

4. Conclusions

We have presented a formal derivation of the pure vibrational classical Hamiltonian in a nonredundant, but not necessarily full set of internal coordinates, which allowed the extension of the quasiclassical trajectory method to constrained, reduced-dimensional systems. The formalism allows one to carry out harmonic vibrational analysis in any reduced set of coordinates, and generate classical states corresponding to a given quantum state of the molecule using normal mode sampling (the “direct problem”), to simulate dynamics and to find normal mode quantum numbers when a classical state is given, for example, at the end of a QCT simulation (the “inverse problem”). For the implementation of the method, the only system specific formulae need to be provided are the definition of the internal coordinates of the chosen model and their connection to a suitable body-fixed Cartesian frame. Once these are provided, the method works as a black box, because, instead of deriving analytical formulae, the inverse mass matrix and the kinetic energy expression of the vibrational Hamiltonian is constructed numerically. The formalism is universal in the sense that it can be applied to both full- and reduced-dimensional models and that the internal and body-fixed Cartesian coordinate systems can be selected arbitrarily, so it can take into account any geometric constraint. The applicability of the method has been demonstrated in a previous study addressing the comparison of RD and FD models in classical simulations.³¹

Semiclassical methods, such as adiabatic switching^43,51,52 (AS) and semiclassical initial value representation^53–55 (SC-IVR), based on classically propagated trajectories can also employ the equations derived here, thus the method can be used for the quantization of rovibrational levels of constrained systems with a better computational scaling than quantum mechanical methods.

Finally, we mention that the RD method can enable one to perform QCT calculations when computational complications arise. One such difficulty is that the high computational cost of electronic structure calculations often does not allow the development of a full-dimensional analytic PES, which is required when long and/or a large number of trajectories need to be simulated in a QCT study. Development of a reduced-dimensional analytic PES combined with the theory presented here can provide a means for the dynamical investigation of such systems. Another computational issue is the undesirable leakage of zero-point energy deposited in each vibrational mode to other modes.

Conflicts of interest

There are no conflicts to declare.

Appendix

A. Orthogonality of translational and rotational basis vectors

Translational basis vectors u^trans,α(x) defined in eqn (14) are orthogonal to each other as their dot products (i.e.u^trans,α,Tu^trans,β, where α = x, y, z and β = x, y, z) give the corresponding (αβ) elements of the unit matrix times the total mass (M).The m_i is the mass of atom i, where i = 1,…,N. δ_αγ is the Kronecker symbol. If α = β the product gives the squared norm.

Rotational basis vectors u^rot,α(x) (see eqn (15)) are orthogonal to each other as their dot products give the corresponding elements of the moment of inertia tensor (Θ_αα) defined in the principal axis (PA) frame, which is by definition a diagonal matrix:

In eqn (A2) vectors e^PA_α are the orthonormal unit vectors along the principal axis α (α = 1, 2, 3) expressed in the body-fixed Cartesian frame and ρ_i := r_i − r_CM are the instantaneous position vector of atom i (i = 1,…,N) relative to the center of mass of the molecule, r_CM. Vector component ρ^PA_iα denotes the α component of the coordinate vector of atom i from the center of mass in PA Cartesian frame. ε_γστ is the Levi-Civita symbol. It is assumed that the three orthogonal e^PA_α (α = 1, 2, 3) vectors form a right-handed system . In the transformations, it was exploited that:

Rotational basis vectors are orthogonal to translational basis vectors regardless of the definition of vectors e^PA_α, because their product is a sum of terms that are proportional to the components of the center of mass position vector.

Acknowledgements

We thank Professors J. M. Bowman for making the potential surface codes available for us. We thank Professor A. G. Császár, Dr V. Szalay and Dr T. Szidarovszky for helpful discussions. This work has been supported by NKFIH, Hungary, (Grant No. K108966 and PD120776 (T. N.)), by the Government of Hungary and the EU via grant VEKOP-2.3.2-16-2017-00013 (G. L.), by COST Action 1401 and by the János Bolyai Research Fellowship (Grant No. BO/00279/16/7 (T.N.)).

Notes and references

1. E. Mátyus, G. Czakó and A. G. Császár, J. Chem. Phys., 2009, 130, 134112.
2. C. Fábri, E. Mátyus and A. G. Császár, J. Chem. Phys., 2011, 134, 74105.
3. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
4. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren and J. Hermans, in Intermolecular Forces. The Jerusalem Symposia on Quantum Chemistry and Biochemistry., ed. B. Pullman, Springer, Dordrecht, 1981, vol. 14, pp. 331–342.
5. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
6. J. M. Bowman and A. F. Wagner, in The Theory of Chemical Reaction Dynamics, ed. D. C. Clary, Springer Netherlands, Dordrecht, 1986, pp. 47–76.
7. Q. Sun, D. L. Yang, N. S. Wang, J. M. Bowman and M. C. Lin, J. Chem. Phys., 1990, 93, 4730–4739.
8. D. C. Clary, J. Phys. Chem., 1994, 98, 10678–10688.
9. J. Palma and D. C. Clary, J. Chem. Phys., 2000, 112, 1859–1867.
10. J. Palma and D. C. Clary, J. Chem. Phys., 2001, 115, 2188–2197.
11. M. Yang, D. H. Zhang and S. Y. Lee, J. Chem. Phys., 2002, 117, 9539–9542.
12. H.-G. Yu and G. Nyman, Chem. Phys. Lett., 1998, 298, 27–35.
13. G. Schiffel and U. Manthe, J. Chem. Phys., 2010, 132, 84103.
14. G. Nyman, in Theory of Chemical Reaction Dynamics, NATO Science Series, ed. A. Laganà and G. Lendvay, Kluwer, Dordrecht, 2004, pp. 253–278.
15. D. H. Zhang, M. Yang, M. A. Collins and S. Lee, in Theory of Chemical Reaction Dynamics, NATO Science Series, ed. A. Laganà and G. Lendvay, Kluwer, Dordrecht, 2004, pp. 279–304.
16. W. Zhang, Y. Zhou, G. Wu, Y. Lu, H. Pan, B. Fu, Q. Shuai, L. Liu, S. Liu, L. Zhang, B. Jiang, D. Dai, S. Lee, Z. Z. Xie, B. J. Braams, J. M. Bowman, M. a. Collins, D. H. Zhang and X. Yang, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 12782–12785.
17. R. Welsch and U. Manthe, J. Chem. Phys., 2015, 142, 64309.
18. E. B. Wilson, J. Chem. Phys., 1939, 7, 1047–1052.
19. E. B. Wilson, J. J. C. Decius and P. C. Cross, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra, Dover Publications, New York, 1980.
20. L. M. Raff and D. L. Thompson, in Theory of Chemical Reaction Dynamics, ed. M. Baer, CRC Press, Boca Raton, Florida, 1985, vol III, pp. 1–122.
21. T. Nagy, A. Vikár and G. Lendvay, J. Chem. Phys., 2016, 144, 14104.
22. D. Lu and W. L. Hase, J. Chem. Phys., 1988, 89, 6723.
23. J.-J. Klossika and R. Schinke, J. Chem. Phys., 1999, 111, 5882–5896.
24. C. L. Stroud and L. M. Raff, J. Chem. Phys., 1980, 72, 5479–5488.
25. N. Sathyamurthy and L. M. Raff, J. Chem. Phys., 1980, 72, 3163–3178.
26. A. Faure, L. Wiesenfeld, M. Wernli and P. Valiron, J. Chem. Phys., 2005, 123, 104309.
27. A. I. Maergoiz, E. E. Nikitin, J. Troe and V. G. Ushakov, J. Chem. Phys., 2002, 117, 4201–4213.
28. A. Faure, C. Rist and P. Valiron, Astron. Astrophys., 1999, 348, 972–977.
29. L. B. Harding, S. J. Klippenstein and Y. Georgievskii, Proc. Combust. Inst., 2005, 30, 985–993.
30. L. B. Harding, Y. Georgievskii and S. J. Klippenstein, J. Phys. Chem. A, 2010, 114, 765–777.
31. A. Vikár, T. Nagy and G. Lendvay, J. Phys. Chem. A, 2016, 120, 5083–5093.
32. G. Fogarasi, X. Zhou, P. W. Taylor and P. Pulay, J. Am. Chem. Soc., 1992, 114, 8191–8201.
33. S. R. Polo, J. Chem. Phys., 1956, 24, 1133.
34. J. K. G. Watson, Mol. Phys., 1968, 15, 479–490.
35. C. S. Sloane and W. L. Hase, J. Chem. Phys., 1977, 66, 1523.
36. W. L. Hase, Encyclopedia of Computational Chemistry, John Wiley & Sons, Ltd,Chichester, UK, 2002.
37. W. H. Miller, N. C. Handy and J. E. Adams, J. Chem. Phys., 1980, 72, 99–112.
38. V. Szalay, J. Chem. Phys., 2014, 140, 234107.
39. P. C. Wilson, E. B. Decius, J. C. Cross, E. B. Wilson, J. C. Decius, P. C. Cross and B. R. Sundheim, Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra, 1955.
40. G. O. Sørensen, Large Amplitude Motion in Molecules II, Springer-Verlag, Berlin/Heidelberg, 1979, pp. 97–175.
41. C. Eckart, Phys. Rev., 1935, 47, 552–558.
42. W. L. Hase, Encyclopedia of Computational Chemistry, John Wiley & Sons, Ltd, Chichester, UK, 2002.
43. T. Nagy and G. Lendvay, J. Phys. Chem. Lett., 2017, 8, 4621–4626.
44. L. Bonnet and J. C. Rayez, Chem. Phys. Lett., 2004, 397, 106–109.
45. G. Czakó, J. Phys. Chem. A, 2012, 116, 7467–7473.
46. G. Schiffel and U. Manthe, J. Chem. Phys., 2010, 133, 174124.
47. R. Welsch and U. Manthe, J. Chem. Phys., 2012, 137, 244106.
48. H.-G. Yu and G. Nyman, J. Chem. Phys., 1999, 110, 7233.
49. H.-G. Yu and G. Nyman, Chem. Phys. Lett., 1999, 312, 585–590.
50. Z. Xie, J. M. Bowman and X. Zhang, J. Chem. Phys., 2006, 125, 133120.
51. R. T. Skodje, F. Borondo and W. P. Reinhardt, J. Chem. Phys., 1985, 82, 4611–4632.
52. C. Qu and J. M. Bowman, J. Phys. Chem. A, 2016, 120, 4988–4993.
53. A. L. Kaledin and W. H. Miller, J. Chem. Phys., 2003, 118, 7174.
54. A. L. Kaledin and W. H. Miller, J. Chem. Phys., 2003, 119, 3078–3084.
55. M. Ceotto, S. Atahan, S. Shim, G. F. Tantardini, A. Aspuru-Guzik, G. E. Scuseria and M. J. Frisch, Phys. Chem. Chem. Phys., 2009, 11, 3861.

---

Footnote

† Electronic supplementary information (ESI) available: Method of Lagrange multipliers for holonomic and scleronomic constraints. An efficient form of constraints allowing simple analytic gradients and Hessians. See DOI: 10.1039/c8cp01600c

---