Molecular Dynamics (MD)

Created in June 01, 2025

2025   ·   physics   MD   simulation   ·   ComputationalPhysics

Basic Idea

Molecular dynamics (MD) is a computer simulation technique that allows one to predict the time evolution of a system of interacting particles (atoms, molecules, granules, etc.). First we have to specify:

  • A set of initial conditions (positions & velocities)
  • Interaction potential for deriving the forces among all the particles.

Second, the evolution in time can be followed by solving a set of classical equations of motion for all particles in the system. Within the framework of classical mechanics, the equations are: \begin{equation} m_i\frac{d^2 \vec r_i}{dt^2}=\vec F_i \end{equation} for the \(i^{\text{th}}\) particle.

If the particles of interest are atoms, the force acting on the \(i^{\text{th}}\) atom at a given time can be obtained from the interatomic potential \(U\) : \begin{equation} \vec F_i = -\vec \nabla_i U(\vec r_1,\vec r_2,\vec r_3,…,\vec r_{N_at}) \end{equation}

Advantages of MD

  • The only input in the model - description of interatomic interaction
  • No assumptions are made about the process
  • Provides a detailed molecular/atomic-level information
  • The main strengths of MD is the ability to study fast non-equilibrium processes with atomic-level resolution.

Limitations of Classic MD

  • Usage of classical forces - ignoring quantum mechanical laws - de Brogile wavelength \(\lambda = \frac{h}{p}\)
  • Quantum effects become important when \(T\) is quite low - Thermal de Brogile wavelength \(\lambda_{th} = \frac{h}{\sqrt{2\pi mk_B T}}\)
  • Time-scale limitations
  • length-scale limitations

MD Algorithm

The basic steps of MD algorithm is as following:

  1. Set the initial conditions: \(\vec r_i(t_0),\vec v_i(t_0,...)\)
  2. Update neighbour list
  3. Get forces \(\vec F_i(t)\)
  4. Solve equations of motion over \(\delta t\)
  5. Perform \(p,T\) control
  6. Update time \(t=t+\delta t\)
  7. Calculate physical quantities
  8. Judge \(t=t_{\text{max}}\)?

1. Setup

Before starting the simulation, several conditoins and parameters should be setup.

Boundary Condition

In practice in most cases the atoms are arranged in an orthogonal simulation cell. There are several tpes of B.C.:

  • Open boundary
    • for a molecule or nanocluster in vacuum
    • not for a continuous medium
  • Fixed boundary
    • fixed boundary atoms
    • completely unphysical
  • Periodic boundary
    • obtaining bulk properties

Initial Velocities

We use the Maxwell-Boltzmann distribution to set the cell temperature to a desired value.

The probability of finding a particle with speed \(v\): \begin{equation} f(v)=\sqrt{(\frac{m}{2\pi k_B T})^3} 4\pi v^2 e^{-\frac{mv^2}{2k_B T}}. \end{equation} Then we can generate random initial atom velocities scaling \(T\) with equipartition theorem: \begin{equation} \frac{3}{2}k_B T=\frac{1}{2}mv^2. \end{equation}

Don’t forget to set the total momentum of the cell to zero

Initial Velocities:

subroutine init
sumv=0
sumv2=0
do i=1,npart
    x(i)=lattice_pos(i)
    v(i)=(ranf()-0.5)
    sumv=sumv+v(i)
    sumv2=sumv2+v(i)**2
enddo
sumv=sumv/npart
sumv2=sumv2/npart
fs=sqrt(3*temp/sumv2)
do i=1,npart
    v(i)=(v(i)-sumv)*fs
    xm(i)=x(i)-v(i)*dt
enddo
return
end

Choosing MD Time Step

  • Too long \(\Delta t\): energy is not conserved
  • \(\Delta r/\Delta t > 1/20\) of the nearest atom distance
  • In practice \(\Delta t \approx 4\text{fs}\)
  • MD is limited to <~100ns

2. Constructing Neighbor

To save CPU time, we can set a cutoff distance \(r_{cut}\) as Calculating force at \(r_{ij}>r_{cut}\) just wastes time. We can achieve this by constructing neighbor list and cell list.

Neighbor List

Looking for atoms within the cutoff distance at every step is quit time-consuming. However, since atoms move within a time step only \(<0.2 \dot{A}\), the local neighbors of a given atom remain the same for many time steps.

Verlet Neighbor List
  • Containing all neighbor atoms within \(r_m\)
  • If updating every \(N_m\) steps, the skin depth \(r_m-r_{cut} > N_m\Delta t V_{typ}\), where \(V_{typ}\) is the typical speed of an atom.
when to update?

A static value of \(N_m\) may lead to problems when there can be energetic process in the system so that \(v_j >> v_{typ}\).

A better way is to have a dynamic way of determing when to update the neighbor list:

  1. Keep track of the two largest atom displacements from the time when the neighbor list was updated the last time.
  2. Update the list when \(r_{\text{max},1}+r_{\text{max},2}>r_m-r_c\).

Cell List

We can divide MD cell into smaller subcells \(n\times n\times n\). The length of subcell \(l\) is chosen s.t. \begin{equation} l=\frac{L}{n}>r_{cut} \end{equation}

By this we can:

  • reduce the atom pairs from \(N(N-1)/2\) to \(27N^2/n^3\),
  • reduce complexity from \(O(N^2)\) to \(O(N)\).

3. Time Integration Algorithm

Among those methods to solve the second-order ODE, we choose Verlet Algorithm due to energy conservation.

Verlet Algorithm

We have the ieration equation: \begin{equation} \vec r (t+\Delta t)=2\vec r (t)-\vec r (t-\Delta t)+\vec a (t) \cdot (\Delta t)^2+o[(\Delta t)^4] \end{equation} where \(\vec a(t)=d^2\vec r /dt^2 = -\nabla V(\vec r(t))/m\). (See Page “ODE”)

Velocity Verlet: ABetter Implementation

In Verlet Algorithm, the error associated to \(\vec v(t)\) is of order \(o(\Delta t^2)\) rather than \(o(\Delta t^4)\). Velocity Verlet scheme, in which the positions, velocities and accelerations at time \(t+\Delta t\) are obtained, is devel

oped to overcome this difficulty.

The detailed steps are:

  1. Compute half next step velocity: \(\vec v(t+\frac{\Delta t}{2})=\vec v(t)+\frac{1}{2}\vec a(t)\Delta t\)
  2. Compute next position: \(\vec r(t+\Delta t) = \vec r (t)+\vec v(t)\cdot \Delta t +\frac{1}{2}\vec a(t)\cdot (\Delta t)^2\)
  3. Compute next acceleration: \(\vec a(t+\Delta t)=-\frac{1}{m} \nabla V(\vec r(t+\Delta t))\)
  4. Compute accurate next velocity: \(\vec v(t+\Delta t)=\vec v(t+\frac{\Delta t}{2})+\frac{1}{2}\vec a(t+\Delta t)\cdot \Delta t\)

4. Force Model

Lennard-Jone Potential

The Lennard-Jones Potential between two atoms has the following expression: \begin{equation} U_{LJ}(r) = 4\epsilon [(\frac{\sigma}{r})^12-(\frac{\sigma}{r}^6)]. \end{equation} And the force: \begin{equation} F_{LJ}(r)=-\frac{dU(r)}{dr}=24\frac{\epsilon}{\sigma}(2(\frac{\sigma}{r})^{13}-(\frac{\sigma}{r})^7). \end{equation}

  • The term \(~1/r_{ij}^{12}\) is the repulsion between atoms when they are brought close to each other, related to the Pauli principle.
  • The term \(~1/r_{ij}^6\), dominating at large distance, constitutes the attractive part and describes the cohesion to the system. The \(1/r^6\) attraction describes van der Waals dispersion forces.

The potential provides a good description of van der Walls interaction in inert gases and molecular systems. It is also often used in simulations when the objective is to model a general class of effects and the only requirement is to have a physically reasonable potential.

how to cutoff?

To make the potential energy become zero at the cutoff distance, we can shift the potential: \begin{equation} U_{\text{shift}}(r_{ij})=U(r_{ij})-U(R_c), r_{ij}<R_c \end{equation}

Most modern potentials for real materials go to zero at \(R_c\) together with several first derivatives of the potential function to avoid physical quantities sudden change.

limitations of pair potentials

Pair potentials work well at systems:

  • Inert gases: L-J model
  • Ionic materials: rigid ion model, shells model
  • Simple metals (s-p bonded)
  • In general for long-range forces: Coulombic, vdW

Pair potentials fail at systems:

  • Transition metals (d orbitals bonding)
  • Semiconductors (covalent bonds)
  • Polymers (covalent bonds)
  • In general for defects, surface configurations, vacancies

Stillinger-Weber Potential

The S-W potential is on of the first potentials for diamond lattices (e.g. Si, GaAs, Ge, C). Directional bounding is introduced in the S-W potential through an explict three-bodu term of the potential energy expansion: \begin{equation} U_2(r_{ij})=A(Br_{ij}^{-p}-r_{ij}^{-q})\text{exp}[\frac{c}{r_{ij}-r_c}] \end{equation} \begin{equation} U_3(r_{ij},r_{ik}) = h(r_{ij},r_{ik},\theta_{jik})+ \cdots \end{equation} where \(h(r_{ij},r_{ik},\theta_{jik})=\lambda \text{exp}[\frac{\gamma}{r_{ij}-r_c}+\frac{\gamma}{r_{ik}-r_c}](\cos{\theta_{jik}}-\beta)^2\)

Embedded Atom Method (EAM)

Pair potentials canot provide an adequate description of metallic systems. An alternative simple but rather realistic approach to the description of bonding in metallic systems is based on the concept of local density that is considered as the key variable: \begin{equation} E_{total}=\sum_i F(\rho_i)+\frac{1}{2}\sum_{i,j(i\neq j)} \phi(r_{ij}), \rho_i=\sum_j f(r_{ij}). \end{equation}

  • The embedding energy \(F\) is the energy associated with placing an atom in the electron environment described by \(\rho\).
  • The pair-potential term \(\phi\) describes electrostatic contributions.
limitations of EAM
  • Bonding is spherical
  • Potential is not unique
  • Limitations in alloys

5. Ensembles

Temperature in MD

The instantaneous value of the temperature is related to the kinetic energy via the particles’ momenta as follows: \begin{equation} \sum_{i=1}^N \frac{|\vec p_i|^2}{2m_i}=\frac{k_B T}{2}(3N-N_c), \end{equation} where \(N_c\) is the number of constraints and so \(3N-N_c\) is the total number of degrees of freedom.

Pressure in MD

Let us consider a system of \(N\) atoms that is evolving in a finite space and let us introduce a function that is called Clausius virial function: \begin{equation} W^{Tot}(\vec r_1,…,\vec r_N)=\sum_{i=1}^N \vec r_i \cdot \vec F_i^{Tot} \end{equation} where \(\vec r_i\) is the position of atom \(i\), \(\vec F_i^{Tot}\) is the total force acting on atom \(i\).

Averaging over the MD trajectory and using Newton’s law, we obtain \begin{equation} \langle W^{Tot} \rangle = \lim_{\tau \to \infty} \frac{1}{\tau} \int_0^\tau \sum_{i=1}^N \vec r_i (t)\cdot m_i\vec{\ddot{r_i}}(t)dt. \end{equation} If the system is localized in a finite region of space, then we have: \begin{equation} \langle W^{Tot}\rangle = -3Nk_B T. \end{equation} As the total force acting on atom \(i\) is composed of internal force and external force from the container walls, the total viral function can be written as: \begin{equation} \langle W_{tot}\rangle = \langle W_{int}\rangle+\langle W_{ext}\rangle = -3Nk_BT, \end{equation} with \(\langle W_{ext}\rangle=-3PV\). Therefore we have: \begin{equation} \langle \sum_{i=1}^N \vec r_i \cdot \vec F_i^{Int}\rangle - 3PV=-3Nk_B T. \end{equation} This equation is known as the viral equation. With this we can calculate \(P\): \begin{equation} P=\frac{Nk_BT}{V}+\frac{1}{3V}\langle \sum_{i=1}^N \vec r_i \cdot \vec F_i^{Int} \rangle. \end{equation}

Thermodynamic Ensembles

(For more see page “Statistical Ensembles”) As in thermodynamics, the ensembles are denoted by letters which indicate which physical quantities are conserved.

  • Microcanonical: \(NVE\) constant
  • Canonical: \(NVT\) constant
  • Isothermal-isobaric: \(NPT\) constant
  • Grand canonical: \(\mu VT\) constant

In MD only the microcanonical ensemble is obtained directly by solving the EOM. For other ones, on has to generate equations of motion which behave according to the desired ensemble \(\rho_{ens}(\Gamma)\)

Thermostats

Velocity Rescaling

Velocity-rescaling method is actually the first method proposed to keep the temperature a fixed value during a simulation without allowing fluctuations of \(T\).

  • Straightforward to implement

  • Good for use in warmup/initialization phase

  • Leads to discontinuities in the momentum part of the phase space trajectory
  • Does not allow the proper temperature fluctuations in canonical ensembles
  • No time reversibility
Berendsen Thermostat

In Berendsen’s method, the velocities are scaled at each time step, such that the rate of change of temperature is proportional to the difference in temperature: \begin{equation} \Delta T(t)=\frac{\delta t}{\tau}(T_0 - T(t)) \end{equation} where \(\tau\) is the coupling parameter.

The scaling factor (on velocity) is \begin{equation} \lambda = \sqrt{1+\frac{\delta t}{\tau}[\frac{T_0}{T(t)}-1]}. \end{equation}

  • It allows the temperature fluctuations, thereby not fixing it to a constant.
  • But the fluctuations are not captured correctly. It does not generate a correct canonical ensemble.