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

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.


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. 1 Despite this fact, it has not found its way into many standard textbooks on classical mechanics.The possibility to consider also nonholonomic 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. 2 and to the monographs by Hoover. 3,4n 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. 5The 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. 6In 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 and Angewandte Mathematik, 1 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 , x ˙, t͒ =0 ͑j =1, ... ,s͒, where x comprises all positions x ␣ and x ˙the corresponding velocities x ˙␣.A set of s linear constraints for the accelerations of the form is obtained by differentiating the constraints h j ͑x , t͒ = 0 twice and the constraints g j ͑x , x ˙, t͒ = 0 once with respect to time.The matrix A in Eq. ͑1͒ has the dimensions s ϫ 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, [7][8][9] which is usually denoted by A + , Here, x ¨ʈ is an undetermined vector in the null space of A, fulfilling Ax ¨ʈ = 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, AA + = 1 s is the s ϫ s unit matrix and A + A = A T ͑AA T ͒ −1 A ϵ P Ќ is a projector on an s-dimensional subspace of R 3N .The latter will be denoted by the symbol V Ќ .The orthogonal complement, V ʈ , is the null space of A and the corresponding projector is P ʈ = 1 − P Ќ .With these preliminaries it follows from relation ͑2͒ by multiplication with P Ќ from the left that Relation ͑2͒ thus represents the decomposition of the 3N-dimensional acceleration vector x ¨into a known component in V Ќ , which is entirely determined by the positions and the velocities of the particles, and an unknown component in V ʈ .The latter must be determined from Newton's equations of motion.Writing x ¨= x ¨ʈ + x ¨Ќ, the latter read Here, f and z contain the external forces and the constraint forces on the particles, respectively, and M is the diagonal 3N ϫ 3N mass matrix

͑5͒
It must be emphasized that Newton's equations of motion constitute a system of 3N linear equations for x ¨ʈ and z.A unique solution for both unknowns exists if one requires that z Ќ x ¨ʈ.͑6͒ The solution of linear systems of equations for two mutually orthogonal unknown vectors has been described by Bott and Duffin 10 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.

͑9͒
Note that the subspace V ʈ ͑t͒ changes with time implicitly through the change of the positions and the velocities.It follows thus from Eqs. ͑6͒, ͑8͒, and ͑9͒ that

͑10͒
where the superscript T denotes a transposition.Condition ͑10͒ is nothing but D'Alembert's principle, which is usually written in the form ͚ ␣ z ␣ T ␦x ␣ = 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 + ␦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 ␦t 2 the latter are given by and one finds that z͑t͒ = M͑x͑t + ␦t͒ − x ͑0͒ ͑t + ␦t͒͒.

͑12͒
Since the unconstrained positions at time t + ␦t are completely determined by the positions, the velocities, and the external forces at time t, there are no degrees of freedom for variations, ␦x ͑0͒ ͑t + ␦t͒ = 0, and one may write

͑14͒
This is precisely Gauss' principle of least constraint.In the literature, one finds more often the equivalent formulation 13,14 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. 5In 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. 6ere, 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 ␦t that are due to the intramolecular dynamics are thus given by the difference between the true, unconstrained positions at time t + ␦t and their virtual counterparts in presence of a rigid-body constraint,

͑16͒
The latter are given by where ␦ is an infinitesimal translational displacement vector, r ␣ = x ␣ − X is the relative position of atom ␣ with respect to the center of the rotation, X, and D is a rotation matrix.As for the translational motion, we consider an infinitesimal rotational motion.In this case, Dr ␣ Ϸ r ␣ + ␦n ∧ r ␣ , where ␦ is an infinitesimal rotation angle and 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 The constraint function depends on ␦ and ␦ ϵ ␦n, where ␦ x , ␦ y , and ␦ z are considered as independent variables.The necessary conditions for a minimum of ,

͑20͒
are precisely the Eckart conditions for the internal displacements ␦u ␣ ͑Refs.15-17͒ and yield six equations for the six unknowns ␦ and ␦.
If the center of rotation is chosen to be the center of mass, such that MX = ͚ ␣ m ␣ x ␣ , with M = ͚ ␣ m ␣ being the total mass, conditions ͑19͒ and ͑20͒ are decoupled and one finds from the translational Eckart condition ͑19͒ X ͑c͒ ͑t + ␦t͒ = X͑t + ␦t͒.

͑21͒
Performing the mass-weighted sum of Eq. ͑17͒ and using the above identity fixes the displacement ␦,

͑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, x ␣ ͑t + ␦t͒ = X͑t + ␦t͒ + r ␣ ͑t + ␦t͒ and x ␣ ͑c͒ ͑t + ␦t͒ = X ͑c͒ ͑t + ␦t͒ + r ␣ ͑c͒ ͑t + ␦t͒.Using now Eq.͑21͒ and that r ͑c͒ ͑t + ␦t͒ = Dr͑t͒, Gauss' principle ͑14͒ amounts to solving the rotational superposition problem

͑23͒
][21][22] 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 ␦ can be easily found from the rotational Eckart condition ͑20͒, which can be written as a linear equation for ␦;

͑24͒
Here, is the tensor if inertia, whose elements are given by It should be noted that the positions r ␣ ͑t + ␦t͒ on the right-hand side of Eq. ͑24͒ can be replaced by the infinitesimal differences r ␣ ͑t + ␦t͒ − r ␣ ͑t͒Ϸr ˙␣͑t͒␦t.Equation ͑24͒ may thus be cast into the form ͑t͒␦ = L͑t͒␦t, ͑25͒ where L = ͚ ␣ m ␣ r ␣ ∧ r ˙␣ is the angular momentum of the molecule with respect to its center of mass.
The rotation matrix D obtained from the minimization problem ͑23͒ relates the so-called Eckart frames of the molecule at time t and t + ␦t, respectively.The Eckart frame of a molecule is a body-fixed, mass-centered, and orthonormal set of vectors ͕f 1 , f 2 , f 3 ͖, which is defined in such a way that the atomic equilibrium positions r ␣ eq related to that basis stay constant with time, i.e., f i T ͑t͒r ␣ eq ͑t͒ = ␣,i , where ␣,i = const.Suppose that the Eckart frame at time t is chosen to coincide with the inertial frame, such that ͑t͒f i ͑t͒ = i f i ͑t͒, where i are the principal moments of inertia of the molecule.The Eckart frame at time t + ␦t is then defined by the set of vectors f i ͑t + ␦t͒ = Df i ͑t͒, which are eigenvectors of ͑t + ␦t͒ only if the molecule is rigid, such that ͑t + ␦t͒ = D͑t͒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 ␦u ␣ for each time frame by fitting the unconstrained positions at time n⌬t, which define a rigid body, onto those at time ͑n +1͒⌬t.The atomic trajectories containing the internal dynamics are then obtained by adding up all internal displacements,

͑26͒
Here, ⌬t is the sampling interval, which is, in general, a multiple of the simulation time step, ␦t.In molecular dynamics simulations the latter is typically of the order of a femtosecond 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 ␦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 ␦t replaced by the finite difference ⌬t, ⌬ = X͑t + ⌬t͒ − X͑t͒.͑27͒ The orthogonal matrix 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 D͑q͒ = q 0 2 + q 1 2 − q 2 2 − q 3 2 2͑− q 0 q 3 + q 1 q 2 ͒ 2͑q 0 q 2 + q 1 q 3 ͒ 2͑q 0 q 3 + q 1 q 2 ͒ q 0 2 + q 2 2 − q 1 2 − q 3 2 2͑− q 0 q 1 + q 2 q 3 ͒ 2͑− q 0 q 2 + q 1 q 3 ͒ 2͑q 0 q 1 + q 2 q 3 ͒ q 0 2 + q 3 2 − q 1 2 − q 2 2 .

͑28͒
Inserting this form into expression ͑23͒ yields a constraint function ͑q͒ to be minimized with respect to the quaternion parameters q.If one observes the normalization condition of the latter, the function may be written as a quadratic form.
Accounting for the normalization condition by the method of Lagrange multipliers, one requires that 21 ͑q͒ = 1 2 q T q − 1 2 ͑q T q − 1͒ = min.͑29͒ Here, q = ͑q 0 , q 1 , q 2 , q 3 ͒ T and, abbreviating r ␣ Ј ϵ r ␣ ͑t + ⌬t͒ and r ␣ ϵ r ␣ ͑t͒, the 4 ϫ 4 matrix can be written as where ␣ has the form The minimization of Eq. ͑29͒ leads then to the eigenvalue problem q = q.͑32͒ Since the target function ͑23͒ fulfills ͑q͒ ജ 0, the matrix is positive semidefinite, leading to four eigenvectors q j and four associated eigenvalues j ͑j =1, ... ,4͒, with j ജ 0. The latter may be ordered by size, such that 1 is the smallest one.For any of the normalized eigenvectors, we have q j T q j = ͑q j ͒ = j , which shows that the eigenvalues are the fit errors.The normalized eigenvector corresponding to the smallest eigenvalue, 1 , is thus the quaternion describing the optimal fit.On account of relation ͑18͒, we have ͑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 where is the angle of the rotation and 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 ␦ introduced in the previous section is obtained from the above eigenvalue problem in the limit of an infinitesimal displacement r ␣ ͑t͒ → r ␣ ͑t + ␦t͒ = r ␣ ͑t͒ + ␦r ␣ .
To show this, one can use a perturbative approach to solve eigenvalue problem ͑32͒.For this purpose, we split = ͑0͒ + ͑1͒ , q = q ͑0͒ + q ͑1͒ and = ͑0͒ + ͑1͒ , where ͑0͒ = 0. Omitting combinations of terms of first order in the perturbation one obtains ͑0͒ q ͑0͒ + ͑0͒ q ͑1͒ + ͑1͒ q ͑0͒ = ͑1͒ q ͑0͒ .͑35͒ Here, ͑0͒ and ͑1͒ contain terms of order zero and one in ␦r ␣ , respectively.On account of relation ͑30͒, one has An infinitesimal displacement of the atoms entails an infinitesimal global rotation, which is characterized by an infinitesimal rotation angle ␦.Developing the form ͑34͒ for the quaternion parameters up to linear terms in shows that q = q ͑0͒ + q ͑1͒ , where q ͑0͒ = ͩ 1 0 ͪ and q ͑1͒ = ͩ 0 ␦/2 ͪ, using again the definition ␦ = ␦n.One observes that the normalization condition q T q Ϸ q ͑0͒T q ͑0͒ +2q ͑0͒T q ͑1͒ = 1 is verified up to first order.Inserting the explicit forms for ͑j͒ and q ͑j͒ ͑j =0,1͒ into Eq.͑35͒ leads to where is the tensor of inertia introduced in the previous section.It follows that ͑1͒ = 0 and that 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 rigidbody 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 1 2 ͑M 1/2 r ¨− M −1/2 f͒ 2 = min, ͑A1͒ which is possible since M is positive definite.The above function is to be minimized with respect to the acceleration vector r ¨, which is subject to constraints of the form ͑1͒.This can be accomplished by the method of Lagrange multipliers.
Introducing the function one looks for the minimum with respect to r ¨and .The latter is a column vector containing s parameters, 1 , ... , s , where s is the number of constraints.The minimization of leads to 3N + s equations for the components of r ¨and , Equation ͑A3͒ shows that the constraint forces are given by z = A T .͑A5͒ The vector of the constraint forces is a linear combination of the rows of A, which span the space V Ќ .The parameters i are obtained by solving Eq. ͑A3͒ for the accelerations and inserting the result in Eq. ͑A4͒, The above set of equations can be solved if the rows in 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 where R is the radius of the sphere.Differentiating the above relation twice leads to x T x ¨=−x ˙2.Here, we have A = x T and b becomes a scalar, b =−x ˙2.With M = m1 Eq. ͑A6͒ takes the form x 2 =−mx ˙2 − x T f and using the constraint, the equation of motion becomes 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, The term m͕...͖ appears here as a "friction constant," which can, however, take positive and negative values.

‫ץ‬ ‫ץ‬r ¨= 0 ⇒
Mr ¨= f + A T , ͑A3͒ ‫ץ‬ ‫ץ‬ = 0 ⇒ Ar ¨= b. ͑A4͒ 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 = mx ˙T and b is again a scalar that is simply zero, b = 0.The equation for the single Lagrange multiplier reads thus mx ˙2 =−x ˙Tf.Using the constraint, we have =−x ˙Tf / ͑3Nk B T͒, and the equations of motion become mx ˙␣ = f ␣ − m ͚ͭ ␤=1 N x ˙␤ T f ␤ 3Nk B T ͮ x ˙␣.