MD/MC simulations of fluids¶

Molecular Dynamics Simulation of a Lennard-Jones Fluid¶

The code in this notebook simulates a collection of particles interacting through the Lennard-Jones (LJ) potential using Molecular Dynamics (MD) in 3D.
The simulation is designed to be efficient, modular, and easy to understand, leveraging NumPy and object-oriented programming (OOP).
Particles are initially arranged in a face-centered cubic (FCC) lattice and evolved over time using the Velocity Verlet algorithm under periodic boundary conditions.
The simulation tracks: Kinetic Energy (KE), Potential Energy (PE), Pressure, Pair correlation function $g(r)$

Physical Model of LJ fluid¶

Each particle experiences a force given by the Lennard-Jones potential:

V(r) = 4 \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)

(1)

$r$ is the distance between two particles.
The potential is cut off at a radius (r_c) to improve efficiency.

The system is simulated at a fixed temperature by rescaling velocities during the early steps (velocity rescaling).

PBC¶

If a particle leaves the box, it re-enters at the opposite side.
What should be the distance between partiles? There is some mbiguity because distance between atoms may now include crossing the system boundary and re-entering. This is identical to introducing copies of the particles around the simulation box. We adopt minimum image convention by choosing the shortest distance possible

figure-md - Unknown Directive

<img src="./figs/pbc1.png" class="bg-primary" width="400px">  

Perioidc Boundary Conditions

figure-md - Unknown Directive

<img src="./figs/pbc2.png" class="bg-primary" width="400px">  

Minimum Image Convention showing that we only track particles with clsoest distance.

MD of LJ fluid¶

Velocity Verlet Algorithm

1. Evaluate the initial force from the current position:

F_t = -\frac{\partial U(r)}{\partial r} \Bigg|_{r(t)}

(2)

2. Update the position:

r_{t+\Delta t} = r_t + v_t \Delta t + \frac{F_t}{2m} \Delta t^2

(3)

3. Partially update the velocity:

v_{t+\Delta t/2} = v_t + \frac{F_t}{2m} \Delta t

(4)

4. Evaluate the force at the new position:

F_{t+\Delta t} = -\frac{\partial U(r)}{\partial r} \Bigg|_{r(t+\Delta t)}

(5)

5. Complete the velocity update:

v_{t+\Delta t} = v_{t+\Delta t/2} + \frac{F_{t+\Delta t}}{2m} \Delta t

(6)

import numpy as np
import matplotlib.pyplot as plt

class LJSystem:
    def __init__(self, rho=0.88, N_cell=3, T=1.0, rcut=2.5):
        # System parameters
        self.N = 4 * N_cell**3
        self.rho = rho
        self.L = (self.N / rho)**(1/3)
        self.T_target = T
        self.rcut = rcut
        self.rcut_sq = rcut**2

        # Create FCC lattice positions
        base = np.array([[0, 0, 0], [0.5, 0.5, 0], [0.5, 0, 0.5], [0, 0.5, 0.5]])
        pos = np.array([b + [i, j, k] for i in range(N_cell) for j in range(N_cell) for k in range(N_cell) for b in base])
        self.pos = pos * (self.L / N_cell)

        # Velocities
        self.vel = np.random.randn(self.N, 3)
        self.vel -= np.mean(self.vel, axis=0)

        # Indices for upper triangular part
        self.I, self.J = np.triu_indices(self.N, k=1)

    def minimum_image(self, r_vec):
        return (r_vec + self.L/2) % self.L - self.L/2

    def compute_forces(self):
        r_vec = self.minimum_image(self.pos[self.I] - self.pos[self.J])
        r_sq = np.sum(r_vec**2, axis=1)

        mask = r_sq < self.rcut_sq
        r_vec = r_vec[mask]
        r_sq = r_sq[mask]

        r2_inv = 1.0 / r_sq
        r6_inv = r2_inv**3

        dU_dr = 48 * r6_inv * (r6_inv - 0.5) * r2_inv
        F_ij = dU_dr[:, None] * r_vec

        # Forces on each particle
        forces = np.zeros_like(self.pos)
        np.add.at(forces, self.I[mask],  F_ij)
        np.add.at(forces, self.J[mask], -F_ij)

        # Potential energy
        U = 4 * np.sum(r6_inv * (r6_inv - 1))

        # Virial for pressure
        Pvir = np.sum(np.sum(F_ij * r_vec, axis=1))

        # Histogram for g(r)
        hist, _ = np.histogram(np.sqrt(r_sq), bins=30, range=(0, self.L/2))

        return forces, U, Pvir, hist

    def integrate(self, dt):
        forces, U, Pvir, hist = self.compute_forces()
        self.vel += 0.5 * dt * forces
        self.pos = (self.pos + dt * self.vel) % self.L
        forces_new, U_new, Pvir_new, hist_new = self.compute_forces()
        self.vel += 0.5 * dt * forces_new
        kin = 0.5 * np.sum(self.vel**2)
        return kin, U_new, Pvir_new, hist_new

    def rescale_velocities(self, kin):
        """Velocity rescaling to maintain target temperature"""
        factor = np.sqrt(3 * self.N * self.T_target / (2 * kin))
        self.vel *= factor

    def simulate(self, dt=0.005, nsteps=5000, freq_out=50):
        kins, pots, Ps, hists = [], [], [], []
        for step in range(nsteps):
            kin, pot, Pvir, hist = self.integrate(dt)

            if step < 1000:  # Temperature stabilization
                self.rescale_velocities(kin)

            if step % freq_out == 0:
                kins.append(kin)
                pots.append(pot)
                Ps.append(Pvir)
                hists.append(hist)

        return np.array(kins), np.array(pots), np.array(Ps), np.mean(hists, axis=0)

Example simulation and analysis¶

# Initialize system
lj = LJSystem(rho=0.88, N_cell=3, T=1.0)

# Run simulation
kins, pots, Ps, hist = lj.simulate(dt=0.005, nsteps=10000, freq_out=50)

# Post-processing
r = np.linspace(0, lj.L/2, 30)
dr = r[1] - r[0]
shell_volumes = 4*np.pi*r**2*dr
g_r = hist / (lj.N * lj.rho * shell_volumes)

# Plot
plt.plot(r, g_r, '-o')
plt.xlabel(r'$r$')
plt.ylabel(r'$g(r)$')
plt.show()

/tmp/ipykernel_3028/3134754071.py:11: RuntimeWarning: invalid value encountered in divide
  g_r = hist / (lj.N * lj.rho * shell_volumes)

MC simulation¶

import numpy as np
import matplotlib.pyplot as plt

class LJ_MC_System:
    def __init__(self, rho=0.88, N_cell=3, T=1.0, rcut=2.5):
        # System parameters
        self.N = 4 * N_cell**3
        self.rho = rho
        self.L = (self.N / rho)**(1/3)
        self.T = T
        self.beta = 1.0 / T
        self.rcut = rcut
        self.rcut_sq = rcut**2

        # Create FCC lattice positions
        base = np.array([[0, 0, 0], [0.5, 0.5, 0], [0.5, 0, 0.5], [0, 0.5, 0.5]])
        pos = np.array([b + [i, j, k] for i in range(N_cell) for j in range(N_cell) for k in range(N_cell) for b in base])
        self.pos = pos * (self.L / N_cell)

        # Pair indices
        self.I, self.J = np.triu_indices(self.N, k=1)

    def minimum_image(self, r_vec):
        return (r_vec + self.L/2) % self.L - self.L/2

    def total_energy(self):
        r_vec = self.minimum_image(self.pos[self.I] - self.pos[self.J])
        r_sq = np.sum(r_vec**2, axis=1)
        mask = r_sq < self.rcut_sq
        r2_inv = 1.0 / r_sq[mask]
        r6_inv = r2_inv**3
        return 4 * np.sum(r6_inv * (r6_inv - 1))

    def particle_energy(self, idx):
        rij = self.minimum_image(self.pos[idx] - self.pos)
        r_sq = np.sum(rij**2, axis=1)
        mask = (r_sq < self.rcut_sq) & (r_sq > 0)  # exclude self-interaction
        r2_inv = 1.0 / r_sq[mask]
        r6_inv = r2_inv**3
        return 4 * np.sum(r6_inv * (r6_inv - 1))

    def mc_move(self, max_disp=0.1):
        idx = np.random.randint(self.N)
        old_pos = self.pos[idx].copy()
        old_energy = self.particle_energy(idx)

        # Propose move
        displacement = (np.random.rand(3) - 0.5) * 2 * max_disp
        self.pos[idx] = (self.pos[idx] + displacement) % self.L

        new_energy = self.particle_energy(idx)
        dE = new_energy - old_energy

        # Metropolis criterion
        if np.random.rand() > np.exp(-self.beta * dE):
            # Reject move
            self.pos[idx] = old_pos

    def sample(self):
        """ Sample pair distances for g(r) """
        r_vec = self.minimum_image(self.pos[self.I] - self.pos[self.J])
        r_sq = np.sum(r_vec**2, axis=1)
        hist, _ = np.histogram(np.sqrt(r_sq), bins=30, range=(0, self.L/2))
        return hist

    def simulate(self, nsteps=50000, freq_out=500):
        hists = []
        energies = []

        for step in range(nsteps):
            self.mc_move()

            if step % freq_out == 0:
                hist = self.sample()
                hists.append(hist)
                energies.append(self.total_energy())

        return np.array(energies), np.mean(hists, axis=0)

Example MC simulation and analysis¶

# ----------------- Main -----------------

# Initialize system
lj_mc = LJ_MC_System(rho=0.88, N_cell=3, T=1.0)

# Run simulation
energies, hist = lj_mc.simulate(nsteps=100000, freq_out=500)

# Post-processing
r = np.linspace(0, lj_mc.L/2, 30)
dr = r[1] - r[0]
shell_volumes = 4*np.pi*r**2*dr
g_r = hist / (lj_mc.N * lj_mc.rho * shell_volumes)

# Plot
plt.plot(r, g_r, '-o')
plt.xlabel(r'$r$')
plt.ylabel(r'$g(r)$')
plt.show()

/tmp/ipykernel_3028/3509366135.py:13: RuntimeWarning: invalid value encountered in divide
  g_r = hist / (lj_mc.N * lj_mc.rho * shell_volumes)

Differences of Monte Carlo vs Molecular Dynamics¶

Monte Carlo (MC)	Molecular Dynamics (MD)
No real dynamics, just configurations	Real-time evolution of particle trajectories
Samples equilibrium ensemble directly	Samples via time evolution
No conservation laws	Energy/momentum conserved after equilibration
Easier to implement	More complex (forces, integration needed)
Faster for static properties like $g(r)$	Necessary for dynamics (diffusion, viscosity)

Problems¶

MD sim of LJ fluid¶

Modify the cut-off radius $r_c$ and observe its effect on pressure and energy.
Remove velocity rescaling after equilibration and check energy conservation.
Implement a thermostat (e.g., Andersen, Berendsen) instead of manual rescaling.
Compute the diffusion coefficient from particle trajectories.
Visualize the particle trajectories over time (animation!).

MC sim of LJ fluid¶

Vary sim parameters

Tune maximum displacement (max_disp) to optimize acceptance rate (~40%-50% is ideal).
Check how $g(r)$ changes when density or temperature changes.
Implement energy cutoffs with shifted potentials to smooth out the energy at $r_c$ .
Try simulating hard spheres (infinite repulsion at contact).

Phases of LJ fluid¶

Run MD or MC simulations of LJ fluid at several temperatures to identify critical temperature.

At each temperature evaluate heat capacity and RDFs.
Plot how energy, heat capacity and RDF change as a function of temperature.
Study dependence on sampling efficiency on magnitude of particle displacement.