Appendix 1
Dynamical systems
1.1 Hamilton’s equations of motion
Newton’s equations of motion are limited to describing the motion of independent point masses in Cartesian coordinates. A more general dynamical postulate is Hamilton’s principle. It can be used to derive equations of motions for arbitrary dynamical variables (e.g. the angle of a pendulum, or in the context of molecular dynamics the internal degrees of freedom of a heat bath; this is useful when discussing constant-temperature molecular dynamics). Hamilton’s principle is the postulate that the action S, an integral quantity of the whole dynamics, has to be stationary. Newton’s equations of motion are contained in this postulate, as will be shown below.
The action S for a system moving between time \(t_{1}\) and \(t_{2}\) is given by \begin {equation} S = \int _{t_{1}}^{t_{2}}{dt\ L\left ( \{ q_i \},\{ {\dot {q}}_i \},t \right )} \end {equation} where \(L=E_\text {kin}-E_\text {pot}\) is called the Lagrangian.1 Hamilton’s principle now states that the system will move along (generalized) coordinates \(\{ q_i( t) \}\) such that the action S is stationary. This means the action does not vary for infinitesimal deviations from this path, it is a maximum or minimum.
To determine this stationary path let us assume we have some path \(\{ q_i(t) \}\) and infinitesimal variations \(\{ \varepsilon _i(t) \}\) around that path. The change of the action functional would then be \begin {equation} \begin {split} \delta S &= \int _{t_{1}}^{t_{2}}{dt\ \left \lbrack L\left ( \{ q_i + \varepsilon _i \},\{ {\dot {q}}_i + {\dot {\varepsilon }}_i \},t \right ) - L\left ( \{ q_i \},\{ {\dot {q}}_i \},t \right ) \right \rbrack } \\ &= \int _{t_{1}}^{t_{2}}{dt\ \sum _i\left \lbrack \varepsilon _i\frac {\partial L}{\partial q_i} + {\dot {\varepsilon }}_i\frac {\partial L}{\partial {\dot {q}}_i} \right \rbrack }. \end {split} \end {equation} Integration by parts yields \begin {equation} \delta S = \sum _i\left \lbrack \varepsilon _i\frac {\partial L}{\partial q_i} \right \rbrack _{t_{1}}^{t_{2}} + \int _{t_{1}}^{t_{2}}{\text {dt\ }\sum _i{\varepsilon _i\left \lbrack \frac {\partial L}{\partial q_i} - \frac {d}{dt}\frac {\partial L}{\partial {\dot {q}}_i} \right \rbrack }}. \label {eq:dS} \end {equation} The first part of Eq. \eqref{eq:dS} vanishes because \(\varepsilon _i\left ( t_{1} \right ) = 0\) and \(\varepsilon _i\left ( t_{2} \right ) = 0\) (by definition; those are the start and endpoints of our path). Stationarity requires \(\delta S = 0\) for any \(\varepsilon _i\). This leads to \begin {equation} \frac {d}{dt}\frac {\partial L}{\partial {\dot {q}}_i} - \frac {\partial L}{\partial q_i} = 0, \label {eq:Lagrange} \end {equation} the Lagrange equation of motion.
Let us reiterate that stationarity of the action is a postulate. By the above derivation, we showed that it is equivalent to the Lagrange equation of motion. The equation of motion can be physically observed; the action itself cannot be observed. An alternative derivation of the Lagrange equation of motion from Newton’s equations of motion employs d’Alemberts principle. This derivation is given in Section 1.2 .
Note: We will now need to compute the Legendre transform of the Lagrangian. The Legendre transform of a function \(F(x)\) is defined as \begin {equation} F^{*}(y) = \min _{x}\left \lbrack xy - F(x) \right \rbrack . \end {equation} In what follows, we will assume that \(F(x)\) is differentiable. Then \(y(x) = \frac {dF}{dx}\) and with this definition for \(y(x)\) and its inverse \(x(y)\) we have \begin {equation} F^{*}(y) = x(y)y - F\left ( x(y) \right ). \label {eq:Legendre} \end {equation} It is instructive to write the total differential of Eq. \eqref{eq:Legendre}, \begin {equation} dF^{*} = xdy + ydx - dF = xdy + \frac {dF}{dx}dx - dF = x(y)dy. \end {equation} Noting that the total differential of the original function \(F\) is \begin {equation} dF = \frac {dF}{dx}dx = y(x)dx \end {equation} we immediately see that we have changed the independent variable from \(x\) to \(y\) when going from \(F\) to \(F^{*}\). The same statement can be expressed as \begin {align} F(x) &= \int \dif F = \int y(x)\dif x \\ F^{*}(y) &= \int \dif F^{*} = \int x(y)\dif y \end {align}
A common way of expressing the relationships above is \begin {equation} \frac {\dif F^{*}}{\dif y} = \left ( \frac {\dif F}{\dif x} \right )^{- 1}(y). \end {equation} The derivatives of the Legendre transformation pairs are their respective inverses.
Starting from the Lagrangian L, we define an additional state function, the Hamiltonian \begin {equation} H( \{ q_i \},\{ p_i \},t ) = \min _{\{ \dot {q}_i \}}\left [ \sum _i \dot {q}_i p_i - L( \{ q_i \},\{ \dot {q}_i \},t) \right ], \label {eq:LegrendeH} \end {equation} from the Legendre transformation of \({\dot {q}}_i\) to \(p_i\). The minimization in Eq. \eqref{eq:LegrendeH} implies that the momenta \(p_i\) are given by \begin {equation} p_i = \frac {\partial L}{\partial {\dot {q}}_i}. \label {eq:momenta} \end {equation} Equation \eqref{eq:momenta} constitutes the definition of what are called generalized momenta. We insert this expression into Lagrange’s equation, Eq. \eqref{eq:Lagrange}, to get \begin {equation} {\dot {p}}_i = \frac {\partial L}{\partial q_i}. \end {equation} Furthermore, Eq. \eqref{eq:LegrendeH} becomes \begin {equation} H\left ( \{ q_i \},\{ p_i \},t \right ) = \sum _i\dot {q}_i p_i - L( \{ q_i \},\{ \dot {q}_i \},t ), \end {equation} and using \(L\left ( \{ q_i \},\{ {\dot {q}}_i \},t \right ) = E_\text {kin}\left ( \{ {\dot {q}}_i \} \right ) - E_\text {pot}(\{ q_i \},t)\) we find \(p_i = \frac {\partial E_\text {kin}}{\partial {\dot {q}}_i}\). Since \(E_\text {kin}\left ( \{ p_i \} \right )\) is the Legendre transformation of \(E_\text {kin}\left ( \{ {\dot {q}}_i \} \right )\), \begin {equation} H\left ( \{ q_i \},\{ p_i \},t \right ) = E_\text {kin}\left ( \{ p_i \} \right ) + E_\text {pot}\left ( \{ q_i \},t \right ). \end {equation} Using Eq. (2.33), the total differential of H is given by \begin {equation} \dif H = \sum _i \dot {q}_i \dif p_i + \sum _ip_id{\dot {q}}_i - \sum _i \frac {\partial L}{\partial q_i} \dif q_i - \sum _i \frac {\partial L}{\partial \dot {q}_i} \dif \dot {q}_i - \frac {\partial L}{\partial t} \dif t \end {equation}
With the definition of the generalized momenta, Eq. (2.31), and their time derivatives, Eq. (2.32), this becomes \begin {equation} \dif H = \sum _i \dot {q}_i \dif p_i - \sum _i\dot {p}_i\dif q_i - \frac {\partial L}{\partial t}\dif t. \end {equation} We can furthermore straightforwardly express the total differential of \(H\left ( \{ q_i \},\{ p_i \},t \right )\) by \begin {equation} dH = \sum _i\frac {\partial H}{\partial p_i}\dif p_i + \sum _i \frac {\partial H}{\partial q_i}\dif q_i + \frac {\partial H}{\partial t}dt. \end {equation} Comparison of Eqs. (2.36) and (2.37) yields \begin {align} {\dot {q}}_i &= \frac {\partial H}{\partial p_i} \label {eq:Hamilton1} \\ {\dot {p}}_i &= - \frac {\partial H}{\partial q_i} \label {eq:Hamilton2} \\ - \frac {\partial L}{\partial t} &= \frac {\partial H}{\partial t} \label {eq:Hamilton3} \end {align}
Equations \eqref{eq:Hamilton1} to \eqref{eq:Hamilton3} are called Hamilton’s equations of motion.
For our molecular systems (without thermostats), the Hamiltonian \(H\) is just the total energy \(E_\text {tot}\) of the system and not explicitly time-dependent. The last Eq. \eqref{eq:Hamilton3} then vanishes and \(H\) becomes a conserved quantity. Hamilton’s equations of motion are two coupled differential equations of first order. They describe the motion of the system in the phase space spanned by the generalized coordinates, \(\v {\Gamma } = \{ q_i \},\{ p_i \}\). In contrast, Lagrange’s equation of motion is a single differential equation of second order.
1.2 D’Alembert’s principle
1.2.1 Constraints and generalized coordinates
An alternative derivation of the Lagrange equations of motion is by explicitly considering constraints. Constraints are sometimes used to approximate some or all interactions in a molecular dynamics calculation. Consider for example a hydrocarbon molecule whose fastest vibrational motion typically occurs at C-H and C-C bonds. One often approximates C-H and C-C bonds as being rigid, i.e. fixing their length. This allows to choose larger time steps for the integration algorithm.
As another example, consider an extended solid body. While the atoms within this body interact via some interaction law, for extended macroscopic objects this motion manifests itself in the elastic deformation of a solid body. It is often desirable to ignore the fact that bodies can deform elastically and model them as rigid bodies. Effectively, all atoms inside this body are then constrained at a certain distance from each other.
We distinguish two types of constraints. Holonomic constraints can be expressed in the form \begin {equation} \phi \left ( \v {r}_{1},\v {r}_{2},\v {r}_{3},\ldots ,\ t \right ) = 0. \end {equation} For example, fixing the distance between particle \(i\) and particle j to \(c_{\text {ij}}\) can be expressed as \(\left ( \v {r}_i - \v {r}_j \right )^{2} - c_{\text {ij}}^{2} = 0\).
If \(k\) holonomic constraints are present, this decreases the total number of degrees of freedom of our \(N\)-particle system from \(3N\) to \(3N - k\). In other words, we can use the \(k\) expressions of type Eq. (2.1) to eliminate the dependency on \(k\) of the initial Cartesian degrees of freedom. This elimination can also be expressed by the introduction of \(3N - k\) independent variables \(q_{1},q_{2},\ldots ,q_{3N - k}\) called generalized coordinates. We can always express the Cartesian degrees of freedom as a function of these generalized coordinates, i.e. \begin {equation} \v {r}_{1} = \v {r}_{1}\left ( q_{1},q_{2},\ldots ,q_{3N - k},t \right ) \end {equation} As a simple example, consider a pendulum where a bead is moving at the end of a stick of length a in a two-dimensional plane. The motion the bead can be expressed by two Cartesian coordinates, \(x\left ( t \right )\) and \(y\left ( t \right )\) with the additional constraints \(x^{2} + y^{2} = a^{2}.\) However, we already know that the only coordinate that changes during the motion of the bead is its angle \(\theta \left ( t \right )\), which the becomes the generalized coordinate. The Cartesian coordinates are then straightforwardly expressed as \(x\left ( t \right ) = a\cos {\theta \left ( t \right )}\) and \(y\left ( t \right ) = a\sin {\theta \left ( t \right )}\).
Constraints that cannot be expressed in the way described above are called nonholonomic. An example would be the wall of a cylinder that cannot be passed by particles inside. This constrained could be expressed as \(r_i^{2} - a^{2} \leq 0\), where \(a\) is the radius of the cylinder. Another example is a cylinder of radius a rolling on a surface. If the center of the cylinder is a position \(\v {r}\), then the distance to the surface cannot become larger than its radius a.
For nonholonomic constraints, we cannot remove a degree of freedom from the description of the problem. Therefore, no general formal solution, such as the introduction of generalized coordinates for holonomic constraints, exists for the imposition of nonholonomic constraints.
1.2.2 D’Alembert’s principle
The existence of a constraint implies that there is a constraint force that acts implicitly to fulfill that constraint. We decompose the force on particle \(i\) into applied force, \(\v {F}_i^{a}\), and constraint force, \(\v {F}_i\), \begin {equation} \v {F}_i = \v {F}_i^{a} + \v {F}_i. \end {equation}
We now consider virtual displacements \(\delta \v {r}_i\) of our particles. Virtual displacements are displacements carried out at some instant in time t under the conditions valid at that instant; i.e. we ignore the fact that during a time interval \(dt\) forces and constraints can change.
The virtual work is carried out by the virtual displacements is given by \begin {equation} \delta \mathcal {W}^\text {s} = \sum _i \v {F}_i \cdot \delta \v {r}_i = \sum _i \v {F}_i^{a} \cdot \delta \v {r}_i + \sum _i\v {f}_i \cdot \delta \v {r}_i. \end {equation} Note that in (static) equilibrium \(\v {F}_i = 0\) and hence the virtual work vanishes, \(\delta \mathcal {W}^\text {s} = 0\).
We now consider systems for which the net virtual work of the forces of constraint, \(\sum _i\v {f}_i \cdot \delta \v {r}_i\) ,is zero. This means the virtual displacements are perpendicular to the constraint forces. For equilibrium we then have, \begin {equation} \delta \mathcal {W}^\text {s} = \sum _i \v {F}_i^{a} \cdot \delta \v {r}_i = 0 \label {eq:virtualwork} \end {equation} Equation \eqref{eq:virtualwork} is called the principle of virtual work, but it only describes equilibrium for static systems.
The corresponding principle for dynamical systems goes back to James Bernoulli and was further developed by Jean le Rond D’Alembert. Newton’s equation of motion can be written as \begin {equation} \v {F}_i - \dot \v {p}_i = 0. \end {equation} The condition for dynamic equilibrium can then be case into the virtual work formulation, \begin {equation} \delta \mathcal {W}^{d} = \sum _i\left ( \v {F}_i - \dot \v {p}_i \right ) \cdot \delta \v {r}_i = 0, \end {equation} or when again considering only virtual displacements perpendicular to the constraint forces, \begin {equation} \delta \mathcal {W}^{d} = \sum _i\left ( \v {F}_i^a - \dot \v {p}_i \right ) \cdot \delta \v {r}_i = 0. \label {eq:dalemberts} \end {equation} Equation \eqref{eq:dalemberts} is often referred to as D’Alemberts principle. In what follows, we will drop the superscript a for brevity.
We now introduce generalized coordinates. Virtual displacements in the Cartesian coordinates can then be expressed as \begin {equation} \delta \v {r}_i = \frac {\partial \v {r}_i}{\partial q_i}\delta q_i, \end {equation} and hence the virtual work becomes \begin {equation} \delta \mathcal {W}^\text {s} = \sum _i \v {F}_i \cdot \delta \v {r}_i = \sum _{ij} \v {F}_i \cdot \frac {\partial \v {r}_i}{\partial q_j}\delta q_j = \sum _j Q_j\delta q_j = 0. \end {equation} The \(Q_j\)’s are called generalized forces.
The second contribution to D’Alembert’s principles, Eq. (2.8), is given by \begin {equation} \delta \mathcal {W}^\text {s} - \delta \mathcal {W}^{d} = \sum _i{\dot {\v {p}_i} \cdot \delta \v {r}_i} = \sum _i{m_i\ddot {\v {r}_i} \cdot \delta \v {r}_i} = \sum _{\text {ij}}{m_i\ddot {\v {r}_i} \cdot \frac {\partial \v {r}_i}{\partial q_j}\delta q_j}, \end {equation} which can be expressed as \begin {equation} \delta \mathcal {W}^\text {s} - \delta \mathcal {W}^{d} = \sum _{\text {ij}}{m_i\left \lbrack \frac {d}{dt}\left ( \dot {\v {r}_i} \cdot \frac {\partial \v {r}_i}{\partial q_j} \right ) - \dot {\v {r}_i} \cdot \frac {d}{dt}\left ( \frac {\partial \v {r}_i}{\partial q_j} \right ) \right \rbrack \delta q_j}. \end {equation} Note that \begin {equation} \dot {\v {r}_i} \equiv \frac {d\v {r}_i}{dt} = \sum _j{\frac {\partial \v {r}_i}{\partial q_j}{\dot {q}}_j} + \frac {\partial \v {r}_i}{\partial t} \end {equation} and hence \begin {equation} \frac {\partial \dot {\v {r}_i}}{\partial \dot {q_j}} = \frac {\partial \v {r}_i}{\partial q_j}. \end {equation} We can therefore rewrite Eq. (2.12) as \begin {equation} \sum _{\text {ij}}{m_i\left \lbrack \frac {d}{dt}\left ( {\v {v}}_i \cdot \frac {\partial {\v {v}}_i}{\partial {\dot {q}}_j} \right ) - {\v {v}}_i \cdot \frac {\partial {\v {v}}_i}{\partial q_j} \right \rbrack \delta q_j} = \sum _{\text {ij}}{\left \{ \frac {d}{dt}\left \lbrack \frac {\partial }{\partial {\dot {q}}_j}\left ( \frac {1}{2}mv_i^{2} \right ) \right \rbrack + \frac {\partial }{\partial q_j}\left ( \frac {1}{2}mv_i^{2} \right ) \right \}\delta q_j}. \end {equation} Since \(T^{*} = \frac {1}{2}mv_i^{2}\), D’Alembert’s principle becomes \begin {equation} \delta \mathcal {W}^{d} = \sum _j{\left \{ \left \lbrack \frac {d}{dt}\left ( \frac {\partial T^{*}}{\partial {\dot {q}}_j} \right ) - \frac {\partial T^{*}}{\partial q_j} \right \rbrack - Q_j \right \}\delta q_j} = 0 \end {equation} Note that Eq. (2.13) must hold for all possible virtual displacements. The individual \(\delta q_j\) are not necessarily independent of each other for general constraints. For holonomic constraints, however, we can find coordinates \(q_j\) that automatically fulfill the constrained, the generalized coordinates. If this is the case, then we can vary each \(\delta q_j\) independently and the constraints are still fulfilled; Eq. (2.15) must hence be fulfilled independently for each of the summands.
If all constraints are implicitly contained in the set of generalized coordinates, then \begin {equation} \left \lbrack \frac {d}{dt}\left ( \frac {\partial T^{*}}{\partial {\dot {q}}_j} \right ) - \frac {\partial T^{*}}{\partial q_j} \right \rbrack - Q_j = 0. \end {equation} Note that if the forces can be derived from a potential, \(\v {F}_i = - \frac {\partial V}{\partial \v {r}_i}\), then \begin {equation} Q_j = \sum _i{\v {F}_i \cdot \frac {\partial \v {r}_i}{\partial q_j}} = - \sum _i{\frac {\partial V}{\partial \v {r}_i} \cdot \frac {\partial \v {r}_i}{\partial q_j}} = - \frac {\partial V}{\partial q_j} \end {equation} and Eq. (2.16) can be written as \begin {equation} \left \lbrack \frac {d}{dt}\left ( \frac {\partial L}{\partial {\dot {q}}_j} \right ) - \frac {\partial L}{\partial q_j} \right \rbrack = 0 \end {equation} with \(L = T^{*} - V\). The function \(L\) is called the Lagrangian and Eq. (2.18) is called the Euler-Lagrange equation or just the Lagrange equation.
Just like the kinetic energy and the potential energy, the Lagrangian is a state function. The generalized coordinates \(q_i\) do not need to have units of length and generalized forces \(Q_i\) do not need to have units of force, but the work \(Q_i\delta q_i\) must have units of energy (work). The utility of the Euler-Lagrange equation is that we get the equations of motion for the generalized coordinates without the need to transform into these coordinates via something like Eq. (2.13). It is important to realize that this is only valid if the potential \(U\) is independent of velocity. A notable difference are electromagnetic forces, see Goldstein chapter 1-5. Because Eq. (2.14) contains terms \({\v {v}}_i \cdot d{\v {v}}_i\) we get the kinetic coenergy \(T^{*}\) rather than the kinetic energy \(T\). Goldstein does not make this distinction but Williams does. Non-conservative forces simply remain as generalized forces in the Euler-Lagrange equation, \begin {equation} \left [ \frac {\dif }{\dif r}\left ( \frac {\partial L}{\partial \dot {q}_j} \right ) - \frac {\partial L}{\partial q_j} \right ] = Q_j^\text {nc} \end {equation} where \(Q_j^{\text {nc}}\) is the non-conservative contribution to the generalized force.