Phase Separation in a Binary Lennard-Jones Mixture

Why does a mixture demix?¶

Two species interacting with identical $\varepsilon$ behave as one — they stay mixed. But if unlike pairs feel a weaker attraction than like pairs,

\varepsilon_{AB} < \frac{1}{2}(\varepsilon_{AA} + \varepsilon_{BB}),

(1)

the system can lower its free energy by separating. Entropy of mixing fights this. Whether the mixture is homogeneous or phase-separated is set by the competition, tuned by temperature and by the interaction asymmetry.

We use 2D MD with velocity-Verlet and periodic boundary conditions. A Berendsen-style thermostat keeps the temperature controlled without ruining the structural properties we care about.

import numpy as np
import matplotlib.pyplot as plt
from numba import njit

rng = np.random.default_rng(0)

A binary LJ force routine¶

Each particle carries a type $t_i \in \{0, 1\}$ (A or B). The pair potential is cut and shifted at $r_c = 2.5\,\sigma$ .

@njit
def forces_energy(pos, types, box, eps_AA, eps_BB, eps_AB, sigma=1.0, rc=2.5):
    N = pos.shape[0]
    F = np.zeros_like(pos)
    U = 0.0
    rc2 = rc * rc
    sig2 = sigma * sigma
    for i in range(N - 1):
        for j in range(i + 1, N):
            dx = pos[i, 0] - pos[j, 0]
            dy = pos[i, 1] - pos[j, 1]
            dx -= box * round(dx / box)
            dy -= box * round(dy / box)
            r2 = dx * dx + dy * dy
            if r2 < rc2 and r2 > 1e-12:
                if types[i] == 0 and types[j] == 0:
                    eps = eps_AA
                elif types[i] == 1 and types[j] == 1:
                    eps = eps_BB
                else:
                    eps = eps_AB
                inv2 = sig2 / r2
                inv6 = inv2 * inv2 * inv2
                inv12 = inv6 * inv6
                U += 4.0 * eps * (inv12 - inv6)
                fmag = 24.0 * eps * (2.0 * inv12 - inv6) / r2
                F[i, 0] += fmag * dx; F[i, 1] += fmag * dy
                F[j, 0] -= fmag * dx; F[j, 1] -= fmag * dy
    return F, U


@njit
def run_md(pos, vel, types, box, T_target=0.8, eps_AA=1.0, eps_BB=1.0,
           eps_AB=0.3, dt=0.005, n_steps=20_000, thermo=0.1):
    N = pos.shape[0]
    F, _ = forces_energy(pos, types, box, eps_AA, eps_BB, eps_AB)
    for step in range(n_steps):
        vel += 0.5 * dt * F
        pos += dt * vel
        pos -= box * np.floor(pos / box)
        F, _ = forces_energy(pos, types, box, eps_AA, eps_BB, eps_AB)
        vel += 0.5 * dt * F
        # Berendsen-style rescale
        KE = 0.5 * np.sum(vel * vel)
        Tnow = KE / N
        lam = np.sqrt(1.0 + thermo * (T_target / Tnow - 1.0))
        vel *= lam
    return pos, vel

Initialize and quench a 50:50 mixture¶

Start from a random mixture on a lattice, give it thermal velocities, and run a short trajectory at $T = 0.8$ with weak A–B attraction $\varepsilon_{AB} = 0.3$ . Like-type particles cluster into domains.

def lattice_init(n_side, box, n_A_fraction=0.5, seed=0):
    r = np.random.default_rng(seed)
    N = n_side * n_side
    xs = np.linspace(0, box, n_side, endpoint=False)
    X, Y = np.meshgrid(xs, xs)
    pos = np.column_stack([X.ravel(), Y.ravel()]) + 0.1 * r.standard_normal((N, 2))
    pos = np.mod(pos, box)
    types = (r.random(N) > n_A_fraction).astype(np.int64)
    vel = r.standard_normal((N, 2)) * np.sqrt(0.8)
    vel -= vel.mean(axis=0)
    return pos, vel, types

n_side, box = 20, 22.0
pos, vel, types = lattice_init(n_side, box)
pos, vel = run_md(pos, vel, types, box, T_target=0.8,
                  eps_AA=1.0, eps_BB=1.0, eps_AB=0.3, n_steps=20_000)

fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(pos[types == 0, 0], pos[types == 0, 1], s=16, color='tab:red', label='A')
ax.scatter(pos[types == 1, 0], pos[types == 1, 1], s=16, color='tab:blue', label='B')
ax.set(xlim=(0, box), ylim=(0, box), aspect='equal',
       title=r'50:50 binary LJ at $T=0.8$, $\varepsilon_{AB}=0.3$')
ax.legend(loc='upper right')
fig.tight_layout()

Problems¶

Phase diagram in $(T, \varepsilon_{AB})$ . Run the simulation on a grid of $T \in \{0.4, 0.6, 0.8, 1.0, 1.5\}$ and $\varepsilon_{AB} \in \{0.1, 0.3, 0.5, 0.8, 1.0\}$ at fixed $\varepsilon_{AA} = \varepsilon_{BB} = 1$ . Classify each run visually as mixed vs. demixed and plot a phase diagram. Where is the critical curve?
Composition dependence. At $T = 0.8$ , $\varepsilon_{AB} = 0.3$ , vary the A:B ratio $\in \{0.2, 0.3, 0.5, 0.7, 0.8\}$ . Does morphology shift from bicontinuous (near 50:50) to isolated droplets (off-critical composition)? Report with snapshots.
Domain size growth. Fix $T = 0.6$ , $\varepsilon_{AB} = 0.2$ , 50:50, and save snapshots at $t \in \{10^3, 5\!\times\!10^3, 2\!\times\!10^4, 5\!\times\!10^4\}$ steps. Estimate the characteristic domain size $L(t)$ (e.g. via the inverse of the first peak position of the A–A structure factor, or by counting same-type neighbors). Is $L(t) \sim t^{1/3}$ as expected for diffusive coarsening?

Reference¶

J. W. Cahn, On spinodal decomposition, Acta Metall. 9, 795 (1961). Cahn (1961)