Molecular dynamics simulations

Molecular dynamics simulations#

Classical molecular dynamics simulations involve numerically integrating Newton’s equations of motion for \(N\) particles interacting through a pair potential \(V(\mathbf{r}_1,\mathbf{r}_2)\):

\[ m_i \ddot{\mathbf{r}}_i = - \sum_{j \neq i} \nabla_i V(\mathbf{r}_i,\mathbf{r}_j)~. \]

If the system is sufficiently large, one can use MD simulations to study its macroscopic (thermodynamic) properties such as the dependence of the pressure on temperature and density – the equation of state.

In many case the pair potential depends only on the distance \(r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|\). In this case the equations of motion read

\[ m_i \ddot{\mathbf{r}}_i = - \sum_{j \neq i} \frac{d V(r_{ij})}{d r_{ij}} \, \frac{\mathbf{r}_i - \mathbf{r}_j}{r_{ij}}~. \]

Lennard-Jones fluid#

Lennard-Jones fluid corresponds to a pair potential

\[ V_{LJ}(r) = 4 \varepsilon \left[ (\sigma/r)^{12} - (\sigma/r)^6 \right]. \]

The Lennard-Jones system has a rich phase diagram containing various phases (gas, liquid, solid) and the associated phase transitions.

LJPD

To evaluate the properties of the LJ fluid, one would run molecular dynamics simulations for sufficiently long time to let the system equilibrate, and calculate its macroscopic properties like pressure as time average. The simulations are typically performed in a box with periodic boundary conditions to minimize finite-size effects.

Velocity Verlet method#

For long molecular dynamics simulations it makes sense to use stable, time-reversal numerical schemes that conserve energy. The leapfrog method is a good choice. When applied to Newton’s equations of motion it is known as Verlet method.

Suppose we are solving Newton’s equations of motion

\[ \frac{d^2 x}{dt^2} = f(x,t). \]

These are rewritten as

\[\begin{align*} \frac{dx}{dt} &= v, \\ \frac{dv}{dt} &= f(x,t). \end{align*}\]

If we apply the leapfrog scheme to this system of equations, it will look like

\[\begin{align*} x(t+h) & = x(t) + h v(t + h/2), \\ v(t+3h/2) &= v(t+h/2) + h f[x(t+h), t+h)], \end{align*}\]

i.e. the coordinates are evaluated at full steps using velocity estimates at half-steps, and vice versa.

This is the essence of the Verlet method. Its feature is that we do not need to keep track of coordinates at half-steps (only the velocities) and of velocities at full steps. In case we are interested in velocity values also at full step, these can be estimated as \(v(t+h) = v(t+h/2) + (h/2) f[x(t+h),t+h]\).

The full algorithm is as follows. Given initial values \(x(t)\) and \(v(t)\) one first computes

\[ v(t+h/2) = v(t) + \frac{h}{2} f[x(t),t]. \]

Then each step of the Verlet algorithm corresponds to

\[\begin{align*} v(t+h/2) & = v(t) + \frac{h}{2} f[x(t),t],\\ x(t+h) & = x(t) + h v(t+h/2), \\ v(t+h) & = v(t+h/2) + \frac{h}{2} f[x(t+h),t+h)]. \end{align*}\]

The algorithm is straightforwardly generalizable to more than one variable.

Let us first write a code for integrating Newton’s equations of motion for arbitrary central potential

import numpy as np

# Simulation parameters (can be changed later)
n_particles = 125      # Number of particles
density     = 0.8      # Density of the fluid
box_length = (n_particles/density)**(1/3.)  # Length of simulation box (determined by the number of particles and density)
temperature0 = 1.0    # Initial (desired) temperature of the system
time_step = 0.004     # Time step for the simulation
num_steps = 1000      # Number of simulation steps


# Initialization
# Initial position at the nodes of the cube
# Helps to avoid particle overlap
def initial_positions():
    ret = np.zeros((n_particles,3))
    Nsingle = np.ceil(n_particles**(1/3.))
    dL = box_length / Nsingle
    for i in range(n_particles):
        ix = i % Nsingle
        iy = (np.trunc(i / Nsingle)) % Nsingle
        iz = np.trunc(i / (Nsingle * Nsingle))
        ret[i][0] = (ix + 0.5) * dL;
        ret[i][1] = (iy + 0.5) * dL;
        ret[i][2] = (iz + 0.5) * dL;
    return ret

# Computes forces for a given vector of positions, interaction potential and its gradient
# Return a tuple: positions, total potential energy, and the virial part of the pressure
def compute_forces(forces, positions, potential, potential_gradient):
    # forces = np.zeros_like(positions)
    forces.fill(0.)
    potential_energy = 0.0
    virial = 0.0
    for i in range(n_particles):
        for j in range(i+1, n_particles):
            # Vector of relative distance
            r_ij = positions[i] - positions[j]
            # Periodic boundary conditions (minimum-image convention)
            r_ij = r_ij - box_length*np.round(r_ij/box_length)
            
            r_sq = np.sum(r_ij**2)
            f_ij = -potential_gradient(r_sq) * r_ij
            
            forces[i] += f_ij
            forces[j] -= f_ij
            potential_energy += potential(r_sq)
            virial += np.dot(f_ij, r_ij)
            
    virial = virial/(3.0*box_length**3)
    return forces, potential_energy, virial


# Compute the kinetic temperature from velocities (kinetic energy) 
def compute_kinetic_temperature(velocities):
    return np.sum(velocities**2) / 3. / n_particles

# Renormalize the velocities to match the desired kinetic temperature 
def renormalize_velocities(velocities, temperature = temperature0):
    Tkin = compute_kinetic_temperature(velocities)
    factor = np.sqrt(temperature / Tkin)
    velocities *= factor


# Apply the velocity Verlet time step
# Returns the tuple of new positions, velocities, accelerations (forces), potential energy and pressure
def velocity_verlet(positions, velocities, accelerations, 
                    time_step, potential, potential_gradient):
    # Update positions
    positions += velocities*time_step + 0.5*accelerations*time_step**2
    positions = positions - box_length*np.floor(positions/box_length)
    # Update velocities
    velocities_half = velocities + 0.5*accelerations*time_step
    # Compute new forces and potential energy
    accelerations, potential_energy, pressure = compute_forces(accelerations, positions, potential, potential_gradient)
    # Update velocities using new accelerations
    velocities = velocities_half + 0.5*accelerations*time_step
    # Add ideal gas contribution to the pressure
    kinetic_temperature = compute_kinetic_temperature(velocities)
    pressure += density * kinetic_temperature
    return positions, velocities, accelerations, potential_energy, pressure

Now the Lennard-Jones potential properties

# Lennard-Jones potential as a function of squared distance
def lj_potential(r_sq):
    r6 = r_sq**3
    r12 = r6**2
    return 4.0*(1./r12 - 1./r6)

# The gradient term dV/dr / r in the rhs of Newton's equations
# for the LJ potential
def lj_potential_gradient(r_sq):
    r6 = r_sq**3
    r12 = r6**2
    return -24.0*(2./r12 - 1./r6) / r_sq

Microcanonical ensemble simulation#

Let us perform the simulation. Given that Newton’s equations of motion conserve energy, we are simulating the microcanonical ensemble.

# Simulation parameters
n_particles = 64       # Number of particles
density     = 0.1      # Density of the fluid
box_length = (n_particles/density)**(1/3.)  # Length of simulation box (determined by the number of particles and density)
temperature0 = 1.0     # Initial (desired) temperature of the system
time_step = 0.01       # Time step for the simulation
num_steps = 1000       # Number of simulation steps

averaging_steps = 100  # Number of steps to compute time average over at the end

# Initial conditions
positions = initial_positions()
velocities = np.random.normal(loc=0.0, scale=np.sqrt(temperature0), size=(n_particles, 3))
accelerations = np.zeros((n_particles, 3))
t = 0

renormalize_velocities(velocities, temperature0)
accelerations, potential_energy, pressure = compute_forces(accelerations, positions, lj_potential, lj_potential_gradient)
kinetic_energy = 0.5*np.sum(velocities**2)
total_energy = potential_energy + kinetic_energy
temperature = temperature0

times = []
kinetic_energies = []
potential_energies = []
total_energies = []
temperatures = []
compressibilities = []

times.append(t)
kinetic_energies.append(kinetic_energy / n_particles)
potential_energies.append(potential_energy / n_particles)
total_energies.append(total_energy / n_particles)
temperatures.append(temperature)
compressibilities.append(pressure / (temperature * density))

# Main loop
for step in range(num_steps):
    # Perform Velocity Verlet integration
    positions, velocities, accelerations, potential_energy, pressure = velocity_verlet(positions, velocities, accelerations, time_step, lj_potential, lj_potential_gradient)
    # Compute kinetic energy
    kinetic_energy = 0.5*np.sum(velocities**2)
    # Compute total energy
    total_energy = potential_energy + kinetic_energy
    # Compute temperature and pressure
    temperature = 2.0*kinetic_energy/(3.0*n_particles)
    # Update time
    t += time_step
    
    times.append(t)
    kinetic_energies.append(kinetic_energy / n_particles)
    potential_energies.append(potential_energy / n_particles)
    total_energies.append(total_energy / n_particles)
    temperatures.append(temperature)
    compressibilities.append(pressure / (temperature * density))

../_images/668cf4fc9f094fed98a9729fc162e2653557998871c740b956a62533754cafe2.png

../_images/8df4482b5bd0007deedc7301865c1417d0d816845dd379b4cf133daa81501167.png

../_images/d4dd2e3fea5b041d93e2ad2a1d4fd367f1d01d7a537c3cfe5e8c7d410c7f10a7.png

Fixing the temperature (canonical ensemble)#

During the equilibration stage, the kinetic energy of the system shows considerable drift. As a result, the equilibrium temperature is not the same as the initial temperature.

To keep the temperature fixed we can renormalize the velocities periodically during the equilibration stage. This effectively models the canonical ensemble. Once equilibration is complete, we can switch to the microcanonical ensemble by removing the velocity renormalization.

# Simulation parameters
n_particles = 64      # Number of particles
density     = 0.1      # Density of the fluid
box_length = (n_particles/density)**(1/3.)  # Length of simulation box (determined by the number of particles and density)
temperature0 = 1.0    # Initial (desired) temperature of the system
time_step = 0.01      # Time step for the simulation
num_steps = 1000       # Number of simulation steps
equilibration_steps = 500 # Number of equilibration steps

# Initial conditions
positions = initial_positions()
velocities = np.random.normal(loc=0.0, scale=np.sqrt(temperature0), size=(n_particles, 3))
accelerations = np.zeros((n_particles, 3))
t = 0

renormalize_velocities(velocities, temperature0)
accelerations, potential_energy, pressure = compute_forces(accelerations, positions, lj_potential, lj_potential_gradient)
kinetic_energy = 0.5*np.sum(velocities**2)
total_energy = potential_energy + kinetic_energy
temperature = temperature0

times = []
kinetic_energies = []
potential_energies = []
total_energies = []
temperatures = []
compressibilities = []

times.append(t)
kinetic_energies.append(kinetic_energy / n_particles)
potential_energies.append(potential_energy / n_particles)
total_energies.append(total_energy / n_particles)
temperatures.append(temperature)
compressibilities.append(pressure / (temperature * density))


# Main loop
for step in range(num_steps):
    # Perform Velocity Verlet integration
    positions, velocities, accelerations, potential_energy, pressure = velocity_verlet(positions, velocities, accelerations, time_step, lj_potential, lj_potential_gradient)
    # Compute kinetic energy
    kinetic_energy = 0.5*np.sum(velocities**2)
    # Compute total energy
    total_energy = potential_energy + kinetic_energy
    # Compute temperature and pressure
    temperature = 2.0*kinetic_energy/(3.0*n_particles)
    # Update time
    t += time_step
    
    times.append(t)
    kinetic_energies.append(kinetic_energy / n_particles)
    potential_energies.append(potential_energy / n_particles)
    total_energies.append(total_energy / n_particles)
    temperatures.append(temperature)
    compressibilities.append(pressure / (temperature * density))
    
    # Renormalize velocities
    if step < equilibration_steps and step % 10 == 0:
        renormalize_velocities(velocities, temperature0)
        kinetic_energy = 0.5*np.sum(velocities**2)
        total_energy = potential_energy + kinetic_energy
        temperature = 2.0*kinetic_energy/(3.0*n_particles)