Compute¶

Trajectory and structure analyses. All available via from molpy.compute import ....

See the Compute overview for the Compute → Result pattern and the Dielectric Spectroscopy guide for the theory behind the dielectric classes.

Quick reference¶

Symbol	Summary	Returns
`Compute`	Base class: configure once, call on data	a `Result` subclass
`DielectricSusceptibility`	\(\varepsilon^*(\omega)\) via Einstein–Helfand / Green–Kubo	`DielectricSusceptibilityResult`
`IonicConductivity`	Ionic conductivity \(\sigma\) (Einstein–Helfand MSD)	`ConductivityResult`
`ACFAnalyzer`	Autocorrelation of selected columns	`ACFResult`
`SpectralAnalyzer`	Windowed time→frequency transform	`SpectralResult`
`MCDCompute`	Mean displacement correlation (diffusion)	`MCDResult`
`PMSDCompute`	Polarization mean-squared displacement	`PMSDResult`
`Onsager`	Onsager transport coefficients \(L_{ij}\) (collective displacement cross-correlation)	`OnsagerResult`
`JACF`	Green–Kubo conductivity \(\sigma\) from the current ACF	`JACFResult`
`Persist`	Pair-survival / residence-time correlation	`PersistResult`
`RDF`	Radial distribution function \(g(r)\)	structural result
`MSD`	Single-particle mean-squared displacement	time series
`StaticStructureFactorDebye`	Static structure factor \(S(k)\) (Debye equation)	structural result
`NeighborList`	Spatial neighbor-pair query within a cutoff	pair list
`LocalDensity`, `GaussianDensity`	Per-particle / grid number density	density field
`Steinhardt`, `Hexatic`, `Nematic`, `SolidLiquid`	Bond-orientational order parameters	per-particle order
`BondOrder`	Bond-orientational \((\theta,\phi)\) diagram	spherical histogram
`PMFTXY`	Potential of mean force and torque on an \((x,y)\) grid	free-energy field
`RadiusOfGyration`, `GyrationTensor`, `InertiaTensor`, `CenterOfMass`	Molecular shape descriptors	per-frame values
`Cluster`, `ClusterCenters`, `Pca`, `KMeans`	Clustering & decomposition	labels / components
`DistanceDistribution`, `AngleDistribution`, `DihedralDistribution`	Geometric distribution functions (distance / ADF / DDF)	distribution
`CombinedDistribution`	Joint multi-axis distribution (CDF)	N-D histogram
`SpatialDistribution`	Spatial distribution function (SDF)	body-fixed density grid
`VanHove`	Van Hove correlation \(G(r,t)\) (self + distinct)	time-resolved \(g(r,t)\)
`LegendreReorientation`	Reorientational TCFs \(C_1(t)\), \(C_2(t)\)	correlation curves
`HBonds`, `HBondCriterion`	Geometric hydrogen-bond detection	per-frame bond lists
`RadicalVoronoi`, `VoronoiIntegration`, `voronoi_domains`, `voronoi_voids`	Radical Voronoi tessellation, domains/voids, charge integration	cells / domains / moments
`PowerSpectrum`, `IRSpectrum`, `RamanSpectrum`, `VcdSpectrum`, `RoaSpectrum`, `ResonanceRamanSpectrum`	Vibrational spectra from ACFs (VDOS / IR / Raman / VCD / ROA)	spectrum
`Workflow`	Directed graph of chained computes	per-node results

Full API¶

Base¶

base ¶

Base class for compute operations.

A :class:Compute is a configurable callable. Construction parameters go to __init__ (stored for serialization via :meth:Compute.dump); data inputs go to __call__. Operators take one or more data inputs directly — there is no single-input restriction:

>>> rdf = RDF(n_bins=100, r_max=10.0)
>>> result = rdf(frames, neighbors)      # two data inputs

The heavy numerics live in molrs; molpy operators are thin, typed shells that forward to the Rust kernels and return the molrs native result unchanged.

Compute ¶

Compute(**config)

Bases: ABC

Abstract base class for compute operations.

Subclasses implement :meth:__call__ with a concrete, fully typed signature (one positional parameter per data input) and pass their construction parameters to super().__init__(**config) so that :meth:dump can round-trip the configuration.

Examples:

>>> class MyCompute(Compute):
...     def __init__(self, scale: float):
...         super().__init__(scale=scale)
...         self.scale = scale
...
...     def __call__(self, frames: Sequence[Frame]) -> MyResult:
...         return MyResult(value=42 * self.scale)
>>>
>>> compute = MyCompute(scale=2.0)
>>> result = compute(frames)
>>> compute.dump()
{'scale': 2.0}

Store construction parameters for serialization.

Parameters:

Name	Type	Description	Default
`**config`	`Any`	Configuration parameters, returned verbatim by :meth:`dump`.	`{}`

dump ¶

dump()

Serialize construction configuration to a dictionary.

Returns:

Type	Description
`dict[str, Any]`	The configuration parameters passed to `__init__`.

Result types¶

result ¶

Result classes for compute operations.

This module defines result types returned by compute operations.

ACFResult `dataclass` ¶

ACFResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    acf=(lambda: np.array([]))(),
    n_lags=0,
)

Bases: TimeSeriesResult

Autocorrelation function result.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time lag values (in ps)
`acf`	`NDArray[float64]`	Autocorrelation values at each time lag, shape (n_lags,)
`n_lags`	`int`	Number of time lags

ConductivityResult `dataclass` ¶

ConductivityResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    msd=(lambda: np.array([]))(),
    sigma=float("nan"),
    slope=float("nan"),
    fit_start=0,
    fit_end=0,
)

Bases: TimeSeriesResult

Einstein-Helfand ionic-conductivity result.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	MSD lag times tau (ps), shape (n_lags,).
`msd`	`NDArray[float64]`	Collective MSD <\|M_J(t+tau) - M_J(t)\|^2> of the ionic translational dipole, (e*A)^2, shape (n_lags,).
`sigma`	`float`	Static ionic conductivity sigma (S/m).
`slope`	`float`	Fitted MSD slope over the diffusive window, (e*A)^2/ps.
`fit_start`	`int`	First lag index used in the linear fit (inclusive).
`fit_end`	`int`	Last lag index used in the linear fit (exclusive).

DebyeFit `dataclass` ¶

DebyeFit(
    tau=float("nan"),
    delta_eps=float("nan"),
    eps_inf=1.0,
    eps_static=float("nan"),
    omega_peak=float("nan"),
)

Single-Debye relaxation parameters fitted from a dielectric spectrum.

Attributes:

Name	Type	Description
`tau`	`float`	Relaxation time (ps).
`delta_eps`	`float`	Relaxation strength epsilon(0) - epsilon_inf (dimensionless).
`eps_inf`	`float`	High-frequency permittivity used in the fit.
`eps_static`	`float`	Static permittivity epsilon(0) (dimensionless).
`omega_peak`	`float`	Angular frequency of the dielectric-loss peak (rad/ps).

epsilon ¶

epsilon(omega)

Evaluate the fitted Debye model at angular frequencies omega.

Returns:

Type	Description
`NDArray`	`(epsilon_real, epsilon_imag)` under the positive-loss convention
`NDArray`	`epsilon* = epsilon' - i epsilon''`.

DielectricResult `dataclass` ¶

DielectricResult(
    meta=dict(),
    frequency=(lambda: np.array([]))(),
    epsilon_real=(lambda: np.array([]))(),
    epsilon_imag=(lambda: np.array([]))(),
    epsilon_static=float("nan"),
    epsilon_inf=1.0,
    route="",
    component="",
    conductivity=None,
)

Bases: Result

Single-route dielectric susceptibility result.

Attributes:

Name	Type	Description
`frequency`	`NDArray[float64]`	Angular frequency grid omega, shape (n_freq,), units rad/ps. Bin 0 is DC; bin 1 is Delta-omega = 2 * pi / (n_pad * dt).
`epsilon_real`	`NDArray[float64]`	Real part epsilon'(omega), shape (n_freq,), dimensionless.
`epsilon_imag`	`NDArray[float64]`	Loss spectrum epsilon''(omega), shape (n_freq,), dimensionless, positive sign convention.
`epsilon_static`	`float`	Static dielectric constant epsilon(0), dimensionless. May be `nan` if the route does not provide a static estimate.
`epsilon_inf`	`float`	High-frequency dielectric constant.
`route`	`str`	Computation route ("einstein-helfand" or "green-kubo").
`component`	`str`	System component ("full", "water", "ion").
`conductivity`	`NDArray[float64] \| None`	Optional conductivity spectrum sigma(omega), shape (n_freq,).

fit_debye ¶

fit_debye()

Fit a single Debye relaxation to this spectrum (NumPy only).

Uses the exact single-Debye identity epsilon''(omega) / (epsilon'(omega) - epsilon_inf) = omega * tau: tau is the least-squares slope through the origin of that ratio versus omega over the low-frequency rising branch (up to the loss peak), with a loss-peak fallback tau = 1 / omega_peak. The relaxation strength is the static limit delta_eps = epsilon(0) - epsilon_inf.

No SciPy: the estimator is closed-form linear regression. For broadened or skewed (Cole-Cole / Havriliak-Negami) line shapes do a nonlinear fit in your analysis script using :meth:DebyeFit.epsilon as the model.

Returns:

Type	Description
`'DebyeFit'`	DebyeFit with tau (ps), delta_eps, eps_inf, eps_static, omega_peak.

DielectricSusceptibilityResult `dataclass` ¶

DielectricSusceptibilityResult(
    meta=dict(), results=dict(), metadata=dict()
)

Bases: Result

Aggregate dielectric susceptibility result.

Attributes:

Name	Type	Description
`results`	`dict[str, DielectricResult]`	Mapping from route-component key to DielectricResult
`metadata`	`dict[str, Any]`	Trajectory parameters and computation info

to_dict ¶

to_dict()

Serialize with nested DielectricResult recursion.

JACFResult `dataclass` ¶

JACFResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    jacf=(lambda: np.array([]))(),
    sigma_running=(lambda: np.array([]))(),
    sigma=float("nan"),
)

Bases: TimeSeriesResult

Results from a Green-Kubo current-autocorrelation conductivity.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time lag values (in ps), shape (n_time_lags,).
`jacf`	`NDArray[float64]`	Current autocorrelation `C(tau) = <J(0).J(tau)>`, (e*A/ps)^2, shape (n_time_lags,).
`sigma_running`	`NDArray[float64]`	Running Green-Kubo conductivity integral `sigma(tau) = 1/(3 V kB T) integral_0^tau C(t) dt` (S/m), shape (n_time_lags,).
`sigma`	`float`	DC ionic conductivity (S/m) — `sigma_running` at the final lag.

MCDResult `dataclass` ¶

MCDResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    correlations=dict(),
)

Bases: TimeSeriesResult

Results from Mean Displacement Correlation calculation.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time lag values (in ps)
`correlations`	`dict[str, NDArray[float64]]`	Dictionary mapping tag names to correlation function arrays (MSD values). Each array has shape (n_time_lags,)

OnsagerResult `dataclass` ¶

OnsagerResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    correlations=dict(),
)

Bases: TimeSeriesResult

Results from an Onsager collective-displacement cross-correlation.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time lag values (in ps), shape (n_time_lags,).
`correlations`	`dict[str, NDArray[float64]]`	Mapping from tag `"i,j"` to the cross-correlation `L_ij(tau) = <DP_i(tau).DP_j(tau)>` of the collective (summed) species displacements, shape (n_time_lags,), units A^2. The diagonal `"i,i"` is the collective MSD of species `i`.

PMSDResult `dataclass` ¶

PMSDResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    pmsd=(lambda: np.array([]))(),
)

Bases: TimeSeriesResult

Results from Polarization Mean Square Displacement calculation.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time lag values (in ps)
`pmsd`	`NDArray[float64]`	Polarization MSD values at each time lag, shape (n_time_lags,)

PersistResult `dataclass` ¶

PersistResult(
    meta=dict(),
    time=(lambda: np.array([]))(),
    correlations=dict(),
)

Bases: TimeSeriesResult

Results from a pair-survival (persistence) correlation.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time lag values (in ps), shape (n_time_lags,).
`correlations`	`dict[str, NDArray[float64]]`	Mapping from tag `"i,j:method:r0[,r1]"` to the persistence correlation `C(tau)` (mean surviving partners per reference particle), shape (n_time_lags,). `C(0)` is the mean coordination number.

Result `dataclass` ¶

Result(meta=dict())

Base class for computation results.

Subclasses should define specific fields for their result data.

to_dict ¶

to_dict()

Convert result to dictionary representation.

SpectralResult `dataclass` ¶

SpectralResult(
    meta=dict(),
    frequency=(lambda: np.array([]))(),
    spectrum=(lambda: np.array([]))(),
)

Bases: Result

Frequency-domain spectrum result.

Attributes:

Name	Type	Description
`frequency`	`NDArray[float64]`	Angular frequency grid omega, shape (n_freq,), units rad/ps.
`spectrum`	`NDArray[float64]`	Spectral density at each frequency, shape (n_freq,).

TimeSeriesResult `dataclass` ¶

TimeSeriesResult(
    meta=dict(), time=(lambda: np.array([]))()
)

Bases: Result

Base class for time-series analysis results.

Attributes:

Name	Type	Description
`time`	`NDArray[float64]`	Time points for the analysis (in ps or frames)

Dielectric¶

dielectric ¶

Dielectric susceptibility Compute classes.

Thin glue layers bridging molpy Trajectory to molrs computational kernels. The Python side does only data extraction (positions, charges) and vectorized NumPy assembly (dipole moment via einsum, minimum-image unwrap via Box.diff_dr); all spectral physics — ACF, windowing, FFT, prefactors — is performed in Rust by the raw computes (molrs.DebyeRelaxation, molrs.GreenKuboConductivity) and the ε(ω) Fits (molrs.EinsteinHelfandSpectrum, molrs.GreenKuboSpectrum), plus the raw molrs.dielectric observables and molrs.signal.

ACFAnalyzer ¶

ACFAnalyzer(
    columns, max_lag, *, unwrap=True, **config_kwargs
)

Bases: Compute

Compute autocorrelation function from trajectory data.

Extracts per-atom columns from each frame, optionally unwraps coordinates via Box.diff_dr, delegates to molrs.signal.acf_fft(), normalizes the ACF (divides by zero-lag value), and returns an ACFResult.

DielectricSusceptibility ¶

DielectricSusceptibility(
    dt,
    temperature,
    max_correlation_time,
    *,
    epsilon_inf=1.0,
    window_type="hann",
    routes=None,
    volume=None,
    **config_kwargs,
)

Bases: Compute

Frequency-dependent dielectric susceptibility from an MD trajectory.

Extracts atomic positions and charges per frame, unwraps coordinates via minimum-image convention, builds the total dipole moment series, and runs one or more spectral routes (Einstein-Helfand and/or Green-Kubo) as the explicit raw-compute + ε(ω)-Fit composition: molrs.DebyeRelaxation (raw fluctuation dipole ACF) + molrs.EinsteinHelfandSpectrum, and molrs.GreenKuboConductivity (raw current ACF) + molrs.GreenKuboSpectrum. No spectra math runs in Python. The static dielectric constant is computed once via Neumann's fluctuation formula and attached to every result.

Parameters:

Name	Type	Description	Default
`dt`	`float`	Frame spacing in ps.	required
`temperature`	`float`	Temperature in K.	required
`max_correlation_time`	`int`	Longest ACF lag in frames (clamped to `n_frames - 1`). Practical choice: ≤ n_frames / 10.	required
`epsilon_inf`	`float`	High-frequency (electronic) permittivity, dimensionless. Use 1.0 for non-polarizable force fields.	`1.0`
`window_type`	`str`	`"hann"` or `"blackman"` window applied to the ACF before FFT.	`'hann'`
`routes`	`list[str] \| None`	Subset of `["einstein-helfand", "green-kubo"]`. Default runs both.	`None`
`volume`	`float \| None`	System volume in Å³. If `None`, uses `frame.box.volume` from the first frame (assumes NVT/NVE).	`None`

Inputs

Each frame's atoms block must contain canonical columns x, y, z (Å) and charge (e). Frames must carry a non-free Box.

IonicConductivity ¶

IonicConductivity(
    dt,
    temperature,
    max_correlation_time,
    *,
    volume=None,
    fit_start_frac=0.1,
    fit_end_frac=0.5,
    **config_kwargs,
)

Bases: Compute

Static ionic conductivity sigma via the Einstein-Helfand relation.

Builds the ionic translational dipole M_J(t) = sum_i q_i r_i(t) from the trajectory (minimum-image unwrapped, same as :class:DielectricSusceptibility), then composes the raw collective-dipole MSD (:class:molrs.EinsteinConductivity) with the diffusive-window slope (:class:molrs.LinearFit) and a slope / (6 V k_B T) S/m prefactor:

sigma = lim_{t->inf} (1 / (6 V k_B T)) d/dt <|M_J(t) - M_J(0)|^2>.

Decomposition is the caller's responsibility and is done with selection, not arithmetic: pass a trajectory whose charge column is non-zero only on the mobile ions (e.g. via a :class:~molpy.Selector over the ion atoms, or by zeroing solvent charges). Including the solvent rotational dipole here would contaminate the translational MSD.

Parameters:

Name	Type	Description	Default
`dt`	`float`	Frame spacing in ps.	required
`temperature`	`float`	Temperature in K.	required
`max_correlation_time`	`int`	Longest MSD lag in frames (clamped to `n_frames - 1`). Practical choice: <= `n_frames / 5`.	required
`volume`	`float \| None`	System volume in A^3. If `None`, uses `frame.box.volume` from the first frame (assumes NVT/NVE).	`None`
`fit_start_frac, fit_end_frac`		Fractions of `max_lag` bounding the linear-fit window over the diffusive regime (default 0.1, 0.5). `sigma` is window-sensitive for few-carrier systems; report a range rather than a single digit.	required

Inputs

Each frame's atoms block must contain x, y, z (A) and charge (e); frames must carry a non-free Box.

SpectralAnalyzer ¶

SpectralAnalyzer(
    dt, *, window_type="hann", **config_kwargs
)

Bases: Compute

Convert time-domain ACF to frequency-domain spectrum.

Applies a window function, generates the frequency grid, and performs the time→frequency conversion. All computation delegated to molrs.signal.

Mean displacement correlation¶

mcd ¶

Mean Displacement Correlation (MCD) computation for diffusion analysis.

This module implements the MCD method for computing self and distinct displacement correlations (MSD) from molecular dynamics trajectories.

Adapted from the tame library (https://github.com/Roy-Kid/tame).

MCDCompute ¶

MCDCompute(tags, max_dt, dt, center_of_mass=None)

Bases: Compute

Compute Mean Displacement Correlations (MSD) for diffusion analysis.

This class implements the MCD method which computes time correlation functions of particle displacements. It returns raw MSD values without fitting. It supports:

Self diffusion: MSD_i = <(r_i(t+dt) - r_i(t))²>
Distinct diffusion: <(r_i(t+dt) - r_i(t)) · (r_j(t+dt) - r_j(t))>

Parameters:

Name	Type	Description	Default
`tags`	`list[str]`	List of atom type specifications. Each tag can be: - Single integer (e.g., "3"): Self-diffusion MSD of type 3 - Two integers separated by comma (e.g., "3,4"): Distinct diffusion correlation between types 3 and 4	required
`max_dt`	`float`	Maximum time lag in ps	required
`dt`	`float`	Timestep in ps	required
`center_of_mass`	`dict[int, float] \| None`	Optional dict mapping element types to masses for COM removal. Format: {element_type: mass}, e.g., {1: 1.008, 6: 12.011}	`None`

Examples:

>>> from molpy.io import read_h5_trajectory
>>> traj = read_h5_trajectory("trajectory.h5")
>>>
>>> # Compute self-diffusion MSD of atom type 3
>>> mcd = MCDCompute(tags=["3"], max_dt=30.0, dt=0.01)
>>> result = mcd(traj)
>>> print(result.correlations["3"])  # MSD values at each time lag
>>>
>>> # Compute distinct diffusion between types 3 and 4
>>> mcd = MCDCompute(tags=["3,4"], max_dt=30.0, dt=0.01)
>>> result = mcd(traj)
>>> print(result.correlations["3,4"])  # Correlation values

Notes (tame-port audit): - Different-species distinct term. This implementation uses the collective cross-correlation <(sum_i Dr_i).(sum_j Dr_j)>. The tame original (tame/recipes/mdc.py) instead uses mean_i for the reference species, i.e. (1/N_i) times this value — an inconsistent normalization relative to its like-species and self terms (both collective). The collective form here is the physically correct Onsager cross term; for explicit, fully-normalized Onsager coefficients use :class:~molpy.compute.Onsager. - Pair-resolved (SSP) tags. tame's mdc also accepts a "i,j:METHOD:r0,r1" form that splits the partner displacement into surviving ("in") vs broken ("out") neighbours via tpairsurvive. That decomposition is not reproduced here; the pair-survival aspect is available separately as :class:~molpy.compute.Persist.

Polarization MSD¶

pmsd ¶

Polarization Mean Square Displacement (PMSD) computation.

This module implements the PMSD method for computing polarization fluctuations in ionic systems from molecular dynamics trajectories.

Adapted from the tame library (https://github.com/Roy-Kid/tame).

PMSDCompute ¶

PMSDCompute(cation_type, anion_type, max_dt, dt)

Bases: Compute

Compute Polarization Mean Square Displacement for ionic systems.

This class computes the PMSD which measures polarization fluctuations in ionic systems. The polarization is defined as:

P(t) = Σ_cations r_i(t) - Σ_anions r_j(t)

And the PMSD is:

PMSD(dt) = <(P(t+dt) - P(t))²>_t

Parameters:

Name	Type	Description	Default
`cation_type`	`int`	Atom type index for cations	required
`anion_type`	`int`	Atom type index for anions	required
`max_dt`	`float`	Maximum time lag in ps	required
`dt`	`float`	Timestep in ps	required

Examples:

>>> from molpy.io import read_h5_trajectory
>>> traj = read_h5_trajectory("ionic_liquid.h5")
>>>
>>> # Compute PMSD for Li+ (type 1) and PF6- (type 2)
>>> pmsd = PMSDCompute(cation_type=1, anion_type=2, max_dt=30.0, dt=0.01)
>>> result = pmsd(traj)
>>>
>>> # Plot results
>>> import matplotlib.pyplot as plt
>>> plt.plot(result.time, result.pmsd)
>>> plt.xlabel("Time lag (ps)")
>>> plt.ylabel("PMSD (Å²)")
>>> plt.show()

Notes (tame-port audit): This returns the PMSD curve only. The tame original (tame/recipes/pmsd.py) additionally fits the long-time slope to a conductivity via the Einstein relation sigma = slope / (6 V kB T). That conductivity is provided here as a dedicated compute, :class:~molpy.compute.IonicConductivity (Einstein-Helfand, S/m), with the equivalent Green-Kubo route in :class:~molpy.compute.JACF.

Onsager coefficients¶

onsager ¶

Onsager transport coefficients from collective mean-displacement correlations.

The Onsager phenomenological coefficients L_ij describe coupled transport of species i and j. They follow from the cross-correlation of the collective (summed) displacements of each species::

L_ij(tau) = <DP_i(tau) . DP_j(tau)>_t ,
    P_s(t)  = sum_{a in species s} r_a(t)   (unwrapped),
    DP_s(tau) = P_s(t+tau) - P_s(t).

The diagonal L_ii is the collective MSD of species i; off-diagonal L_ij captures the cross-correlated drift distinguishing the Onsager picture from the bare Nernst-Einstein sum. A long-time linear fit of L_ij yields the transport coefficient (left to the caller).

The numerically heavy windowed (all-time-origins) cross-correlation runs in Rust (molrs.transport.Onsager); this wrapper extracts coordinates, performs the minimum-image unwrap via :meth:molrs.Box.delta, and builds the per-species collective coordinates.

Adapted from the tame library (https://github.com/Roy-Kid/tame), tame/recipes/onsager.py.

Onsager ¶

Onsager(tags, max_dt, dt, center_of_mass=None)

Bases: Compute

Compute Onsager collective-displacement cross-correlations.

Parameters:

Name	Type	Description	Default
`tags`	`list[str]`	Species pairs `"i,j"` (e.g. `["1,1", "1,2", "2,2"]`). Each yields the cross-correlation of the collective displacements of species `i` and `j`.	required
`max_dt`	`float`	Maximum time lag in ps.	required
`dt`	`float`	Timestep in ps.	required
`center_of_mass`	`dict[int, float] \| None`	Optional `{type: mass}` mapping; when given, the system center of mass is removed each frame before forming the collective coordinates.	`None`

Examples:

>>> from molpy.io import read_h5_trajectory
>>> traj = read_h5_trajectory("electrolyte.h5")
>>> ons = Onsager(tags=["1,1", "1,2", "2,2"], max_dt=20.0, dt=0.01)
>>> result = ons(traj)
>>> result.correlations["1,2"]  # L_12(tau), shape (n_cache,)

Current-ACF conductivity (Green–Kubo)¶

jacf ¶

Ionic conductivity from the charge-current autocorrelation (Green-Kubo).

The DC ionic conductivity follows from the Green-Kubo relation for the collective charge current J(t) = sum_a q_a v_a(t)::

sigma = 1 / (3 V kB T) * integral_0^inf <J(0).J(t)> dt .

This wrapper assembles the collective current J = sum v_cation - sum v_anion (unit charges +/-1) from per-atom velocities and delegates the current autocorrelation + trapezoidal Green-Kubo integral to Rust via the explicit raw-compute + fit composition (molrs.GreenKuboConductivity for the raw current ACF, then molrs.RunningIntegral for the cumulative integral). The only Python-side step is multiplying by the fixed Green-Kubo unit prefactor e^2 A^2 ps^-1 / (3 A^3 kB) to convert the integral to SI S/m — no ACF / integration math is reimplemented in Python.

Units (LAMMPS real, matching :mod:molpy.compute.dielectric): velocities in A/ps (so J is e*A/ps), dt in ps, volume in A^3, temperature in K; the output sigma is SI S/m. Pass velocities already in A/ps; conductivity scales linearly so other velocity units can be rescaled afterwards.

Adapted from the tame library (https://github.com/Roy-Kid/tame), tame/recipes/jacf.py (whose published version never evaluates the autocorrelation before integrating — corrected here).

JACF ¶

JACF(
    cation_type,
    anion_type,
    max_dt,
    dt,
    temperature,
    volume=None,
)

Bases: Compute

Green-Kubo ionic conductivity from the charge-current autocorrelation.

Parameters:

Name	Type	Description	Default
`cation_type`	`int`	Atom type index for cations (charge +1).	required
`anion_type`	`int`	Atom type index for anions (charge -1).	required
`max_dt`	`float`	Maximum correlation time in ps.	required
`dt`	`float`	Timestep in ps.	required
`temperature`	`float`	Temperature in K.	required
`volume`	`float \| None`	System volume in A^3. If `None`, the mean box volume over the trajectory is used.	`None`

Examples:

>>> from molpy.io import read_h5_trajectory
>>> traj = read_h5_trajectory("ionic_liquid.h5")
>>> jacf = JACF(cation_type=1, anion_type=2, max_dt=30.0, dt=0.01,
...             temperature=300.0)
>>> result = jacf(traj)
>>> result.sigma          # DC ionic conductivity, S/m
>>> result.jacf           # <J(0).J(t)>, shape (n_cache,)

Pair persistence¶

persist ¶

Pair-survival (persistence) time-correlation functions.

Measures how long pairs of particles remain within a distance cutoff as a function of time lag — residence-time / hydrogen-bond-dynamics analysis. For a reference species i and partner species j::

C(tau) = < (1/N_i) sum_i sum_j S_ij(t, t+tau) >_t ,

where S_ij in {0,1} is the survival indicator for the pair born at t and observed at t+tau. C(0) is the mean coordination number.

Three survival criteria (see :class:molrs.transport.Persist):

continuous — within the survival cutoff at every frame since birth.
intermittent — within the cutoff at t+tau (re-formation allowed).
ssp — stable-state picture: born within inner cutoff r0, continuously within outer cutoff r1 (r1 >= r0) since.

The per-pair, per-frame survival accounting runs in Rust (molrs.transport.Persist); this wrapper extracts per-species coordinates and per-frame orthorhombic box edge lengths.

Adapted from the tame library (https://github.com/Roy-Kid/tame), tame/recipes/persist.py / tame/ops/time.py (tpairsurvive). The published persist.py recipe is non-functional (undefined names); this port implements the intended correlation with explicit survival criteria.

Persist ¶

Persist(tags, max_dt, dt)

Bases: Compute

Compute pair-survival (persistence) time-correlation functions.

Parameters:

Name	Type	Description	Default
`tags`	`list[str]`	Pair specifications `"t1,t2:method:r0[,r1]"` — e.g. `"3,4:ssp:3.0,4.0"` (cation-anion stable-state pairs born within 3 A, surviving while within 4 A) or `"1,1:continuous:3.5"` (like-species, single cutoff). `method` is one of `continuous` / `intermittent` / `ssp`.	required
`max_dt`	`float`	Maximum time lag in ps.	required
`dt`	`float`	Timestep in ps.	required

Examples:

>>> from molpy.io import read_h5_trajectory
>>> traj = read_h5_trajectory("electrolyte.h5")
>>> p = Persist(tags=["3,4:ssp:3.0,4.0"], max_dt=30.0, dt=0.01)
>>> result = p(traj)
>>> result.correlations["3,4:ssp:3.0,4.0"]  # C(tau), shape (n_cache,)

Radial distribution¶

rdf ¶

Radial distribution function g(r) — molrs-backed.

Returns molrs.compute.density.RDFResult directly (no molpy wrapper). The result is finalized eagerly inside RDF.compute, so result.rdf is the normalized g(r) and not the raw histogram.

References¶

V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, S. C. Glotzer, Comput. Phys. Commun. 254, 107275 (2020) — the freud library, on which this kernel is modelled.
M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, 2nd ed., Oxford (2017) — g(r) and coordination numbers.

RDF ¶

RDF(n_bins, r_max, r_min=0.0)

Bases: Compute

Histogram pair distances into g(r) over one or more frames.

Parameters¶

n_bins : int Number of histogram bins. r_max : float Upper edge of the last bin in Angstroms. r_min : float, default 0.0 Lower edge of bin 0 in Angstroms.

Notes¶

RDF takes two data inputs (frames + neighbor lists).

Mean-squared displacement¶

msd ¶

Mean Squared Displacement — molrs-backed.

Returns molrs.compute.msd.MSDTimeSeries directly. Frame[0] is the reference; series.mean[i] is ⟨|r(i) - r(0)|²⟩ averaged over particles.

MSD ¶

MSD()

Bases: Compute

Mean squared displacement against frame[0].

Examples¶

series = MSD()(trajectory_frames) series.mean.shape # (n_frames,)

Static structure factor¶

diffraction ¶

Diffraction observables — molrs-backed.

StaticStructureFactorDebye computes the static structure factor

.. math::

S(k) = \frac{1}{N} \left\langle \sum_{i,j}
        \frac{\sin(k r_{ij})}{k r_{ij}} \right\rangle

via the Debye scattering equation. Thin shell over molrs.compute.diffraction.StaticStructureFactorDebye; takes (frames).

References¶

P. Debye, Ann. Phys. 351, 809 (1915) — the Debye scattering equation.
V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, on which this kernel is modelled.

StaticStructureFactorDebye ¶

StaticStructureFactorDebye(k_values)

Bases: Compute

Static structure factor S(k) via the Debye equation.

Parameters¶

k_values : array-like Scattering wavenumbers (1/Angstrom) at which to evaluate S(k).

Neighbor list¶

neighborlist ¶

Spatial neighbor list — molrs-backed.

Returns molrs.NeighborList directly (no molpy wrapper). Coordinates are stacked once via frame["atoms"][["x", "y", "z"]] (the only unavoidable copy, internal to Block.__getitem__(list)); from that point through to the returned indices/distances the path is zero-copy borrowed views into the molrs Rust buffers.

References¶

V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, whose locality neighbor query this mirrors.

NeighborList ¶

NeighborList(cutoff)

Bases: Compute

Spatial neighbor-pair query within a cutoff radius.

Parameters¶

cutoff : float Cutoff radius in Angstroms.

Examples¶

nlist = NeighborList(cutoff=2.5)(frame) nlist.n_pairs 1834

Local & grid density¶

density ¶

Density-field operators — molrs-backed.

Thin Compute shells over molrs.compute.density.*, returning the molrs native results unchanged. LocalDensity takes (frames, nlists) (like RDF); GaussianDensity takes (frames) only.

References¶

V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, on which these density kernels are modelled.

GaussianDensity ¶

GaussianDensity(nx, ny, nz, sigma)

Bases: Compute

Gaussian-smeared number density on a 3-D grid.

Parameters¶

nx, ny, nz : int Grid resolution along each axis. sigma : float Gaussian smearing width (Angstrom).

LocalDensity ¶

LocalDensity(r_max, diameter=0.0)

Bases: Compute

Per-particle local number density within a cutoff sphere.

Parameters¶

r_max : float Cutoff radius (Angstrom) of the counting sphere. diameter : float Particle diameter correction; 0.0 counts particle centres only.

Bond-orientational order parameters¶

order ¶

Bond-orientational order operators — molrs-backed.

Thin Compute shells over molrs.compute.order.*. Each forwards verbatim to the Rust kernel and returns the molrs native result; molpy adds no wrapping and copies nothing. Like compute.RDF, these take two data inputs.

References¶

Steinhardt, Nelson & Ronchetti, Phys. Rev. B 28, 784 (1983).
Nelson & Halperin, Phys. Rev. B 19, 2457 (1979) (hexatic).
ten Wolde, Ruiz-Montero & Frenkel, J. Chem. Phys. 104, 9932 (1996) (solid-liquid).
de Gennes & Prost, The Physics of Liquid Crystals, 2nd ed. (1993) (nematic).
W. Lechner, C. Dellago, J. Chem. Phys. 129, 114707 (2008) (averaged variant).
V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, on which these order parameters are modelled.

Hexatic ¶

Hexatic(k)

Bases: Compute

Two-dimensional :math:\psi_k hexatic bond-orientational order.

Parameters¶

k : int Symmetry order (6 for the hexatic order parameter).

Nematic ¶

Nematic()

Bases: Compute

Nematic order parameter and Q-tensor from per-particle directors.

The second call argument is a (N, 3) array of per-particle orientation vectors. Returns (order, eigenvalues, director, q_tensor).

SolidLiquid ¶

SolidLiquid(l, q_threshold=0.7, n_threshold=6)

Bases: Compute

Solid-liquid classification via :math:q_\ell bond correlations.

Parameters¶

l : int Spherical-harmonic degree used for the bond correlation. q_threshold : float Minimum dot-product for a bond to count as solid-like. n_threshold : int Minimum number of solid-like bonds for a particle to be solid.

Steinhardt ¶

Steinhardt(l, average=False, wl=False, wl_normalize=False)

Bases: Compute

Steinhardt :math:q_\ell / :math:w_\ell bond-orientational order.

Parameters¶

l : Sequence[int] Spherical-harmonic degrees, e.g. [6] or [4, 6]. average : bool Use the locally averaged variant. wl : bool Compute third-order :math:w_\ell invariants instead of :math:q_\ell. wl_normalize : bool Normalize the :math:w_\ell invariants.

Bond-orientational environment¶

environment ¶

Local-environment observables — molrs-backed.

BondOrder histograms neighbor bond directions on a (theta, phi) grid. Thin shell over molrs.compute.environment.BondOrder; takes (frames, nlists) like RDF.

References¶

V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, whose environment.BondOrder this mirrors.

BondOrder ¶

BondOrder(n_theta, n_phi)

Bases: Compute

Bond-orientational order diagram on a spherical (theta, phi) grid.

Parameters¶

n_theta : int Number of polar-angle bins. n_phi : int Number of azimuthal-angle bins.

Potential of mean force & torque¶

pmft ¶

Potential of mean force and torque — molrs-backed.

PMFTXY accumulates the 2-D potential of mean force in the local (x, y) frame of each particle. Thin shell over molrs.compute.pmft.PMFTXY; takes (frames, nlists, orientations=None).

References¶

G. van Anders, D. Klotsa, N. K. Ahmed, M. Engel, S. C. Glotzer, ACS Nano 8, 931 (2014) — potential of mean force and torque.
V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, on which this kernel is modelled.

PMFTXY ¶

PMFTXY(x_max, y_max, n_x, n_y)

Bases: Compute

2-D potential of mean force and torque on an (x, y) grid.

Parameters¶

x_max, y_max : float Half-extent of the integration window along each local axis. n_x, n_y : int Grid resolution along each axis.

Shape descriptors¶

shape ¶

Per-cluster shape descriptors — molrs-backed.

Thin wrappers around molrs.compute.cluster:

CenterOfMass (frames, clusters) → mass-weighted centers
GyrationTensor (frames, clusters, centers) → 3×3 tensors per cluster
InertiaTensor (frames, clusters, com) → 3×3 tensors per cluster
RadiusOfGyration (frames, clusters, com) → Rg per cluster

References¶

V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, on which these per-cluster reductions are modelled.
D. N. Theodorou, U. W. Suter, Macromolecules 18, 1206 (1985) — gyration tensor and shape descriptors.

CenterOfMass ¶

CenterOfMass(masses=None)

Bases: Compute

Mass-weighted center per cluster.

Parameters¶

masses : ndarray | None Per-particle masses (length N). None falls back to unit mass.

GyrationTensor ¶

GyrationTensor()

Bases: Compute

Gyration tensor per cluster (unweighted).

InertiaTensor ¶

InertiaTensor(masses=None)

Bases: Compute

Inertia tensor per cluster.

Parameters¶

masses : ndarray | None Per-particle masses (length N). None falls back to unit mass.

RadiusOfGyration ¶

RadiusOfGyration(masses=None)

Bases: Compute

Radius of gyration per cluster.

Parameters¶

masses : ndarray | None Per-particle masses (length N). None falls back to unit mass.

Decomposition¶

decomposition ¶

Dimensionality reduction + clustering ML primitives — molrs-backed.

Pca (alias for molrs Pca2): 2-component PCA over a list of DescriptorRow objects.
KMeans: k-means clustering over a PcaResult.

DescriptorRow is re-exported as the input wrapper: each row is a 1-D float ndarray passed through DescriptorRow(row).

KMeans ¶

KMeans(k, max_iter=100, seed=0)

Bases: Compute

k-means over a PcaResult.

Parameters¶

k : int Number of clusters. max_iter : int, default 100 Lloyd-iteration cap. seed : int, default 0 RNG seed for centroid initialization.

Pca ¶

Pca()

Bases: Compute

Two-component PCA. Input: list of DescriptorRow.

Clustering¶

cluster ¶

Distance-based clustering — molrs-backed.

Cluster takes (frames, nlists) and returns one ClusterResult per frame. ClusterCenters takes (frames, clusters) and returns the geometric centers per cluster.

Both wrappers take two data inputs — mirroring the RDF pattern.

References¶

V. Ramasubramani et al., Comput. Phys. Commun. 254, 107275 (2020) — the freud library, whose cluster module this mirrors.

Cluster ¶

Cluster(min_cluster_size)

Bases: Compute

Group particles into clusters via a neighbor-list connectivity graph.

Parameters¶

min_cluster_size : int Minimum size for a connected component to count as a cluster.

ClusterCenters ¶

ClusterCenters()

Bases: Compute

Geometric centers per cluster (unweighted).

ClusterProperties ¶

ClusterProperties()

Bases: Compute

Per-cluster properties (sizes, centers, masses, gyration tensors, Rg).

Takes (frames, clusters) where clusters is the sequence of ClusterResult returned by :class:Cluster. Returns one dict per frame.

Distribution functions (ADF / DDF / distance / combined)¶

distribution ¶

Geometric distribution functions (ADF / DDF / distance-DF, CDF) — molrs-backed.

Thin Compute shells over molrs TRAVIS-parity geometric distributions. Each forwards verbatim to the Rust kernel and returns the molrs native result. All take (frames, groups) where groups is an integer index array selecting the atom tuples to histogram (pairs for distances, triplets for angles, quadruplets for dihedrals).

References¶

M. Brehm, B. Kirchner, J. Chem. Inf. Model. 51, 2007 (2011) — TRAVIS; radial/angular/dihedral and combined distribution functions.
M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, J. Chem. Phys. 152, 164105 (2020) — TRAVIS, current feature set.

AngleDistribution ¶

AngleDistribution(n_bins, min=0.0, max=180.0)

Bases: Compute

Angular distribution function (ADF) over (i, j, k) triplets.

The angle is taken at the middle atom j. The result's density_sin_corrected removes the trivial sin(theta) solid-angle weighting so a structureless distribution is flat.

Parameters¶

n_bins : int Number of angular bins. min, max : float, default 0.0 / 180.0 Angle range in degrees.

CombinedDistribution ¶

CombinedDistribution(axes)

Bases: Compute

Joint (combined) distribution over several geometric observables — the TRAVIS combined distribution function (CDF).

Parameters¶

axes : Sequence[tuple[str, int, float, float, bool]] One (kind, n_bins, min, max, sin_weight) per axis, where kind is "distance", "angle", or "dihedral".

Notes¶

Called as compute(frames, groups) with one index array per axis.

DihedralDistribution ¶

DihedralDistribution(n_bins, min=-180.0, max=180.0)

Bases: Compute

Dihedral distribution function (DDF) over (i, j, k, l) quadruplets.

Parameters¶

n_bins : int Number of angular bins. min, max : float, default -180.0 / 180.0 Dihedral range in degrees.

DistanceDistribution ¶

DistanceDistribution(n_bins, min, max)

Bases: Compute

Distance distribution function over (i, j) atom pairs.

Parameters¶

n_bins : int Number of histogram bins. min, max : float Lower and upper distance edges (Angstrom).

Spatial distribution function¶

spatial ¶

Spatial distribution function (SDF) — molrs-backed.

SpatialDistribution accumulates the density of target atoms on a body-fixed 3-D grid, after Kabsch-aligning each reference molecule to a template geometry. This is the orientation-resolved generalization of the radial distribution function: it reveals where around a molecule its neighbours sit, not just how far. Thin shell over the molrs TRAVIS-parity kernel; takes (frames).

References¶

I. M. Svishchev, P. G. Kusalik, J. Chem. Phys. 99, 3049 (1993); P. G. Kusalik, I. M. Svishchev, Science 265, 1219 (1994) — spatial distribution functions.
M. Brehm, B. Kirchner, J. Chem. Inf. Model. 51, 2007 (2011) — TRAVIS SDF.

SpatialDistribution ¶

SpatialDistribution(
    reference,
    template,
    target,
    n,
    extent,
    bulk_density=None,
    orientation_pairs=None,
)

Bases: Compute

Target-atom density on a molecule body-fixed grid (SDF).

Parameters¶

reference : Sequence[int] Atom indices whose instances define the body-fixed frame each step. template : ndarray (n_ref, 3) reference geometry the per-frame reference atoms are Kabsch-aligned to. target : Sequence[int] Atom indices whose density is accumulated on the grid. n : tuple[int, int, int] Grid resolution along each body-fixed axis. extent : tuple[float, float, float] Half-extent (Angstrom) of the grid along each axis. bulk_density : float, optional If given, the result also exposes g_sdf (the density normalized by the bulk number density). orientation_pairs : ndarray, optional Atom-index pairs whose mean orientation is accumulated per voxel.

Van Hove correlation¶

van_hove ¶

Van Hove correlation function G(r, t) — molrs-backed.

VanHove computes the self and distinct parts of the Van Hove correlation function, the time-resolved generalization of the radial distribution function: G_s(r, t) is the probability that a particle has moved a distance r in time t (it integrates to the self-diffusion picture), while G_d(r, t) tracks how the structure around a particle decorrelates. Thin shell over the molrs TRAVIS-parity kernel; takes (frames).

References¶

L. Van Hove, Phys. Rev. 95, 249 (1954) — the correlation function G(r, t).
M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, J. Chem. Phys. 152, 164105 (2020) — TRAVIS.

VanHove ¶

VanHove(n_rbins, r_max, lags, stride=1)

Bases: Compute

Van Hove correlation function G(r, t) (self + distinct parts).

Parameters¶

n_rbins : int Number of radial bins. r_max : float Upper edge of the radial grid (Angstrom). lags : Sequence[int] Time lags (in frames) at which to evaluate G(r, t). stride : int, default 1 Stride between time origins.

Notes¶

The result exposes r_centers, lags, g_self, and g_distinct (the latter only when has_distinct is true).

Reorientational correlations¶

reorientation ¶

Legendre reorientational correlation functions — molrs-backed.

LegendreReorientation computes the first- and second-order Legendre reorientational time-correlation functions of bond (or molecular) vectors,

.. math::

C_\ell(t) = \big\langle P_\ell\big(\mathbf u(0)\cdot\mathbf u(t)\big)\big\rangle,

with P_1(x) = x and P_2(x) = (3x^2 - 1)/2. C_2(t) is the quantity probed by NMR and dielectric relaxation; its decay time is the reorientational correlation time. Thin shell over the molrs TRAVIS-parity kernel; takes (frames, pairs) where pairs selects the vector endpoints.

References¶

B. J. Berne, R. Pecora, Dynamic Light Scattering, Wiley (1976) — reorientational correlation functions.
M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, J. Chem. Phys. 152, 164105 (2020) — TRAVIS.

LegendreReorientation ¶

LegendreReorientation(max_lag, stride=1)

Bases: Compute

First/second Legendre reorientational TCFs C_1(t), C_2(t).

Parameters¶

max_lag : int Longest lag (in frames) to evaluate. stride : int, default 1 Stride between time origins.

Notes¶

Called as compute(frames, pairs) where pairs is an integer array of (i, j) endpoints defining each tracked vector. The result exposes lags, c1, and c2.

Hydrogen bonds¶

hbond ¶

Hydrogen-bond detection — molrs-backed.

HBonds detects geometric hydrogen bonds per frame from explicit donor (D, H) pairs and acceptor atoms, using the :class:HBondCriterion (donor-acceptor distance and D-H...A angle). Pair these counts with the :class:~molpy.compute.persist.Persist survival analysis for hydrogen-bond lifetime dynamics. Thin shell over the molrs TRAVIS-parity kernel; takes (frames).

References¶

A. Luzar, D. Chandler, Nature 379, 55 (1996); Phys. Rev. Lett. 76, 928 (1996) — geometric hydrogen-bond criterion and bond-lifetime kinetics.
M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, J. Chem. Phys. 152, 164105 (2020) — TRAVIS.

HBonds ¶

HBonds(donors, acceptors, criterion=None)

Bases: Compute

Detect hydrogen bonds per frame from explicit donors and acceptors.

Parameters¶

donors : ndarray (n_donor, 2) integer array of (D, H) atom-index pairs. acceptors : ndarray Integer array of acceptor atom indices. criterion : HBondCriterion, optional Geometric criterion; defaults to Luzar-Chandler values (3.5 Angstrom donor-acceptor cutoff, 150 deg angle cutoff).

Notes¶

The result exposes per_frame (lists of (D, H, A, distance, angle)) and counts (hydrogen bonds per frame).

Radical Voronoi¶

voronoi ¶

Radical (Laguerre) Voronoi tessellation, domains, voids & integration — molrs-backed.

The radical (power / Laguerre) Voronoi tessellation partitions space by radius-weighted planes, so atoms of different size get cells proportional to their radii — the physically correct partition for polydisperse systems. On top of the tessellation:

:func:voronoi_domains merges cells sharing a label into connected domains (e.g. polar vs. apolar nanostructuring in ionic liquids).
:func:voronoi_voids aggregates the cells flagged as empty into void volumes.
:class:VoronoiIntegration integrates an electron density over the cells to yield per-molecule charges and dipoles (Voronoi/atomic-charge partitioning), the basis for predicting infrared spectra from ab initio MD.

Thin shells over the molrs TRAVIS-parity kernels.

References¶

B. J. Gellatly, J. L. Finney, J. Non-Cryst. Solids 50, 313 (1982) — radical (power) Voronoi tessellation.
M. Thomas, M. Brehm, B. Kirchner, Phys. Chem. Chem. Phys. 17, 3207 (2015) — Voronoi integration of the electron density for molecular dipoles.
M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, J. Chem. Phys. 152, 164105 (2020) — TRAVIS; domain and void analysis.

RadicalVoronoi ¶

RadicalVoronoi()

Bases: Compute

Radical (Laguerre / power) Voronoi tessellation under periodic boundaries.

Notes¶

Called as compute(positions, radii, box) and returns :class:VoronoiCells (per-cell volumes, faces, and neighbours).

VoronoiIntegration ¶

VoronoiIntegration()

Bases: Compute

Integrate an electron density over radical-Voronoi cells.

Aggregates a volumetric electron density into per-molecule charges and dipole moments, the Voronoi/atomic-charge partition used to derive infrared intensities from ab initio MD trajectories.

Notes¶

Called as compute(positions, radii, atomic_numbers, atom_to_mol, n_mol, grid, box) and returns the per-molecule moments.

Vibrational spectra¶

spectra ¶

Vibrational spectra from autocorrelation functions — molrs-backed.

These operators are spectral transforms: each takes a raw, precomputed autocorrelation function (ACF) sampled at dt_fs femtoseconds and returns the corresponding spectrum (Fourier transform with the appropriate quantum/temperature prefactor). They are the time-correlation route to vibrational spectroscopy from molecular-dynamics trajectories:

:class:PowerSpectrum — vibrational density of states from a velocity ACF.
:class:IRSpectrum — infrared absorption from a dipole-derivative (flux) ACF.
:class:RamanSpectrum — Raman from isotropic + anisotropic polarizability ACFs.
:class:VcdSpectrum — vibrational circular dichroism.
:class:RoaSpectrum — Raman optical activity.
:class:ResonanceRamanSpectrum — resonance Raman.

Thin shells over the molrs TRAVIS-parity kernels; called as compute(acf, dt_fs) (or compute(acf_iso, acf_aniso, dt_fs) for the polarizability-based spectra).

References¶

D. A. McQuarrie, Statistical Mechanics, Harper & Row (1976) — time-correlation functions and spectral densities.
M. Thomas, M. Brehm, R. Fligg, P. Vöhringer, B. Kirchner, Phys. Chem. Chem. Phys. 15, 6608 (2013) — IR and Raman spectra from AIMD via TCFs.
M. Brehm, M. Thomas, J. Phys. Chem. Lett. 8, 3409 (2017) — VCD, ROA and resonance Raman from MD (TRAVIS).
M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, J. Chem. Phys. 152, 164105 (2020) — TRAVIS.

IRSpectrum ¶

IRSpectrum()

Bases: Compute

Infrared absorption spectrum from a dipole-flux ACF.

Called as compute(acf, dt_fs).

PowerSpectrum ¶

PowerSpectrum()

Bases: Compute

Vibrational density of states (VDOS) from a velocity ACF.

Called as compute(acf, dt_fs).

RamanSpectrum ¶

RamanSpectrum(
    incident_frequency_cm1=0.0,
    temperature_k=0.0,
    averaged=False,
)

Bases: Compute

Raman spectrum from isotropic + anisotropic polarizability ACFs.

Parameters¶

incident_frequency_cm1 : float, default 0.0 Laser frequency for the Raman prefactor (0 disables it). temperature_k : float, default 0.0 Temperature for the Bose-Einstein prefactor (0 disables it). averaged : bool, default False Apply orientational averaging.

Called as compute(acf_iso, acf_aniso, dt_fs).

ResonanceRamanSpectrum ¶

ResonanceRamanSpectrum(
    incident_frequency_cm1=0.0,
    temperature_k=0.0,
    averaged=False,
)

Bases: Compute

Resonance-Raman spectrum from resonant isotropic + anisotropic ACFs.

Parameters mirror :class:RamanSpectrum. Called as compute(acf_iso, acf_aniso, dt_fs).

RoaSpectrum ¶

RoaSpectrum(
    incident_frequency_cm1=0.0,
    temperature_k=0.0,
    averaged=False,
)

Bases: Compute

Raman optical activity from isotropic + anisotropic ROA ACFs.

Parameters mirror :class:RamanSpectrum. Called as compute(acf_iso, acf_aniso, dt_fs).

VcdSpectrum ¶

VcdSpectrum()

Bases: Compute

Vibrational circular dichroism from an electric/magnetic cross-ACF.

Called as compute(acf, dt_fs).

Time series¶

time_series ¶

Time-series analysis operations for trajectory data.

This module provides utilities for computing time-correlation functions, mean squared displacements, and other time-series statistics commonly used in molecular dynamics trajectory analysis.

Adapted from the tame library (https://github.com/Roy-Kid/tame).

TimeAverage ¶

TimeAverage(shape, dtype=np.float64, dropnan='partial')

Compute running time average with NaN handling.

This class accumulates data over time and computes the average, with options for handling NaN values.

Parameters:

Name	Type	Description	Default
`shape`	`tuple[int, ...]`	Shape of data arrays to average	required
`dtype`	`dtype \| type`	Data type for accumulated arrays	`float64`
`dropnan`	`Literal['none', 'partial', 'all']`	How to handle NaN values: - 'none': Include NaN values in average (result may be NaN) - 'partial': Ignore individual NaN entries - 'all': Skip entire frame if any NaN is present	`'partial'`

Examples:

>>> avg = TimeAverage(shape=(10,), dropnan='partial')
>>> avg.update(np.array([1.0, 2.0, np.nan, 4.0]))
>>> avg.update(np.array([2.0, 3.0, 3.0, 5.0]))
>>> result = avg.get()  # [1.5, 2.5, 3.0, 4.5]

get ¶

get()

Get current time-averaged value.

Returns:

Type	Description
`NDArray`	Time-averaged data array

reset ¶

reset()

Reset accumulator to initial state.

update ¶

update(new_data)

Add new data to running average.

Parameters:

Name	Type	Description	Default
`new_data`	`NDArray`	New data array to include in average	required

TimeCache ¶

TimeCache(
    cache_size, shape, dtype=np.float64, default_val=np.nan
)

Cache previous N frames of trajectory data for correlation calculations.

Uses an in-place ring buffer for O(1) per update (no array allocation).

Parameters:

Name	Type	Description	Default
`cache_size`	`int`	Number of frames to cache (maximum time lag)	required
`shape`	`tuple[int, ...]`	Shape of data arrays to cache (e.g., (n_atoms, 3) for coordinates)	required
`dtype`	`dtype \| type`	Data type for cached arrays	`float64`
`default_val`	`float`	Default value to fill cache initially (default: NaN)	`nan`

Examples:

>>> cache = TimeCache(cache_size=100, shape=(10, 3))
>>> coords = np.random.randn(10, 3)
>>> cache.update(coords)
>>> cached_data = cache.get()  # Shape: (100, 10, 3)

cache `property` ¶

cache

Return data ordered newest-first (for backward compatibility).

get ¶

get()

Get cached data array, ordered newest-first.

Returns:

Type	Description
`NDArray`	Cached data with shape (cache_size, *data_shape)

reset ¶

reset()

Reset cache to initial state.

update ¶

update(new_data)

Add new frame to cache (O(1) in-place write).

Parameters:

Name	Type	Description	Default
`new_data`	`NDArray`	New data array to add (shape must match self.shape)	required

compute_acf ¶

compute_acf(data, cache_size, dropnan='partial')

Compute autocorrelation function over trajectory.

Calculates: _{i,t}

The particle dimension is averaged, and the time dimension is accumulated using a rolling cache to compute correlations at different time lags.

Parameters:

Name	Type	Description	Default
`data`	`NDArray`	Trajectory data with shape (n_frames, n_particles, n_dim)	required
`cache_size`	`int`	Maximum time lag (dt) to compute, in frames	required
`dropnan`	`Literal['none', 'partial', 'all']`	How to handle NaN values in averaging	`'partial'`

Returns:

Type	Description
`NDArray`	ACF array with shape (cache_size,) containing ACF at each time lag

Examples:

>>> # Velocity autocorrelation
>>> n_frames, n_particles = 1000, 100
>>> velocities = np.random.randn(n_frames, n_particles, 3)
>>> acf = compute_acf(velocities, cache_size=100)

compute_msd ¶

compute_msd(data, cache_size, dropnan='partial')

Compute mean squared displacement over trajectory.

Calculates: <(r_i(t+dt) - r_i(t))^2>_{i,t}

The particle dimension is averaged, and the time dimension is accumulated using a rolling cache to compute correlations at different time lags.

Parameters:

Name	Type	Description	Default
`data`	`NDArray`	Trajectory data with shape (n_frames, n_particles, n_dim)	required
`cache_size`	`int`	Maximum time lag (dt) to compute, in frames	required
`dropnan`	`Literal['none', 'partial', 'all']`	How to handle NaN values in averaging	`'partial'`

Returns:

Type	Description
`NDArray`	MSD array with shape (cache_size,) containing MSD at each time lag

Examples:

>>> # Simple 1D random walk
>>> n_frames, n_particles = 1000, 100
>>> positions = np.cumsum(np.random.randn(n_frames, n_particles, 1), axis=0)
>>> msd = compute_msd(positions, cache_size=100)
>>> # MSD should grow linearly with time for random walk

Workflow¶

workflow ¶

Lightweight DAG orchestration for Compute nodes.

Zero non-stdlib dependencies (uses graphlib.TopologicalSorter).

Workflow ¶

Workflow()

Compose Compute nodes into a DAG and execute them in topological order.

Parameters are bound by name: each node is called as node(**resolved) where resolved maps its parameter names to upstream results or externally-supplied values. The Workflow never inspects node signatures.

external_inputs `property` ¶

external_inputs

All source names that are not registered nodes.

nodes `property` ¶

nodes

Registered node names in insertion order.

add ¶

add(name, compute, inputs=None)

Register a compute node.

Parameters:

Name	Type	Description	Default
`name`	`str`	Unique node name.	required
`compute`		Callable object called as `compute(**resolved)`.	required
`inputs`	`dict[str, str] \| None`	Mapping from the node's parameter names to source names. Each source is either a registered node name or an external input name.	`None`

Returns:

Type	Description
`str`	name, for fluent chaining.

Raises:

Type	Description
`WorkflowDuplicateNodeError`	name already registered.
`WorkflowCycleError`	Adding this node creates a cycle.

predecessors ¶

predecessors(name)

Return the set of node-name predecessors for name.

External inputs are excluded.

run ¶

run(**externals)

Execute all nodes in topological order.

Parameters:

Name	Type	Description	Default
`**externals`	`Any`	Values for every external input.	`{}`

Returns:

Type	Description
`dict[str, Any]`	`{node_name: result}` for every registered node.

Raises:

Type	Description
`WorkflowMissingInputError`	One or more external inputs are absent.

topological_order ¶

topological_order()

Return nodes in topological execution order.

WorkflowCycleError ¶

Bases: WorkflowError

Adding this edge would create a cycle in the DAG.

WorkflowDuplicateNodeError ¶

Bases: WorkflowError

A node with the same name is already registered.

WorkflowError ¶

Bases: Exception

Base exception for Workflow errors.

WorkflowMissingInputError ¶

WorkflowMissingInputError(missing)

Bases: WorkflowError

Required external inputs were not provided to run().

Compute¶

Quick reference¶

Full API¶

Base¶

base ¶

Compute ¶

dump ¶

Result types¶

result ¶

ACFResult dataclass ¶

ConductivityResult dataclass ¶

DebyeFit dataclass ¶

epsilon ¶

DielectricResult dataclass ¶

fit_debye ¶

DielectricSusceptibilityResult dataclass ¶

to_dict ¶

JACFResult dataclass ¶

MCDResult dataclass ¶

OnsagerResult dataclass ¶

PMSDResult dataclass ¶

PersistResult dataclass ¶

Result dataclass ¶

to_dict ¶

SpectralResult dataclass ¶

TimeSeriesResult dataclass ¶

Dielectric¶

dielectric ¶

ACFAnalyzer ¶

DielectricSusceptibility ¶

IonicConductivity ¶

SpectralAnalyzer ¶

Mean displacement correlation¶

mcd ¶

MCDCompute ¶

Polarization MSD¶

pmsd ¶

PMSDCompute ¶

Onsager coefficients¶

onsager ¶

Onsager ¶

Current-ACF conductivity (Green–Kubo)¶

jacf ¶

JACF ¶

Pair persistence¶

persist ¶

Persist ¶

Radial distribution¶

rdf ¶

References¶

RDF ¶

Parameters¶

Notes¶

Mean-squared displacement¶

msd ¶

MSD ¶

Examples¶

Static structure factor¶

diffraction ¶

References¶

StaticStructureFactorDebye ¶

Parameters¶

Neighbor list¶

neighborlist ¶

References¶

NeighborList ¶

Parameters¶

Examples¶

Local & grid density¶

density ¶

References¶

GaussianDensity ¶

Parameters¶

LocalDensity ¶

Parameters¶

Bond-orientational order parameters¶

order ¶

References¶

Hexatic ¶

Parameters¶

ACFResult `dataclass` ¶

ConductivityResult `dataclass` ¶

DebyeFit `dataclass` ¶

DielectricResult `dataclass` ¶

DielectricSusceptibilityResult `dataclass` ¶

JACFResult `dataclass` ¶

MCDResult `dataclass` ¶

OnsagerResult `dataclass` ¶

PMSDResult `dataclass` ¶

PersistResult `dataclass` ¶

Result `dataclass` ¶

SpectralResult `dataclass` ¶

TimeSeriesResult `dataclass` ¶