Lab: beyond, by hand

The setup

Three independent parts, each with its own tiny system:

Everything is checkable: the hand numbers of Lab 1, the exact Morse curve, and the exact barrier of 5.

Step 1

Lab 1 in LAMMPS: the same numbers

Why: A real MD code does exactly what we did by hand, only for millions of atoms. Before trusting it with a research system, it is worth checking that it reproduces a calculation we understand completely. This also shows how a code reports its results.

What the numbers say: All six digits agree, because LAMMPS fix nve is the velocity Verlet we did by hand in Lab 1. The one trap: in units lj LAMMPS divides energies by the number of atoms in its thermo output by default, which gives exactly half of our numbers for 2 atoms. Use thermo_modify norm no when comparing totals. Checking a code on a tiny system with a known answer is the cheapest bug test there is, and it is the first step in the Research section on modifying LAMMPS.

Step 2

Units: the same run in real and metal units

Why: Research inputs use physical units: real (kcal/mol, Å, fs; for biomolecules and liquids) or metal (eV, Å, ps; for materials and most ML potentials). A unit slip in ε or in the time step silently gives a different simulation, so let us convert once by hand and check that LAMMPS agrees.

What the numbers say: Three unit systems, one trajectory: dividing by ε gives back the reduced numbers exactly. The classic mistakes are all visible here: the time step is in fs for real but in ps for metal (writing 1 instead of 0.001 gives a 1000 times larger step), energies are in kcal/mol vs eV (a factor of 23.06), and k_B differs between them. Always check the units line of an input before reading any number from it.

Step 3

A machine-learned potential in miniature: least squares by hand

Why: An ML potential is a flexible function fitted to energies (and forces) from quantum chemistry. The simplest version is a linear model: energy = a sum of fixed functions of the geometry, each with an unknown coefficient. The coefficients come from least squares, a 2 × 2 system here.

What the numbers say: The whole "training" is two sums of squares and one 2 × 2 solve; real ML potentials do the same with thousands of coefficients and millions of energies and forces (MACE, NequIP and ACE are built on this idea, with smarter functions of the geometry). But this model has only the Lennard-Jones shape: with two coefficients it cannot follow the Morse curve, and it misses the third point by 0.35. A model that is too rigid cannot be fixed by fitting.

Step 4

Interpolation, extrapolation and active learning

Why: A fitted potential is only tested where it was trained. During a simulation, atoms visit geometries outside the training set, for example two atoms pushed closer than any training point. What happens there?

What the numbers say: Outside the training data the error explodes: at r = 0.90 the fit says +4.2 where the truth is −0.69, a fake wall or, with another model, a fake hole the atoms fall into. Adding the geometries where the model is worst ("active learning": run, find uncertain configurations, compute them with DFT, refit) reduces the worst error six-fold. The remaining error is the rigid model again: improving one region worsens another. Limitation for ML potentials: they are reliable only inside the space covered by their training set, so the training data must include the high-temperature, compressed and reactive geometries the simulation will visit.

Step 5

Flexible models need data

Why: The cure for a rigid model is more flexibility. But a flexible model can pass through every training point and still be wrong in between. We test both effects with two flexible models.

What the numbers say: Zero training error means nothing by itself: the parabola fits all three points exactly and then rises to +6.7 at r = 3, where the true energy goes to 0. With only 3 points, even 8 basis functions leave a test error of 0.26; with 20 points the same model is accurate to 0.004 everywhere in the range. The lesson behind every ML potential: always judge a model on test data it has not seen, the amount and spread of training data matters more than the model, and physics built into the model (the energy must go to zero at large r, smoothness, symmetry) is what makes it safe outside the data.

Step 6

Metadynamics: three hills by hand

Why: A barrier of 5 k_B T is crossed only rarely, so a plain simulation spends almost all its time in one well. Metadynamics adds small Gaussian hills at the places the system visits, along a chosen collective variable s. The hills fill the current well until the system has to leave it, and their sum, the bias, records the free energy.

What the numbers say: Three hills already lower the effective barrier from 5.0 to 3.53, which makes escape about 4 times more likely (0.029 vs 0.0067). Each new hill is added where the system is, so the well fills from the bottom up; once it is full, the walker spills over and starts filling the other well. At the end, the bias is a mirror image of the free energy: F(s) ≈ −V(s) + constant. The hills are only as good as the variable s: if s misses a slow motion (a "hidden" variable), filling along s does not help, and choosing the collective variable is the real difficulty of metadynamics.

Step 7

A long run: standard vs well-tempered

Why: Does the bias really converge to the free energy, and how fast? We run a Langevin walker in the double well for 10⁶ steps, adding a hill every 200 steps, and read the barrier from the bias as the run goes on. Well-tempered metadynamics makes each new hill smaller where the bias is already high, w × exp(−V/ΔT).

What the numbers say: Both runs cross the barrier 10–16 times more often than the unbiased walker (5 crossings), which is the whole point. But standard metadynamics never settles: with full-height hills forever, the bias keeps overshooting, and the barrier estimate swings between 1.8 and 7.6. Well-tempered hills shrink as the bias grows, so the estimate settles at 5.00 ± 0.3. Limitations: hill height, width, stride and ΔT must be chosen (too big gives a rough, overfilled surface; too small converges slowly); the dynamics are biased, so rates need extra theory; and everything depends on a good collective variable. In practice this is done with PLUMED, coupled to LAMMPS or GROMACS.

Summary