Relaxing an MD Structure: Three Techniques Before You Hit Run

Before you run a single step of production MD, you need a reasonable starting structure. Skip this and your simulation becomes unstable, takes forever to equilibrate, or simply refuses to start. This post walks through three techniques for relaxing an initial structure — gradient descent, NVE-limit, and a soft-core pushoff — and shows each one running live in Jungular.

Download the files

Why relaxation matters

The first step of any MD simulation is building the geometry and topology. It's worth setting yourself up for success here: give the engine a sensible initial structure and relax the geometry properly. If you don't, the consequences are unstable dynamics, painfully slow equilibration, or a run that dies on step one.

Here's a classic bad starting point — two molecules placed too close together, with overlapping atoms. Hit run on this and the engine loses atoms immediately:

All three methods below use LAMMPS as the underlying engine, but the same ideas carry over to most MD software.

Method 1: Gradient descent (energy minimization)

If you've taken a calculus class, this one clicks immediately. We know the total potential energy of the system — the sum of the Lennard-Jones potential for non-bonded atoms plus bond, angle, dihedral, and (sometimes) improper terms:

pe = lj_potential + bond_energy + angle_energy + dihedral + improper

The force on each atom is the negative gradient of that total potential energy. From there we can pick a minimization style. The most straightforward is steepest descent — just follow the downhill direction given by the force vector. LAMMPS also offers conjugate gradient (its default), which is more involved but converges faster. If you want the details, the LAMMPS documentation covers both.

In the demo, we nudge some atoms far from their equilibrium positions, then click Relax → Minimize. You can inspect the exact LAMMPS command it runs, hit go, and watch the atoms snap back to their equilibrium configuration.

Gradient descent is perfect for a structure that's already reasonable but has some strain to work out. It won't, however, rescue atoms that are sitting on top of each other — for that, we need dynamics.

Method 2: NVE-limit

The second method is NVE-limit. As the name suggests, it's very similar to running dynamics in the NVE ensemble — total energy is held constant — with one crucial addition: we cap how far any atom can move in a single timestep. Set that limit to, say, 0.1 Å, and no atom moves more than 0.1 Å per step.

That cap is exactly what makes it useful for bad geometries. If atoms are overlapping all over the place, the total energy (and the forces) are enormous. But because NVE-limit bounds each atom's displacement, the system doesn't explode — it stabilizes over time.

Here's a genuinely bad structure: a polymer packed into the box 20 times with overlaps allowed. The red regions are where atoms clash.

Relax it under NVE-limit — the command is just fix relax all nve/limit 0.1, where 0.1 is the per-step displacement cap — and the red overlap regions disappear. After a few moments you get a clean configuration with all the overlaps resolved, and the system runs fine under NVT afterward (bring the temperature down to 300 K and you're set).

One important caveat: NVE-limit can stabilize a system over time, but the forces still have to be computable. Two atoms can't be at the exact same position — if r = 0, you get a divide-by-zero, and no amount of NVE-limit can fix that.

Method 3: Soft push

NVE-limit handles maybe 95% of cases. But every so often you hit something like this: two molecules stuck in a weird, glued-together configuration. They look like one blob, but they're actually two separate molecules trapped in a local energy minimum — the potential energy is absurdly high for a two-molecule system, and NVE-limit can't get them out.

To see why overlaps are so dangerous, look at the Lennard-Jones potential:

V = 4·epsilon·[ (sigma/r)^12 − (sigma/r)^6 ]

When r is very small, the (sigma/r)^12 term dominates — it's an enormous number — so V blows up and the whole system's energy goes with it.

Soft push sidesteps this by replacing the Lennard-Jones pair potential with a soft potential. Here's the LAMMPS definition:

The soft potential is E = A·[1 + cos(π·r/rc)], with three parameters: A is a prefactor, r is the interatomic distance, and rc is a cutoff. The key property: when r = 0, E is not infinity. cos(0) = 1, so E = 2A — a perfectly normal number. That's why it's called soft: even when two atoms overlap, the energy doesn't blow up.

Better still, A can be varied over time. In practice we set A = 0 at t = 0, so overlapping atoms literally can't feel each other, then ramp it up so the interactions gradually harden and the atoms separate gently. It only replaces the non-bonded pair term — bond, angle, and dihedral terms stay intact, so molecules keep their shape.

Back in the simulation, click Relax → Unstick on those two stuck molecules and they separate cleanly, escaping the unphysical state. Run dynamics afterward and everything looks normal — the total energy is sensible and the temperature settles around the 300 K we set for the NVT run.

Try it yourself

Three tools, in order of escalation: gradient descent for a structure that just needs strain relieved, NVE-limit for overlapping atoms (your workhorse), and soft push for the stubborn cases stuck in a local minimum. Grab the template LAMMPS scripts below and try them on your own systems.

Get the relaxation scripts