autoFF

Automated free-energy simulation and force-field parameter fitting with Tinker and Tinker9.

One configuration file describes an entire study: any number of solutes whose hydration free energies you want, any number of neat liquids whose densities you want, and any number of dimer geometries. autoFF either reports every one of those properties (single-point) or fits a shared parameter file so they all match reference values at once (optimize).

• Install
• Quick start
• How it works
• Configuration reference
• Job types
• Directory layout
• Command reference
• Examples
• Notes and limitations

git clone https://github.com/leucinw/autoFF.git
cd autoFF
pip install -e .

Dependencies are numpy, scipy and ruamel.yaml; environment.yml provides a conda environment. Installing puts two commands on your path: autoff and autoff-submit.

You also need working Tinker8 and Tinker9 builds. Their locations come from a small environment file listing $DYNAMIC8/9, $BAR8/9, $ANALYZE8/9 and $MINIMIZE8/9. The bundled default is autoff/data/tinker.env; point shared.tinker_env at your own copy:

export TINKER8=/path/to/tinker
export  DYNAMIC8="$TINKER8/dynamic"
export  ANALYZE8="$TINKER8/analyze"
export      BAR8="$TINKER8/bar"
export MINIMIZE8="$TINKER8/minimize"

export tk9home=/path/to/tinker9/build
export  DYNAMIC9="$tk9home/dynamic9"
export  ANALYZE9="$tk9home/analyze9"
export      BAR9="$tk9home/bar9"
export MINIMIZE9="$tk9home/minimize9"

cd examples/Phenol-HFE

autoff run config.yaml --dry-run   # generate every input file, submit nothing
autoff run config.yaml             # submit, wait, and report
autoff check config.yaml           # progress of a run already under way

--dry-run is the fastest way to see what a config will actually do: it writes every key file and job script and records the intended submissions in results/submitted_jobs.txt, without touching the cluster.

Hydration free energy. A solute is decoupled from water in stages: first its electrostatics are switched off, then its van der Waals interactions. Each stage is a lambda window with its own MD trajectory. Consecutive windows are combined with Bennett Acceptance Ratio (BAR), and the per-window free energies sum to the transfer free energy. The same decoupling is repeated in the gas phase and subtracted, which removes the solute's intramolecular contribution. Solutes with fewer than five atoms have no such contribution, so their gas leg is skipped automatically.

Neat-liquid density. One NPT trajectory per temperature; the density is the mean over production frames, after discarding an equilibration segment.

Dimer energies. E_int = E_dimer - E_mon1 - E_mon2 from single-point energy evaluations, optionally after letting the trial parameters relax the dimer geometry.

Fitting. Every property with a reference value becomes one entry in a single residual vector:

residual = weight * (calculated - reference) / denominator

The denominator normalizes across units, so a density in kg/m³ and a free energy in kcal/mol pull on the fit comparably. It defaults to the spread (standard deviation) of the reference values when several exist, and to sqrt(|value|) for a lone one — the ForceBalance convention. scipy's least_squares (trust-region reflective, soft-L1 loss) minimizes the sum of squares subject to the bounds you give.

Why fitting is affordable. Derivatives never rerun dynamics. To evaluate a perturbed parameter set, autoFF adds an extra window at a fictitious lambda above 1.0 that reweights the existing fully-coupled trajectory — so a trial parameter set costs one BAR evaluation instead of a full simulation. Density derivatives use the fluctuation formula (Eq. 4 of Wang et al., J. Chem. Theory Comput. 2013):

d⟨ρ⟩/dλ = -β ( ⟨ρ·dE/dλ⟩ - ⟨ρ⟩⟨dE/dλ⟩ )

where dE/dλ comes from re-analyzing the stored production frames.

MD blocks take total_time (ns), time_step (fs), write_freq (ps), temperature (K); liquid also takes pressure (atm) and ensemble (NPT/NVT). Setting a gas total_time of 0 disables the gas leg.

job:
  type: single-point            # or: optimize
  optimize:
    opt_params:   ["vdw-36 3.4050 0.1100"]
    params_range: ["0.20   0.05"]
    diff_step: 0.0001           # finite-difference step
    ftol: 0.0001                # convergence tolerances
    gtol: 0.0001
    xtol: 0.0001

Each opt_params entry is "<term_key> <value1> <value2> ..." with a matching range string. The term key uses hyphens where the Tinker line has fields, so vdw-36 writes vdw 36 ... and vdwpair-401-402 writes vdwpair 401 402 .... Bounds are value ± range; a range of 0 pins that value — it is written unchanged but hidden from the optimizer.

Optimized parameters are written as override lines appended to a pristine copy of your .prm. Tinker takes the last definition of a term, so the appended line supersedes the original. The result lands in prm/<name>.prm.final.

Runs every configured simulation and reports every derived property, with the deviation from each reference value you supplied. Output goes to results/singlepoint.txt (human-readable), results/singlepoint.yaml (machine-readable), and one results/hfe_<name>.txt per solute with its per-window BAR breakdown.

Fits the shared parameter file to all targets jointly, reusing the single-point machinery to evaluate the objective. Each step logs a table of every target with its current value, reference, difference, and weighted residual.

Systems without a reference value still run — they are simply reported rather than fitted.

Everything a run generates lives under workdir, separated per system so nothing collides:

<workdir>/
├── config.yaml
├── prm/
│   ├── amoeba09.prm            # working copy, rewritten each optimizer step
│   ├── amoeba09.prm.orig       # pristine snapshot; your input is never touched
│   ├── amoeba09.prm_01 …       # perturbation sidecars
│   └── amoeba09.prm.final      # optimized output
├── systems/
│   ├── phenol/{gas,liquid}/    # lambda windows, trajectories, BAR output, FEP_NN/
│   ├── sodium/liquid/
│   └── water_neat/             # one trajectory per temperature
└── results/
    ├── singlepoint.txt / .yaml
    ├── hfe_<name>.txt
    ├── autoff.log
    └── submitted_jobs.txt      # dry-run manifest

Your input files are only ever read. Coordinates are copied into the system directory before minimization, and the parameter file is snapshotted, so rerunning never mutates what you supplied.

autoff run    config.yaml [--dry-run] [-s/--skip-check] [-v N]
autoff setup  config.yaml [--dry-run]
autoff check  config.yaml
autoff report config.yaml [-s/--skip-check]

• run — execute the job named by job.type, end to end.
• setup — generate directories, keys and scripts; submit nothing.
• check — print per-system completeness without changing anything.
• report — collect results from output files already on disk.
• --skip-check — assume .arc/.ene files are complete. Skips the expensive scans; only use it when you are sure.

Jobs are dispatched over SSH by autoff-submit, which polls nproc/top for CPU nodes and nvidia-smi for free GPUs and retries until every job is placed. It can also be used standalone:

autoff-submit -x run_md.sh -t GPU -nodes node103 node104

Submission is fire-and-forget — there is no queue system and no job IDs — so completion is detected by polling output files. All systems are submitted together and advanced by one loop, so a fast system reaches BAR while a slow one is still running dynamics; the wall time is set by the slowest system rather than by their sum.

Run any of them with autoff run config.yaml --dry-run first to inspect what would be submitted.

The dimer geometry and reference energy in Phenol-HFE-Dimer are illustrative placeholders, not QM-derived values. Replace them before drawing conclusions from a fit.

• BAR equilibration. The first 20% of every trajectory is discarded. A window therefore needs at least 6 snapshots, which is validated at load time.
• Density source. Neat-liquid densities are read from the box line of each .arc frame, not from the MD log. The archive is the record of what was simulated, and it is also what the reweighted energies of the density derivative are computed over, so both series come from one file and stay aligned frame for frame. Volume is the full triclinic expression, so non-cubic cells are handled; the run must be NPT, since an NVT archive carries no box line to take a volume from. Scanning a multi-GB archive for its box lines costs a few seconds per temperature per step.
• Production window. The neat-liquid production sample is the last production_time worth of frames, not everything past equil_time. A trajectory rarely lands on exactly the planned length — a resumed run can overshoot — and counting back from the end keeps the sample the intended length and the intended age instead of folding the surplus into the average. The window is floored at equil_time, so a short trajectory averages fewer frames (with a warning) rather than reaching back into equilibration. The same window feeds both the densities and the trimmed -prod.arc the reweighted energies are computed over, so the two cannot drift apart.
• Parameter count. Perturbation sidecars are numbered _01.., and a gradient needs two per free parameter plus one, so at most 49 parameters can be fitted at once. This is checked when the config loads.
• Monomer caching. The dimer_opt target caches relaxed monomer energies once, which assumes the fitted terms do not change intramolecular energy — true for a van der Waals class whose 1-2 interactions are excluded, but not in general.
• Saturated clusters. Submission blocks until every job is placed, so a full cluster stalls the polling loop.
• Resuming. A partial trajectory with a .dyn checkpoint is resumed for only its outstanding steps; a complete one is reused. Stale .bar/.ene files from a run with different sampling are detected and regenerated.
• Dead jobs. Nothing resubmits, so an MD that ends early has to be noticed or the poll loop waits forever. A blow-up is read out of the log or the .err dump; a job killed with its node leaves neither, and is caught only by md_stall_timeout — its .arc and .log simply stop being written. Both raise, naming the temperature or window to restart. The timeout has to sit well above the gap between frames: these runs write one every few seconds, so the hour-long default is far outside the noise, but a much slower protocol should raise it.

pytest

The suite runs without a cluster or a Tinker build: parsers are checked against real BAR output from a completed run, generated scripts are asserted against their exact expected command lines, and the optimizer's residual and Jacobian assembly is exercised with the simulation layer stubbed out.

Chengwen Liu — liuchw2010@gmail.com