# Eckart axis conditions, Gauss' principle of least constraint, and the optimal superposition of molecular structures

Gerald R. Kneller

► To cite this version:

Gerald R. Kneller. Eckart axis conditions, Gauss' principle of least constraint, and the optimal superposition of molecular structures. The Journal of Chemical Physics, 2008, 128 (19), pp.194101. 10.1063/1.2902290. hal-01378837

HAL Id: hal-01378837

https://hal.science/hal-01378837

Submitted on 10 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Eckart axis conditions, Gauss' principle of least constraint, and the optimal superposition of molecular structures

Gerald R. Kneller

Centre de Biophysique Moléculaire, CNRS, Rue Charles Sadron, 45071 Orléans, France, Université d'Orléans, Chateau de la Source-Av. du Parc Floral, 45067 Orléans, France, and Synchrotron Soleil, L'Orme de Merisiers, B.P. 48, 91192 Gif-sur-Yvette, France

The relation of the Eckart axis conditions for polyatomic vibrating molecules to the problem of optimal superposition of molecular structures has been pointed out recently [J. Chem. Phys. 122, 224105 (2005)]. Here, it is shown that both problems are intimately related to Gauss' principle of least constraint, for which a concise derivation is presented. In the context of this article, Gauss' principle leads to a rotational superposition problem of the unconstrained atomic displacements and the corresponding displacements due to a molecular rigid-body motion. The Eckart axis conditions appear here as necessary conditions for a minimum of the constraint function. The importance of Eckart's problem for extracting the internal motions of macromolecules from simulated molecular dynamics trajectories is pointed out, and it is shown how the case of coarse-grained sampled trajectories can be treated.

## I. INTRODUCTION

Although classical point mechanics can be considered as the oldest branch of theoretical physics, there are still applications to discover that concern less well known but very useful formulations of mechanical problems. An example is Gauss' principle of least constraint, an elegant reformulation of D'Alembert's principle that allows one to establish the equations of motion of systems of point particles under constraints in a very intuitive and general way.¹ Despite this fact, it has not found its way into many standard textbooks on classical mechanics. The possibility to consider also non-holonomic constraints in a simple way, which involve explicitly the velocities of the particles, was the reason for a rediscovery of Gauss' formulation of mechanics for molecular dynamics simulations in the 1980's. Its application allowed one to consider nonstandard situations, such as systems with constant kinetic energy or with a given shear flow, and promoted the development of nonequilibrium molecular dynamics simulations. For more details, the reader is referred to an article by Evans et al.² and to the monographs by Hoover.³,⁴

In this article, I will show an application of Gauss' principle that concerns Eckart's problem of following the internal motions in a polyatomic molecule in the gas phase.⁵ The idea is to subtract the rigid-body dynamics from the total dynamics of the molecule by using Gauss' least-squares principle for the atomic displacements of a constrained dynamical system. The method leads immediately to the Eckart conditions for the intramolecular atomic displacements and, in particular, to a superposition problem of molecular structures. This has been recently pointed out by Kudin and Dymarsky.⁶ In this paper, I will also briefly comment on the relation of Gauss' principle with the more familiar formula-

tion of mechanics by D'Alembert and on the justification of the latter from the point of view of modern linear algebra.

## II. GAUSS' MECHANICS

As Gauss pointed out in his famous article in the Journal für Reine und Angewandte Mathematik,¹ he did not claim having found new formulation of classical mechanics of constrained systems, but a useful reformulation of D'Alembert's principle. In the following, I give a derivation of Gauss' principle that discusses also a subtle point in the application of D'Alembert's principle for dynamical problems, which is apparently not mentioned in the literature. Consider a system of N pointlike particles with positions $x_α$ and masses $m_α$, which are subjected to s constraints of the form $h_j(x,t)=0$ or $g_j(x, \dot{x}, t)=0$ ($j=1, ..., s$), where x comprises all positions $x_α$ and $\dot{x}$ the corresponding velocities $\dot{x}_α$. A set of s linear constraints for the accelerations of the form

$$A\ddot{x} = b \qquad (1)$$

is obtained by differentiating the constraints $h_j(x,t)=0$ twice and the constraints $g_j(x, \dot{x}, t)=0$ once with respect to time. The matrix A in Eq. (1) has the dimensions $s \times 3N$ and b is an s-dimensional column vector. Both A and b depend, in general, on the positions and velocities of the particles.

Equation (1) may be considered as an underdetermined set of equations for the accelerations of the particles and a solution can be given in terms of the generalized inverse of A,⁷⁻⁹ which is usually denoted by $A^+$,

$$\ddot{x} = A^+b + \ddot{x}_{||}. \qquad (2)$$

Here, $\ddot{x}_{||}$ is an undetermined vector in the null space of A, fulfilling $A\ddot{x}_{||}=0$. The most important features of generalized inverses are that $AA^+$ and $A^+A$ are projectors on the spaces spanned by the columns and rows of A, respectively. If all
rows in **A** are linearly independent, $\mathbf{A}\mathbf{A}^+ = \mathbf{1}_s$ is the $s \times s$ unit matrix and $\mathbf{A}^+\mathbf{A} = \mathbf{A}^T(\mathbf{A}\mathbf{A}^T)^{-1}\mathbf{A} = \mathbf{P}_\perp$ is a projector on an $s$-dimensional subspace of $\mathbb{R}_{3N}$. The latter will be denoted by the symbol $\mathcal{V}_\perp$. The orthogonal complement, $\mathcal{V}_\parallel$, is the null space of **A** and the corresponding projector is $\mathbf{P}_\parallel = \mathbf{1} - \mathbf{P}_\perp$. With these preliminaries it follows from relation (2) by multiplication with $\mathbf{P}_\perp$ from the left that

$$ \mathbf{A}^{+}\mathbf{b} = \ddot{\mathbf{x}}_{\perp} \in \mathcal{V}_{\perp}. \qquad (3) $$

Relation (2) thus represents the decomposition of the $3N$-dimensional acceleration vector $\ddot{\mathbf{x}}$ into a known component in $\mathcal{V}_\perp$, which is entirely determined by the positions and the velocities of the particles, and an unknown component in $\mathcal{V}_\parallel$. The latter must be determined from Newton's equations of motion. Writing $\ddot{\mathbf{x}} = \ddot{\mathbf{x}}_\parallel + \ddot{\mathbf{x}}_\perp$, the latter read

$$ \mathbf{M}(\ddot{\mathbf{x}}_{\parallel} + \ddot{\mathbf{x}}_{\perp}) = \mathbf{f} + \mathbf{z}. \qquad (4) $$

Here, **f** and **z** contain the external forces and the constraint forces on the particles, respectively, and **M** is the diagonal $3N \times 3N$ mass matrix

$$ \mathbf{M} = \begin{pmatrix} m_1 & 1 & 0 & \dots & 0 \\ 0 & m_2 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & m_N & 1 \end{pmatrix}. \qquad (5) $$

It must be emphasized that Newton's equations of motion constitute a system of $3N$ linear equations for $\ddot{\mathbf{x}}_\parallel$ and $\mathbf{z}$. A unique solution for both unknowns exists if one requires that

$$ \mathbf{z} \perp \ddot{\mathbf{x}}_{\parallel}. \qquad (6) $$

The solution of linear systems of equations for two mutually orthogonal unknown vectors has been described by Bott and Duffin¹⁰ and, in a recent paper, it has been shown how the Bott-Duffin inverse can be used to establish equations of motions for constrained dynamical systems.¹¹,¹² It will now be demonstrated that condition (6) implies D'Alembert's and Gauss' formulation of constrained mechanics.

Let us consider all possible positions of the system at time $t + \delta t$ that are compatible with the given constraints. It follows by the Taylor expansion up to order $\delta t^2$ that

$$ \mathbf{x}(t + \delta t) = \mathbf{x}(t) + \delta t \dot{\mathbf{x}}(t) + \frac{\delta t^2}{2} \ddot{\mathbf{x}}(t). \qquad (7) $$

Since the positions and velocities at time $t$ are given, the possible variations of the positions at time $t + \delta t$ are determined by the possible variations of the accelerations. The latter are obtained by varying the component $\ddot{\mathbf{x}}_\parallel$ that is left undetermined by the constraints,

$$ \delta\ddot{\mathbf{x}}(t) = \delta\ddot{\mathbf{x}}_\parallel(t) \in \mathcal{V}_\parallel(t), \qquad (8) $$

and one may write

$$ \delta\mathbf{x}(t + \delta t) = \frac{\delta t^2}{2} \delta\ddot{\mathbf{x}}_\parallel(t). \qquad (9) $$

Note that the subspace $\mathcal{V}_\parallel(t)$ changes with time implicitly through the change of the positions and the velocities. It follows thus from Eqs. (6), (8), and (9) that

$$ \mathbf{z}^T(t) \delta\mathbf{x}(t + \delta t) = 0, \qquad (10) $$

where the superscript $T$ denotes a transposition. Condition (10) is nothing but D'Alembert's principle, which is usually written in the form $\sum_{\alpha} \mathbf{z}_{\alpha}^T \delta\mathbf{x}_{\alpha} = 0$. This notation masks, however, the fact that the constraint forces and the virtual displacements are, strictly speaking, not to be considered at the same time.

The idea of Gauss was to turn D'Alembert's principle into a true minimization problem, in which the positions of the particles at time $t + \delta t$ are the variables to be optimized. To derive the function to be minimized, one compares the evolution of the constrained positions with the unconstrained ones. Up to order $\delta t^2$ the latter are given by

$$ \mathbf{x}^{(0)}(t + \delta t) = \mathbf{x}(t) + \delta t \dot{\mathbf{x}}(t) + \frac{\delta t^2}{2} \mathbf{M}^{-1} \mathbf{f}(t), \qquad (11) $$

and one finds that

$$ \mathbf{z}(t) = \mathbf{M}(\mathbf{x}(t + \delta t) - \mathbf{x}^{(0)}(t + \delta t)). \qquad (12) $$

Since the unconstrained positions at time $t + \delta t$ are completely determined by the positions, the velocities, and the external forces at time $t$, there are no degrees of freedom for variations, $\delta\mathbf{x}^{(0)}(t + \delta t) = 0$, and one may write

$$ \delta\mathbf{x}(t + \delta t) = \delta(\mathbf{x}(t + \delta t) - \mathbf{x}^{(0)}(t + \delta t)). \qquad (13) $$

D'Alembert's principle (10) appears thus as the necessary condition for

$$ \frac{1}{2} \sum_{\alpha=1}^{N} m_{\alpha} (\mathbf{x}_{\alpha}(t + \delta t) - \mathbf{x}_{\alpha}^{(0)}(t + \delta t))^2 = \min. \qquad (14) $$

This is precisely Gauss' principle of least constraint. In the literature, one finds more often the equivalent formulation¹³,¹⁴

$$ \sum_{\alpha=1}^{N} \frac{1}{2m_{\alpha}} (m_{\alpha} \ddot{\mathbf{x}}_{\alpha} - \mathbf{f}_{\alpha})^2 = \min, \qquad (15) $$

where the accelerations are the variables to be optimized and the time argument can be omitted since the forces and accelerations are considered at the same time. The equivalence of Eqs. (14) and (15) follows immediately from the Taylor expansions (7) and (11). It must be emphasized that Eq. (14) is the original formulation of the Gauss' principle, but Eq. (15) is the one that is predominantly used in the literature. The reason is that most problems in classical mechanics are concerned with the treatment of kinematic conditions that lead straightforwardly to linear acceleration constraints of the form (1). For completeness, I give a short résumé of the use of Gauss' principle for such cases in the Appendix.

### III. ECKART AXIS CONDITIONS FROM GAUSS' PRINCIPLE

We consider now the problem of Eckart that consists in defining a reference frame for the internal motions of a flexible polyatomic molecule in a liquid or a gas.⁵ In a recent paper by Kudin and Dymarsky, it has been pointed out that the so-called Eckart axis conditions are closely related to the problem of optimal superposition of molecular structures.⁶
Here, it will be shown that both the conditions for the Eckart frame and the superposition problem are direct consequences of Gauss' principle of least constraint.

The essential point in Eckart's problem is to separate the conformational changes from arbitrary displacements of a flexible molecule, which contain also global translations and rotations. Gauss' principle offers a very simple route to perform this task. By definition, a rigid-body displacement of a molecule does not lead to an internal deformation—and thus not to a change of its internal energy. Given the atomic positions at time t, the atomic displacements within an infinitesimal time interval $\delta t$ that are due to the intramolecular dynamics are thus given by the difference between the true, unconstrained positions at time $t+\delta t$ and their virtual counterparts in presence of a rigid-body constraint,

$$ \delta\mathbf{u}_{\alpha} = \mathbf{x}_{\alpha}(t + \delta t) - \mathbf{x}_{\alpha}^{(c)}(t + \delta t). \quad (16) $$

The latter are given by

$$ \mathbf{x}_{\alpha}^{(c)}(t + \delta t) = \mathbf{x}_{\alpha}(t) + \delta\tau + \mathbf{D}\mathbf{r}_{\alpha}(t), \quad (17) $$

where $\delta\tau$ is an infinitesimal translational displacement vector, $\mathbf{r}_{\alpha}=\mathbf{x}_{\alpha}-\mathbf{X}$ is the relative position of atom $\alpha$ with respect to the center of the rotation, $\mathbf{X}$, and $\mathbf{D}$ is a rotation matrix. As for the translational motion, we consider an infinitesimal rotational motion. In this case, $\mathbf{D}\mathbf{r}_{\alpha} \approx \mathbf{r}_{\alpha} + \delta\mathbf{\phi}\mathbf{n} \wedge \mathbf{r}_{\alpha}$, where $\delta\mathbf{\phi}$ is an infinitesimal rotation angle and $\mathbf{n}$ is a unit vector defining the axis of the rotation. The rigid-body displacements can be obtained from Gauss' principle [Eq. (14)] that takes the form here

$$ \zeta = \sum_{\alpha=1}^{N} m_{\alpha} \delta\mathbf{u}_{\alpha}^{2} = \min. \quad (18) $$

The constraint function $\zeta$ depends on $\delta\tau$ and $\delta\phi = \delta\phi_{\mathbf{n}}$, where $\delta\phi_{\mathbf{x}}$, $\delta\phi_{\mathbf{y}}$, and $\delta\phi_{\mathbf{z}}$ are considered as independent variables. The necessary conditions for a minimum of $\zeta$,

$$ \frac{\partial \zeta}{\partial \delta\tau} = \sum_{\alpha=1}^{N} m_{\alpha} \delta\mathbf{u}_{\alpha} = 0, \quad (19) $$

$$ \frac{\partial \zeta}{\partial \delta\phi} = \sum_{\alpha=1}^{N} m_{\alpha} \mathbf{r}_{\alpha} \wedge \delta\mathbf{u}_{\alpha} = 0, \quad (20) $$

are precisely the Eckart conditions for the internal displacements $\delta\mathbf{u}_{\alpha}$ (Refs. 15–17) and yield six equations for the six unknowns $\delta\tau$ and $\delta\phi$.

If the center of rotation is chosen to be the center of mass, such that $M\mathbf{X}=\sum_{\alpha}m_{\alpha}\mathbf{x}_{\alpha}$, with $M=\sum_{\alpha}m_{\alpha}$ being the total mass, conditions (19) and (20) are decoupled and one finds from the translational Eckart condition (19)

$$ \mathbf{X}^{(c)}(t + \delta t) = \mathbf{X}(t + \delta t). \quad (21) $$

Performing the mass-weighted sum of Eq. (17) and using the above identity fixes the displacement $\delta\tau$,

$$ \delta\tau = \mathbf{X}(t + \delta t) - \mathbf{X}(t). \quad (22) $$

To determine the optimal infinitesimal rotation the unconstrained and the constrained atomic positions are decomposed into a component due to the center of

mass motion and a component describing the motion relative to the center of mass, $\mathbf{x}_{\alpha}(t+\delta t)=\mathbf{X}(t+\delta t)+\mathbf{r}_{\alpha}(t+\delta t)$ and $\mathbf{x}_{\alpha}^{(c)}(t+\delta t)=\mathbf{X}^{(c)}(t+\delta t)+\mathbf{r}_{\alpha}^{(c)}(t+\delta t)$. Using now Eq. (21) and that $\mathbf{r}^{(c)}(t+\delta t)=\mathbf{D}\mathbf{r}(t)$, Gauss' principle (14) amounts to solving the rotational superposition problem

$$ \zeta = \sum_{\alpha=1}^{N} m_{\alpha} (\mathbf{D}\mathbf{r}_{\alpha}(t) - \mathbf{r}_{\alpha}(t + \delta t))^{2} = \min. \quad (23) $$

From a mathematical point of view, problem (23) can be solved for any two sets of vectors, and essentially two methods are routinely used for this purpose: (a) The method by Kabsch,¹⁸ which optimizes directly the elements of the rotation matrix $\mathbf{D}$, and (b) quaternion-based algorithms, which have been proposed and discussed by different authors.¹⁹⁻²² The latter method will be briefly discussed in the next section, when optimal finite rotations are considered. Here, one looks for an optimal infinitesimal rotation, which superposes constrained and unconstrained atomic displacements in presence of physical forces. The optimal rotation vector $\delta\mathbf{\phi}$ can be easily found from the rotational Eckart condition (20), which can be written as a linear equation for $\delta\mathbf{\phi}$;

$$ \theta(t) \delta\mathbf{\phi} = \sum_{\alpha=1}^{N} m_{\alpha} \mathbf{r}_{\alpha}(t) \wedge \mathbf{r}_{\alpha}(t + \delta t) \quad (|\delta\mathbf{\phi}| \ll 1). \quad (24) $$

Here, $\theta$ is the tensor if inertia, whose elements are given by $\theta_{ij}=\sum_{\alpha=1}^{N}m_{\alpha}(r_{ij}^{2}\delta_{ij}-r_{\alpha,i}r_{\alpha,j})$, with $r_{\alpha}=|\mathbf{r}_{\alpha}|$ being the modulus of $\mathbf{r}_{\alpha}=\mathbf{r}_{\alpha}(t)$. It should be noted that the positions $\mathbf{r}_{\alpha}(t+\delta t)$ on the right-hand side of Eq. (24) can be replaced by the infinitesimal differences $\mathbf{r}_{\alpha}(t+\delta t)-\mathbf{r}_{\alpha}(t) \approx \mathbf{r}_{\alpha}(t)\delta t$. Equation (24) may thus be cast into the form

$$ \theta(t) \delta\mathbf{\phi} = \mathbf{L}(t) \delta t, \quad (25) $$

where $\mathbf{L}=\sum_{\alpha}m_{\alpha}\mathbf{r}_{\alpha} \wedge \mathbf{r}_{\alpha}$ is the angular momentum of the molecule with respect to its center of mass.

The rotation matrix $\mathbf{D}$ obtained from the minimization problem (23) relates the so-called Eckart frames of the molecule at time $t$ and $t+\delta t$, respectively. The Eckart frame of a molecule is a body-fixed, mass-centered, and orthonormal set of vectors ${\mathbf{f}_1, \mathbf{f}_2, \mathbf{f}_3}$, which is defined in such a way that the atomic equilibrium positions $\mathbf{r}_{\alpha}^{\text{eq}}$ related to that basis stay constant with time, i.e., $\mathbf{f}_i^T(t)\mathbf{r}_{\alpha}^{\text{eq}}(t)=\rho_{\alpha,i}$, where $\rho_{\alpha,i}=\text{const}$. Suppose that the Eckart frame at time $t$ is chosen to coincide with the inertial frame, such that $\theta(t)\mathbf{f}_i(t)=\theta_i\mathbf{f}_i(t)$, where $i$ are the principal moments of inertia of the molecule. The Eckart frame at time $t+\delta t$ is then defined by the set of vectors $\mathbf{f}_i(t+\delta t)=\mathbf{D}\mathbf{f}_i(t)$, which are eigenvectors of $\theta(t+\delta t)$ only if the molecule is rigid, such that $\theta(t+\delta t)=\mathbf{D}\theta(t)\mathbf{D}^T$.

## IV. CREATING TRAJECTORIES FOR INTERNAL MOLECULAR MOTIONS

In the light of computer simulations of macromolecules, in particular, of proteins, the solution of Eckart's problem is of great practical importance since it allows to analyze internal molecular motions independently from the global ones, which is often not possible in experiments. To generate a trajectory for the (simulated) internal dynamics of a macromolecule one starts from a trajectory containing the total
dynamics and calculates the internal displacements $\delta u_{\alpha}$ for each time frame by fitting the unconstrained positions at time $n\Delta t$, which define a rigid body, onto those at time $(n+1)\Delta t$. The atomic trajectories containing the internal dynamics are then obtained by adding up all internal displacements,

$$ \mathbf{x}_{\alpha}^{\text{int}}(n\Delta t) = \mathbf{x}_{\alpha}(0) + \sum_{k=1}^{n} \delta\mathbf{u}_{\alpha}(k\Delta t). \qquad (26) $$

Here, $\Delta t$ is the sampling interval, which is, in general, a multiple of the simulation time step, $\delta t$. In molecular dynamics simulations the latter is typically of the order of a femto-second and is to be considered as a “differential.” Relation (25) may be safely used to subtract the global rotation of a flexible molecule if the trajectory of the total dynamics is available for each simulation step. In many situations, the

sampling interval is, however, by a factor of 10–100 larger than $\delta t$, in order to save storage space. One is thus led to consider rigid-body displacements of the form (17) involving a finite translation vector and a rotation matrix for a finite rotation. If one requires again that the center of rotation is the center of mass, the translation vector has the same form (22) as for an infinitesimal displacement, the time differential $\delta t$ replaced by the finite difference $\Delta t$,

$$ \Delta\tau = \mathbf{X}(t + \Delta t) - \mathbf{X}(t). \qquad (27) $$

The orthogonal matrix $\mathbf{D}$ describing a finite rotation is conveniently expressed in terms of four quaternion parameters, $q=\{q_0, q_1, q_2, q_3\}$, which are normalized such that $q_0^2+q_1^2+q_2^2+q_3^2=1$. In terms of these parameters, the rotation matrix takes the form$^{23}$

$$ \mathbf{D}(q) = \begin{pmatrix} q_0^2 + q_1^2 - q_2^2 - q_3^2 & 2(-q_0q_3 + q_1q_2) & 2(q_0q_2 + q_1q_3) \\ 2(q_0q_3 + q_1q_2) & q_0^2 + q_2^2 - q_1^2 - q_3^2 & 2(-q_0q_1 + q_2q_3) \\ 2(-q_0q_2 + q_1q_3) & 2(q_0q_1 + q_2q_3) & q_0^2 + q_3^2 - q_1^2 - q_2^2 \end{pmatrix}. \qquad (28) $$

Inserting this form into expression (23) yields a constraint function $\zeta(q)$ to be minimized with respect to the quaternion parameters $q$. If one observes the normalization condition of the latter, the function $\zeta$ may be written as a quadratic form. Accounting for the normalization condition by the method of Lagrange multipliers, one requires that$^{21}$

$$ \zeta(q) = \frac{1}{2} \mathbf{q}^T \boldsymbol{\mu} \mathbf{q} - \frac{1}{2} \lambda (\mathbf{q}^T \mathbf{q} - 1) = \min. \qquad (29) $$

Here, $\mathbf{q}=(q_0, q_1, q_2, q_3)^T$ and, abbreviating $\mathbf{r}'_{\alpha} \equiv \mathbf{r}_{\alpha}(t+\Delta t)$ and $\mathbf{r}_{\alpha} \equiv \mathbf{r}_{\alpha}(t)$, the $4 \times 4$ matrix $\boldsymbol{\mu}$ can be written as

$$ \boldsymbol{\mu} = \sum_{\alpha=1}^{N} m_{\alpha} \boldsymbol{\mu}_{\alpha}, \qquad (30) $$

where $\boldsymbol{\mu}_{\alpha}$ has the form

$$ \boldsymbol{\mu}_{\alpha} = \begin{pmatrix} (\mathbf{r}_{\alpha} - \mathbf{r}'_{\alpha})^2 & 2(\mathbf{r}'_{\alpha} \wedge \mathbf{r}_{\alpha})^T \\ 2(\mathbf{r}'_{\alpha} \wedge \mathbf{r}_{\alpha}) & (\mathbf{r}_{\alpha} + \mathbf{r}'_{\alpha})^2 \mathbf{1} - 2(\mathbf{r}'_{\alpha} \mathbf{r}_{\alpha}^T + \mathbf{r}_{\alpha} \mathbf{r}'_{\alpha}^T) \end{pmatrix}. \qquad (31) $$

The minimization of Eq. (29) leads then to the eigenvalue problem

$$ \boldsymbol{\mu} \mathbf{q} = \lambda \mathbf{q}. \qquad (32) $$

Since the target function (23) fulfills $\zeta(q) \ge 0$, the matrix $\boldsymbol{\mu}$ is positive semidefinite, leading to four eigenvectors $\mathbf{q}_j$ and four associated eigenvalues $\lambda_j$ ($j=1, \dots, 4$), with $\lambda_j \ge 0$. The latter may be ordered by size, such that $\lambda_1$ is the smallest one. For any of the normalized eigenvectors, we have $\mathbf{q}_j^T \boldsymbol{\mu} \mathbf{q}_j = \zeta(q_j) = \lambda_j$, which shows that the eigenvalues are the fit errors. The normalized eigenvector corresponding to the smallest eigenvalue, $\lambda_1$, is thus the quaternion describing the optimal fit. On account of relation (18), we have

$$ \lambda_1 = \sum_{\alpha=1}^{N} m_{\alpha} \delta\mathbf{u}_{\alpha}^2. \qquad (33) $$

The quaternion parameters resulting from the solution of the eigenvalue problem (32) can be related to more familiar parameters describing a rotation. One has

$$ \mathbf{q} = \begin{pmatrix} \cos(\phi/2) \\ \sin(\phi/2) \mathbf{n} \end{pmatrix}, \qquad (34) $$

where $\phi$ is the angle of the rotation and $\mathbf{n}$ is the unit vector in the direction of the rotation axis.

It is not obvious that the infinitesimal rotation described by the rotation vector $\delta\boldsymbol{\phi}$ introduced in the previous section is obtained from the above eigenvalue problem in the limit of an infinitesimal displacement $\mathbf{r}_{\alpha}(t) \rightarrow \mathbf{r}_{\alpha}(t+\Delta t) = \mathbf{r}_{\alpha}(t) + \delta\mathbf{r}_{\alpha}$. To show this, one can use a perturbative approach to solve eigenvalue problem (32). For this purpose, we split $\boldsymbol{\mu} = \boldsymbol{\mu}^{(0)} + \boldsymbol{\mu}^{(1)}$, $\mathbf{q} = \mathbf{q}^{(0)} + \mathbf{q}^{(1)}$ and $\lambda = \lambda^{(0)} + \lambda^{(1)}$, where $\lambda^{(0)}=0$. Omitting combinations of terms of first order in the perturbation one obtains

$$ \boldsymbol{\mu}^{(0)} \mathbf{q}^{(0)} + \boldsymbol{\mu}^{(1)} \mathbf{q}^{(1)} + \boldsymbol{\mu}^{(1)} \mathbf{q}^{(0)} = \lambda^{(1)} \mathbf{q}^{(0)}. \qquad (35) $$

Here, $\boldsymbol{\mu}^{(0)}$ and $\boldsymbol{\mu}^{(1)}$ contain terms of order zero and one in $\delta\mathbf{r}_{\alpha}$, respectively. On account of relation (30), one has $\boldsymbol{\mu} = \sum_{\alpha=1}^{N} m_{\alpha} (\boldsymbol{\mu}_{\alpha}^{(0)} + \boldsymbol{\mu}_{\alpha}^{(1)})$, where

$$ \boldsymbol{\mu}_{\alpha}^{(0)} = \begin{pmatrix} 0 & 0^T \\ 0^T & 4\mathbf{r}_{\alpha}^2 \mathbf{1} - 4\mathbf{r}_{\alpha} \mathbf{r}_{\alpha}^T \end{pmatrix}, $$
$$ \boldsymbol{\mu}_{\alpha}^{(1)} = \begin{pmatrix} 0 & 2(\boldsymbol{\delta}\boldsymbol{r}_{\alpha} \wedge \boldsymbol{r}_{\alpha})^{\mathrm{T}} \\ 2(\boldsymbol{\delta}\boldsymbol{r}_{\alpha} \wedge \boldsymbol{r}_{\alpha}) & 4\boldsymbol{\delta}\boldsymbol{r}_{\alpha}^{\mathrm{T}}\boldsymbol{r}_{\alpha} - 1 - 2(\boldsymbol{\delta}\boldsymbol{r}_{\alpha}\boldsymbol{r}_{\alpha}^{\mathrm{T}} + \boldsymbol{r}_{\alpha}\boldsymbol{\delta}\boldsymbol{r}_{\alpha}^{\mathrm{T}}) \end{pmatrix}. $$

An infinitesimal displacement of the atoms entails an infinitesimal global rotation, which is characterized by an infinitesimal rotation angle $\delta\phi$. Developing the form (34) for the quaternion parameters up to linear terms in $\phi$ shows that $\mathbf{q} = \mathbf{q}^{(0)} + \mathbf{q}^{(1)}$, where

$$ \mathbf{q}^{(0)} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \quad \text{and} \quad \mathbf{q}^{(1)} = \begin{pmatrix} 0 \\ \delta\phi/2 \end{pmatrix}, $$

using again the definition $\delta\phi = \delta\phi\mathbf{n}$. One observes that the normalization condition $\mathbf{q}^{\mathrm{T}}\mathbf{q} \approx \mathbf{q}^{(0)\mathrm{T}}\mathbf{q}^{(0)} + 2\mathbf{q}^{(0)\mathrm{T}}\mathbf{q}^{(1)} = 1$ is verified up to first order. Inserting the explicit forms for $\boldsymbol{\mu}^{(j)}$ and $\mathbf{q}^{(j)}$ ($j=0,1$) into Eq. (35) leads to

$$ \begin{pmatrix} 0 & \mathbf{0}^{\mathrm{T}} \\ \mathbf{0} & \boldsymbol{\theta} \end{pmatrix} \begin{pmatrix} 0 \\ \boldsymbol{\delta}\phi \end{pmatrix} + \left( \sum_{\alpha=1}^{N} m_{\alpha} (\boldsymbol{\delta}\boldsymbol{r}_{\alpha} \wedge \boldsymbol{r}_{\alpha}) \right) = \begin{pmatrix} \mathbf{0} \\ \boldsymbol{\lambda}^{(1)} \end{pmatrix}, $$

where $\boldsymbol{\theta}$ is the tensor of inertia introduced in the previous section. It follows that $\boldsymbol{\lambda}^{(1)}=0$ and that

$$ \boldsymbol{\theta}\boldsymbol{\delta}\boldsymbol{\phi} = \sum_{\alpha=1}^{N} m_{\alpha}(\boldsymbol{r}_{\alpha} \wedge \boldsymbol{\delta}\boldsymbol{r}_{\alpha}), $$

confirming thus Eq. (25).

## V. CONCLUSION

It has been shown that Gauss' principle of least constraint allows one to derive the so-called Eckart axis conditions for a vibrating polyatomic molecule in a particularly simple way. The basic problem here is to separate the internal motions of the molecule from the global ones, which include rigid-body translations and rotations. Since rigid-body constraints are a special class of kinematic conditions, Gauss' principle can be applied to determine their contribution to the total motion of a molecule. For this task, the original formulation by Gauss is particularly useful, in which the evolution of the particle positions in a constrained system is formulated as a least-squares problem. A detailed derivation based on concepts of modern linear algebra was given in a separate section. In the context of this article, Gauss' principle leads to a rotational superposition problem of constrained and unconstrained atomic positions with respect to the molecular center of mass. The resulting optimal rotation is a priori an infinitesimal rotation. It was shown that Eckart's problem can also be solved in cases where the consecutive atomic positions are not separated by differentials of positions, but by finite differences. This point is of importance when trajectories for the internal dynamics of flexible molecules are to be constructed from coarse-grained molecular dynamics trajectories. In the case of finite displacements, the rotational superposition problem can be solved by a quaternion-based method, which leads for each molecule to the solution of an eigenvalue problem for a positive 4×4 matrix, and the solution for infinitesimal displacements is retrieved by a perturbative solution of the eigenvalue problem.

## APPENDIX: GAUSS' PRINCIPLE FOR ACCELERATIONS—TWO EXAMPLES

To treat mechanical problems subjected to constraints leading to linear acceleration constraints of the form (1) it is useful to write Gauss' principle (15) in matrix form

$$ \frac{1}{2}(\mathbf{M}^{1/2}\ddot{\mathbf{r}} - \mathbf{M}^{-1/2}\mathbf{f})^2 = \min, \quad (A1) $$

which is possible since $\mathbf{M}$ is positive definite. The above function is to be minimized with respect to the acceleration vector $\ddot{\mathbf{r}}$, which is subject to constraints of the form (1). This can be accomplished by the method of Lagrange multipliers. Introducing the function

$$ \zeta(\ddot{\mathbf{r}}, \boldsymbol{\lambda}) = \frac{1}{2}(\mathbf{M}^{1/2}\ddot{\mathbf{r}} - \mathbf{M}^{-1/2}\mathbf{f})^2 + (\mathbf{A}\ddot{\mathbf{r}} - \mathbf{b})^T\boldsymbol{\lambda}, \quad (A2) $$

one looks for the minimum with respect to $\ddot{\mathbf{r}}$ and $\boldsymbol{\lambda}$. The latter is a column vector containing $s$ parameters, $\lambda_1, \dots, \lambda_s$, where $s$ is the number of constraints. The minimization of $\zeta$ leads to $3N+s$ equations for the components of $\ddot{\mathbf{r}}$ and $\boldsymbol{\lambda}$,

$$ \frac{\partial \zeta}{\partial \ddot{\mathbf{r}}} = \mathbf{0} \Rightarrow \mathbf{M}\ddot{\mathbf{r}} = \mathbf{f} + \mathbf{A}^T\boldsymbol{\lambda}, \quad (A3) $$

$$ \frac{\partial \zeta}{\partial \boldsymbol{\lambda}} = \mathbf{0} \Rightarrow \mathbf{A}\ddot{\mathbf{r}} = \mathbf{b}. \quad (A4) $$

Equation (A3) shows that the constraint forces are given by

$$ \mathbf{z} = \mathbf{A}^T\boldsymbol{\lambda}. \quad (A5) $$

The vector of the constraint forces is a linear combination of the rows of $\mathbf{A}$, which span the space $V_{\perp}$. The parameters $\lambda_i$ are obtained by solving Eq. (A3) for the accelerations and inserting the result in Eq. (A4),

$$ \mathbf{A}\mathbf{M}^{-1}\mathbf{A}^T\boldsymbol{\lambda} = \mathbf{b} - \mathbf{A}\mathbf{M}^{-1}\mathbf{f}. \quad (A6) $$

The above set of equations can be solved if the rows in $\mathbf{A}$ are linearly independent.

Two simple examples may illustrate the use of the above equations. Consider first the motion of a single particle on a sphere. The constraint is here

$$ \mathbf{x}^2 = R^2, $$

where $R$ is the radius of the sphere. Differentiating the above relation twice leads to $\mathbf{x}^T\ddot{\mathbf{x}} = -\dot{\mathbf{x}}^2$. Here, we have $\mathbf{A} = \mathbf{x}^T$ and $\mathbf{b}$ becomes a scalar, $\mathbf{b} = -\dot{\mathbf{x}}^2$. With $\mathbf{M} = m\mathbf{1}$ Eq. (A6) takes the form $\mathbf{x}^2\boldsymbol{\lambda} = -m\dot{\mathbf{x}}^2 - \mathbf{x}^T\mathbf{f}$ and using the constraint, the equation of motion becomes

$$ m\ddot{\mathbf{x}} = \left( \mathbf{1} - \frac{\mathbf{x}\mathbf{x}^T}{R^2} \right) \mathbf{f} - m\frac{\dot{\mathbf{x}}^2}{R^2}. $$

One realizes that only the tangential component of the external force acts and that the second term on the right-hand side is the centripetal force.

The second example treats a nonholonomic constraint. We consider a system at constant kinetic energy,
$$ \frac{1}{2} \sum_{\alpha=1}^{N} m \dot{\mathbf{x}}_{\alpha}^{2} = \frac{3}{2} N k_{B} T, $$

where all particles are supposed to have the same mass. As usual, $T$ denotes the absolute temperature and $k_B$ is the Boltzmann constant. The matrix $A$ here takes the form $A=m\dot{\mathbf{x}}^T$ and $b$ is again a scalar that is simply zero, $b=0$. The equation for the single Lagrange multiplier $\lambda$ reads thus $m\dot{\mathbf{x}}^2\lambda=-\dot{\mathbf{x}}^T\mathbf{f}$. Using the constraint, we have $\lambda=-\dot{\mathbf{x}}^T\mathbf{f}/(3Nk_B T)$, and the equations of motion become

$$ m \dot{\mathbf{x}}_{\alpha} = \mathbf{f}_{\alpha} - m \left\{ \sum_{\beta=1}^{N} \frac{\dot{\mathbf{x}}_{\beta}^{T} \mathbf{f}_{\beta}}{3Nk_{B}T} \right\} \dot{\mathbf{x}}_{\alpha}. $$

The term $m\{...\}$ appears here as a “friction constant,” which can, however, take positive and negative values.

<footnote n="1">¹C. Gauss, J. Reine Angew. Math. 4, 232 (1829).</footnote>

<footnote n="2">²D. Evans, W. Hoover, B. Failor, B. Moran, and A. Ladd, Phys. Rev. A 28, 1016 (1983).</footnote>

<footnote n="3">³W. Hoover, *Time Reversibility, Computer Simulation, and Chaos*, Advanced Series in Nonlinear Dynamics (World Scientific, Singapore, 1999), Vol. 13.</footnote>

<footnote n="4">⁴W. Hoover, *Molecular Dynamics*, Lecture Notes in Physics (Springer, Berlin, 1986).</footnote>

<footnote n="5">⁵C. Eckart, Phys. Rev. 47, 552 (1935).</footnote>

<footnote n="6">⁶K. Kudin and A. Dymarsky, J. Chem. Phys. 122, 224105 (2005).</footnote>

<footnote n="7">⁷R. Penrose, Proc. Cambridge Philos. Soc. 51, 406 (1955).</footnote>

<footnote n="8">⁸A. Ben-Israel and T. Greville, *Generalized Inverses: Theory and Applications* (Wiley, New York, 1974).</footnote>

<footnote n="9">⁹S. Campbell and C. Meyer, Jr., *Generalized Inverses of Linear Transformations* (Pitman, New York, 1979).</footnote>

<footnote n="10">¹⁰R. Bott and R. Duffin, Trans. Am. Math. Soc. 74, 99 (1953).</footnote>

<footnote n="11">¹¹G. Kneller, J. Chem. Phys. 125, 114107 (2006).</footnote>

<footnote n="12">¹²G. Kneller, J. Chem. Phys. 127, 164114 (2007).</footnote>

<footnote n="13">¹³L. Pars, *A Treatise on Analytical Dynamics* (Heinemann, London, 1965).</footnote>

<footnote n="14">¹⁴C. Lanczos, *The Variational Principles of Mechanics*, 4th ed. (University of Toronto Press, Toronto, 1974).</footnote>

<footnote n="15">¹⁵L. Biedenham and J. Louck, in *Encyclopedia of Mathematics and its Applications*, edited by G.-C. Rota (Addison-Wesley, Reading, 1981), Vol. 8.</footnote>

<footnote n="16">¹⁶D. Papoušek and M. Aliev, *Molecular Vibrational-Rotational Spectra* (Elsevier, New York, 1982).</footnote>

<footnote n="17">¹⁷S. Califano, *Vibrational States* (Wiley, New York, 1976).</footnote>

<footnote n="18">¹⁸W. Kabsch, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr. 32, 922 (1976).</footnote>

<footnote n="19">¹⁹J. L. Tietze, J. Astronaut. Sci. 30, 171 (1982).</footnote>

<footnote n="20">²⁰R. Diamond, Acta Crystallogr., Sect. A: Found. Crystallogr. 44, 211 (1988).</footnote>

<footnote n="21">²¹G. Kneller, Mol. Simul. 7, 113 (1991).</footnote>

<footnote n="22">²²G. Kneller, J. Comput. Chem. 26, 1660 (2005).</footnote>

<footnote n="23">²³S. Altmann, *Rotations, Quaternions, and Double Groups* (Clarendon, Oxford, 1986).</footnote>