Geometry Optimization

Geometry optimization searches for a stationary structure on the potential-energy surface defined by a particular calculator. Starting from an imported or constructed geometry, the optimizer repeatedly evaluates energy and forces, constructs a local model of the surface, chooses a bounded displacement, and updates that model until the residual force satisfies the declared threshold.

The word optimization can be misleading: the procedure does not improve a molecule in a general chemical sense and does not search every possible structure. It finds a nearby stationary point under one Hamiltonian or interatomic model, within the released coordinates and constraints. Changing the calculator, charge state, periodic cell, or constrained degrees of freedom changes the surface being optimized.

The potential energy surface

Within the Born–Oppenheimer approximation, the potential energy $E(\mathbf{R})$ is a function of the $3N$ nuclear coordinates $\mathbf{R} = (\mathbf{R}_1, \dots, \mathbf{R}_N)$. The force on atom $i$ is the negative gradient of that surface:

A stationary point is any geometry where all forces vanish, $\nabla E = \mathbf{0}$. Stationary points are classified by the eigenvalues of the Hessian $\mathbf{H}_{ij} = \partial^2 E / \partial R_i \partial R_j$:

• Local minimum — all eigenvalues in the internal or vibrational subspace are positive; overall translations and rotations of an isolated molecule supply zero modes. This is what optimization targets: a stable or metastable structure.
• First-order saddle point — exactly one negative eigenvalue. This is a transition state; you find these with NEB or the dimer method, not with a plain optimization.
• Higher-order saddle — several negative eigenvalues. Usually an artifact of symmetry or a bad starting guess.

The classification diagram separates numerical stopping from curvature classification. Optimization can reach a stationary point, but only the Hessian or a corresponding vibration/phonon calculation distinguishes a minimum from a saddle:

An optimizer only walks downhill into the nearest basin. A starting structure in the wrong basin — a different conformer, a different polymorph — will not cross a barrier into the right one. Global searches are a sampling problem (molecular dynamics, manual perturbation), not a local optimization problem.

The local quadratic model

Near a geometry $\mathbf R_k$, the potential energy may be expanded as

where $\mathbf g_k=\nabla E(\mathbf R_k)$ is the gradient and $\mathbf H_k$ is the Hessian. If this quadratic model were exact and $\mathbf H_k$ were positive definite, its stationary step would be

This Newton step explains both the power and difficulty of geometry optimization. The gradient says which way is downhill; the Hessian describes curvature and therefore how far the optimizer can move in each collective direction. Soft torsions require larger displacements than stiff bond stretches for the same force. A steepest-descent method ignores that distinction. A Newton method uses it directly but would require an expensive, sufficiently accurate Hessian.

Quasi-Newton methods occupy the useful middle ground. They begin with an approximate curvature model and update it from successive position and gradient differences. The approximation improves along directions the optimization has actually explored. It is therefore normal for early steps to be cautious and later steps near a smooth minimum to converge much faster.

The quadratic model is local. Far from the current structure it can predict a displacement that crosses into a region where its curvature estimate is invalid. Geometry optimizers consequently limit the step through a trust radius or maximum displacement. Tako exposes maxstepAngstrom, an atomic displacement cap for fixed-cell runs and an augmented-row cap for variable-cell runs.

The convergence criterion

For a fixed cell, Tako judges convergence on the largest per-atom force magnitude,

and declares the run converged when $f_{\max} < f_{\text{tol}}$ with $f_{\text{tol}}$ given by the fmaxEvPerAngstrom setting (default 0.01 eV/Å). Note this is the norm of the whole force vector on each atom, not the largest single component. In variable-cell optimization, the same maximum-row-norm test is applied to the augmented $N+3$ array containing atomic and scaled stress-derived cell rows; the reported value is then a generalized filter diagnostic rather than necessarily an atomic force. The check runs before the first step too: an input that already satisfies the applicable threshold returns converged: true with steps: 0.

There is no separate energy-change or displacement criterion — forces alone decide. Because forces are the gradient, a small $f_{\max}$ implies the energy is stationary to first order, which is the property that downstream calculations (vibrations, phonons, barriers) actually rely on.

Coordinate representation

An optimizer may express displacement in Cartesian coordinates, internal coordinates, or a mixture of atomic and cell variables. The representation changes the shape of the numerical problem even though the underlying potential energy is the same.

Cartesian coordinates use the $3N$ atomic components directly. They are unambiguous and work for molecules, crystals, fragments, changing connectivity, and machine-learned potentials. However, a molecular bond stretch, angle bend, and torsion are spread across several Cartesian components and can have very different characteristic curvature.

Redundant internal coordinates represent chemically recognizable motions—bond lengths, angles, dihedrals, and related primitives—then project out their algebraic redundancy. This construction relies on the Wilson $\mathbf B$ matrix,

to transform between Cartesian coordinates $\mathbf x$ and internal coordinates $\mathbf q$. This can make molecular optimization more efficient because the coordinates align better with the surface, at the cost of building, projecting, and sometimes repairing the internal-coordinate system.

Tako’s exposed structure optimizers operate in Cartesian atomic coordinates. They do not construct redundant internal coordinates, build model internal-coordinate Hessians, take rational-function steps, or fall back automatically to a different coordinate system. In Tako, difficult coordinate coupling is instead addressed through the optimizer choice, initial geometry, constraints, displacement cap, or a staged relaxation.

Variable-cell optimization augments the atomic variables with cell-deformation variables. That representation and its stress-derived generalized forces are developed in Cell relaxation.

Optimization algorithms

Tako’s engine, a Rust port of ASE optimizer conventions running as WebAssembly, exposes BFGS, FIRE, and MDMin for structure relaxation. In a fixed-cell run, all three cap the largest per-atom displacement; if any atom would move farther, the whole Cartesian step is scaled down uniformly. In a variable-cell run the same limiter acts on the augmented $N+3$ optimizer rows, so the limiting row can be a scaled cell variable rather than an atom. A limited-memory L-BFGS implementation exists in the engine but is not exposed by the current public operation.

BFGS (default)

BFGS is a quasi-Newton method: it maintains an approximate Hessian $\mathbf{H}$, built up from the observed force changes along the trajectory, and takes Newton-like steps toward the quadratic minimum. Tako’s implementation stores the full dense $3N \times 3N$ Hessian, initialized as a scaled identity,

so the very first step is plain steepest descent with stiffness $\alpha$. Larger $\alpha$ means more cautious early steps. After each step the Hessian is updated from the displacement $\Delta\mathbf{x} = \mathbf{x}_k - \mathbf{x}_{k-1}$ and the force change $\Delta\mathbf{f} = \mathbf{f}_k - \mathbf{f}_{k-1}$:

which is the classic BFGS update of Broyden, Fletcher, Goldfarb, and Shanno written in terms of forces rather than gradients. The update is skipped when the step was numerically negligible (below $10^{-7}$ bohr), which protects the Hessian from noise-driven corruption.

The step itself diagonalizes $\mathbf{H} = \sum_j \omega_j \mathbf{v}_j \mathbf{v}_j^{\mathsf T}$ and takes

with $\mathbf{g} = \nabla E$. Using $\lvert \omega_j \rvert$ instead of $\omega_j$ means that along directions of negative curvature the step still points downhill in the gradient, so the optimizer walks away from saddle points instead of converging onto them — the standard ASE-style safeguard for minimization.

The cost is one dense eigendecomposition per step, which scales as $O((3N)^3)$ with $O((3N)^2)$ memory. For the system sizes that run comfortably in a browser this is negligible next to the force evaluation itself, and BFGS’s superlinear convergence near a smooth minimum makes it the best default.

FIRE

FIRE (Fast Inertial Relaxation Engine, Bitzek et al. 2006) treats optimization as damped fictitious dynamics. Velocities $\mathbf{v}$ are optimizer state, not physical velocities. Each step computes the power

and adapts the dynamics:

• If $P > 0$ (moving downhill): the velocity is steered toward the force direction, $\mathbf{v} \leftarrow (1-\alpha)\,\mathbf{v} + \alpha\, \lvert\mathbf{v}\rvert\, \hat{\mathbf{F}}$. After more than $N_{\min}$ consecutive downhill steps the timestep grows, $\Delta t \leftarrow \min(f_{\text{inc}}\,\Delta t, \Delta t_{\max})$, and the mixing weakens, $\alpha \leftarrow f_\alpha\,\alpha$ — the system accelerates while the going is good.
• If $P \le 0$ (overshot uphill): the velocity is zeroed, $\Delta t \leftarrow f_{\text{dec}}\,\Delta t$, and $\alpha$ resets to $\alpha_{\text{start}}$ — the system stops and restarts cautiously.

Coordinates then advance by semi-implicit Euler: $\mathbf{v} \leftarrow \mathbf{v} + \Delta t\,\mathbf{F}$, $\Delta\mathbf{x} = \Delta t\,\mathbf{v}$, capped at maxstepAngstrom. Here $\mathbf{x}$ is Cartesian for fixed-cell work and the augmented filter array for variable-cell work. Tako’s FIRE uses the internal parameter set $\Delta t_0 = 0.1$, $\Delta t_{\max} = 1.0$, $N_{\min} = 5$, $f_{\text{inc}} = 1.1$, $f_{\text{dec}} = 0.5$, $\alpha_{\text{start}} = 0.1$, $f_\alpha = 0.99$ (these are engine-internal and not exposed as settings), plus a downhill check on top of the published algorithm: if a step raises the energy, the displacement is halved repeatedly (up to 20 times) until the energy drops, and the inertia is reset. FIRE needs no Hessian, tolerates terrible starting geometries, and degrades gracefully when forces are noisy — but it converges linearly, so the last stretch to a tight fmaxEvPerAngstrom takes more steps than BFGS.

MDMin

MDMin is the simplest of the three: velocity-projection dynamics. Each step the velocity is replaced by its projection onto the current force — $\mathbf{v} \leftarrow \mathbf{F}\,(\mathbf{v}\cdot\mathbf{F})/(\mathbf{F}\cdot\mathbf{F})$ when $\mathbf{v}\cdot\mathbf{F} > 0$, zeroed otherwise — then advanced by $\mathbf{v} \leftarrow \mathbf{v} + \Delta t\,\mathbf{F}$ with a fixed pseudo-timestep ($\Delta t = 0.2$, internal). It has no adaptive machinery at all, which makes it predictable and occasionally useful as a fallback when both BFGS and FIRE misbehave, but it is rarely the fastest choice.

Which one to use

For NEB reaction paths the same names appear as reactionOptimizer, but the trade-offs differ — there, stick with bfgs; see Nudged Elastic Band Method.

Implementation scope in Tako

Tako evaluates geometry steps with either the GFN2-family electronic backend or a selected machine-learned interatomic potential. The optimizer consumes energy and Cartesian forces in a shared engine representation; variable-cell work additionally consumes stress. Engine-facing atomic units are converted at the operation boundary, while reader-facing thresholds and traces are expressed in eV, ångström, and eV/Å.

The exposed minimizers are BFGS, FIRE, and MDMin. Their mathematical behavior is described above. Tako does not provide a redundant-internal coordinate optimizer, model-Hessian selection, rational-function optimization, or exact/recalculated Hessian controls. The analogous Tako decisions concern starting geometry, Cartesian constraints, optimizer family, maximum displacement, force threshold, and cell mask.

Each accepted step can contribute a stored trajectory frame and energy/force trace. Storage cadence affects diagnostic resolution and workspace size, not the optimization equations. Exact settings, defaults, promise behavior, and result fields belong to the separate Geometry Optimization with Tako Script and generated settings reference.

Cell relaxation

For a crystal, the atomic positions are only half the story — the lattice vectors are degrees of freedom too. Tako implements variable-cell optimization with a unit-cell filter in the spirit of Pfrommer et al. (1997) and ASE’s UnitCellFilter: the optimizer sees an augmented system of $N + 3$ coordinate triples. The first $N$ rows are the atomic positions expressed in the undeformed frame; the last three rows hold the deformation gradient $\mathbf{F}_{\mathrm{def}}$ (the matrix relating the current cell to the original cell, $\mathbf{h} = \mathbf{F}_{\mathrm{def}}\,\mathbf{h}_0$), scaled by a cell factor $c_f$ so cell strain and atomic displacement live on comparable numeric scales. By default $c_f = N$, the number of atoms.

The generalized force on the three cell rows comes from the stress tensor $\boldsymbol{\sigma}$ through the virial:

where $V$ is the cell volume and $\mathbf{M}$ is a strain mask that switches individual strain components on or off. Atomic forces are mapped into the same frame by $\mathbf{F}_{\mathrm{def}}^{\mathsf T}$. At convergence both the atomic forces and the masked stress have dropped below the threshold — the structure is simultaneously at (numerically) zero force and, for the released components, zero stress.

The mask is what the cellRelaxation modes control, in Voigt order $(xx, yy, zz, yz, xz, xy)$:

The exact cellRelaxation fields, defaults, and accepted values are generated in the optimization settings reference.

Three practical consequences of the implementation:

• Stress is required. Cell relaxation needs "stress" in properties and a calculator that provides cell derivatives. MLIPs generally do; for GFN2 this restricts you to the small periodic cells the backend supports. Without a cell at all, the run errors immediately.
• Convergence and traces include the cell rows. The optimizer tests $f_{\max}$ over all $N+3$ generalized rows, so a run can keep iterating after the atomic forces are below threshold while a stress-derived cell row remains large. The reported final_fmax_ev_per_angstrom, per-frame fmax, and per-frame maxstep come from that augmented optimizer array. In a variable-cell run they are generalized-coordinate diagnostics in the filter’s scaling convention, not guaranteed atomic force norms or atomic displacements.
• Fixed-cell and variable-cell energies answer different questions. A fixed-cell relaxation gives the energy at an imposed lattice; a variable-cell relaxation gives the equilibrium lattice of the model. Do not mix the two in one energy comparison without saying so.

Constraints: selective dynamics

Tako structures can carry VASP-style per-axis selective-dynamics flags. A true flag means that the corresponding Cartesian component is free to move; a false flag fixes it. When all three components are fixed the atom is immobile. These flags are part of the calculation input and must be retained with the optimized structure’s provenance.

In the engine, fully frozen atoms map to a FixAtoms constraint and partially frozen atoms to FixCartesian groups. During optimization the constraint acts through positions: after every trial step the frozen coordinates are restored, so constrained atoms cannot move no matter what forces act on them. One implementation detail to know: the optimizer’s convergence test uses the unmasked forces, so the residual force on a frozen atom still counts toward $f_{\max}$. A constrained system whose anchors carry more force than fmaxEvPerAngstrom can report converged: false even after every free atom has stopped moving — if the energy trace has flattened and the trajectory is static while fmax plateaus, that plateau is the anchor force; loosen the threshold or reconsider which atoms you froze. In molecular dynamics, NEB, dimer, vibration, and phonon runs the engine additionally zeroes force components on constrained coordinates, so frozen anchors there are fully inert. Freezing anchor atoms in reaction paths is discussed in Nudged Elastic Band Method.

Remember what a constrained result means physically: the reported minimum is a minimum of the constrained surface. Forces on the frozen atoms are generally nonzero there, and releasing them changes the energy.

The starting structure defines the search

Local optimization is conditional on its starting point. If the initial geometry is denoted $\mathbf R_0$, the algorithm normally follows a connected sequence of locally acceptable steps into one basin of attraction. A different conformer, molecular orientation, proton position, magnetic state, adsorbate site, crystal ordering, or cell shape can lead to a different stationary structure even when every numerical setting is identical.

This dependence is the mathematical consequence of searching a high-dimensional surface locally. A single relaxed structure establishes that one basin was reached from one preparation protocol; it does not establish that the structure is the global minimum or the thermodynamically dominant state.

Initial structures should satisfy three conditions before a force evaluation is trusted:

1. atom identities, ordering, cell, and connectivity represent the intended system;
2. no unintended interatomic overlap places the model deep on a repulsive wall;
3. charge, spin, model, constraints, and periodicity describe the intended state.

A severe overlap can produce forces so large that the local quadratic assumption is useless. The displacement cap can prevent one catastrophic step but cannot make the initial electronic state meaningful. Repair the geometry before optimization rather than relying on a very small step limit.

For conformational, adsorption, polymorph, or defect problems, generate an ensemble of chemically distinct starting states. Relax every member under one declared protocol, deduplicate structures using an appropriate geometric or symmetry criterion, and compare only energies that share composition and state conventions. The optimizer is one stage inside that sampling design.

Step quality and displacement control

Every proposed displacement is a prediction from a local model of the surface, and its quality can be described by comparing predicted with actual energy change — the quadratic model, the trust ratio $\rho$, and their consequences are developed below in Curvature, trust, and acceptance of a step. Tako’s exposed BFGS path uses a simple maximum-row limiter rather than an adaptive trust-radius interface: for a fixed cell the rows are atomic displacements; for a variable cell they include three scaled cell rows. maxstepAngstrom therefore bounds the optimizer coordinate update but does not report a model-quality ratio to the user.

A smaller displacement cap is useful when the start is rough, forces vary sharply, a flexible fragment can swing into another fragment, or the calculator becomes unreliable for distorted structures. An unnecessarily small cap increases the number of expensive force evaluations and can make a healthy run appear stalled. A larger cap is appropriate only when the local model and calculator remain smooth over that distance.

Energy need not decrease monotonically at every accepted step for every optimizer. FIRE and damped dynamics carry velocity-like state; a quasi-Newton curvature update can also propose a step whose finite-model error produces a small rise. The decisive evidence is whether the energy/force history stabilizes toward the intended basin without repeated large oscillation, model failure, or structural damage.

Force noise and attainable convergence

The force threshold should be stricter than the scientific resolution required downstream but looser than the irreducible numerical or model noise. Write a computed force as

where $\boldsymbol\eta$ collects SCF error, finite precision, model approximation noise, discontinuities in neighbor lists or cutoffs, and other numerical effects. If the norm of $\boldsymbol\eta$ is comparable to the requested threshold, the optimizer may alternate around a noise floor rather than converge smoothly.

For a self-consistent electronic model, force quality depends on electronic convergence. A loose or unstable SCF can give an apparently finite energy but inconsistent gradients. Tightening the geometry threshold without ensuring sufficiently stable electronic convergence is internally inconsistent. Conversely, excessively tight SCF and force thresholds increase cost without improving a screening-level model beyond its physical accuracy.

For an MLIP, smooth inference is not proof of accuracy. Near-equilibrium force errors can displace the predicted minimum, and an off-domain distorted geometry can produce confident but unphysical forces. Validate representative energies and forces against an appropriate reference. When a target force threshold is below the model’s validated force error, report the model limitation rather than implying that numerical convergence creates higher physical accuracy.

An optimizer trace that plateaus above the threshold should be examined before adding steps. Distinguish at least four patterns:

Maximum force versus other convergence tests

Optimization programs use different combinations of maximum gradient, root-mean-square gradient, displacement, and energy change. If Cartesian force components are $F_a$ over $3N$ components, common measures include

whereas Tako’s fixed-cell criterion is the maximum norm over per-atom three-vectors; its variable-cell criterion extends that maximum to the three scaled cell rows. These measures have different numerical values. A threshold copied from another optimizer measures a different quantity even when the unit is eV/Å.

Energy-change convergence can be misleading on a flat surface: a tiny energy difference between steps does not guarantee a small gradient. Displacement convergence can be misleading when a restrictive step cap prevents motion. Force convergence directly tests first-order stationarity, but it still does not classify the Hessian or establish the intended basin.

Report the exact measure and threshold. “Tightly optimized” is not reproducible without the norm definition, allowed coordinate space, calculator convergence, and whether cell degrees of freedom were active.

Translation, rotation, and soft modes

For an isolated molecule, exact translation and rotation leave the energy unchanged. Numerical forces can contain small components along these nearly zero-curvature directions. Cartesian optimization may permit the molecule to drift or rotate without changing its internal structure; this is normally harmless but can complicate trajectory comparison.

Soft torsions and intermolecular motions are physically different: they can connect shallow conformational basins separated by small barriers. A low force along a soft coordinate can coexist with a visually appreciable change in geometry because the curvature is small. Conversely, a stiff bond can have a substantial force for a tiny displacement.

Alignment should therefore be applied only for analysis and visualization, not by silently rewriting the optimization state. Compare structures after removing overall translation and rotation when appropriate, but preserve the raw coordinates, cell, and atom mapping in the calculation record.

Symmetry and constrained subspaces

An exactly symmetric starting structure can retain symmetry under symmetric forces and remain at a stationary point that is unstable to a symmetry-breaking displacement. Plain minimization has no downhill signal if the gradient is exactly zero. A Hessian calculation can reveal the negative-curvature mode; a small physically motivated displacement along it then supplies a new starting geometry.

Tako does not impose a general molecular or crystallographic symmetry constraint during the exposed optimization. Numerical noise may break symmetry, or a structure may relax into a lower-symmetry basin. If symmetry must be retained as a scientific approximation, express the allowed degrees of freedom through an explicit structure/constraint protocol and document it. Do not infer that a symmetric final geometry proves stability.

Selective dynamics restricts the coordinates that survive a trial update, but Tako’s fixed-cell optimizer and convergence test still receive the unmasked calculator force described in Constraints: selective dynamics. The intended physical object is a stationary structure within the released coordinate subspace; the numerical result can nevertheless remain unconverged when a frozen component carries a large raw force. This is often appropriate for deep substrate layers, embedded fragments, or controlled scans, but the mobility mask and constrained force pattern must be inspected to ensure that the approximation is not carrying the phenomenon of interest.

Cell optimization as a thermodynamic model

Atomic relaxation at fixed cell minimizes $E(\mathbf R;\mathbf h)$ with respect to $\mathbf R$. Variable-cell relaxation changes the question by adding strain coordinates. At imposed external stress $\boldsymbol\sigma^{\mathrm{ext}}$, the generalized target is related to an enthalpy-like potential. Under hydrostatic pressure $p$, the familiar scalar form is

At zero external pressure, a fully relaxed equilibrium cell should have small residual stress under the calculator’s convention. A force-only convergence flag does not by itself prove stress convergence unless the engine’s combined atom/cell generalized-force criterion is understood. Inspect the final tensor and the cell history.

The five cell modes tabulated in Cell relaxation define different constrained thermodynamic models, not interchangeable numerical options. A lengths-only relaxation cannot discover a shear instability. A full relaxation can change crystal symmetry, orientation, or cell representation. For a slab, relaxing the vacuum direction usually has no intended physical meaning; mask it. For a molecule in an artificial periodic box, variable-cell relaxation is normally meaningless.

Cell scaling also creates conditioning challenges because atomic displacements are measured in ångström while strain is dimensionless. The cell filter maps strain-like variables into an optimizer coordinate scale. Its cellFactor changes numerical conditioning, not the underlying energy. An extreme value can make the optimizer move atoms while barely changing the cell, or overemphasize cell response relative to atomic relaxation.

Staged optimization protocols

A staged protocol is justified when each stage changes a declared numerical approximation while preserving the scientific state needed for the final result. Common patterns include:

1. repair or pre-relax a rough geometry with a conservative displacement cap and loose force threshold;
2. continue with BFGS and the target threshold on the same calculator;
3. evaluate a final single point with the requested analysis channels;
4. perform vibrations or phonons when stationary-point classification is required.

Changing calculator between stages can also be useful, but it changes the surface. The final stage must reoptimize sufficiently under the reported final method; a high-level single point on a low-level minimum is a composite fixed-geometry protocol, not a high-level optimized structure.

When a run reaches its step budget, continuing from the last geometry starts a new optimizer state. The structure is warm, but a full-memory BFGS curvature estimate is not by itself a portable checkpoint across independent calls. The second run can take several steps to rebuild useful curvature. Record both workspaces or combine their trajectories explicitly in the analysis rather than presenting them as one uninterrupted optimizer history.

Curvature, trust, and acceptance of a step

An optimizer does not know the potential-energy surface away from the points already evaluated. Every proposed displacement is therefore a prediction made from a local model. For a quasi-Newton method, write that model around iteration $k$ as

where $\mathbf{s}$ is the trial optimizer-coordinate displacement, $\mathbf{g}_k$ is the corresponding gradient, and $\mathbf{B}_k$ is the current Hessian approximation. Minimizing the quadratic without a displacement restriction gives $\mathbf{s}=-\mathbf{B}_k^{-1}\mathbf{g}_k$. That expression is useful only while the quadratic represents the sampled basin. Far from a minimum, an inaccurate curvature model can predict a step that is too large, points across a ridge, or drives two atoms through one another. Max step A is therefore not merely a performance control. It restricts extrapolation in the Cartesian atomic array for fixed-cell work and in the augmented atomic/cell array for variable-cell work.

The quality of a trial step can be described by the ratio

The numerator is the energy reduction actually observed; the denominator is the reduction predicted by the local model. A value near one means that the model predicted the energy change well. A negative value means that the proposed downhill step increased the energy. Tako does not expose $\rho_k$ as a result field and its optimizers do not present a user-level trust-radius history, so it would be incorrect to claim that a displayed energy trace is an explicit trust-region diagnostic. The ratio nevertheless explains why conservative displacement limits help poor starting structures and why a curvature history built in one basin should not be assumed valid after a discontinuous rearrangement.

BFGS updates its approximate curvature from the displacement $\mathbf{s}_k$ and gradient change $\mathbf{y}_k=\mathbf{g}_{k+1}-\mathbf{g}_k$. The secant condition $\mathbf{B}_{k+1}\mathbf{s}_k=\mathbf{y}_k$ asks the new matrix to reproduce the curvature observed along the last displacement. When $\mathbf{s}_k^{\mathrm T}\mathbf{y}_k$ is small or negative, the sampled curvature is uninformative or inconsistent with a positive-definite minimum model. This occurs near inflection regions, after abrupt bonding changes, or when force noise is comparable with the gradient change. A nominally sophisticated Hessian update cannot create reliable curvature information from such data; shorter steps, a better initial structure, or a dynamics-based relaxation stage may be more defensible than repeatedly restarting BFGS with a larger budget.

FIRE and MDMin treat the force as a dynamical driving term rather than constructing an explicit inverse Hessian. They can dissipate motion through rough parts of the surface, but their trajectories are still numerical optimization paths rather than physical molecular dynamics. The fictitious velocities, damping behavior, and accepted steps must not be interpreted as time-dependent molecular motion or used to infer kinetic barriers. Their advantage is algorithmic: directional persistence can move through an initially poor region without requiring a trustworthy curvature matrix.

In an ideal projected constrained optimization, a projector $\mathbf{P}$ would define the mobile subspace, with projected gradient $\mathbf{P}^{\mathrm T}\mathbf{g}$ and curvature $\mathbf{P}^{\mathrm T}\mathbf{H}\mathbf{P}$. This mathematics explains stationarity only with respect to admitted motions, but it is not Tako’s fixed-cell force convention. Tako restores constrained coordinates after a trial position update while the atomic optimizer and convergence test still receive the unmasked calculator forces. A frozen anchor can therefore dominate fmax, prevent converged: true, and influence the numerical update even though its coordinate is restored. Interpret the plateau using the trajectory and mobility mask; do not relabel the reported force as projected. Removing the constraint creates a different optimization problem and may release the stored force.

Cell relaxation adds another conditioning problem because atomic displacements have length units while strain is dimensionless. A cell factor supplies a numerical scale between these blocks. If cell variables are scaled too aggressively, the lattice can change faster than atoms accommodate the new environment; if scaled too weakly, atomic forces converge while stress remains far from the intended condition. The final decision should therefore inspect both geometry and lattice behavior, not only a single threshold. Tako records the normalized cell mode and masks in result.settings.cell_relaxation, but cell_factor is the requested optional value: it is null when omitted, and null or a nonpositive request makes the engine use its atom-count default without writing that resolved value back into the result. Reproduction must therefore retain both the request and the documented default rule.

Cost and scaling

The dominant cost is normally the number of calculator evaluations multiplied by the cost per energy/force evaluation. If $n_{\mathrm{eval}}$ is the number of force calls,

BFGS additionally stores a dense matrix for the released Cartesian variables and diagonalizes it during a step, producing quadratic memory and cubic algebraic scaling in coordinate count. For small browser-scale systems the calculator often dominates; for large structures, the dense optimizer can become material. FIRE and MDMin carry less curvature state and can be attractive for large or rough relaxations even when they require more force evaluations.

Variable-cell calculations require stress and can increase calculator cost. Tight electronic thresholds, implicit solvent, large periodic cells, and expensive MLIP weights affect each evaluation. Saving every trajectory frame increases workspace storage but not the model evaluation count. Choose the frame interval according to diagnostic needs without discarding the final structure or convergence trace.

Reading results

An optimization result contains the numerical convergence decision, accepted-step count, final energy, final convergence diagnostic, stored trajectory frames, and final atomic coordinates. That diagnostic is a maximum per-atom force norm for fixed-cell runs and a maximum over atomic plus scaled cell-force rows for variable-cell runs. Periodic calculations additionally retain the cell and can retain the final stress when it was requested and supported. Exact field names and units are defined in the generated optimization reference, not repeated here.

The raw result must be read as a connected record. A final energy without the convergence flag and optimizer trace is insufficient evidence. During variable-cell relaxation, inspect the final stress tensor and cell history together with the reported generalized-row maximum; do not label that maximum as a purely atomic force. A sparse stored trajectory can still establish the final state, but it may be inadequate for diagnosing an early basin change or force spike.

Establishing that the result is a minimum

Force convergence establishes first-order stationarity within the released coordinate space. Minimum assignment requires curvature evidence. For a nonlinear isolated molecule, the mass-weighted Hessian should yield $3N-6$ internal modes after translational and rotational modes are treated; a linear molecule has $3N-5$. A physically meaningful imaginary mode reveals negative curvature and therefore contradicts a local-minimum assignment.

Very small imaginary or low-frequency modes require judgment. They can arise from incomplete force convergence, finite-difference noise, weak intermolecular motion, approximate rotational/translational removal, periodic acoustic modes, or genuine shallow negative curvature. Tighten the reference geometry and numerical differentiation consistently before declaring them artifacts. Inspect the displacement pattern: an internal bond rearrangement is different from uniform translation or a nearly free fragment rotation.

For a periodic crystal, phonon stability at one wave vector is not stability throughout the Brillouin zone. A unit-cell Hessian or zone-center calculation can miss a supercell distortion. The phonon chapter defines the necessary reciprocal-space analysis. Geometry optimization prepares the stationary reference but does not replace dynamical-stability evaluation.

Repeated relaxation from perturbed structures provides complementary evidence about the basin. Small random or mode-directed displacements should return to the same minimum when the basin is robust under the model. If perturbations relax reproducibly into distinct minima, report the ensemble and their relative energies rather than selecting one run as uniquely correct.

What to report

The shared record skeleton and reporting standard are defined in Reproducibility and Reporting. Beyond that skeleton, a geometry-optimization account specifically identifies:

• the exact force-norm definition and threshold (fixed-cell atomic maximum versus variable-cell generalized rows);
• fixed-cell or variable-cell policy, released strain components, cell factor, and imposed external stress when applicable;
• optimizer, maximum displacement, step budget, and any staging policy;
• enough trajectory/trace information to establish basin continuity;
• the curvature, perturbation, or comparison evidence used to assign the stationary point.

Comparing optimized structures

Energy comparison between separately optimized structures is a comparison of local minima only when each run uses the same method and compatible state, composition, constraints, and cell policy. If one candidate is constrained and another is free, or one cell is fixed and another relaxed, their energies answer different variational questions.

Structural comparison requires atom mapping and an alignment convention. Root-mean-square displacement after optimal rigid alignment is useful for equal molecular compositions but can hide torsional or permutational differences. Crystals additionally require symmetry/cell standardization and periodic atom matching. Report chemically meaningful internal coordinates, coordination changes, volume or strain, and symmetry alongside any scalar RMSD.

An optimized energy difference remains an electronic or model energy difference unless vibrational, thermal, entropic, standard-state, pressure, or reservoir terms are added explicitly. Geometry relaxation changes the evaluated coordinates but adds no enthalpy or Gibbs free-energy terms.

Worked optimization protocol

Consider a molecular conformer intended for a subsequent vibrational calculation. The starting geometry is first checked for atom identity, charge/spin state, and close contacts. A GFN2-family model with the declared dispersion treatment defines the surface. The first relaxation uses a conservative displacement cap and a moderate force threshold; its purpose is to leave the repulsive or poorly conditioned region without allowing a large structural jump.

If the trajectory shows smooth relaxation into the intended conformational basin, a second BFGS stage begins from the first stage’s final geometry and applies the target force threshold. The second workspace is a new numerical run: it rebuilds its curvature approximation even though its coordinates are warm. Its result is accepted only if the declared maximum per-atom force norm is below threshold and the geometry preserves the intended constitution and conformer.

Finally, a vibrational calculation classifies local curvature. No meaningful imaginary mode supports assignment as a local minimum under the combined numerical protocol; an imaginary torsion or rearrangement mode indicates that the force-converged geometry is not the desired minimum. The complete record retains both optimization trajectories, the final structure, calculator state, convergence measures, and vibrational modes.

For a crystal, the same logic includes an explicit cell policy. A fixed-cell atomic relaxation answers whether atoms are stationary at the imposed lattice. A lengths-only stage tests lattice constants without shear. A full-cell stage releases all strain components and can change symmetry or cell representation. The selected protocol must match the physical boundary conditions rather than the number of exposed variables.

This example is intentionally method-level rather than executable. Exact Script construction and output writing are isolated in the Script chapter, while the interface sequence is isolated in Optimize a Structure.

Diagnostics and failure modes

Sketched geometries can diverge catastrophically under GFN2 — the cramped-cell failure mode. A hand-built prototype with interatomic distances below roughly 1.5 Å puts atoms deep on the repulsive wall. Instead of relaxing, the first steps launch atoms outward, the SCF struggles, and the “optimization” wanders off to a nonsensical geometry — sometimes even reporting large negative energies along the way. This is a verified failure mode in real Tako runs. Before optimizing a sketched structure, inspect chemically relevant distances and separate unintended overlaps before the first force evaluation. Reducing maxstepAngstrom limits the damage per step but does not repair a bad input.

Choose the threshold for the job. The 0.01 eV/Å default is a good general-purpose target. Loosen to 0.02–0.05 eV/Å for screening large candidate sets; tighten only when a downstream calculation demands it (phonons and vibrations are the classic case — residual forces show up as spurious imaginary modes). Beware that very tight thresholds can chase model noise: an MLIP’s force error near a minimum puts a floor under the reachable $f_{\max}$, and FIRE or more steps will not push through it.

Noisy or rising energy trace. Open energy_trace.json. A healthy run shows a broadly downward trend that stabilizes; individual accepted steps need not be monotonic. Repeated large oscillation usually means steps are too aggressive — lower maxstepAngstrom (e.g. 0.02), or raise alpha to stiffen early BFGS steps. If BFGS keeps overshooting on a rough start, switch to fire, converge coarsely (0.05 eV/Å), then finish with bfgs at the final threshold.

Converged does not mean minimum. $f_{\max}$ is small at any stationary point, including saddles left by a symmetric starting guess. The BFGS $\lvert\omega\rvert$ safeguard pushes away from saddles it can detect, but the definitive check is a vibration or phonon calculation (molecules or crystals, respectively): imaginary frequencies mean you are not at a minimum — displace along the offending mode and re-optimize.

Pick the backend by chemistry. GFN2 for organic molecules, lighter elements, and small periodic cells (Bannwarth et al. 2019 parametrizes $Z \le 86$); MLIP for dense inorganic crystals, alloys, surfaces, and battery phases where the model’s training domain covers your elements. Energies from different backends are never comparable — see Calculation Overview.

Cell relaxation fails or stalls. Confirm that the structure has a meaningful periodic cell, the calculator supplies stress, the released strain components match the physical boundary conditions, and the final stress tensor is interpreted under the backend convention.

References

1. C. G. Broyden, The convergence of a class of double-rank minimization algorithms, J. Inst. Math. Appl. 6, 76 (1970); R. Fletcher, Comput. J. 13, 317 (1970); D. Goldfarb, Math. Comp. 24, 23 (1970); D. F. Shanno, Math. Comp. 24, 647 (1970). — The BFGS update.
2. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97, 170201 (2006). — FIRE.
3. B. G. Pfrommer, M. Côté, S. G. Louie, and M. L. Cohen, Relaxation of crystals with the quasi-Newton method, J. Comput. Phys. 131, 233 (1997). — Combined atom and cell relaxation.
4. H. B. Schlegel, Geometry optimization, WIREs Comput. Mol. Sci. 1, 790 (2011). — Survey of optimization methods on potential energy surfaces.
5. A. Hjorth Larsen et al., The atomic simulation environment — a Python library for working with atoms, J. Phys.: Condens. Matter 29, 273002 (2017). — The optimizer conventions Tako’s engine follows.
6. C. Bannwarth, S. Ehlert, and S. Grimme, GFN2-xTB — an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions, J. Chem. Theory Comput. 15, 1652 (2019). — The GFN2 backend.