Refitting tabulated EAM potentials

tadah refit converts a tabulated EAM potential (a setfl file in the LAMMPS eam/alloy / eam/fs format) into a native Tadah! potential. The three tabulated functions — pair \(\phi(r)\), electron density \(\rho(r)\) and embedding \(F(\rho)\) — are re-expressed as small basis-function expansions and fitted, producing a native pot.tadah. Its descriptors follow REFIT_BASIS: a D2_* pair descriptor (e.g. D2_Knot5, D2_Blip) together with the DM_REAM many-body descriptor. The result is usable with tadah predict, the Tadah!LAMMPS plugin, and nested fitting. refit also writes a round-trip setfl (for direct LAMMPS comparison) and a tadah hpo starter target.

Typical workflow

  1. Write a refit.cfg naming the input EAM and the basis (see Quick start).

  2. Run refit on a setfl EAM: tadah refit -c refit.cfg -v 2.

  3. Read the crystal self-test at the end of refit_opt.log (energy–volume accuracy of the fitted model vs the original table).

  4. Use pot.tadah directly, or hand refit.hpotarget to tadah hpo for nested fitting.

Quick start

Minimal config (all other settings take their defaults):

EAM             Ti1.eam.fs
REFIT_OBJECTIVE curve
REFIT_BASIS     Knot5

Run it:

tadah refit -c refit.cfg -v 2

A fuller curve-mode config (pinned window edges, explicit basis ceiling):

EAM             Ti1.eam.fs
OUTFILE         pot.tadah
SETFL_OUT       refit.eam
REFIT_OBJECTIVE curve
REFIT_BASIS     Knot5            # Knot3 suits some elements (see REFIT_BASIS)
REFIT_COMPACT_NMAX 35            # one value, or "phi rho F" e.g. 35 12 35

INIT
  PHI XMIN PIN
  RHO XMIN PIN
  F   XMIN PIN
  F   XMAX PIN
ENDINIT

The pipeline stages EXPLORE / REFINE / FUSE use built-in defaults when omitted; see How the fit works for what each does and Stage blocks to override them.

How the fit works

Each function is fitted by a basis-count ladder: start from a small seed, then add one basis function per rung at the largest residual, until the per-function ceiling REFIT_COMPACT_NMAX or a stall. Each rung is optimised in stages:

  • EXPLORE — a global seed then a per-rung local search that place the basis-function centres and widths. The optimiser at each stage is configurable (see Stage blocks).

  • REFINE — a high-precision local polish of centres, widths and coefficients (Ceres LM).

  • FUSE — a final refinement: in curve mode a per-function polish of the weighted curve, in energy mode a joint fit of all three functions to the crystal energy–volume curve.

REFIT_OBJECTIVE selects what every stage optimises. curve fits the tabulated functions directly and concentrates accuracy near the equilibrium structure (see REFIT_CURVE_WEIGHTING); energy fits the crystal energy–volume curve.

Reading the result

The self-test builds an ideal crystal from the lattice type and constant in the EAM header (the <Z> <mass> <a0> <lattice> line). fcc, bcc and hcp are supported (case-insensitive); hcp uses the ideal \(c/a=\sqrt{8/3}\). If the header declares no supported lattice or a non-positive constant (for example ... 0.0 NONE), refit stops with an error — edit that line to the actual ground-state lattice and lattice constant and rerun. The test compares the original and fitted models on the same geometry, so the reported errors are exact even when the header constant is not the potential’s own energy minimum.

The end of refit_opt.log reports a crystal self-test over \(a/a_0 \in [0.75, 1.10]\):

ambient    RMSE=0.05 meV/atom  max=0.10
compressed RMSE=0.20 meV/atom  max=0.40
cohesive: original=-5.34 refit=-5.34 delta=0.05 meV/atom

ambient is \(a/a_0\,0.95\!-\!1.10\), compressed is \(0.75\!-\!0.95\); REFIT_TOL_AMBIENT / REFIT_TOL_COMPRESSED are the pass thresholds. cohesive delta is the equilibrium-energy error.

Configuration reference

General and output

Key

Default

Effect

EAM (path)

—

Input setfl EAM file to refit (required).

OUTFILE (path)

pot.tadah

Output native Tadah! potential.

SETFL_OUT (path)

refit.eam

Round-trip setfl written from the fitted model for LAMMPS comparison.

HPOFILEOUT (path)

refit.hpotarget

Generated tadah hpo starter target.

REFIT_LOG (path)

refit_opt.log

Live, line-flushed optimisation log; the self-test is at its end.

OUTPUT_GRID_SCALE (float)

1.0

Output setfl grid density relative to the input table.

OPTIM_PCT (float)

5.0

±% half-width of the OPTIM bounds written into the HPO target.

Mode and basis

Key

Default

Effect

REFIT_MODE

compact

compact fits a small optimised basis via the ladder; dense fits fixed grids (see Dense mode).

REFIT_BASIS

Knot5

Basis family (compact): Knot5 quintic knots, Knot3 cubic knots, KnotMix per-knot cubic/quintic, Blip, Blip5. Knot5 is smoothest per knot; Knot3 captures sharper features. Best family is element-dependent.

Objective and curve weighting

Key

Default

Effect

REFIT_OBJECTIVE

energy

curve fits the tabulated functions (value + derivatives) with pressure weighting; energy fits the crystal energy–volume curve. Use curve for a faithful crystal from a small basis.

REFIT_CURVE_WEIGHTING

true

(curve) Weights \(\phi,\rho\) fully between the equilibrium nearest-neighbour distance and the cutoff, and \(F\) from zero to the equilibrium density, relaxing toward compression — tighter equilibrium energies for a given basis size.

REFIT_CURVE_RELAX

0.05

(curve, weighting on) Floor weight in the compressed region, in (0, 1]; smaller focuses the fit harder on the equilibrium region.

REFIT_CURVE_DERIV_WEIGHT

0.0

(curve) Weight of the 1st/2nd-derivative terms relative to value. 0 fits values only (best crystal energies); raise to fit slope/curvature shape.

Compact ladder

Key

Default

Effect

REFIT_COMPACT_NMAX

40

Per-function basis ceiling. One value (all functions) or three phi rho F (e.g. 35 12 35). Higher is more flexible and slower.

REFIT_COMPACT_SAMPLES

800

Sample points per function for the fit and RMSE (raised to at least 4×NMAX). More is finer and slower.

REFIT_COMPACT_OPT_EVERY

1

Run the full per-rung optimisation every K-th rung; intermediate rungs grow by insertion only. Higher is faster and coarser.

REFIT_COMPACT_STALL

12

Stop a function’s ladder after this many rungs without improvement (Blip families: 6).

REFIT_COMPACT_INSERT_TOPK

3

Try this many top residual peaks per rung and keep the best insertion (Blip families: 1).

REFIT_COMPACT_JITTERS

2

Randomised restarts per rung (jitters+1 candidates); more explores the landscape harder.

REFIT_COMPACT_SHARE

0.5

Per-function share of the tolerance budget used by the rung-selection metric.

REFIT_COMPACT_SEED

12345

RNG seed for the jittered restarts; fixes the run for reproducibility.

REFIT_COMPACT_MERGE

auto

Knot coalescence guard: penalises near-coincident knot pairs during EXPLORE and merges them after REFINE, re-inserting a fresh knot at the worst residual. auto is off for curve and on for energy; true/false force it. Knot families only.

Window and self-test

Key

Default

Effect

REFIT_TOL_AMBIENT (meV/atom)

0.1

Pass threshold for the ambient band (\(a/a_0\,0.95\!-\!1.10\)); also weights the energy-mode metric there.

REFIT_TOL_COMPRESSED (meV/atom)

1.0

Pass threshold for the compressed band (\(a/a_0\,0.75\!-\!0.95\)).

REFIT_REACH_SMIN

0.75

Smallest \(a/a_0\) the fit must stay faithful at; sets the reachable \(r\)/\(\rho\) window.

REFIT_REACH_MARGIN

1.25

Safety margin on the reachable \(F(\rho)\) window.

REFIT_R_LO_FRAC

0.20

Inner radial fit bound as a fraction of the cutoff; also the lower edge of the \(\phi/\rho\) weight ramp.

REFIT_PLACEMENT

window

window places the basis on the reachable window; uniform spreads it over the full range.

Stage blocks

Four optional multi-line blocks tune the pipeline; omit any to use defaults.

INIT declares basis placements (no optimisation), one per line:

INIT
  <PHI|RHO|F>  <placement>  [PIN]
ENDINIT

placement is MANUAL <mu> [<eta>] (fixed centre, plus width for Blip families), XMIN / XMAX (window edges; XMAX applies to F only), YMIN / YMAX (target minimum/maximum location) or SMAX / CMAX (maximum slope/curvature location). PIN freezes that centre through all later stages (its coefficient is still fitted). Pinning the window edges sets the extrapolation behaviour.

EXPLORE / REFINE / FUSE each accept ENABLED true|false (REFINE/FUSE; ENABLED false skips the stage) and a nested OPTIMIZER block. EXPLORE also accepts MAX_ADD (one value or three phi rho F) capping the basis functions added in that stage.

EXPLORE
  MAX_ADD 45
  OPTIMIZER
    LIB NLOPT  ALGO GN_MLSL_LDS  MAXEVAL 1000000
    INNER
      LIB NLOPT  ALGO LN_SBPLX  MAXEVAL 50000  FTOL_REL 1e-7
    ENDINNER
  ENDOPTIMIZER
ENDEXPLORE
REFINE
  OPTIMIZER
    LIB CERES  ALGO LM  MAXEVAL 1000  FTOL_REL 1e-14  GTOL 1e-14
  ENDOPTIMIZER
ENDREFINE
FUSE
  OPTIMIZER
    LIB CERES  ALGO DOGLEG  MAXEVAL 2000  FTOL_REL 1e-14
  ENDOPTIMIZER
ENDFUSE

Legal optimisers and the stop-criteria each stage consumes:

Stage

LIB / ALGO

Stop-criteria used

EXPLORE outer

NLOPT GN_MLSL_LDS | GN_DIRECT_L | GN_CRS2_LM; DLIB GLOBAL; UNIFORM

MAXEVAL

EXPLORE INNER

NLOPT LN_SBPLX | LN_NELDERMEAD | LN_BOBYQA | LN_PRAXIS | LD_LBFGS; DLIB LBFGS; NONE

MAXEVAL, FTOL_REL, XTOL_REL

REFINE

CERES LM | DOGLEG; NLOPT LD_LBFGS; DLIB LBFGS; NONE

MAXEVAL (alias MAXITER), FTOL_REL, GTOL

FUSE

CERES LM | DOGLEG; NONE

MAXEVAL (alias MAXITER), FTOL_REL

MAXEVAL is the NLopt evaluation budget / Ceres iteration cap; FTOL_REL on the Ceres stages is Ceres’ relative function tolerance.

Output files

  • pot.tadah — fitted native potential.

  • refit.eam — round-trip setfl for LAMMPS.

  • refit.hpotarget — tadah hpo starter target.

  • refit_opt.log — optimisation log and crystal self-test.

  • refit_basis_info.dat — final basis: function, centre, width, coefficient, degree.

  • refit_curve_{phi,rho,F}.dat — fitted vs tabulated curves and derivatives.

  • refit_compare.dat, refit_range_report.dat, refit_ladder_trace.dat — energy–volume comparison, reachable window, per-rung ladder trace.

Dense mode

REFIT_MODE dense fits fixed basis grids instead of the ladder; set the grid sizes directly. Large dense grids reach very high accuracy, but every basis function becomes a search dimension in nested fitting — a high-count dense potential makes the downstream tadah hpo search much larger. Prefer compact mode (a small optimised basis) when the refit will be tuned further.

Key

Default

Effect

REFIT_LEVEL

balanced

Preset grid density: coarse | balanced | fine | accurate.

N2B

—

\(\phi(r)\) basis count.

NRHO_BASIS

—

\(\rho(r)\) basis count.

NF_BASIS

—

\(F(\rho)\) basis count.

REFIT_MAX_ESCALATIONS

4

Automatic basis-escalation rounds while above tolerance.

REFIT_CONSTRAIN_F0

false

Pin \(F(0)\) to the table value.

Dense REFIT_BASIS families are Blip | Blip5 | Gaussian.

Editing a refitted potential (--edit)

tadah refit --edit adds or removes ONE basis function of one curve of an existing refit potential — no EAM input, no new fit:

# insert a pair-function basis function at r = 2.75 A
tadah refit --edit pot.tadah --func phi --add 2.75

# remove the embedding basis function nearest rho = 0.9
tadah refit --edit pot.tadah --func F --remove 0.9

--func selects the curve: phi (pair, CGRID2B), rho (density, CGRIDMB / AMPGRIDMB) or F (embedding, CEMBFUNC). Exactly one of --add <centre> / --remove <centre> is required. Both knot (Knot5 / Knot3) and blip (Blip / Blip5) family pair descriptors are supported; D2_mJoin and Gaussian potentials are not, and BIAS / NORM must be false (their columns are WEIGHTS-aligned and cannot survive an edit).

The many-body term is optional:

  • pure two-body potentials (TYPE2B only — e.g. hand-built knot pair potentials) are fully supported through --func phi;

  • DM_REAM potentials (as written by tadah refit) expose all three curves;

  • any other many-body descriptor restricts editing to --func phi — the many-body keys and their WEIGHTS block ride along untouched (the pair block is always the WEIGHTS prefix).

Adding inserts the function with coefficient 0, so the edited potential predicts bit-identically to the input — value, every derivative, and the extrapolation region. The new coefficient only adds capacity: it appears as an extra OPTIM FIXWEIGHT line in the regenerated hpotarget, ready for a tadah hpo search. Knot families take an optional --deg 3|5 (default quintic; 3 materialises the per-knot KDEG* degree arrays) and reject --eta (the knot width is a pure gauge). Blip families take an optional --eta (default: the mean width of the two nearest neighbours); a blip whose support lies entirely beyond RCUT is rejected.

Removing deletes the basis function whose centre is nearest the given value (the choice is reported). Unless its coefficient is exactly 0, the two remaining nearest-in-centre coefficients are least-squares adjusted against the removed contribution’s value, first and second derivative over the region it controlled (per-derivative-order normalisation keeps the orders balanced). Removal is inherently lossy — the log reports the residual before and after compensation; expect observables to change and re-validate the potential. Removing the outermost radial knot also shrinks the effective cutoff and the uncovered tail cannot be compensated (warned).

Outputs default to pot_edited.tadah and refit_edited.hpotarget — the input file is never overwritten silently (an explicit -o / --hpofileout wins; writing onto the input is refused). All bookkeeping (TYPE2B / TYPEMB count tokens, WEIGHTS, FIXINDEX/ FIXWEIGHT, KDEG*, RCUT2B/RCUTMB) is regenerated consistently.

Rewriting HPO search bounds (--retarget)

The hpotarget written by refit uses ±OPTIM_PCT (default 5 %) bounds. --retarget rewrites them without re-running anything:

# widen every bound to +-10%, coefficients to +-2%
tadah refit --retarget refit.hpotarget --pct 10 FIXWEIGHT=2

# absolute bounds for one parameter, edit the file in place
tadah refit --retarget refit.hpotarget --set 'CGRID2B(3)=2.1,2.4' --in-place

Every OPTIM line is rewritten — active and commented alike, preserving the comment state — around the parameter values recovered from the existing bounds; everything else in the file is preserved byte-for-byte. --pct takes a bare number (the global percentage) and/or KEY=N per-key overrides; --set 'KEY(i)=lo,hi' (repeatable) is applied last. Output goes to <input>_edited.hpotarget by default, -o <name>, or back into the input with --in-place.

Retargeting a freshly-emitted file at the same percentage reproduces it byte-for-byte, and any retargeted file is a fixed point of further retargeting at the same settings. Bound recovery is exact at the printed precision (8 significant digits); a hand-pruned group (deleted lines) shifts the group’s scale floor relative to a fresh emission (warned).

Cutoffs: knot channels are derived, blip channels are fixed

The two basis families treat RCUT2B / RCUTMB differently:

Family

Effective cutoff

RCUT keys

Knot (Knot5/Knot3)

The largest knot position — each knot function \((t - r)_+^k\) vanishes \(C^4\)-smoothly at its own knot, so nothing exists beyond the largest one.

Derived. Every consumer (refit, --edit, --retarget, train, predict, hpo) recomputes them from the largest knot; a differing stated value is corrected with a warning. Moving RCUT anywhere at or beyond the largest knot is prediction-exact.

Blip (Blip/Blip5)

The stated RCUT — blips near the cutoff overhang it and are hard-truncated there by design (the truncation was present during the fit).

Authoritative. Never auto-changed; warnings only (a centre beyond RCUT, or OPTIM bounds that let a blip’s support cross it).

Practical consequences for knot potentials:

  • To change the effective cutoff, move the CGRID2B / CGRIDMB knot positions (OPTIM CGRID2B in the hpotarget, or --edit --add beyond the current largest knot). Do not search OPTIM RCUT2B / RCUTMB — the keys are re-derived on every HPO evaluation, so that dimension is wasted (warned at startup).

  • During tadah hpo the neighbour list is sized once for the largest cutoff the search can reach (including the upper bounds of OPTIM’d knot positions), and the RCUT keys follow the moving knots automatically; the written best_pot.tadah is always self-consistent.

  • Hand-editing knot positions in a pot.tadah is safe: the stale RCUT is detected and corrected (with a warning) by whichever command reads the file next.