Molecular Simulation/Molecular Dynamics

Molecular Dynamics Simulation of Argon

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 {dr_{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:

The Verlet Algorithm $r_{i}(t+\Delta t)=2r_{i}(t)-r_{i}(t-\Delta t)+{\frac {F_{i}(t)}{m_{i}}}\Delta t^{2}$

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:

The Velocity Verlet Algorithm $v_{i}(t+\Delta t)=v_{i}(t)+{\frac {\Delta t}{2m_{i}}}\left[F_{i}(t)+F_{i}(t+\Delta t)\right]$

See Also

Wikipedia:Verlet integration Chemwiki: Lennard-Jones

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)

Conservation of the total kinetic and potential energy in a molecular dynamics simulation of liquid argon using the Verlet equation

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

Molecular Dynamics Simulation of POPC bilayer

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.

See Also

Wikipedia:Molecular Dynamics

Tuckerman Notes: The Hamiltonian formulation of classical mechanics