Introduction

alchemrs is a Rust-first toolkit for alchemical free energy analysis.

The project is organized as one package containing a library crate and a single CLI binary. Inside the library, the code is split into focused modules for:

parsing AMBER output into typed Rust data structures
preprocessing time series by trimming and decorrelation
estimating free energies with TI, BAR, MBAR, EXP, and DEXP
computing overlap diagnostics
exposing the workflow through a CLI and a clean top-level Rust API

The project is intentionally Rust-native. The alchemrs library crate is the canonical implementation surface, and a CLI is included for convenience. Python bindings will come later.

At the moment, the practical workflow is:

Parse AMBER outputs into DhdlSeries or UNkMatrix.
Trim and decorrelate those data if needed.
Apply an estimator.
Inspect uncertainties and overlap diagnostics.

This can be done from the CLI or the library API.

This book documents the current codebase rather than a future roadmap.

Getting Started

Prerequisites

Rust 1.85 or newer
AMBER output files
- Currently, the parser only support AMBER output files. More are coming.

The repository is a normal Cargo package and can be built and tested with standard Rust tooling.

Build the package

cargo build

Run tests

cargo test

Build the CLI binary

cargo build --release

Run the top-level examples

cargo run --example amber_ti -- 300 ./fixtures/amber/acetamide_tiny/\*/acetamide.prod.out

cargo run --example amber_mbar -- 300 ./fixtures/amber/acetamide_tiny/\*/acetamide.prod.out

Build this book

This documentation uses mdBook.

Install it once:

cargo install mdbook

Preview the docs locally:

mdbook serve docs

Build static HTML:

mdbook build docs

Repository Layout

The repository is centered on a single Cargo package, alchemrs.

That package contains:

the alchemrs library crate
the alchemrs command-line binary

`alchemrs` library layout

The alchemrs crate is organized into modules rather than separate library crates:

data: domain types and scientific data-model structures
error: shared crate-level error handling
parse: engine-specific parsers; today this means AMBER parsing
prep: time-series preparation such as duplicate cleanup, sorting, equilibration detection, and decorrelation
estimators: TI, BAR, MBAR, EXP, and DEXP implementations
analysis: overlap diagnostics built on top of MBAR log weights

The crate root also re-exports the most common types and functions directly for a flatter API.

CLI binary

The alchemrs binary wraps parsing, preprocessing, estimation, and output formatting into a command-line workflow. The scientific logic lives in the library; the CLI is a thin consumer of that API.

Layering

The intended dependency direction inside the library is:

parse -> data + error
prep -> data + error
estimators -> data + error
analysis -> data + error + estimators

The binary depends on the library through the same package source tree.

Data Model

The alchemrs::data module defines the typed structures that the rest of the library operates on.

`StatePoint`

StatePoint stores:

one or more lambda values
temperature in Kelvin

Validation rules:

temperature must be finite and positive
every lambda value must be finite

Current estimators require one-dimensional lambda states. The type itself can store multiple lambda dimensions, but the estimator layer currently rejects multidimensional states.

`DhdlSeries`

DhdlSeries represents a time-indexed dH/dlambda series for one state.

It stores:

the associated StatePoint
time_ps
scalar values

Validation rules:

time_ps.len() == values.len()
all times must be finite
all values must be finite
time must be monotonic nondecreasing

This is the primary input type for TI.

`UNkMatrix`

UNkMatrix stores reduced energies evaluated across states.

Conceptually:

rows are samples
columns are evaluated states
sampled_state identifies the state the trajectory was sampled from
evaluated_states identifies the column states

Validation rules:

data.len() == n_samples * n_states
evaluated_states.len() == n_states
time_ps.len() == n_samples
time must be monotonic nondecreasing
u_nk values may be finite or positive infinity, but not NaN or negative infinity

That last point matters for preprocessing:

some u_nk-derived observables fail on positive infinity
observable-based preprocessing using EPtot can still work if the external observable is finite

`FreeEnergyEstimate`

Represents a scalar free energy result with:

delta_f
optional uncertainty
from_state
to_state

This is used by TI.

`DeltaFMatrix`

Represents pairwise free energy differences between all states.

It stores:

a dense n_states x n_states matrix of values
optional uncertainties
the ordered state list

This is used by BAR, MBAR, EXP, and DEXP.

`OverlapMatrix`

Represents the overlap matrix between states and is used by the diagnostics layer.

Error model

The crate uses the shared CoreError type from alchemrs::error for:

shape mismatches
invalid state metadata
invalid time ordering
non-finite numerical values
parse failures converted into the shared error model
convergence failures
unsupported inputs

This keeps error handling consistent across parsing, preprocessing, estimation, and analysis.

Parsing AMBER Outputs

The current parser implementation is AMBER-specific and lives in alchemrs::parse::amber.

Available entry points

extract_dhdl(path, temperature_k)
extract_u_nk(path, temperature_k)
extract_u_nk_with_potential(path, temperature_k)

`extract_dhdl`

This parser extracts:

temp0
clambda
dt
ntpr
the begin time coordinate
DV/DL values from NSTEP blocks

The extracted gradients are converted to reduced units by multiplying by beta = 1 / (k_B T).

Returned value:

DhdlSeries

Common failure modes:

missing metadata fields such as temp0, clambda, dt, or ntpr
no DV/DL samples found
temperature mismatch between the file and the requested temperature

`extract_u_nk`

This parser reads AMBER MBAR output blocks and constructs a UNkMatrix.

Key behavior:

the parser requires a valid MBAR lambda list in the file header
clambda must be present in the evaluated-state grid
the returned matrix keeps sample rows with positive infinity values
rows that are entirely unusable can still cause parsing to fail if no finite MBAR samples remain

Returned value:

UNkMatrix

`extract_u_nk_with_potential`

This parser does everything extract_u_nk does, and additionally extracts EPtot from NSTEP blocks.

Returned value:

(UNkMatrix, Vec<f64>)

The parser requires the number of EPtot samples to match the number of retained MBAR samples exactly.

This function exists because EPtot is often useful as a finite external observable for:

auto-equilibration
decorrelation

especially when the u_nk matrix contains positive infinity values that make DE-based preprocessing invalid.

Temperature validation

The parser checks the temperature in the AMBER file against the temperature passed by the caller and fails if they differ by more than a small tolerance.

This is deliberate: parser output is temperature-dependent because reduced energies depend on beta.

Current scope

The public parser surface in this repository is AMBER-only. The architecture leaves room for additional engines later, but the implemented and tested parser today is the AMBER parser.

Preprocessing

Preprocessing lives in the alchemrs::prep module and is where most of the scientifically important workflow details live.

Core concepts

The prep crate supports:

duplicate-time cleanup
sorting by time
optional time slicing
equilibration detection
decorrelation / subsampling

The main option struct is DecorrelationOptions.

Default values are:

drop_duplicates = true
sort = true
conservative = true
remove_burnin = false
fast = false
nskip = 1
lower = None
upper = None
step = None

Time cleanup

Before decorrelation or equilibration detection, the prep crate can:

drop duplicate time values
sort by time

This matters because the timeseries logic assumes one contiguous series ordered in time.

For UNkMatrix, cleanup preserves row/state alignment by selecting or reordering entire rows, not individual values.

Time slicing

lower, upper, and step allow simple time-domain slicing before decorrelation.

These options are applied to:

the series itself for DhdlSeries
both the selected scalar observable and the aligned u_nk rows for UNkMatrix

Statistical inefficiency and `g`

Decorrelation is based on an estimate of statistical inefficiency g.

The implementation follows the same basic pymbar.timeseries logic:

center the series
estimate autocorrelation contributions at increasing lag
stop once the correlation function becomes non-positive after a minimum lag
clamp g to at least 1

`fast`

fast changes how g is estimated.

fast = false evaluates lag contributions one lag at a time
fast = true increases the lag increment as it goes, which is cheaper but less accurate

This affects:

plain decorrelation
equilibration detection

because equilibration detection repeatedly estimates g on suffixes of the series.

`conservative`

conservative changes how the code turns g into retained sample indices.

conservative = true uses a uniform stride of ceil(g)
conservative = false uses rounded multiples of fractional g

The conservative mode keeps fewer samples and is intentionally cautious.

Equilibration detection

Equilibration detection returns:

t0: the chosen start index for equilibrated data
g: statistical inefficiency for the retained suffix
neff_max: the estimated number of effectively uncorrelated samples

The heuristic scans possible start positions and chooses the one that maximizes:

Neff = (N - t0) / g

This is the same style of automated equilibration detection used by pymbar, but alchemrs counts the actual retained suffix length N - t0 when reporting Neff_max. That means Neff_max can differ by 1 from pymbar/alchemlyb, while the broader workflow remains the same. For further reading on the topic, see Chodera 2016.

Note: when nskip > 1, alchemrs maximizes Neff only over the sampled candidate time origins (0, nskip, 2*nskip, ...). This intentionally differs from pymbar's current implementation, which pre-fills arrays for all indices before taking the maximum, because the documented meaning of nskip is to try only every nskip-th sample as a potential origin.

`remove_burnin`

In the prep crate, remove_burnin = true means:

run equilibration detection
discard all data before t0
then decorrelate the remaining suffix using the g estimated for that suffix

This is different from the CLI flag --remove-burnin <N>, which is a fixed-count trim done before any automated detection.

`u_nk` observables

u_nk decorrelation needs a scalar time series. The prep crate supports two native derived observables and one external-observable path.

`UNkSeriesMethod::DE`

For each sample row:

identify the sampled-state column
use the next evaluated-state column if it exists
otherwise use the previous column for the last state
compute:

DE_t = u_nk[t, other] - u_nk[t, sampled]

This matches the alchemlyb dE convention.

Important detail:

“adjacent” means adjacent in evaluated-state order, not nearest by numeric lambda distance

`UNkSeriesMethod::All`

For each sample row, sum all evaluated-state reduced energies in that row.

This is a generic matrix-derived scalar, but it is usually less targeted than DE.

`decorrelate_u_nk_with_observable`

This path accepts an external scalar observable, most commonly EPtot.

The observable determines which sample indices are retained, and those retained indices are then applied back to the full u_nk rows.

This is useful when:

the matrix contains infinite energy values that make DE invalid

Non-finite `u_nk` behavior

decorrelate_u_nk rejects non-finite derived scalar series.

In practice:

DE or All can fail if the derived scalar is not finite
decorrelate_u_nk_with_observable can still succeed if the supplied observable itself is finite

That is why the CLI exposes epot as an observable choice instead of only supporting de.

CLI preprocessing order

For CLI commands, preprocessing order is:

fixed-count --remove-burnin <N>
--auto-equilibrate
--decorrelate

TI uses dH/dlambda as the scalar series.

BAR, MBAR, EXP, and DEXP use the observable chosen by --u-nk-observable <de|all|epot>.

CLI auto-equilibrate overrides

When --auto-equilibrate is enabled in the CLI, the effective preprocessing flags become:

fast = true
conservative = false

even if the user supplied different --fast or --conservative values.

This is intentional and mirrors the alchemlyb / pymbar automated-equilibration workflow:

detect_equilibration(..., fast=True)
then subsample_correlated_data(..., conservative=False)

Those effective values are still recorded in CLI provenance so the output remains auditable.

Estimators

Estimator implementations live in the alchemrs::estimators module.

Shared assumptions

For UNkMatrix-based estimators, the current implementation assumes:

one-dimensional lambda states
consistent evaluated-state grids across windows
sampled_state is set and present in evaluated_states
state temperatures are consistent across all windows

If those invariants are violated, the estimators fail early.

TI

Types:

TiEstimator
TiOptions
IntegrationMethod::{Trapezoidal, Simpson}

Input:

slice of DhdlSeries

Behavior:

sorts windows by lambda
computes the mean dH/dlambda for each window
integrates across lambda using trapezoidal or Simpson integration

Current uncertainty behavior:

trapezoidal mode reports an uncertainty estimate
Simpson mode currently returns no uncertainty

BAR

Types:

BarEstimator
BarOptions
BarMethod::{FalsePosition, SelfConsistentIteration, Bisection}

Input:

slice of UNkMatrix

Behavior:

expects windows mapped to sampled states on one evaluated-state grid
computes work values between adjacent states
fills a full pairwise DeltaFMatrix

Important limitation:

uncertainty is only computed for adjacent windows
non-adjacent uncertainties are reported as NaN

This matches the current CLI documentation and the reference-comparison tests in the repository.

MBAR

Types:

MbarEstimator
MbarOptions

Important options:

max_iterations
tolerance
initial_f_k
compute_uncertainty
parallel

Input:

slice of UNkMatrix

Behavior:

combines windows into a shared reduced-energy representation
solves for free energies iteratively
returns a full DeltaFMatrix

The analysis crate uses MBAR-derived log weights to compute overlap diagnostics.

EXP and DEXP

Types:

ExpEstimator
ExpOptions

Input:

slice of UNkMatrix

Behavior:

computes exponential averaging from each sampled-state window to all evaluated states
returns a full pairwise DeltaFMatrix

CLI direction conventions:

exp reports the forward direction
dexp reports the reverse direction

DEXP is not a separate estimator type in the library; it is a CLI convention over the same estimator results.

Parallel execution

Some estimators support a parallel option and use Rayon-backed parallel loops internally.

Parallelism affects performance, not the conceptual API.

Overlap Diagnostics

Overlap diagnostics live in the alchemrs::analysis module.

Available functions:

overlap_matrix
overlap_eigenvalues
overlap_scalar

`overlap_matrix`

This function:

computes MBAR log weights from the supplied windows
constructs the overlap matrix from those weights and state sample counts
returns an OverlapMatrix

Input:

slice of UNkMatrix
optional MbarOptions

`overlap_eigenvalues`

Computes the eigenvalues of the overlap matrix and returns them sorted in descending order.

The implementation rejects eigenvalues with significant imaginary components.

`overlap_scalar`

Computes:

1 - second_largest_eigenvalue

This is the scalar overlap diagnostic currently surfaced by the CLI.

CLI usage

bar, mbar, exp, and dexp can include overlap diagnostics with:

--overlap-summary

The CLI then reports:

overlap scalar
overlap eigenvalues

in text, JSON, or CSV output.

CLI Guide

The CLI is the alchemrs binary in the same package as the alchemrs library crate.

Commands:

ti
bar
mbar
exp
dexp

Build

cargo build --release

Common workflow

Typical CLI usage is:

pass one AMBER output per lambda window
optionally trim burn-in
optionally run auto-equilibration
optionally decorrelate
fit the requested estimator
emit text, JSON, or CSV

Shared flags

--temperature <K>
--remove-burnin <N>
--auto-equilibrate
--decorrelate
--fast
--conservative[=true|false]
--nskip <N>
--output-units <kt|kcal|kj>
--output-format <text|json|csv>
--output <PATH>
--parallel

For bar, mbar, exp, and dexp:

--u-nk-observable <de|all|epot>

For overlap-aware commands:

--overlap-summary

Observable selection

For u_nk-based estimators:

de is the default and matches the alchemlyb-style adjacent-state dE observable
all sums the full u_nk row
epot uses EPtot parsed from the AMBER output

Use epot when:

you want the CLI's external-observable path
the u_nk matrix contains positive infinity values that make de invalid

TI-specific behavior

TI uses dH/dlambda, not u_nk.

The CLI accepts --u-nk-observable on ti only to provide a more helpful error message. If supplied, the command fails with a domain-specific explanation instead of a generic unknown-flag parse error.

Command examples

TI

cargo run --release -- ti \
  --temperature 300 \
  --method trapezoidal \
  --decorrelate \
  path/to/*/prod.out

BAR

cargo run --release -- bar \
  --temperature 300 \
  --decorrelate \
  --u-nk-observable de \
  --overlap-summary \
  path/to/*/prod.out

MBAR with EPtot fallback

cargo run --release -- mbar \
  --temperature 300 \
  --auto-equilibrate \
  --decorrelate \
  --u-nk-observable epot \
  --output-format json \
  path/to/*/prod.out

EXP / DEXP

cargo run --release -- exp --temperature 300 path/to/*/prod.out

cargo run --release -- dexp --temperature 300 path/to/*/prod.out

Effective settings

When --auto-equilibrate is enabled, the CLI reports the effective preprocessing policy in provenance:

fast = true
conservative = false

That override is deliberate and mirrors the alchemlyb / pymbar equilibration-detection workflow.

Outputs and Provenance

The CLI supports three output formats:

text
JSON
CSV

Units

Output units can be selected with:

--output-units <kt|kcal|kj>

The CLI converts the reported free energy and uncertainty from reduced units into the selected output unit using the requested temperature.

Text output

Text output is intended for direct reading in the terminal.

It includes:

delta_f
uncertainty
lambda endpoints
sample counts
u_nk_observable when relevant
overlap summary when requested

JSON output

JSON output is suitable for shell pipelines and downstream tooling.

The payload contains:

the scalar result fields
optional overlap summary
a provenance object

CSV output

CSV output appends provenance fields after the result columns so the row remains self-describing.

Provenance

In this project, “provenance” means the analysis settings and sample-count metadata that explain how a reported result was produced.

Current provenance fields include:

estimator name
temperature
whether decorrelation ran
fixed-count burn-in trimming
whether auto-equilibration ran
the effective fast value
the effective conservative value
nskip
u_nk_observable when applicable
number of windows
number of samples before preprocessing
number of samples after burn-in trimming / auto-equilibration
number of retained samples after decorrelation

This matters because preprocessing choices directly change the data used by the estimators.

Example JSON payload

{
  "delta_f": -113.58,
  "uncertainty": 1.17,
  "from_lambda": 0.0,
  "to_lambda": 1.0,
  "units": "kT",
  "overlap": {
    "scalar": 0.0183,
    "eigenvalues": [1.0, 0.9817]
  },
  "provenance": {
    "estimator": "mbar",
    "temperature_k": 300.0,
    "decorrelate": true,
    "remove_burnin": 0,
    "auto_equilibrate": false,
    "fast": false,
    "conservative": true,
    "nskip": 1,
    "u_nk_observable": "de",
    "windows": 15,
    "samples_in": 300,
    "samples_after_burnin": 300,
    "samples_kept": 126
  }
}

Library Guide

The top-level alchemrs crate re-exports the main types and functions needed for direct Rust use.

Common imports

#![allow(unused)]
fn main() {
use alchemrs::{
    DecorrelationOptions, MbarEstimator, MbarOptions, TiEstimator, TiOptions,
    decorrelate_dhdl, decorrelate_u_nk_with_observable, extract_dhdl,
    extract_u_nk_with_potential, overlap_matrix, overlap_scalar,
};
}

TI example

#![allow(unused)]
fn main() {
use alchemrs::estimators::{TiEstimator, TiOptions};
use alchemrs::parse::amber::extract_dhdl;

fn run_ti(paths: &[String], temperature_k: f64) -> Result<(), Box<dyn std::error::Error>> {
    let mut series = Vec::new();
    for path in paths {
        series.push(extract_dhdl(path, temperature_k)?);
    }

    let estimator = TiEstimator::new(TiOptions::default());
    let result = estimator.fit(&series)?;
    println!("TI dG = {:.6}", result.delta_f());
    Ok(())
}
}

MBAR example with `EPtot` decorrelation

#![allow(unused)]
fn main() {
use alchemrs::{
    DecorrelationOptions, MbarEstimator, MbarOptions, UNkMatrix,
    decorrelate_u_nk_with_observable, extract_u_nk_with_potential,
    overlap_matrix, overlap_scalar,
};

fn load_windows(
    paths: &[String],
    temperature_k: f64,
) -> Result<Vec<UNkMatrix>, Box<dyn std::error::Error>> {
    let mut windows = Vec::new();
    for path in paths {
        let (u_nk, epot) = extract_u_nk_with_potential(path, temperature_k)?;
        let decorrelated =
            decorrelate_u_nk_with_observable(&u_nk, &epot, &DecorrelationOptions::default())?;
        windows.push(decorrelated);
    }
    Ok(windows)
}

fn run_mbar(windows: &[UNkMatrix]) -> Result<(), Box<dyn std::error::Error>> {
    let estimator = MbarEstimator::new(MbarOptions::default());
    let result = estimator.fit(windows)?;
    let delta_index = result.n_states() - 1;
    println!("MBAR dG = {:.6}", result.values()[delta_index]);

    let overlap = overlap_matrix(windows, Some(MbarOptions::default()))?;
    println!("overlap = {:.6}", overlap_scalar(&overlap)?);
    Ok(())
}
}

Choosing the right preprocessing entry point

Use decorrelate_dhdl for TI inputs.
Use decorrelate_u_nk(..., UNkSeriesMethod::DE, ...) when you want an internal u_nk-derived scalar series.
Use decorrelate_u_nk_with_observable when you already have a finite external observable such as EPtot.

Top-level re-exports

The root crate re-exports:

analysis
data
error
estimators
parse
prep

You can use either:

the direct re-exported functions and types
the crate namespaces for more explicit imports

depending on your style preference.

Compatibility and Limitations

This section describes the current scope of the codebase and where it intentionally matches external tools.

Current compatibility targets

The repository includes fixtures and tests intended to compare behavior against:

alchemlyb
pymbar

Important matches include:

u_nk DE decorrelation semantics
pymbar-style statistical inefficiency estimation
pymbar-style automated equilibration detection
BAR reporting NaN uncertainties for non-adjacent pairs

Important intentional difference:

alchemrs reports equilibration Neff_max using the actual suffix length N - t0, so this scalar can differ by 1 from pymbar / alchemlyb

de and all observable choices can fail during decorrelation if the derived scalar series is non-finite
epot can remain usable if the external observable is finite

Documentation scope

This book documents the current implementation, not a frozen public-stability promise. The project is still early enough that:

APIs may still evolve
defaults may still be refined
more formal hosted documentation infrastructure may still change

What is not covered here

This repository does not currently provide:

plotting in the Rust library
production-grade multi-engine parser support
Python bindings in the current package
versioned public API guarantees

Those may come later, but they are not part of the implemented surface documented in this book.

Keyboard shortcuts

alchemrs