Classical mechanics

Newton’s equations

You are all familiar with this formulation of classical mechanics, so we will go over it very briefly. In his masterpiece “Principia”, first published in 1687(!), Newton formulated his three celebrated laws of motion. The first two can be succinctly summarized in the equation,

(1)\[ m \ddot{ \textbf{r} } = \textbf{F}, \]

where \(m\) and \( \textbf{r}(t) = [x(t), y(t), z(t)] \) are the mass and the Cartesian position vector of the particle, respectively. The single and double dot represent the first and second time-derivatives, respectively, and \(\textbf{F}\) is the force acting on the particle. We will denote vectors by bold letters throughout the course. Eq. (1) is a second-order differential equation. To solve it, we require two initial conditions, namely, the position and velocity at \(t=0\). The solution of the equation specifying the position of a particle in time (which also gives the velocity by differentiation) is called a trajectory.

Eq. (1) shows that the acceleration produced by the force \(\textbf{F}\) acting on the particle is equal to the force divided by the particle’s mass, also known as Newton’s second law or Newton’s equation of motion. It also embodies the first law which states that in the absence of force, a body will remain at rest or move in a straight line with constant velocity. This can be seen by setting \(\textbf{F} = 0\) in Eq. (1) and integrating twice to get \(\dot{ \textbf{r} } = \textbf{v}\) and \(\textbf{r} = \textbf{r}(0) + \textbf{v}t\). Newton’s third law states that “for every action, there is an equal and opposite reaction”, i.e., that if \(\textbf{F}_{12}\) is the force that particle 1 exerts on particle 2, then

(2)\[ \textbf{F}_{21} = -\textbf{F}_{12}. \]

Active learning

Solve Newton’s equation of motion for a single particle of mass \(m\) in a one-dimensional harmonic trap with frequency \(\omega_0\). Use the initial conditions \(q(0) = q_0 \) and \( v(0) = v_0\).

For a collection for \(N\) particles, Newton’s equations of motion are

(3)\[ m \ddot{\textbf{r}}_j = \textbf{F}_j, \]

where \(j=1,...,N\) and the force acting on particle \(j\) depends on the coordinates of all particles, \(\textbf{F}_j(\textbf{r}_1,...,\textbf{r}_N)\). We define conservative forces as those that are the gradient some potential energy,

(4)\[ \textbf{F}_j = -\nabla U(\textbf{r}_1,...,\textbf{r}_N). \]

The gradient \(\nabla_j = [\frac{\partial}{\partial x_j}, \frac{\partial}{\partial y_j}, \frac{\partial}{\partial z_j}]\) is sometimes also denoted as \(\frac{\partial}{\partial \textbf{r}_j}\). The name conservative forces comes from the fact that the total energy is conserved throughout the classical trajectory for such forces.

Active learning

Prove that Newton’s equations conserve the total energy \(E = \frac{1}{2} \sum_{j=1}^N m_j \dot{ \textbf{r}} _j^2 + U(\textbf{r}_1,...,\textbf{r}_N)\). Hint: Show that its time derivative is zero.

We require the initial positions and velocities of all \(N\) particles to solve Eq. (3). In practice, however, the exact many-body forces are complicated, nonlinear functions and the resulting dynamics are highly complex. Any attempt to find an analytic solution to the equations of motion of a many-body system is therefore futile. Moreover, if we give up trying to obtain an analytic solution and are content with solving Eq. (3) numerically for a macroscopic system with \(\sim 10^{23}\) particles, it would require more computational resources than are currently (or in the foreseeable future) available. This infamous problem is called the “n-body problem” when the interaction between the particles is gravitational and it is still an area of active research in Astronomy. So, is all hope lost? No! The essence of MD simulations is to instead solve Eq. (3) numerically for a large, but much smaller, number of particles, typically ranging between \(100-10^6\), typically in a box of given volume under some specified boundary conditions (representing a solid, for example, more about this later). It turns out that this approach works remarkably well, and many properties converge to their bulk values for this much smaller number of particles. Before describing MD simulations in detail, we first introduce the concept of phase space, and then look into alternative formulations of classical mechanics which will be of use in the following lessons.

Phase space

The solution to Newton’s equations of motion specify the positions \(\textbf{r}(t)\) and velocities \(\textbf{v}(t)\) of the particles at every time, given we know their value at some instance (usually, \(t=0\)). In statistical mechanics and quantum mechanics it turns out to be more convenient to work with the momenta,

(5)\[ \textbf{p}_j(t) = m_j \textbf{v}_j(t) = m_j \dot{\textbf{r}}_j(t). \]

Eq. (3) can be written in terms of the momenta as well,

(6)\[ \dot{\textbf{p}}_j = \textbf{F}_j. \]

This equation already tells us something very interesting. When there is no total force acting on the particles,\(\sum_{j=1}^N \dot{\textbf{p}}_j = \sum_{j=1}^N \textbf{F}_j = 0\) and the total linear momentum \(\textbf{P} = \sum_{j=1}^N \textbf{p}_j\) is conserved.

At any moment in time, the microscopic state (“microstate”) of the system is specified by \(6N\) numbers - the components of the position and momentum vectors of all particles. These numbers can be considered as a point in a \(6N\) dimensional space, called phase space. We will denote a point in phase space as

(7)\[ \textbf{x}(t) = [\textbf{r}_1(t), ..., \textbf{r}_N(t),\textbf{p}_1(t),...,\textbf{p}_N(t)]. \]

A trajectory obtained from the solution of Newton’s equations of motion defines a parametric path in phase space. Visualizing it in a \(6N\) dimensional space is inconvenient, but we can do it for a single particle moving in one dimension. Let us look at the phase space path for a single particle of mass \(m\) in a one-dimensional harmonic trap with frequency \(\omega_0\). We will show that the phase space trajectory is an ellipse, whose size is determined by the initial conditions which define the conserved energy. The shape of the ellipse is determined by the mass of the particle and the trap frequency.

Active learning

Show that the trajectory for a single particle of mass \(m\) in a one-dimensional harmonic trap with frequency \(\omega_0\) is an ellipse in phase space. Use the initial conditions \(q(0)=q_0\) and \(v(0)=0\). Hint: Use energy conservation.

Now that we have mastered Newtonian Mechanics and the phase space representation of trajectories, we will move on to two equivalent formulations of classical mechanics which are more convenient for molecular simulations.

Lagrangian mechanics

We will first present the formalism of Lagrangian mechanics and only later show how to derive it, and all other formulations, from one general axiom called Hamilton’s principle.

We define the Lagrangian as the difference between the kinetic and potential energy of the system, expressed in terms of the Cartesian positions and velocities,

(8)\[ \mathcal{L} (\textbf{r}_1,...,\textbf{r}_N,\dot{\textbf{r}}_1,...,\dot{\textbf{r}}_N) = K(\dot{\textbf{r}}_1,...,\dot{\textbf{r}}_N) - U(\textbf{r}_1,...,\textbf{r}_N), \]

where the kinetic energy is given by

(9)\[ K(\dot{\textbf{r}}_1,...,\dot{\textbf{r}}_N) = \frac{1}{2} \sum_{j=1}^N m_j \dot{\textbf{r}}^2_j. \]

The classical equations of motion can be recovered from \(\mathcal{L}\) through the celebrated Euler-Lagrange equations,

(10)\[ \frac{\mathrm{d}}{\mathrm{d}t} \frac{\partial \mathcal{L}}{\partial \dot{\textbf{r}}_j} = \frac{\partial \mathcal{L}}{\partial \textbf{r}_j}. \]

Active learning

Show that for the Lagrangian of Eq. (8) and the Kinetic energy defined in Eq. (9), Euler-Lagrange equations are equivalent to Newton’s equations of motion for conservative forces.

So, why do we need Lagrangian Mechanics? In what way are they more useful than Newtonian Mechanics? The short (and insufficient) answer is that they allow us to write equations of motions easily when non-Cartesian coordinates are natural to the problem. The new coordinates are often called generalized coordinates. This is because the Euler-Lagrange equations, Eq. (10), do not change under such point transformations. Lagrangian mechanics also allow us to handle constraints in a much more straightforward manner than Newtonian dynamics. Finally, Lagrangian mechanics are convenient to prove a remarkable result - that conservation laws are related to symmetry! This is a topic of some depth, and we will barely scratch the surface with a simple example.

Consider a pendulum composed of a mass \(m\) attached to a string of fixed length \(l\) hanging off a surface by an angle \(\phi\). In Newtonian Mechanics, we would have first considered the different forces: gravity applies a force \(m\textbf{g}\) directed along the \(y\) axis and the string applies a radial force \(\textbf{T}\), which is unknown, to compensate for the radial component of the gravitational pull. In other words, the tension force is required to make sure the mass doesn’t move in the radial direction, i.e., that the constraint on the length of the spring is fulfilled. Then, we would write down the equations of motions for each force component and after some algebraic manipulations, we would get the equations for the dynamics of the angle \(\phi\) and an expression for the tension force. Let us see how simple this process gets with the Lagrangian formalism instead.

We define two generalized coordinates, \( x = l\sin{\phi}\) and \( y= l \cos{\phi}.\) The kinetic energy is then given by,

\[ K = \frac{1}{2} m (\dot{x}^2 + \dot{y}^2) = \frac{1}{2} m l^2 \dot{\phi}^2 (\cos^2{\phi} + \sin^2{\phi}) = \frac{1}{2} m l^2 \dot{\phi}^2, \]

and the potential energy by,

\[ U = m g (l-y) = m g l(1-\cos{\phi}). \]

The Lagrangian is then,

\[ \mathcal{L} = \frac{1}{2} m l^2 \dot{\phi}^2 - m g l(1-\cos{\phi}), \]

which results in the Euler-Lagrange equation for \(\phi\),

\[ \ddot{\phi} = - \frac{g}{l} \sin{\phi}. \]

This is a non-linear differential equation which is generally very hard to solve. What you did in your first-year mechanics course is to assume that \(\phi\) is small so that \(\sin{\phi} \approx \phi\), which leads to simple harmonic motion with frequency \(\omega = \sqrt{\frac{g}{l}}\). Look at how easy it was to derive it using the Lagrangian formulation! We now consider another formulation of classical mechanics, due to Hamilton, which you already encountered in quantum mechanics but probably were not aware of its utility in the classical context.

Hamiltonian mechanics

While Lagrangian mechanics is based on a description of the dynamics of a system in terms of the positions and velocities, the Hamiltonian formulation is intimately connected to the phase space representation we introduced in Section Phase space. We previously defined the Cartesian momentum, \(\textbf{p}_j = m_j \dot{\textbf{r}}_j\) in Eq. (5). A wider definition is to consider a generalized momentum,

(11)\[\begin{equation} \textbf{p}_j = \frac{\partial \mathcal{L}}{\partial \dot{\textbf{q}}_j}, \end{equation}\]

where \(\textbf{q}_j\) is some generalized coordinate. When we write the Lagrangian in Cartesian coordinates, \(\textbf{r}_j\), this definition reduces to Eq. (5), but it holds for any choice of coordinates! (Remember, this was a major strength also of the Lagrangian formulation).

The Hamiltonian is then defined1 as

(12)\[\begin{equation} \label{Hamiltonian_def} \mathcal{H}(\textbf{q}_1,...,\textbf{q}_N,\textbf{p}_1,...,\textbf{p}_N) = \sum_{j=1}^N \textbf{p}_j \cdot \dot{\textbf{q}}_j(\textbf{p}_j) - \mathcal{L} \left (\textbf{q}_1,...,\textbf{q}_N,\dot{\textbf{q}}_1(\textbf{p}_1),...,\dot{\textbf{q}}_N(\textbf{p}_N) \right ), \end{equation}\]

and the equations of motion are called Hamilton’s equations,

(13)\[ \dot{\textbf{q}}_j = \frac{\partial \mathcal{H}}{\partial \textbf{p}_j}, \dot{\textbf{p}}_j = - \frac{\partial \mathcal{H}}{\partial \textbf{q}_j} . \]

If we choose to work in Cartesian coordinates, \(\dot{\textbf{r}}_j(\textbf{p}_j) = \textbf{p}_j / m_j \), according to Eq. (5), we get

(14)\[ \mathcal{H}(\textbf{r}_1,...,\textbf{r}_N,\textbf{p}_1,...,\textbf{p}_N) = \sum_{i=1}^N \frac{\textbf{p}_j^2}{2 m_j} + U(\textbf{r}_1,...,\textbf{r}_N), \]

which shows that the Hamiltonian is just the total energy of the system in terms of the Cartesian positions and momenta,

(15)\[\begin{equation} \mathcal{H}(\textbf{r}_1,...,\textbf{r}_N,\textbf{p}_1,...,\textbf{p}_N) = K(\dot{\textbf{r}}_1,...,\dot{\textbf{r}}_N) + U(\textbf{r}_1,...,\textbf{r}_N). \end{equation}\]

Hamilton’s equations in this case also reduce to Newton’s equations of motion.

Active learning

Show that for the Hamiltonian of Eq. (14), Hamilton’s equations are equivalent to Newton’s equations of motion for conservative forces.

We are now in a position to appreciate a very beautiful result of the Hamiltonian formalism. It allows us to write a specific condition for when some quantity, \(f(\textbf{q}_1,...,\textbf{q}_N,\textbf{p}_1,...,\textbf{p}_N)\), a function of the phase space variables, is conserved. We write the time derivative of \(f\) using the chain rule as

\[ \frac{\mathrm{d} f}{\mathrm{d} t} = \sum_{j=1}^N \left[ \frac{\partial f}{\partial \textbf{q}_j} \cdot \dot{\textbf{q}}_j + \frac{\partial f}{\partial \textbf{p}_j} \cdot \dot{\textbf{p}}_j \right]. \]

Using Hamilton’s equations, we can write

(16)\[ \frac{\mathrm{d} f}{\mathrm{d} t} = \sum_{j=1}^N \left[ \frac{\partial f}{\partial \textbf{q}_j} \cdot \frac{\partial \mathcal{H}}{\partial \textbf{p}_j} - \frac{\partial f}{\partial \textbf{p}_j} \cdot \frac{\partial \mathcal{H}}{\partial \textbf{q}_j} \right] \equiv \left \{ f,\mathcal{H} \right \}, \]

where \(\left \{ f, g \right \}\) are termed the Poisson bracket between \(f\) and \(g\). Eq. (16) shows that a quantity \(f\) is conserved when its Poisson bracket with the Hamiltonian is zero. If you are scratching your head thinking, “Hey, this looks familiar!”, you are right. Indeed, the commutator between operators in quantum mechanics plays a similar role to that of the Poisson bracket between functions of phase space variables in classical mechanics. Unfortunately, we can’t discuss this in more detail. Since the Hamiltonian of Eq. (14) commutes with itself (and assuming it is time-independent, as we have been doing), this constitutes another proof of the conservation of energy along classical trajectories for conservative forces.

Active learning

Prove that for a system of \(N\) particles with no external forces acting on them, \(\sum_{j=1}^N \textbf{F}_j = 0\), the total momentum is conserved.

The action

We have so far presented the classical equations of motions, be it in the Hamiltonian or the Lagrangian formalism, without deriving them. The question arises - can we derive the Euler-Lagrange equations, and from them Hamilton’s and Newton’s equations, from some underlying general principle? The answer is positive.

To see how, we need to define the Action integral,

(17)\[ S\left[ \textbf{q}_1,...,\textbf{q}_N, \dot{\textbf{q}}_1,...,\dot{\textbf{q}}_N \right] = \int_{t_1}^{t_2} \mathcal{L} (\textbf{q}_1(t),...,\textbf{q}_N(t), \dot{\textbf{q}}_1(t),...,\dot{\textbf{q}}_N(t)) \mathrm{d}t. \]

For brevity, we will denote collectively \(\textbf{q}_1(t),...,\textbf{q}_N(t)\) as \(\textbf{Q}(t) \) and similarly for the time derivatives. Note that while the Lagrangian and Hamiltonian are functions of the coordinates and velocities or momenta, the action is a functional. A functional is a mathematical object that takes a function as input and maps it uniquely to some number, similarly to how a function takes one number as input and maps it uniquely to another number. Every \(N\)-particles trajectory, \(\textbf{Q}(t), \dot{\textbf{Q}}(t)\), that fulfills the boundary conditions, (i.e., starts from \(\textbf{Q}(t_1)\) with \(\dot{\textbf{Q}}(t_1)\) and ends at \(\textbf{Q}(t_2)\) with \(\dot{\textbf{Q}}(t_2))\), uniquely defines the time evolution of the Lagrangian which in turn translates to a specific value of the action. Alternatively, the difference between a functional and a function is that the action depends on all values of the trajectory between two points while the Lagrangian and Hamiltonian depend only on the value of the trajectory in some specific instance in time.

It turns out that the Euler-Lagrange equations can be derived by postulating that the classical trajectory makes the first variation of the action integral vanish. The first variation of a functional is an extension of the definition of derivatives to the world of functionals. We consider how the action changes when the trajectories undergo an arbitrary variation

(18)\[\begin{align} \textbf{q}_j(t,\alpha) &= \textbf{q}_j(t,0) + \alpha \textbf{n}_j(t), & \dot{\textbf{q}}_j(t,\alpha) &= \dot{\textbf{q}}_j(t,0) + \alpha \dot{\textbf{n}}_j(t), \end{align}\]

and denote the action for such trajectories as

(19)\[\begin{equation} S(\alpha) = \int_{t_1}^{t_2} \mathcal{L} (\textbf{q}_1(t,\alpha),...,\textbf{q}_N(t,\alpha), \dot{\textbf{q}}_1(t,\alpha),...,\dot{\textbf{q}}_N(t,\alpha)) \mathrm{d}t. \end{equation}\]

The first variation of \(S\) is defined as

(20)\[ \delta S = \frac{\mathrm{d} S}{\mathrm{d} \alpha} \Bigg|_{\alpha=0} \mathrm{d} \alpha = \int_{t_1}^{t_2} \frac{\mathrm{d \mathcal{L}}}{\mathrm{d} \alpha} \Bigg|_{\alpha=0} \mathrm{d} \alpha \mathrm{d}t, \]

and we search for the conditions that satisfy \(\delta S = 0.\)

Taking the derivative of the Lagrangian we get

\[ \frac{\mathrm{d} \mathcal{L}}{\mathrm{d} \alpha} = \sum_{j=1}^N \left[ \frac{\partial \mathcal{L}}{\partial \textbf{q}_j} \cdot \frac{\partial \textbf{q}_j}{\partial \alpha} + \frac{\partial \mathcal{L}}{\partial \dot{\textbf{q}}_j} \cdot \frac{\partial \dot{\textbf{q}}_j}{\partial \alpha} \right]. \]

The derivative of \(S\) is just the time-integral of the derivative of \(\mathcal{L}\). If we realize that \(\frac{\partial \dot{\textbf{q}}_j}{\partial \alpha} = \frac{\partial^2 \textbf{q}_j}{\partial t \partial \alpha}\), and integrate the second term by parts, we get

\[ \frac{\mathrm{d} S}{\mathrm{d} \alpha} = \frac{\partial \mathcal{L}}{\partial \dot{\textbf{q}}_j} \cdot \frac{\partial \textbf{q}_j}{\partial \alpha} \Bigg|_{t_1}^{t_2} + \sum_{j=1}^N \int_{t_1}^{t_2} \left[ \frac{\partial \mathcal{L}}{\partial \textbf{q}_j} - \frac{\mathrm{d}}{\mathrm{d} t} \frac{\partial \mathcal{L}}{\partial \dot{\textbf{q}}_j} \right] \cdot \frac{\partial \textbf{q}_j}{\partial \alpha} \mathrm{d}t. \]

The first term vanishes because the path variations at the endpoints have to be zero to fulfill the boundary conditions. Plugging the second term into Eq. (20) gives

(21)\[ \delta S = \sum_{j=1}^N \int_{t_1}^{t_2} \left[ \frac{\partial \mathcal{L}}{\partial \textbf{q}_j} - \frac{\mathrm{d}}{\mathrm{d} t} \frac{\partial \mathcal{L}}{\partial \dot{\textbf{q}}_j} \right] \cdot \delta \textbf{q}_j \mathrm{d}t = 0, \]

where \( \delta \textbf{q}_j = \left( \frac{\partial \textbf{q}_j}{\partial \alpha} \right)_{\alpha=0} \mathrm{d} \alpha\). But since the path variations \(\delta \textbf{q}_j\) are independent, Eq. (21) has to vanish for each \(j\) separately. This is achieved by having the term in square brackets equal to zero which gives the Euler-Lagrange equations (10).

This result tells us that out of all the possible paths that connect \(\textbf{Q}(t_1), \dot{\textbf{Q}}(t_1)\) and \(\textbf{Q}(t_2), \dot{\textbf{Q}}(t_2)\) the classical path is the one for which the action integral is stationary (also known as Hamilton’s principle).

Is this postulate involving the first variation of the action more fundamental than the three laws of Newton or the Euler-Lagrange equations? In the last project of this course, some of you will work on incorporating quantum effects into your classical MD simulations. You will then encounter an ingenious formulation of quantum mechanics, due to Nobel Laureate Richard P. Feynman, which beautifully connects classical and quantum mechanics on the same framework involving the action integral. Its core observation is that while in classical mechanics there is one path that connects the initial and final points, the one that makes the action stationary, in quantum mechanics we have to sum over all possible paths - called a path integral.


1

We will not go into the details, but the Hamiltonian is defined as a Legendre transform of the Lagrangian with respect to the velocities.