Molecular Dynamics

We will now introduce a method that will allow us to sample the phase space. The method goes under the name of molecular dynamics (MD).

An isolated system of $N$ interacting classical particles obeys the Newton's equations of motion:

$$M \ddot{{\bf r}}_i = - \frac{\partial U(\{{\bf r})\}}{\partial {\bf r}_i} = {\bf f}_i.$$ (7.87)

Integrating  over time one obtains a trajectory $\{{\bf r}(t)\}$, represented by a collection of positions of all particles in the system as function of time. Since the forces are conservative, the total energy of the system (kinetic plus potential) is constant and the trajectory samples the microcanonical ensemble. If, after a sufficiently long time, the trajectory passes within any arbitrary distance from any state of the ensemble, then the system is said to be ergodic and time averages are equivalent to ensemble averages. Ergodicity is not always satisfied, and so care must be exercised when assessing the suitability of this method for sampling the phase space. In particular, a harmonic system is a clear example of a system that is not ergodic, as energy cannot be transferred between different normal modes. We shall come back to this point later.

We now discuss how to generate such a trajectory in a computer, which means that we have to discretise . Let us consider a Taylor expansion of the position ${\bf r}_i(t)$:

$${{\bf r}}_i(t+ \delta t) = {\bf r}_i(t) + \dot{\bf r}_i(t) \delta... ...r}_i(t) \delta t^2 + \frac{1}{3!} \dddot{\bf r}_i(t) \delta t^3 + o(\delta t^4)$$ (7.88)

and

$${{\bf r}}_i(t- \delta t) = {\bf r}_i(t) - \dot{\bf r}_i(t) \delta... ...}_i(t) \delta t^2 - \frac{1}{3!} \dddot{\bf r}_i(t) \delta t^3 + o(\delta t^4),$$ (7.89)

where $\delta t$ is a small time interval. If we sum  and  the odd powers of $\delta t$ cancel out and we are left with:

$${{\bf r}}_i(t + \delta t) + {{\bf r}}_i(t- \delta t) = 2 {\bf r}_i(t) + \ddot {\bf r}_i(t) \delta t^2 + o(\delta t^4),$$ (7.90)

which we rearrange as:

$${{\bf r}}_i(t + \delta t) = 2 {\bf r}_i(t) - {{\bf r}}_i(t- \delta t) + \frac{1}{m} {\bf f}_i(t) \delta t^2 + o(\delta t^4),$$ (7.91)

where we have replaced the acceleration $\ddot{{\bf r}}_i$ with the force divided by the mass, $\frac{1}{M} {\bf f}_i(t)$. The algorithm described in  is known as the Verlet algorightm. We see that we can obtain the positions at time $t+\delta t$ from knowledge of the positions at time $t$ and $t-\delta t$, and from the forces at time $t$. Eq.  can therefore be used to generate a trajectory in time, with an error that is proportional to $\delta t^4$. For a simulation of length $\tau$ divided in $L$ time steps the total error is proportional to $L \delta t^4 = \tau \delta t^3$ and can therefore easily be reduced by reducing the time step (at an increased computer cost). This error will mostly show up in a non perfect conservation of the total energy and, for this reason, it is common practice to monitor the behaviour of the total energy to assess the quality of the simulations. A rule of thumb for choosing a reasonable value of the time step is to take it to be $\approx 1/30^{th}$ of the the shortest period in the system. It should be pointed out, however, that if the sole purpose of the simulation is to sample phase space, a faithful conservation of the constant of motion is not necessarily a stringent requirement.

We can re-express Eq.  in terms of the velocities:

$${\bf v}_i(t) = \frac{{\bf r}_i(t+\delta t) - {\bf r}_i(t-\delta t)}{2 \delta t},$$ (7.92)

and obtain:

$${{\bf r}}_i(t + \delta t) = {\bf r}_i(t) + {\bf v}_i(t) \delta t + \frac{1}{2M} {\bf f}_i(t) \delta t^2 + o(\delta t^4),$$ (7.93)

which gives access to the positions at time $t+\delta t$ with knowledge only of the positions, velocities and forces at time $t$. This expression is particularly useful at the beginning of the simulation, where only the initial positions are available, and shows that to begin a simulation we also need to provide the initial velocities.

Because of the equipartition theorem, the temperature of the system can be obtained form the ensemble average of the kinetic energy, given by Eq. :

$$\langle E_k \rangle = \frac{3N}{2} k_{\rm B}T,$$ (7.94)

where $\langle E_k \rangle$ is the time average of the instantaneous kinetic energy $E_k(t)$, given by:

$$E_k(t) = \sum_i \frac{1}{2} M v^2_i(t),$$ (7.95)