Analysis

Analysis objects sample the system at a given frequency during a simulation and optionally write results to output files. They are defined in the analysis section of the YAML input:

analysis:
  - !Energy
    file: energy.csv.gz
    frequency: !Every 100
  - !RadialDistribution
    selections: ["atomtype Na", "atomtype Cl"]
    file: rdf.dat
    resolution: 0.1
    frequency: !End

Sampling frequency

Each analysis requires a frequency field controlling when it is evaluated.

Frequency	Description
!Every N	Every N steps
!Once N	Once at step N
!End	Once after the last step

In output.yaml every analysis reports num_samples, the number of frames it sampled. Analyses that loop over molecules also report how many molecules they accumulated: num_groups_sampled for Polymer Shape, num_blocks for Widom Rotation. Those counts, not num_samples, sit behind the reported means and error bars.

Output file formats

Most analysis files are column data. For those files, the format is inferred from the file extension:

Extension	Format
.dat	Space-separated; header prefixed with #
.csv/.tsv	Comma/Tab-separated; plain header row (column labels)

Both formats support transparent gzip compression by appending .gz (e.g. energy.dat.gz).

Grid-based analyses document their own file format.

Loading CSV output in Python:

import numpy as np
data = np.loadtxt("data.csv.gz", delimiter=",", skiprows=1)

import pandas as pd
data = pd.read_csv("data.csv.gz").to_numpy()

---

Collective Variable

Monitors a collective variable (CV) over the course of a simulation, recording the instantaneous value and a running average at each sampled step. This is useful for verifying that a CV remains within the expected range, for convergence checking, and for post-processing of time/step series data.

If file is given, each sampled step writes a line with columns step, value, and running_average. If no file is given, only the mean and RMS are written to output.yaml.

Example

analysis:
  - !CollectiveVariable
    property: mass_center_position
    selection: "molecule protein"
    projection: z
    file: cv.dat
    frequency: !Every 100

Options

Key	Required	Default	Description
property	yes		CV type (see table below)
frequency	yes		Sample frequency, e.g. !Every 100
projection	no	xyz	Axis projection (x, y, z, xy, …); alias: dimension
selection	depends		Selection expression for one atom or group
selection2	depends		Second selection (for two-group properties)
resolution	no		Bin width (only used by Penalty)
file	no		Output file path (see Output file formats); omit to only track the mean

Supported properties

Property	Selection	Description
volume	none	Cell measure via dimension: volume (xyz), area (xy), or length (z). Note: volume uses dimension, not projection
atom_position	one atom	Signed component for single axis (x,y,z); Euclidean norm for multi-axis (xy etc.). Selection resolved live each evaluation — works with speciation/GCMC where the matching atom changes. Returns NaN if selection matches ≠ 1 atom. Use atomtype (not name) for atoms defined via explicit atoms: lists
count	atoms or groups	Molecule instances for Molecular groups; atom count for Atomic/Reservoir groups
molarity	atoms or groups	Molar concentration (mol/L); molecule-based for Molecular groups, atom-based for Atomic/Reservoir
charge	atoms or groups	Sum of charges of active atoms matching selection
size	one group	Number of active particles in a group
end_to_end	one group	End-to-end distance of a molecular group
gyration_radius	one group	Radius of gyration (default xyz = full Rg; single axis gives spread along it)
dipole_moment	one group	Electric dipole moment (default xyz = magnitude; single axis gives signed component)
mass_center_position	one group	Mass center position along projection
mass_center_separation	two groups	Distance between two group mass centers
dipole_product	two groups	Normalized dipole dot product μ̂₁·μ̂₂ = cos(θ) (default xyz = full 3D; projection filters dipoles before comparing)

---

Polymer Shape

Computes size and shape anisotropy descriptors from the mass-weighted gyration tensor (doi:10/d6ff). For each matching group the 3×3 symmetric gyration tensor is built as S_ij = (1/M) Σ m_k r_i r_j, where r is the PBC-aware displacement from the center of mass. Eigenvalue decomposition yields the principal moments λ₁ ≤ λ₂ ≤ λ₃ from which the following descriptors are derived.

Following IUPAC 2014 nomenclature, the per-configuration squared radius of gyration is s² = Tr(S) (Def. 1.14), and the reported Rg = √⟨s²⟩ is the root-mean-square radius of gyration (Def. 1.16). Note that ⟨s⟩ ≠ √⟨s²⟩ in general. The reported ratio ⟨Re²⟩/⟨Rg²⟩ equals 6 for ideal Gaussian chains (Wall & Erpenbeck 1959).

Descriptor	Formula	Reference
Radius of gyration squared	Rg² = λ₁ + λ₂ + λ₃	IUPAC 2014 Def. 1.14
Asphericity	b = λ₃ − (λ₁ + λ₂)/2	Aronovitz & Nelson 1986
Acylindricity	c = λ₂ − λ₁	Aronovitz & Nelson 1986
Relative shape anisotropy	κ² = (b² + ¾c²) / Rg⁴ ∈ [0, 1]	Rudnick & Gaspari 1986
Prolateness	S = 27(λ₁−λ̄)(λ₂−λ̄)(λ₃−λ̄) / Rg⁶ ∈ [−0.25, 2]	Theodorou & Suter 1985
Westin linear (rod-like)	Cl = (λ₃ − λ₂) / Rg²	Westin 2002
Westin planar (disc-like)	Cp = 2(λ₂ − λ₁) / Rg²	Westin 2002
Westin spherical	Cs = 3λ₁ / Rg² (Cl + Cp + Cs = 1)	Westin 2002

Mean gyration tensor components (Sxx, Sxy, Sxz, Syy, Syz, Szz) are reported in the YAML output. The ratio ⟨Szz⟩ / ⟨Rg²⟩ indicates orientation relative to the z-axis: ~0 for flat objects lying in the xy-plane, ~1/3 for isotropic, ~1 for rod-like objects aligned with z.

If file is given (single-molecule selection only), each sampled step writes a line with columns step, Rg, and the upper triangle of the gyration tensor (Sxx Sxy Sxz Syy Syz Szz).

Example

analysis:
  - !PolymerShape
    selection: "molecule polymer"
    file: shape.csv.gz
    frequency: !Every 100

Options

Key	Required	Default	Description
selection	yes		Selection expression for molecule group(s)
file	no		Streaming output file, single molecule only (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 100

---

Multipole

Per-group charge, dipole, and quadrupole analysis, averaged over all groups matching a selection. Useful for tracking titration state and charge fluctuations.

Reports:

• net charge \(\langle Z\rangle \pm \sigma\) — mean and standard deviation
• capacitance \(C = \langle Z^2\rangle - \langle Z\rangle^2\) — charge variance
• dipole magnitude \(\langle|\boldsymbol{\mu}|\rangle \pm \sigma\) (eÅ)
• quadrupole magnitude \(\langle\|\boldsymbol{\Theta}\|\rangle \pm \sigma\) — Frobenius norm (eŲ)
• quadrupole tensor — mean and standard deviation \(\sigma\) of the six independent components
• per-atom ⟨q⟩ and ⟨q²⟩-⟨q⟩² — for atoms with fluctuating charge (e.g. from titration or atom swaps)

The dipole and quadrupole moments are computed relative to each group's center of mass with periodic boundary conditions applied. Handles atom-type swaps (titration) and GCMC (only active groups contribute).

The reported quadrupole is the traceless form (Buckingham convention)

\[\Theta_{\alpha\beta} = \tfrac{1}{2}\sum_i q_i\left(3 d_{i\alpha} d_{i\beta} - d_i^2\,\delta_{\alpha\beta}\right),\]

where \(\boldsymbol{d}_i\) points from the group's center of mass to atom \(i\). The isotropic part dropped here does not affect the far-field potential, and removing it makes the moment vanish for isotropic charge distributions, so the values compare directly with experiment and literature.

Selection resolves group-wise

The selection picks whole groups that contain at least one active matching atom — not just the matched atoms themselves. So an atom-level selection like atomtype CA pulls in every group that has an active CA atom and accumulates the full group's charge and dipole moment. This is useful when a molecule changes name during speciation, for example across protonation states, but still contains a stable marker atom:

- !Multipole
  selection: "atomtype INO"   # picks every phytate, regardless of protonation state
  frequency: !Every 10

By contrast, !CollectiveVariable property: charge sums the charge of only the matched atoms, not their parent groups — see Supported properties.

Per-atom charge statistics are reported in output.yaml as an atoms list only when the selection resolves to a single molecular molecule kind. Selections matching multiple molecule kinds are rejected, and atomic/reservoir groups do not emit per-atom output. The list includes only atoms whose charge variance is nonzero:

multipole:
  selection: molecule MOL1
  num_samples: 2000
  charge: {mean: 4.7785, error: 2.2373}
  capacitance: 5.005
  dipole_moment: {mean: 141.6532, error: 23.8816}
  quadrupole_moment: {mean: 312.4187, error: 45.1902}
  quadrupole_tensor:            # order [xx, xy, xz, yy, yz, zz], loads straight into NumPy
    order: [xx, xy, xz, yy, yz, zz]
    values: [180.21, 12.03, -4.11, -95.44, 8.90, -84.77]
    errors: [22.14, 3.02, 1.88, 15.30, 2.71, 11.05]
  atoms:
    - {index: 1, name: NP, ⟨q⟩: -0.52, ⟨q²⟩-⟨q⟩²: 0.2496}
    - {index: 6, name: NP, ⟨q⟩: -0.48, ⟨q²⟩-⟨q⟩²: 0.2496}

Example

analysis:
  - !Multipole
    selection: "molecule peptide"
    frequency: !Every 10

Options

Key	Required	Default	Description
selection	yes		Selection expression for molecule group(s)
frequency	yes		Sample frequency

---

Multipole Distribution

This analysis decomposes the electrostatic interaction between two molecular species into ion–ion, ion–dipole, dipole–dipole, and ion–quadrupole terms, each averaged as a function of the center-of-mass separation R. Comparing their sum with the exact Coulomb energy shows how far a multipole expansion reproduces the true interaction, and how the molecules' mutual orientation changes with distance. Each species may contain many molecules: every frame the analysis averages over all molecule pairs — all pairs of distinct molecules when the two selections name the same species, or all cross pairs when they name different species — so it characterises a whole solution of N macroions, not just an isolated pair.

Each molecule is reduced to a net charge, a dipole moment, and a traceless quadrupole moment about its center of mass, with periodic boundaries applied. The pair energy at separation \(\boldsymbol{R} = \boldsymbol{r}_a - \boldsymbol{r}_b\) is

\[ u = \frac{q_a q_b}{R} + \frac{q_a(\boldsymbol{\mu}_b\cdot\boldsymbol{R}) - q_b(\boldsymbol{\mu}_a\cdot\boldsymbol{R})}{R^3} + \frac{\boldsymbol{\mu}_a\cdot\boldsymbol{\mu}_b}{R^3} - \frac{3(\boldsymbol{\mu}_a\cdot\boldsymbol{R})(\boldsymbol{\mu}_b\cdot\boldsymbol{R})}{R^5} + \frac{q_a\,\boldsymbol{R}^{\!\top}\boldsymbol{\Theta}_b\,\boldsymbol{R} + q_b\,\boldsymbol{R}^{\!\top}\boldsymbol{\Theta}_a\,\boldsymbol{R}}{R^5} \]

Each term is thermally averaged over orientations and co-solute configurations. Energies are reported in kJ/mol: the geometric factors above carry the Bjerrum length and RT of the medium. The exact energy is the explicit sum over atomic charges, \(u_\text{exact} = \sum_i^a\sum_j^b q_i q_j / |\boldsymbol{r}_i-\boldsymbol{r}_j|\).

Alongside the energies, five orientational measures are accumulated (column name in parentheses):

• dipole correlation \(\langle\hat{\boldsymbol{\mu}}_a\cdot\hat{\boldsymbol{\mu}}_b\rangle\) (mucorr), between \(-1\) and \(1\) — mean cosine between the dipole directions: \(+1\) parallel, \(-1\) antiparallel, \(0\) uncorrelated.
• nematic order \(\langle P_2\rangle = \langle(3\cos^2\theta - 1)/2\rangle\) (p2), between \(-0.5\) and \(1\) — alignment blind to dipole sign: \(+1\) collinear (parallel or antiparallel), \(0\) random, \(-0.5\) mutually perpendicular.
• longitudinal projection \(\langle(\hat{\boldsymbol{\mu}}_a\cdot\hat{\boldsymbol{R}})(\hat{\boldsymbol{\mu}}_b\cdot\hat{\boldsymbol{R}})\rangle\) (long), between \(-1\) and \(1\) — orientation relative to the separation axis: near \(+1\) both point along \(\boldsymbol{R}\) (head-to-tail, end-on), near \(0\) both lie perpendicular to it (side-by-side, broadside).
• Kirkwood factor \(g_K(R) = 1 + \frac{1}{N}\sum_{i \neq j,\, r_{ij} \leq R}\langle\hat{\boldsymbol{\mu}}_i\cdot\hat{\boldsymbol{\mu}}_j\rangle\) (gK) — the running dipole-correlation integral, where \(N\) counts the reference molecules that carry a dipole (apolar molecules are excluded from the normalization): \(g_K > 1\) signals net parallel ordering of the surrounding dipoles, \(g_K < 1\) net antiparallel ordering. It applies when both selections are the same species; for two different species the bare cumulative sum is reported, without the \(+1\).
• quadrupole correlation — the tensor overlap \(\langle\boldsymbol{\Theta}_a : \boldsymbol{\Theta}_b\rangle\) (quadcorr) and its normalized form \(\langle\hat{\boldsymbol{\Theta}}_a : \hat{\boldsymbol{\Theta}}_b\rangle\) (quadcorr_norm), between \(-1\) and \(1\): \(+1\) identically oriented quadrupoles, \(-1\) maximally anti-aligned, \(0\) uncorrelated.

Both selections must resolve to molecules with a center of mass. Results are written to a CSV file, one row per distance bin, with columns R,exact,tot,ii,id,dd,iq,mucorr,p2,long,gK,quadcorr,quadcorr_norm.

Example

analysis:
  - !MultipoleDistribution
    selections: ["molecule protein", "molecule protein"]
    resolution: 2.0
    frequency: !Every 100

The output defaults to multipole_dist.csv; set file to override it.

Options

Key	Required	Default	Description
selections	yes		Two selection expressions, a and b
file	no	multipole_dist.csv	Output CSV file
resolution	no	1.0	Distance resolution along R (Å)
max_r	no	half box	Maximum separation (Å)
frequency	yes		Sample frequency

A medium is required for the Bjerrum length.

---

Energy

Streams energy values to a file at each sampled step. Two modes are supported:

• Total (default): writes every Hamiltonian term plus the total. Output columns: step term1 term2 ... total.
• Partial: writes the nonbonded energy between two sets of atoms selected with VMD-like expressions. Output columns: step energy running_average.

Examples

analysis:
  # Total energy with per-term breakdown
  - !Energy
    file: energy.csv.gz
    frequency: !Every 100

  # Nonbonded energy between two molecules
  - !Energy
    file: mol1_mol2_energy.csv.gz
    frequency: !Every 100
    selections: ["molecule MOL1", "molecule MOL2"]

  # Nonbonded energy between hydrophobic atoms in two molecules
  - !Energy
    file: hydrophobic_energy.csv.gz
    frequency: !Every 100
    selections: ["hydrophobic and molecule MOL1", "hydrophobic and molecule MOL2"]

Options

Key	Required	Default	Description
file	yes		Output file path (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 100
selections	no		Pair of selection expressions for partial nonbonded energy

When selections is omitted, total mode is used. When given, only nonbonded energy terms contribute; other terms (bonded, constraints, etc.) are skipped. Selections resolve at the atom level, so atom-type filters (e.g. hydrophobic and molecule MOL1) work correctly. When both selections resolve to the same atoms, self-pairs and duplicates are automatically excluded.

---

Trajectory

Writes particle coordinates to a trajectory file at each sampled step. The format is auto-detected from the file extension (.xyz or .xtc).

On finalization a companion PSF topology file and a VMD scene script are written alongside the trajectory. The PSF contains atoms, bonds, angles, and dihedrals; the Tcl script loads the PSF and trajectory, sets VDW radii from sigma values, colors by charge, and draws the periodic box.

For speciation/GCMC/titration simulations, two additional companion files enable per-frame visualization updates in VMD:

• .sizes.dat — per-frame group active counts (written when any group has inactive atoms). The VMD script hides inactive atoms (radius = 0).
• .charges.dat — per-frame atom charges (only when save_charges: true). The VMD script then updates charges each frame so that atom-type swaps from titration and speciation are reflected in the charge coloring. Without it, VMD colors by the static charges stored in the PSF.

Example

analysis:
  - !Trajectory
    file: traj.xtc
    frequency: !Every 100
    save_frame_state: true
  - !Trajectory
    file: protein.xyz
    frequency: !Every 100
    selection: "protein"

After the simulation, visualize with:

vmd -e traj.tcl

Options

Key	Required	Default	Description
file	yes		Output file path (.xyz, .xtc, etc.)
frequency	yes		Sample frequency, e.g. !Every 100 or !End
selection	no		VMD-like molecule selection, e.g. "molecule protein". Writes only matching groups.
save_frame_state	no	false	Write a binary .aux file alongside the trajectory (see Rerun). Cannot be combined with selection.
save_charges	no	false	Write a per-frame .charges.dat file for VMD charge coloring of titration/speciation swaps. One float per atom per frame, so large for big systems.

Output files

Given file: traj.xtc, the following files are produced:

File	Description
traj.xtc	Trajectory (coordinates per frame)
traj.psf	X-PLOR PSF topology (atoms, bonds, angles, dihedrals, charges, masses)
traj.tcl	VMD scene script (vmd -e traj.tcl loads everything)
traj.charges.dat	Per-frame atom charges (only when save_charges: true)
traj.sizes.dat	Per-frame group active counts (only when groups have inactive atoms)
traj.aux	Frame state file (only when save_frame_state: true)

Frame state file (.aux)

XTC stores only atom positions. Rigid-body simulations also need quaternions, mass centers, and group sizes (for GC moves); swap moves change atom_id. When save_frame_state: true, a binary .aux file is written alongside each XTC frame, storing this per-frame microstate data.

The .aux file is required by faunus rerun (see Rerun) to fully reconstruct the simulation state from each trajectory frame.

---

Radial Distribution Function

Computes the radial distribution function g(r) for pairs of particles or molecule centers of mass. Distances use the minimum image convention for periodic boundary conditions. The output file is updated at each sample with columns r and g(r), normalized using the average volume (NPT-compatible).

Two modes are supported:

• Atom-atom (default): pairwise distances between individual atoms. Intramolecular pairs (atoms in the same molecule) are excluded by default.
• COM-COM (use_com: true): pairwise distances between molecular centers of mass.

The dimension option controls which spatial components are used for the pair distance and shell normalization:

• xyz (default): full 3D distance, normalized with spherical shells (4πr²).
• xy, xz, yz: projected distance in a plane, normalized with cylindrical shells (2πr).
• x, y, z: projected distance along an axis, normalized with linear shells.

This is useful for systems where particles are constrained to a plane or a line.

When both selections are identical, self-pairs and duplicates are automatically excluded.

Examples

analysis:
  # Atom-atom RDF between Na and Cl
  - !RadialDistribution
    selections: ["atomtype Na", "atomtype Cl"]
    file: rdf_nacl.dat
    resolution: 0.1
    frequency: !Every 100

  # COM-COM RDF for polymer molecules
  - !RadialDistribution
    selections: ["molecule polymer", "molecule polymer"]
    use_com: true
    file: rdf_com.dat
    resolution: 0.5
    max_r: 30.0
    frequency: !Every 100

Options

Key	Required	Default	Description
selections	yes		Pair of selection expressions, e.g. ["atomtype Na", "atomtype Cl"]
file	yes		Output file path (see Output file formats)
resolution	yes		Bin width in distance units
frequency	yes		Sample frequency, e.g. !Every 100
max_r	no	half shortest box dim	Maximum distance for histogram
use_com	no	false	Use center-of-mass distances instead of atom-atom
exclude_intramolecular	no	true (atom-atom)	Skip pairs within the same molecule (atom-atom only)
dimension	no	xyz	Axes for distance projection and normalization (x, y, z, xy, …)

---

Spatial Distribution Function

Computes the spatial distribution function (SDF) of an atom selection around a rigid molecular reference selection. The result is a scalar grid in the reference molecule body frame, averaged over all matching reference molecules and all sampled frames.

This is useful for mapping ions or solvent around rigid macromolecules in a periodic box. Target positions are measured from each reference center using minimum-image periodic displacements and then rotated by the inverse reference quaternion before binning.

The reference selection must match rigid molecular groups. Flexible reference molecules are not supported by this analysis.

Example

analysis:
  - !SpatialDistribution
    reference: "molecule Macro"
    selection: "atomtype Na"
    frequency: !Every 100

Options

Key	Required	Default	Description
reference	yes		Molecular group selection defining the body-fixed frame
selection	yes		Atom selection accumulated on the grid
frequency	yes		Sample frequency, e.g. !Every 100
file	no	spatial.dx	OpenDX output grid
reference_file	no		Optional XYZ structure of the reference molecule, for visualizing the density around it
resolution	no	1.0	Cubic grid spacing in Å
padding	no	8.0	Margin in Å added on every side of the reference molecule's bounding box
bulk_normalize	no	true	Normalize by bulk density to produce dimensionless relative density
exclude_reference	no	true	Skip target atoms belonging to the current reference group

The grid is the axis-aligned bounding box of the active reference molecule(s) in the body frame, rounded to the grid spacing and enlarged by padding on every side. Set reference_file to also write the reference molecule itself (in the same body frame) so it can be overlaid on the density. Output is written once at the end of the run.

Normalization

With the default bulk_normalize: true, grid values are normalized by the instantaneous bulk density of the target selection. A homogeneous ideal gas therefore gives SDF ≈ 1. The normalization uses the instantaneous target count, reference count, and cell volume at each sample, so it works with GCMC where particle numbers fluctuate.

Set bulk_normalize: false to write molar concentration in mol/L instead of relative density.

Output

The output file is an OpenDX scalar grid suitable for visualization in VMD and PyMOL. Grid coordinates are in the reference body frame, not the lab frame.

---

Widom Insertion

Measures the excess chemical potential of a single ion species using the scaled Widom method (Svensson & Woodward, 1988). A ghost particle is inserted at random positions; charge scaling maintains electroneutrality in the finite periodic box, correcting the Coulombic size error inherent in naive single-ion Widom insertion.

The excess chemical potential is decomposed into short-range and electrostatic contributions. The electrostatic part is evaluated by numerical integration over a charging parameter λ ∈ [0, 1]. Results are block-averaged with standard error of the mean reported.

Pair interactions for the ghost particle are defined directly in the analysis block using the same syntax as energy.nonbonded.default, allowing arbitrary short-range potentials (WCA, LJ, HardSphere, etc.) beyond the primitive model.

Example

analysis:
  - !ScaledWidomInsertion
    atom: Na
    insertions: 20
    lambda_points: 11
    frequency: !Every 10
    default:
      - !Coulomb { cutoff: 1000.0 }
      - !WCA { mixing: arithmetic }

Options

Key	Required	Default	Description
atom	yes		Ghost atom type (must exist in topology)
insertions	no	10	Number of ghost insertions per sample
lambda_points	no	11	Quadrature points for charge scaling (odd recommended)
default	yes		Pair interactions (same syntax as nonbonded.default)
frequency	yes		Sample frequency, e.g. !Every 10

---

Virtual Translate

Performs a virtual displacement of a single molecule and measures the mean force by Widom perturbation (Widom, 1963):

\[f = \frac{k_BT \ln\langle e^{-\Delta U / k_BT}\rangle}{\delta L}\]

where \(\Delta U\) is the energy change due to the displacement \(\delta L\). If file is given, each sampled step writes columns step, dL, dU/kT, and <force>/kT/Å.

Example

analysis:
  - !VirtualTranslate
    selection: "molecule protein"
    dL: 0.01
    directions: !z
    file: force.dat
    frequency: !Every 10

Options

Key	Required	Default	Description
selection	yes		Selection matching exactly one molecule
dL	yes		Displacement magnitude (Å)
directions	no	z	Displacement direction (x, y, z, xy, …)
file	no		Output file path (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 10

---

Virtual Volume Move

Performs a virtual volume perturbation and measures the excess pressure by Widom perturbation (doi:10.1063/1.472721):

\[P_{\text{ex}} = \frac{k_BT \ln\langle e^{-\Delta U / k_BT}\rangle}{\delta V}\]

where \(\Delta U\) is the energy change due to the volume displacement \(\delta V\). All particle positions are scaled according to the chosen method. If file is given, each sampled step writes columns step, dV/ų, dU/kT, and <Pex>/kT/ų.

Example

analysis:
  - !VirtualVolumeMove
    volume_displacement: 0.2
    method: Isotropic
    file: pressure.csv
    frequency: !Every 10

Options

Key	Required	Default	Description
volume_displacement	yes		Volume displacement (ų)
method	no	Isotropic	Scaling policy: Isotropic, ScaleZ, ScaleXY
file	no		Output file path (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 10

---

Double Layer Pressure

Osmotic (disjoining) pressure between two uniformly charged planes with explicit point-charge counterions, by the midplane method of Guldbrand et al. (doi:10.1063/1.446912):

\[P_\text{osm} = k_BT\sum_i C_i(0) \;+\; \frac{F_z^{AB}}{\text{area}}\]

The first term is the entropic (repulsive) contribution from the ion concentration at the midplane; the second is the cross-midplane electrostatic force. Its ion–ion correlation part is the attractive, van-der-Waals-like force that dominates for divalent ions and that mean-field Poisson–Boltzmann misses. The surface charge density \(\sigma\) is set by electroneutrality from the counterion charges — you do not specify it.

The analysis is hardcoded to a Slit cell (periodic in xy, walls at \(z=\pm L_z/2\), so the midplane is \(z=0\)) and to electrostatics only; it suits the primitive-model double layer with point-charge counterions. Valency mixtures need nothing special — select all the counterions (e.g. "atomtype Na Ca") and report. Results are block-averaged and reported as mean ± standard error.

Example

analysis:
  - !DoubleLayerPressure
    selection: "atomtype Na Ca"   # mono- + divalent counterions
    file: pressure.csv
    frequency: !Every 10

Options

Key	Required	Default	Description
selection	yes		Mobile counterions; sets the electroneutral surface charge \(\sigma\) (the Coulomb prefactor comes from the medium)
midplane_halfwidth	no	gap/12	Half-width (Å) of the midplane density sampling window. The default scales with the plate separation and is usually fine; treat it as a convergence parameter
density_bins	no	50	Resolution of the internal charge profile used by the long-range correction; increase and re-check if results drift
file	no		Output file path (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 10

A long-range correction for the finite lateral box, reported as F_iPB, is computed and added to the pressure automatically — it needs no configuration and handles valency mixtures. For reliable numbers, equilibrate first and start production from a state file (-s), and confirm the result is stable against density_bins and the lateral box size.

output.yaml reports {mean, error} for the midplane density rho_mid/Å⁻³ and the pressures p_ideal/mM (entropic), p_corr/mM (electrostatic), and p_osm/mM (total — the number you usually want), plus F_iPB/mM and p_osm/Pa; the p_corr and p_osm means include F_iPB. With file set, each sampled step writes step, rho_mid/Å⁻³, p_ideal/mM, p_corr/mM, p_osm/mM — there the per-step p_corr/p_osm columns exclude F_iPB (a run-averaged constant folded into the final means only).

---

Electric Potential Profile

Mean electric potential \(\varphi(z)\) along the \(z\)-axis of a slab, for an electrolyte with implicit salt (Debye–Hückel / Yukawa screening). A uniformly charged plane in a screened electrolyte produces a potential that decays exponentially away from it, so the laterally-averaged potential is a screened sum over the charge in each thin slab:

\[\varphi(z) = \frac{2\pi\,l_B}{\kappa}\sum_{z'} \sigma(z')\,e^{-\kappa\,|z-z'|}\]

where \(\sigma(z')\) is the average charge per unit area in the slab at \(z'\), \(l_B\) the Bjerrum length, and \(1/\kappa\) the Debye length — both taken from the medium. Because screening keeps the sum convergent, no infinite-plane correction term is needed. The walls are assumed neutral; only the explicit ions contribute.

Without salt (no Debye length) the kernel falls back automatically to bare Coulomb, the unscreened limit \(\kappa\to0\). The plane potential is then the one-dimensional Poisson Green's function

\[\varphi(z) = -2\pi\,l_B\sum_{z'}\sigma(z')\,|z-z'|,\]

which is finite and physical only for an electroneutral slab (the linear background cancels when \(\sum_{z'}\sigma(z')=0\)).

The cell geometry is detected automatically: a cuboid or slit must have a square base (\(L_x = L_y\)), and a cylinder uses its circular cross-section. The exponential treatment assumes each charged plane is effectively infinite, which holds only when the lateral box size is much larger than the Debye length — a warning is printed otherwise (screened kernel only).

Finite-box correction (optional)

For a box that is not much larger than the Debye length, set finite_box_correction: true. The infinite-plane potential is then replaced by that of the finite minimum-image cross-section, valid at any box size:

\[\varphi_\text{box}(z) = \varphi_\infty(z) - \varphi_\text{ext}(z), \qquad \varphi_\infty(z) = \frac{2\pi\,l_B}{\kappa}\,e^{-\kappa|z|},\]

where \(\varphi_\text{ext}\) is the contribution of the charge outside the cross-section. This is the screened (Yukawa) analogue of the finite-box correction of Greberg et al., doi:10/dhb9mj; screening makes every term finite, so the divergent and linear pieces of the bare-Coulomb construction cancel before quadrature. For a square base of half-width \(a\) (\(L_x = 2a\)),

\[\varphi_\text{ext}(z) = \frac{8\,l_B}{\kappa}\int_0^{\pi/4} \exp\!\Big(\!-\kappa\sqrt{a^2/\cos^2\theta + z^2}\,\Big)\,\mathrm{d}\theta,\]

a smooth integral evaluated by quadrature; for a disk of radius \(R\) (cylinder) it has the closed form

\[\varphi_\text{ext}(z) = \frac{2\pi\,l_B}{\kappa}\,e^{-\kappa\sqrt{R^2 + z^2}}.\]

Both vanish as the cross-section grows (\(\varphi_\text{ext}\to0\) for \(\kappa a\gg1\)), so the correction matters only for thin boxes. Enable it only when the simulation itself does not already apply such an external correction, otherwise the far field is subtracted twice.

For the unscreened (bare-Coulomb) kernel the correction reduces to Greberg's original square-base form,

\[\varphi_\text{ext}(z) = -2\pi\,l_B\,|z| - l_B\,u_\text{box}(z), \qquad u_\text{box}(z) = 8a\,\ln\!\frac{\sqrt{2a^2+z^2}+a}{\sqrt{a^2+z^2}} - 2z\left(\frac{\pi}{2} + \arcsin\frac{a^4 - z^4 - 2a^2z^2}{(a^2+z^2)^2}\right),\]

recovered from the screened form as \(\kappa\to0\) once the regularizing constant \(2\pi l_B/\kappa\) is removed. Because Greberg's construction models a square minimum-image box, the unscreened correction is only defined for a square base; enabling it for a cylinder is an error.

ζ-potential (optional)

Setting zeta_threshold reports a ζ-potential — the potential at the hydrodynamic shear plane — directly comparable to electrokinetic measurements. For a charged wall or brush the shear plane sits at the outer edge of the fixed charge, where \(\sigma(z)\) has decayed into the diffuse layer. The edge is located as the slab nearest the solution where

\[|\sigma(z_\text{shear})| > \texttt{zeta\_threshold}\cdot\max_z|\sigma(z)|,\]

with the charge-weighted centroid of \(|\sigma(z)|\) selecting which wall the fixed charge sits on. The reported ζ is the mean potential \(\varphi(z_\text{shear})\) with its statistical error, alongside the slab position \(z_\text{shear}\). A cutoff of 0.02 (2 % of the peak) is a reasonable starting point, but the right value is system-dependent, so inspect \(\sigma(z)\) in potential.csv and adjust. The estimate is only meaningful where the fixed charge is localized against a wall and decays into solution before the midplane; a slit charged symmetrically on both walls reports the equivalent value at one of them, and it is skipped when no charge is sampled. Note that the shear plane is here defined from the charge density, not the solvent (matter) density that sets the true no-slip surface, so the magnitude depends on this choice while the sign and trends are robust.

Example

analysis:
  - !ElectricPotentialProfile
    frequency: !Every 10
    # selection, resolution and file are optional:
    # selection: "all"      # atoms contributing charge (default: all)
    # resolution: 0.5       # slab thickness Δz in Å (default: 0.5)
    # zeta_threshold: 0.02  # also report a ζ-potential at the fixed-charge edge (off by default)
    # file: potential.csv   # output file (default: potential.csv)

Options

Key	Required	Default	Description
selection	no	all	Atoms whose charge contributes to the profile
resolution	no	0.5	Slab thickness Δz (Å) along the \(z\)-axis
finite_box_correction	no	false	Report the finite-box (Greberg) potential instead of the infinite-plane one
zeta_threshold	no	(off)	Report a ζ-potential at the fixed-charge edge, as a fraction of peak $
file	no	potential.csv	Output file path (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 10

The output file contains, per slab: the position z/Å, the slab charge density (per area e·Å⁻² and per volume e·Å⁻³), the potential potential/mV with its statistical error potential_error/mV, and the electric field field/mV·Å⁻¹. output.yaml additionally reports the potential at each wall and at the midplane (with the wall-to-midplane drops) as mean ± error, and, when zeta_threshold is set, the ζ-potential zeta/mV and the shear-plane position z_shear/Å. For reliable numbers, equilibrate first and start production from a state file (-s).

---

Density Profile

Average density of a chosen species along the \(z\)-axis. The cell is cut into slabs of equal thickness \(\Delta z\), and the density in the slab centred at \(z\) is the mean number of particles it holds divided by its volume,

\[\rho(z) = \frac{\langle N(z)\rangle}{A\,\Delta z},\]

where \(A\) is the cross-sectional area of the cell. Use it to measure how a species distributes itself between a wall and the bulk: ion layering at a charged surface, adsorption in a slit pore, or the density of a fluid along a cylindrical channel.

By default each selected atom counts towards the slab holding it. Set use_com: true to count each selected molecule in the slab holding its centre of mass instead, which profiles whole molecules rather than their individual atoms. Molecules declared atomic have no centre of mass, so selecting one with use_com is an error.

Every configuration contributes its own count, which is then averaged. The profile therefore stays correct when particles are created and destroyed: \(\langle N(z)\rangle\) is the mean occupancy of the slab in the grand canonical ensemble. Particles absent from a configuration do not count.

The cell must have the same cross-section at every height — a cuboid, slit, cylinder, or hexagonal prism — since a slab volume is otherwise undefined. Spherical and unbounded cells are rejected. The slabs are laid out once at the start of the run and keep their size, so the profile is meaningful only at constant volume; rescaling the cell later triggers a warning.

Example

Sodium and chloride profiles across a slit, plus the profile of whole water molecules:

analysis:
  - !DensityProfile
    selection: "atomtype Na"
    file: density_na.csv
    frequency: !Every 10
  - !DensityProfile
    selection: "atomtype Cl"
    file: density_cl.csv
    frequency: !Every 10
  - !DensityProfile
    selection: "molecule water"
    file: density_water.csv
    use_com: true          # bin mass centres, not individual atoms
    resolution: 0.5        # slab thickness Δz in Å (default: 1.0)
    frequency: !Every 10

Each block profiles one species and writes one file; add a block per species.

Options

Key	Required	Default	Description
selection	yes		Atoms, or molecules when use_com is set, whose density is profiled
file	yes		Output file path (see Output file formats)
use_com	no	false	Bin the centre of mass of each selected molecule instead of each atom
resolution	no	1.0	Slab thickness Δz (Å) along the \(z\)-axis
frequency	yes		Sample frequency, e.g. !Every 10

The output file gives, per slab, its centre z/Å followed by three equivalent measures of the same density, each with its statistical error: the number density density/Å⁻³, the molar concentration molarity/mol·L⁻¹, and the mass density mass_density/g·mL⁻¹. output.yaml additionally reports the slab layout and the mean total number of selected particles, mean_count, as mean ± error — a useful check, since summing \(\rho(z)\,A\,\Delta z\) over all slabs must return it.

---

Mean Along Coordinate

Computes the average of one collective variable (CV1) binned along another (CV2). CV2 is discretised into uniform bins of width resolution; no range is required since bins are created on demand using a sorted map. The output file contains columns for the bin center, running mean of CV1, and sample count per bin. The file is rewritten at each sample so partial results survive crashes.

Example

Average charge of GLU residue as a function of distance from a protein:

analysis:
  - !MeanAlongCoordinate
    property: charge
    selection: "molecule GLU"
    coordinate:
      property: mass_center_separation
      selection: "molecule GLU"
      selection2: "molecule protein"
      resolution: 0.5
    file: charge_vs_dist.dat
    frequency: !Every 100

Options

Key	Required	Default	Description
property	yes		CV1 type to average (see Supported properties)
selection	depends		Selection for CV1
coordinate	yes		CV2 block (must include resolution)
file	yes		Output file path (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 100

The coordinate block accepts all collective variable fields (property, selection, dimension, etc.) plus a required resolution for the bin width.

---

Rotational Diffusion

Estimates the anisotropic rotational diffusion tensor from the quaternion covariance matrix \(\tilde{Q}(\tau)\) of molecular orientations (Favro 1960; Holtbrügge & Schäfer 2025).

For each lag \(\tau\) (in snapshot units), the body-frame reorientation quaternion \(q(t,\tau) = q(t)\cdot q^{-1}(t+\tau)\) is computed, and its vector components form the \(3\times3\) covariance matrix \(\tilde{Q}_{ij}(\tau) = \langle q_i\cdot q_j\rangle\). In the principal coordinate system, \(\tilde{Q}\) is diagonal and follows the Favro model

\[\tilde{Q}_{ii}(\tau) = \tfrac{1}{4}\left(1 + e^{-(D_j+D_k)\tau} - e^{-(D_i+D_j)\tau} - e^{-(D_i+D_k)\tau}\right),\]

where \(D_x\), \(D_y\) and \(D_z\) are the principal rotational diffusion coefficients (rad²/snapshot). Correlations are averaged over all molecules matching the selection, exploiting ensemble averaging.

The selection must resolve to molecular groups (not atomic) that carry rigid-body quaternion state.

Example

analysis:
  - !RotationalDiffusion
    selection: "molecule Water"
    file: rotdiff.dat.gz
    frequency: !Every 100
    max_lag: 1000

Options

Key	Required	Default	Description
selection	yes		Selection for molecular group(s) to track
file	no		Streaming output file (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 100
max_lag	no	1000	Max lag in snapshots; actual time window = max_lag × frequency

The max_lag should be large enough for \(\tilde{Q}_{ii}(\tau)\) to approach its plateau at \(\tfrac{1}{4}\). A warning is emitted if convergence is not reached.

Output

The YAML output includes the covariance matrix at log-spaced lags, time-dependent diffusion coefficients \(D_x(\tau)\), \(D_y(\tau)\), \(D_z(\tau)\) from eigenvalue decomposition, and an isotropic estimate from the trace.

---

Widom Rotation

Measures how the surroundings restrict a molecule's orientation. For each snapshot a molecule is held at its center of mass and rigidly turned to M trial orientations spread evenly over all rotations (Alexa 2022). Only its interaction with the rest of the system is evaluated, using the run's own energy function. Turning a molecule about its center of mass leaves its internal geometry unchanged, so its internal energy plays no part.

From the trial energies u(Ω) the analysis reports, by Widom perturbation (Widom 1963):

\[W = -k_BT \ln\left\langle \frac{1}{M}\sum_k e^{-u_k/k_BT}\right\rangle\]

the orientation-averaged free energy W, the mean interaction, and the effective number of accessible orientations N_eff. The analysis also reports how strongly the surroundings restrict the molecule, as an entropy relative to free rotation:

\[S_\text{orient}/R = -\sum_k w_k \ln (M w_k)\]

where w_k is the Boltzmann weight of trial orientation k. The entropy follows directly from the sampled orientations and assumes nothing about the shape of the orientational well. It is dimensionless and never positive; zero means the molecule rotates freely. That restricting a molecule's orientation carries a free energy cost is the idea of Akke et al. 1993. Together with W and the mean interaction, the entropy splits the free energy into energetic and entropic parts.

For each requested molecular vector the analysis also reports the Lipari–Szabo order parameter (Lipari & Szabo 1982):

\[S^2 = \frac{3}{2}\sum_{\alpha\beta}\langle v_\alpha v_\beta\rangle^2 - \frac{1}{2}\]

where S²=1 means the vector is fully locked and S²=0 means it is free. A vector is given as a pair of atoms (!pair [i, j]), a gyration principal axis (!axis 0|1|2), or an explicit direction in the molecule's reference frame (!body [x, y, z]). Optionally, the analysis also reports the mean torque and the librational stiffness of the orientational cage.

The analysis reports one set of numbers per run; for a spatial profile, run separate umbrella windows and combine them afterwards. With implicit solvent the reported energies are free energies relative to pure solvent, not mechanical energies.

The selection must resolve to molecular groups (not atomic) of a single kind.

Example

analysis:
  - !WidomRotation
    selection: "molecule phytate"
    orientations: 500
    vectors:
      - !pair [0, 5]
      - !axis 2
      - !body [0.0, 0.0, 1.0]
    torque: true
    stiffness: true
    file: widomrot.csv.gz
    frequency: !Every 100

Options

Key	Required	Default	Description
selection	yes		Selection for molecular group(s) of one kind
orientations	yes		Number of trial orientations M (typically 100–1000)
vectors	no	[]	Vectors for S²: !pair [i,j], !axis 0\|1\|2, or !body [x,y,z]
torque	no	false	Also measure the mean torque
dtheta	no	0.01	Rotation step for the torque (rad)
stiffness	no	false	Also measure the librational stiffness
file	no		Streaming output file (see Output file formats)
frequency	yes		Sample frequency, e.g. !Every 100

Output

The YAML output includes W, the mean interaction, the orientational entropy, N_eff, and the per-vector and mean S², each with a statistical error. Torque and stiffness are added when enabled. If file is given, each sampled step writes columns step, W/kJ/mol, mean_S2, and N_eff.

W, the mean interaction, and the stiffness are in kJ/mol, the torque is in kT per radian, and the entropy and S² are dimensionless. Each output key names its own unit.

---

Rerun

The rerun subcommand replays a trajectory through a (possibly different) Hamiltonian, running the analysis pipeline on each frame. This decouples analysis from propagation, enabling e.g. comparison of explicit nonbonded energies against tabulated 6D potentials for the same configurations.

The input YAML provides the Hamiltonian and analysis configuration; the propagate: section is ignored. All analysis frequencies are overridden to sample every frame.

Fluctuating cells

Each frame is evaluated in the cell it was generated in, read from the box stored in the trajectory, so a constant-pressure trajectory reruns at its own volumes rather than at the volume declared in the input. A trajectory records an orthorhombic box, which determines a Cuboid, a Slit, a Sphere and a Cylinder. A HexagonalPrism is written to a trajectory as an expanded orthorhombic supercell and cannot be rerun. An Endless cell has no box, so it is left as the input declares it.

Reweighting biased trajectories

If the Hamiltonian contains a penalty term (from a converged Wang-Landau run), rerun automatically reweights all analyses by \(w = \exp(-\ln g(\text{bin}))\), recovering correct ensemble averages from the biased trajectory. This is logged at startup:

Reweighting enabled: penalty bias detected (Δg=X.X kT)

No special configuration is needed — include the same energy.penalty section used during the biased simulation.

Usage

faunus rerun -i input.yaml --traj traj.xtc [--aux traj.aux]

Flag	Required	Default	Description
-i	yes		Input YAML with Hamiltonian + analysis config
--traj	yes		XTC trajectory file
--aux	no	traj.aux (from --traj)	Frame state file

Requirements

The trajectory must have been produced with save_frame_state: true on the Trajectory analysis, so that a matching .aux file exists. The .aux file header must match the topology in the rerun input (same number of groups, particles, and molecule types).

Example workflow

1. Run the original simulation with frame state output:

   analysis:
     - !Trajectory
       file: traj.xtc
       frequency: !Every 100
       save_frame_state: true
     - !RadialDistribution
       selections: ["molecule A", "molecule B"]
       file: rdf_explicit.dat
       resolution: 0.1
       frequency: !Every 100
   

2. Rerun with a different Hamiltonian (e.g. 6D tabulated potential):

   faunus rerun -i input_6dtable.yaml --traj traj.xtc -o output_rerun.yaml
   

   where input_6dtable.yaml uses the same topology but a different energy section, and defines the desired analysis objects (e.g. RDF, energy time series).

3. Compare the RDFs from the original and rerun outputs.