Molecular Simulation/Molecular Dynamics

Molecular dynamics simulation is a tool of simulating motions of the atoms of a many-body system. Molecular dynamics is used to compute the equilibrium and transport properties (viscosity, thermal conductivity, diffusion, reaction rate, protein folding time, structure and surface coating) of a classical system. To simulate the dynamics of molecules, classical Newtonian mechanics can be used. In practice, molecular dynamics simulations are performed by moving in small increments of time. It is the solution of the classical equations of motion for atoms and molecules to obtain the evolution of a system. Because molecular systems typically consist of a vast number of particles, it is impossible to determine the properties of such complex systems analytically, so MD is applied to these systems to solve this problem by using numerical methods.

Verlet Algorithm
One popular algorithm that predict positions of atoms at a later point in time is Verlet algorithm. This algorithm approximates the position of particles at time $$t + \Delta t$$ (one time step in future from time t) by using a Taylor series approximation.

A Taylor series approximates a function at point x based on the values and derivatives of that function at another point a.

$$ f(x) = f(a) + \frac{\textrm{d} f(a)}{\textrm{d}t} \left(x-a\right) + \frac{1}{2!} \frac{\textrm{d}^2 f(a)}{\textrm{d}t^2} \left(x-a\right)^2 + \frac{1}{3!} \frac{\textrm{d}^3 f(a)}{\textrm{d}t^3} \left(x-a\right)^3 + \ldots $$

Derivation of the Verlet Algorithm
A Taylor series approximates an equation for the position of the molecule in $$t + \Delta t$$.

$$ r_i(t + \Delta t) \approx r_i(t) + \frac{\textrm{d} r_i(t)}{dt}(t + \Delta t - t) + \frac{1}{2} \frac{\textrm{d}^2 r_i(t)}{\textrm{d}t^2}(t + \Delta t -t)^2 $$

$$ r_i(t + \Delta t) \approx r_i(t) + \frac{\textrm{d} r_i(t)}{\textrm{d}t} \Delta t + \frac{1}{2} \frac{\textrm{d}^2 r_i(t)}{\textrm{d}t^2} \Delta t^2 $$

First and second derivation of position correspond to velocity and acceleration, and acceleration is related to force through Newton’s second law.

$$ \frac{\textrm{d} r_i(t)}{\textrm{d}t} = v_i(t) $$

$$ \frac{\textrm{d}^2 r_i(t)}{\textrm{d}t^2} = a_i(t) = \frac{F_i(t)}{m_i} $$

$$ r_i(t + \Delta t) \approx r_i(t) + v_i(t) \Delta t + \frac{1}{2} \frac{F_i(t)}{m_i} \Delta t^2 $$

Considering one time step backward in time:

$$ r_i(t - \Delta t) \approx r_i(t) + \frac{d r_i(t)}{dt}(t - \Delta t - t) + \frac{1}{2} \frac{d^2 r_i(t)}{dt^2}(t - \Delta t -t)^2 $$

$$ r_i(t - \Delta t) \approx r_i(t) - v_i(t) \Delta t + \frac{1}{2} \frac{F_i(t)}{m_i} \Delta t^2 $$

Taking sum of positions at previous and next step:

$$ r_i(t + \Delta t) + r_i(t - \Delta t) = r_i(t) + v_i(t) \Delta t + \frac{1}{2} \frac{F_i(t)}{m_i} \Delta t^2 + r_i(t) - v_i(t) \Delta t + \frac{1}{2} \frac{F_i(t)}{m_i} \Delta t^2 $$

$$ r_i(t + \Delta t) + r_i(t - \Delta t) = 2r_i(t) + \frac{F_i(t)}{m_i} \Delta t^2 $$

Rearranges to:

$$

Derivation of the Velocity Verlet Algorithm
An alternative form of Verlet equation that involves velocities. By using Taylor series for later time step from current time:

$$ r_i(t + \Delta t) \approx r_i(t) + v_i(t) \Delta t + \frac{1}{2} \frac{F_i(t)}{m_i} \Delta t^2 $$

Using current time step from the next time step:

$$ r_i(t) \approx r_i(t + \Delta t) - v_i(t + \Delta t) \Delta t + \frac{1}{2} \frac{F_i(t + \Delta t)}{m_i} \Delta t^2 $$

By substitution:

$$

Time step
The time step for the most of molecular dynamics simulations is on the femtosecond scale which is the scale of chemical bond vibrations. Molecular dynamics simulations are limited by the highest frequency vibration and normally the time step should be ten times lower that the highest frequency. Generally, molecular dynamics simulations of organic molecules use a timestep of 1-2 fs ($${10^{-15}}$$ s) because of the fast vibrations of bonds containing hydrogen.

$$ \text{Number of time steps} = \frac{\text{length of simulation}}{\text{time step length}} $$

The first problem is that some processes, such as protein folding, require relatively long simulation times to observe in their entirety. To observe these processes in an MD simulation we need at least simulation with microseconds length.

$$ \text{Number of time steps} = \frac{10^{-6}}{10^{-15}} = 1 \text{ billion time steps} $$

This is a very big number of time steps to do, even for a computer. To solve this problem some constraints are applied. Many molecular dynamics simulations constrain all hydrogen containing bonds to be at equilibrium value. Because hydrogen is a very light atom and has very fast bond vibrations. By constraining these hydrogen containing bonds to be fixed at equilibration value, fast vibrations will be neglected.

Another problem is that the vibrational energy levels of the vibration of covalent bonds have large spaces between quantum energy levels so it is poorly described using classical mechanics, it is better to assume that bonds are rigid.

Time reversibility
An isolated system that is propagating through Newton's laws of motion should be time-reversible. This is to say, if you. In many cases molecular dynamics algorithms are time reversible. It means that by running the trajectory forward and then backward by the same number of steps, get back to the starting point by the same sequence of steps. So two MD simulation with the same initial conditions and time steps will follow exactly the same trajectories.

A second consequence of this is that the system is said to be deterministic. If you know the positions and velocities of the particles at a given point of time, and can calculate the forces acting on them in any position, in principle, you could integrate the laws of motion to determine the state of the system at any point in the future. In this way, we can determine the state of the system (i.e., positions and velocities) in the future based on the state we are now.

Microcanonical ensemble (NVE)


A molecular dynamics simulation with an integrator that approximates Newtonian dynamics (e.g., the Verlet algorithm) will conserve the total energy of the system (i.e., the sum of kinetic and potential energies remains constant at each time step. As a result, these simulations will sample the microcanonical ensemble (NVE), because the total energy of the system (E) is conserved.

The Verlet equations and all other integrators have some error associated with the truncation of higher order terms in the Taylor series expansion. This error is more significant when larger time steps are used. Further, there is numerical error that results from the rounding error by performing these calculations on finite-precision computers (e.g., 32 or 64 bit registers). As a result, the total energy of a molecular dynamics simulation can drift from its initial value over the course of a simulation.

Sampling Different Ensembles
To sample a different ensemble, the equations of motion must be modified. This is commonly described as performing a molecular dynamics simulation with a thermostat or barostat.

Canonical ensemble (NVT)
In the canonical ensemble, amount of substance (N), volume (V) and temperature (T) are constant. In this ensemble, the energy of endothermic and exothermic processes is exchanged with a thermostat. Several methods have been developed that modify the equations of motion of a molecular dynamics simulation so that they sample the canonical ensemble, such as the wikipedia::Andersen thermostat and wikipedia::Nosé–Hoover thermostat.

Isothermal–isobaric ensemble (NPT)
In the isothermal–isobaric ensemble, amount of substance (N), pressure (P) and temperature (T) are conserved. Configuration of particles are Boltzmann-weighted this ensemble. Also, a barostat is needed in this ensemble.

Visualizing the MD Simulations
Time evolution of positions and velocities of the particles of the system is called a trajectory. This can be rendered into an animation that shows the motion of the particles over the course of the simulation. VMD is a popular program to visualize molecular dynamics simulations.