Structure optimization

A relaxation moves atoms downhill until the largest force left on any atom falls below a threshold you pick. The API follows ASE: build an Atoms, attach a calculator, hand it to an optimizer, call run(fmax=...).

Benzene relaxation

Benzene with every atom displaced by a seeded random vector of 0.18 Å, relaxed with BFGS and GFN2-xTB. Drag to orbit, and use the playback bar to step through the relaxation.

C6H6, 12 atoms, every atom displaced by a seeded random vector · 30 BFGS steps · one frame per step

import numpy as np

from atomli.build import molecule
from atomli.calculators.xtb import XTB
from atomli.optimize import BFGS

rng = np.random.default_rng(20260807)

atoms = molecule("C6H6")
atoms.set_positions(atoms.get_positions() + rng.normal(0.0, 0.18, (len(atoms), 3)))
atoms.calc = XTB(method="gfn2")

opt = BFGS(atoms, trajectory="benzene.traj", logfile="-")
opt.run(fmax=0.02)

That run took 30 steps and dropped the energy by 14.21 eV, from a starting fmax of 29.5 eV/Å to 0.0147 eV/Å.

Benzene is the test case because the answer is known in advance. The scrambled ring starts with C-C bonds spanning 0.580 Å and sitting 0.116 Å out of plane. After the relaxation all six C-C bonds agree to 0.17 mÅ at a mean of 1.3847 Å, and the ring is flat to 0.8 mÅ RMS. The experimental C-C distance in benzene is 1.39 Å. That agreement is the check worth running on your own system: relax something whose geometry you already know, and confirm you get it back.

Force tolerance

fmax is the largest force vector magnitude over atoms, in eV/Å, not the largest force component. The distinction is a factor of up to √3, so a component-wise check reports convergence that has not happened:

import numpy as np

fmax = np.linalg.norm(atoms.get_forces(), axis=1).max()  # what the optimizer tests
wrong = np.abs(atoms.get_forces()).max()                 # smaller, and not it

Choosing a value is a question about energy, not force, so here is the same benzene relaxation stopped at six different thresholds. The last column is how much energy the geometry is still carrying relative to the tightest run:

fmax (eV/Å)	Steps	C-C spread	Out of plane	Energy left
0.5	8	7.3 mÅ	46.7 mÅ	68.21 meV
0.2	12	2.9 mÅ	32.3 mÅ	13.60 meV
0.1	15	1.6 mÅ	25.7 mÅ	5.76 meV
0.05	20	1.1 mÅ	13.0 mÅ	1.81 meV
0.02	30	0.2 mÅ	0.8 mÅ	0.03 meV
0.01	32	0.1 mÅ	0.3 mÅ	0.00 meV

Read the ends against each other. Stopping at 0.5 costs 8 steps and leaves 68 meV on the table, about 2.7 times room-temperature kT, and the ring is still visibly buckled at 47 mÅ. Stopping at 0.02 costs 30 steps and leaves 0.026 meV, which is below the noise of any method you would run it with. The last factor of two after that is 2 more steps for nothing.

So: 0.05 eV/Å for a structure you are going to feed to something else, 0.01 to 0.02 eV/Å if the energy difference itself is the answer, and tighter than that only for a vibrational analysis, where residual forces contaminate the low frequencies. Do not pick a threshold from habit. Run the ladder once on a system like yours and pick the point where the energy stops moving.

Energy and force convergence

The single most misleading thing about a relaxation log is that fmax does not fall monotonically. It is a maximum over atoms, and it jumps every time a different atom becomes the worst one. In the benzene run above, fmax rose between consecutive steps 7 times out of 30, the first time at step 3, where it went from 4.18 to 7.02 eV/Å.

The energy did not. It fell on every single step of both runs on this page. That is the quantity to watch:

opt.attach(lambda: print(atoms.get_potential_energy()), interval=1)

A run where the energy also rises, or oscillates, is a real problem: the timestep or trust radius is too large, or the calculator’s forces are not the gradient of its energy. A run where only fmax bounces around is healthy.

Optimizer selection

The same two systems, the same starting geometry, the same threshold, all seven optimizers atomli ships:

Optimizer	Benzene, fmax 0.02	Ar adatom, fmax 0.005
BFGS	30 steps	112 steps
LBFGS	30 steps	111 steps
LBFGSLineSearch	27 steps	54 steps
FIRE	107 steps	281 steps
FIRE2	85 steps	147 steps
ABCFIRE	71 steps	111 steps
MDMin	88 steps	379 steps

Every one of those converged, and on benzene they agree on the final energy to within 0.52 meV. On that one molecule the step counts still differ by a factor of 4.0.

The ranking is not the same on both systems, which is the point: it is a property of the potential energy surface, not of the optimizer.

• BFGS builds an approximate inverse Hessian and is the default choice for small systems. Its memory is O(N²) in the number of atoms.
• LBFGS keeps only the last few updates instead of the full matrix. Same behaviour on a small system, tractable on a large one. Use it above a few hundred atoms.
• LBFGSLineSearch adds a line search along each direction. It won both systems here, at 27 and 54 steps, because each step is longer and better chosen. Each step costs more than one force evaluation, so compare wall time and not just step count when the calculator is expensive.
• FIRE, FIRE2, ABCFIRE are damped molecular dynamics with an adaptive timestep, not quasi-Newton methods. They build no Hessian, so they are hard to destabilise from a very bad starting geometry and they cost nothing per step in memory. They also take 5.2× more steps here. FIRE2 is the revised algorithm and ABCFIRE adds a bias correction to it.
• MDMin is the simplest of the lot, a velocity quench. It was the slowest on both systems. Reach for it when everything else diverges.

Start with BFGS, move to LBFGS when the system gets big, and switch to a FIRE variant when a quasi-Newton run blows up rather than converges.

Constraint validation

A surface calculation freezes the deep layers of the slab. Those layers stand in for bulk that is not in the cell, and letting them relax would let the whole slab drift and contaminate the surface energy.

Here is an argon adatom starting 3 Å above an fcc(100) slab, off every mirror plane, with the bottom two layers held by FixAtoms. Scrub it and watch the adatom slide sideways across the surface into the fourfold hollow while the bottom of the slab stays exactly where it started.

Ar fcc(100), 4 layers, 3x3 cells, 72 slab atoms + 1 adatom. 36 atoms fixed, 37 free, 112 BFGS steps.

from atomli.atoms import Atoms
from atomli.build import add_adsorbate, bulk
from atomli.calculators import LennardJones
from atomli.constraints import FixAtoms
from atomli.optimize import BFGS

slab = bulk("Ar", "fcc", a=5.26, cubic=True).repeat((3, 3, 2))
slab.set_pbc([True, True, False])
slab.center(vacuum=8.0, axis=[False, False, True])
add_adsorbate(slab, Atoms("Ar", positions=[[0, 0, 0]]), height=3.0, position=(2.4, 0.5))

z = slab.get_positions()[:, 2]
slab.set_constraint(FixAtoms(indices=[i for i, zi in enumerate(z) if zi < z.min() + 3.0]))
slab.calc = LennardJones(epsilon=0.0103, sigma=3.4, rc=5.0, shift=True)

BFGS(slab).run(fmax=0.005)

The constraint is attached to the atoms, not to the optimizer, so every optimizer respects it and so does molecular dynamics. Over the whole run the 36 fixed atoms moved by 0.0 Å, exactly zero, while the adatom travelled 2.06 Å across the surface. Check it rather than trusting it:

before = slab.get_positions()[fixed].copy()
BFGS(slab).run(fmax=0.005)
assert np.abs(slab.get_positions()[fixed] - before).max() == 0.0

get_forces() returns zero on a constrained atom, which is what makes fmax the right convergence measure here: it is a maximum over the degrees of freedom that are actually free. An unconstrained fmax would include the large forces holding the frozen layers in place and would never fall. See Constraints for the full set.

The periodic cell brings its own trap. This one is 15.78 Å across, so the Lennard-Jones cutoff is 5.0 Å. A cutoff past half the shortest cell vector makes an atom interact with two images of the same neighbour, and the forces stop being the gradient of the energy. Assert it, because nothing else will tell you:

assert calc_rc < min(np.linalg.norm(slab.get_cell()[:2], axis=1)) / 2

Minima and saddle points

An optimizer stops when the force vanishes. The force also vanishes at saddle points, and a symmetric starting geometry walks straight into one, because symmetry keeps the net force zero in the direction that would take it out.

The same adatom, started above four different points on the same surface. All four converged:

Start	Steps	Final fmax	Height	Energy above lowest
hollow	19	0.0037	2.96 Å	0.0 meV
bridge	28	0.0043	3.59 Å	13.1 meV
on top	35	0.0048	3.93 Å	20.2 meV
off-site	112	0.0048	3.03 Å	0.3 meV

The bridge and on-top runs finished with fmax below the threshold, in 28 and 35 steps, and are 13 and 20 meV above the hollow. They did not move laterally at all: their final positions are their starting positions to machine precision. They are diffusion barriers, not adsorption sites, and nothing in the optimizer’s output says so.

Two habits fix this:

1. Never start an adsorbate exactly on a symmetry point. The off-site run started at (2.4, 0.5) Å, on no mirror plane, and found the hollow on its own in 112 steps.
2. Rattle the converged structure and relax again. If it comes back to the same geometry it was a minimum; if it falls somewhere lower it was not. That is also exactly what the benzene demo does in reverse.

The off-site and hollow runs both land in the same site and still differ by 0.31 meV, because fmax 0.005 on a surface this flat leaves the adatom about 0.2 Å of slack. That is the fmax ladder again, in a system where it matters more.

Convergence troubleshooting

There are two distinct failures and they look nothing alike.

Out of steps. run() returns whether or not it converged, so the return value alone is not a check. Ask the optimizer. The same benzene relaxation given only 5 steps:

opt = BFGS(atoms)
opt.run(fmax=0.02, steps=5)

opt.converged()   # False
opt.nsteps        # 5

It stopped at fmax 1.90 eV/Å, a factor of 95 short, with 0.47 eV still to give up. Always test opt.converged() before using the geometry for anything.

Plateaued. The slab run needed 112 steps but was within a factor of two of the threshold by step 30. Across the 82 steps that followed, fmax rose and fell without ever clearing the bar, and the energy moved by a median of 0.12 meV per step for a total of 18.6 meV. The adatom was crossing an almost flat corrugation. A shallow direction on the energy surface looks exactly like a stall.

The bare slab with no adatom converges in 11 steps. The step count is a property of the softest mode in the system, not of the number of atoms.

To tell the two apart, look at the energy:

• Energy still falling, fmax flat: a soft mode. Give it more steps, or try LBFGSLineSearch, which took 54 instead of 112 on exactly this problem.
• Energy flat and fmax flat well above the threshold: genuinely stuck. The starting geometry has atoms on top of each other, the calculator’s forces disagree with its energy, or a constraint is fighting the relaxation.
• Energy rising: the step size is too large. Switch to a FIRE variant, which controls its own step length.

opt = BFGS(atoms)
opt.run(fmax=0.02, steps=200)
if not opt.converged():
    raise RuntimeError(f"stopped at {opt.nsteps} steps, fmax still too high")