Introduction

Methane in nitrogen is one of the most economically important gas mixtures on Earth. Natural gas wells routinely produce CH₄ contaminated with N₂ at concentrations that push the fuel below pipeline heating-value specifications; upstream processors spend billions each year on nitrogen rejection units, and the design of those units depends on how the two gases mix, separate, and diffuse under high pressure. On the cryogenic side, the boiling points of CH₄ (111 K) and N₂ (77 K) differ by only 34 K, so the same mixture is the working system for LNG production and for many air-separation cascades. Understanding how these two species interpenetrate at the molecular scale — how fast, in which regime, and under what thermal driving force — is a materials-science question dressed up as an engineering one.

Molecular dynamics is uniquely suited to answering it. Below the continuum limit (roughly < 100 nm channels or < 10 ns timescales) Fick's law starts to fray, and any predictive statement about diffusion coefficients, thermal separation, or interfacial thickness must come from an atomistic simulation. This post walks through the physics of CH₄–N₂ mixing and shows how a compact Lennard-Jones LAMMPS setup can capture density profiles, temperature gradients, and transport signatures across a non-isothermal binary interface.

Self-diffusion ≠ mutual diffusion – The rate at which a labelled CH₄ molecule wanders through a bath of identical CH₄ molecules (self-diffusion, D_s) is a different quantity from the rate at which CH₄ smears into a body of N₂ (mutual, or interdiffusion, D_AB). Only the latter has an interface. Mixing them up is one of the classic beginner MD errors.

1 Fundamental Physics

The Lennard-Jones potential itself is deceptively simple:

$$U_\text{LJ}(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$$

For methane, the united-atom parameters commonly used are ε/k_B ≈ 148 K and σ ≈ 3.73 Å. For diatomic nitrogen (two-site model), ε/k_B ≈ 36 K per site, σ ≈ 3.31 Å with a bond length of 1.10 Å. The 12–6 form is not fundamental — it is a computationally convenient approximation to the true intermolecular potential — but for these two small, weakly-polar molecules it reproduces vapour–liquid equilibria, viscosity, and diffusion coefficients to within a few percent of experiment.

2 Why Modeling Diffusion Is Hard

1. Timescale gap – Real diffusion coefficients are ~10⁻⁵ cm²/s. Traversing a 10 nm interface takes ~1 ns in MD, ~100 μs in experiment. The atomistic run is fine; extrapolating to macroscopic timescales requires care.
2. Finite-size effects – Yeh–Hummer corrections show that D_MD systematically under-predicts D_∞ by O(k_BT/(6πηL)). Box sides < 5 nm can be 10–20% off.
3. Cutoff artefacts – Truncating the LJ potential at r_c without a tail correction shifts pressure, phase boundaries, and interfacial tension. Use pair_modify tail yes or explicit long-range dispersion when accuracy matters.
4. Thermostat coupling – Applying a global Berendsen or Langevin thermostat to a mixture with two nominal temperatures homogenises them within tens of picoseconds. If you want a steady thermal gradient, thermostat each species (or each spatial region) separately.
5. Interface definition – The instantaneous interface between two miscible gases fluctuates at the ~1 nm scale. Averaging density over 100+ ps and 20+ chunks along the mixing axis is standard.

3 MD Workflow in LAMMPS

3.1 Choose an Appropriate Force Field

For pure mixing at ambient conditions the united-atom LJ (or the two-site N₂ model) is the correct choice. All-atom + Coulomb (lj/cut/coul/long with PPPM electrostatics) becomes necessary only when explicit hydrogens or polarisation effects matter — for most transport-property studies the united-atom flavour is faster and equally accurate.

3.2 Build the Initial Configuration

1. Two-region initial condition – place CH₄ molecules in one slab (say 0 < x < 100 Å) and N₂ in the adjacent slab (100 < x < 200 Å). LAMMPS region ... block + create_atoms + delete_atoms porosity is the most compact way; packmol is friendlier for larger systems.
2. Charge neutrality – if using lj/cut/coul/long, verify each molecule's partial charges sum to zero.
3. Minimise before dynamics – always: 200 steps of conjugate-gradient prevents blow-ups from overlapping placements.

3.3 Assign Ensembles and Boundary Conditions

• Equilibration – fix press/berendsen iso 20.0 20.0 1000.0 at target pressure, plus per-species temp/berendsen at the target temperatures. Do NOT couple pressure and temperature with a single Nose–Hoover chain at this stage; species-specific thermostats let you hold hot CH₄ and cold N₂ simultaneously without artificially equilibrating them.
• Production – switch to fix nvt temp <T> <T> 100 per species. Berendsen is a rescaling thermostat and generates wrong canonical distributions; NVT (Nose–Hoover) is correct for equilibrium sampling.
• Periodic boundaries – full 3D periodicity is standard. If the interface is normal to x, ensure the box is long enough in x that the interfacial region and both bulk reservoirs are well-separated (typically ≥ 5σ each).

4 Sample Applications

5 See It In Motion

The animation below shows a slab of methane (hot, 400 K) meeting a slab of nitrogen (cold, 300 K) at 20 bar. Watch the interface widen as mutual diffusion carries CH₄ molecules into the N₂ bath, while species-specific thermostats hold each reservoir at its target temperature so a persistent thermal gradient develops across the mixing zone.

The subtle asymmetry in the mixing front — CH₄ advances slightly faster into the cold N₂ than the reverse — is the fingerprint of the Soret effect predicted by the Onsager coupling in Section 1.

6 Common Pitfalls & Remedies

1. Both species drift to the same temperature within picoseconds – you accidentally used a global thermostat. Split into per-group Nose–Hoover chains (fix nvt temp per group).
2. Pressure ratchets up during equilibration – Berendsen barostat time constant too small; increase to ≥ 1000 fs and lengthen the equilibration run.
3. Density profile looks like a straight line right from step 1 – the initial condition is already mixed. Check that your CH₄ and N₂ groups were placed in disjoint x-regions and that delete_atoms overlap didn't blend them.
4. PPPM warnings about slab geometry – if you have vacuum in z, add kspace_modify slab 3.0 or switch to 2D-Ewald.
5. D_AB disagrees with experiment by 30% – check the Yeh–Hummer finite-size correction and box size, then verify your force-field parameters (TraPPE-UA methane ≠ OPLS-AA methane).

Stay Connected

Questions? Email support@simuxai.com or join our community forum.

Happy Simulating!