seamm_md_analysis package#

Submodules#

seamm_md_analysis.correlate module#

Correlation functions of time series, computed with FFTs.

seamm_md_analysis.correlate.TENSOR_COMPONENTS = ('xx', 'yy', 'zz', 'xy', 'xz', 'yz')#: The order of the components of a symmetric tensor built from a vector.

seamm_md_analysis.correlate.acf_err(acf)[source]#

The standard error of an autocorrelation function (Bartlett’s formula).

The same as the estimate in statsmodels.tsa.stattools.acf.

Parameters:: acf (numpy.ndarray(n)) – The autocorrelation function.
Returns:: The standard error at each lag.
Return type:: numpy.ndarray(n)

seamm_md_analysis.correlate.autocorrelation_functions(X, m=None, demean=True)[source]#

The autocorrelation function of each series and its running integral.

This is the Green–Kubo input for independent series, such as the five components of the shear stress used for the viscosity (seamm_md_analysis.transport.shear_stress_components()).

Parameters:

X (numpy.ndarray(k, n)) – The k series, each of n samples.
m (int) – The number of lags to keep; all of them by default.
demean (bool = True) – Subtract the mean of each series first. When the true mean is known to be zero, as for the shear stress of a fluid, do not: the sample mean of a series of length T lowers every lag by about 2 sigma^2 tau / T, so the running integral falls linearly and the transport coefficient is underestimated.

Returns:

numpy.ndarray(k, m) – The autocorrelation functions, <X_i(t) X_i(0)>.
numpy.ndarray(k, m) – Their running integrals, in units of the sampling interval.

seamm_md_analysis.correlate.ccf_err(acf1, acf2)[source]#

The standard error of a cross-correlation function.

Parameters:: acf1, acf2 (numpy.ndarray(n)) – The autocorrelation functions of the two series.
Returns:: The standard error at each lag.
Return type:: numpy.ndarray(n)

seamm_md_analysis.correlate.ccovf(x, y, adjusted=True, demean=True)[source]#

The cross-covariance function of two series.

The same definition as statsmodels.tsa.stattools.ccovf:

\[c(k) = \frac{1}{n - k} \sum_t (x_{t+k} - \bar{x}) (y_t - \bar{y})\]

Parameters:

x, y (numpy.ndarray(n)) – The two series, of the same length.
adjusted (bool = True) – Divide each lag by the number of terms, n - k, rather than by n.
demean (bool = True) – Subtract the means first.

Returns:

The cross-covariance for lags 0 .. n-1.

Return type:

numpy.ndarray(n)

seamm_md_analysis.correlate.moving_average(a, n)[source]#

The moving average of a series over a window of n points.

Parameters:

a (numpy.ndarray) – The series.
n (int) – The width of the window.

Returns:

The averages.

Return type:

numpy.ndarray(len(a) - n + 1)

seamm_md_analysis.correlate.tensor_correlation_functions(J, m=None)[source]#

The correlation functions of a vector series and their running integrals.

This is the Green–Kubo input for a vector flux such as the heat flux: the six distinct components <J_a(t) J_b(0)> of the correlation tensor, in the order of TENSOR_COMPONENTS, and their cumulative trapezoidal integrals in units of the sampling interval.

Parameters:

J (numpy.ndarray(3, n)) – The x, y and z components of the series.
m (int) – The number of lags to keep; all of them by default.

Returns:

numpy.ndarray(6, m) – The correlation functions.
numpy.ndarray(6, m) – Their running integrals.

seamm_md_analysis.electric module#

Electrical properties from equilibrium MD: the static dielectric constant and the ionic conductivity.

Both are computed from the trajectory – unwrapped positions, so molecules stay whole, and fixed (forcefield) charges:

the dielectric constant from the fluctuations of the dipole of the cell, for Ewald or PPPM electrostatics with conducting (“tin-foil”) boundary conditions: ε = 1 + (⟨M²⟩ − ⟨M⟩²) / (3 V ε0 kB T). For an electrolyte, M is the dipole of the neutral molecules only: the ions carry a current, not a bounded dipole.
the ionic conductivity from the collective displacement of the ions’ charges (Einstein–Helfand), σ = lim ⟨‖P(t) − P(0)‖²⟩ / (6 V kB t), P = Σ qᵢ rᵢ over the ions; and the Nernst–Einstein estimate from the ions’ self-diffusion, σ_NE = e² / (V kB T) Σ N_s z_s² D_s. Their ratio is the ionicity.

Positions are in Å, charges in e, volumes in Å³, temperatures in K and times in ps; results are in SI units (S/m for conductivities, m²/s for diffusion coefficients).

seamm_md_analysis.electric.E = 1.602176634e-19#: The elementary charge, C

seamm_md_analysis.electric.EPS0 = 8.8541878128e-12#: The vacuum permittivity, F/m

seamm_md_analysis.electric.KB = 1.380649e-23#: The Boltzmann constant, J/K

seamm_md_analysis.electric.cell_dipole(positions, atoms, charges)[source]#

The dipole Σ qᵢ rᵢ of a set of atoms, e·Å.

Parameters:

positions (numpy.ndarray(n_atoms, 3)) – All the (unwrapped) positions of one frame, Å.
atoms (numpy.ndarray of int) – The atoms to include – whole molecules, so that the dipole of neutral molecules does not depend on the origin.
charges (numpy.ndarray) – The charges of those atoms, e.

seamm_md_analysis.electric.classify_molecules(charges, molecules, tolerance=0.05, scale=1.0)[source]#

Sort the molecules into neutral ones, ions, and others by their total charge.

Parameters:

charges (array_like) – The charge of each atom, e.
molecules ([[int]]) – The atoms of each molecule.
tolerance (float) – How far from an integer a total charge may be, in units of the formal charge.
scale (float) – The factor the charges of the ions were scaled by in the MD, as in the scaled-charge models of electrolytes (the parameter ion_charge_scale of md.h5md): an ion’s total charge is then scale times its formal charge.

Returns:

neutral, ions and other: lists of molecule indices; charge: the total charge of each molecule; z: the formal charge of each, the nearest integer to its charge divided by scale.

Return type:

dict

seamm_md_analysis.electric.collective_msd(P)[source]#

The mean square displacement of the collective charge displacement.

Parameters:: P (numpy.ndarray(n, 3)) – Σ qᵢ rᵢ over the ions in each frame, e·Å (unwrapped positions).
Returns:: ⟨‖P(t) − P(0)‖²⟩ over time origins, e² Å², for lags 0 .. n-1.
Return type:: numpy.ndarray(n)

seamm_md_analysis.electric.conductivity_from_msd(msd, dt, volume, temperature, start=0.1, end=0.5)[source]#

The Einstein–Helfand conductivity from the collective MSD.

The slope is fitted between start and end of the curve (fractions): the late part has few time origins and is noisy.

Returns:: sigma and stderr (S/m, the error of the fit), slope (e² Å² / ps) and the fitted line (xs, ys).
Return type:: dict

seamm_md_analysis.electric.dielectric_constant(M, volume, temperature, n_blocks=5, n_running=100)[source]#

The static dielectric constant from the dipole of the cell over time.

Parameters:

M (numpy.ndarray(n, 3)) – The dipole of the cell in each frame, e·Å.
volume (float) – The (mean) volume of the cell, Å³.
temperature (float) – The mean temperature, K.
n_blocks (int) – Consecutive blocks for the standard error; with fewer than 2, none.
n_running (int) – The number of points of the running estimate.

Returns:

epsilon and stderr (None without blocks), variance (e² Å²), mean_dipole (e·Å), running – the estimate from the first k frames, as (frame counts, values) – and n_blocks.

Return type:

dict

seamm_md_analysis.electric.dipole_correlation(M, dt, m=None)[source]#

The normalized autocorrelation of the dipole of the cell and its time.

Parameters:

M (numpy.ndarray(n, 3)) – The dipole in each frame.
dt (float) – The time between frames, ps.
m (int) – The lags to keep; a quarter of the frames by default.

Returns:

times (ps), acf (the average over x, y and z, normalized) and tau, its integral up to the first zero (ps).

Return type:

dict

seamm_md_analysis.electric.nernst_einstein(diffusion, counts, charges, volume, temperature)[source]#

The Nernst–Einstein conductivity, S/m.

Parameters:

diffusion ({str: float}) – The self-diffusion coefficient of each ionic species, m²/s.
counts ({str: int}) – The number of each.
charges ({str: int}) – The charge of each, in e.

seamm_md_analysis.h5md module#

Writing SEAMM’s H5MD files.

An MD step writes one H5MD file per run, md.h5md, from its code’s native output when the run ends. H5MDWriter holds the layout, so that every MD step writes the same thing; the layout is described in the documentation (“The H5MD file”).

The file is written as <name>.tmp and renamed when it is complete, so an interrupted conversion never leaves a damaged file under the real name.

class seamm_md_analysis.h5md.H5MDWriter(path, creator='SEAMM', version='', author=None, compress_particles=False)[source]#

Bases: object

Write one MD run as a SEAMM H5MD file.

Use it as a context manager. The file is renamed to its final name only if the block finishes without an exception; otherwise the temporary file is removed.

>>> with H5MDWriter("md.h5md", creator="SEAMM lammps_step", version="1") as w:
...     w.particles(ids=..., species=..., masses=..., cell=..., periodic=True)
...     w.observable("temperature", steps, times, T, "K", "fs")

Parameters:

path (str or pathlib.Path) – The final name of the file.
creator (str) – The program writing the file, e.g. “SEAMM lammps_step”.
version (str) – Its version.
author (str = None) – The user, by default from the environment.
compress_particles (bool = False) – Compress positions, velocities and forces with gzip. It gains little on coordinates, so it is off by default.

abort()[source]#: Discard the file.

box_appender(group='all', length_unit='Angstrom', time_unit='fs')[source]#

An appender for a box that changes with time, as in NPT.

Returns:: Call append(step, time, cell) for each frame, with the cell vectors as the rows of a 3×3 array.
Return type:: object

close()[source]#: Finish the file and give it its final name.

observable(name, steps, times, values, unit, time_unit='fs', average_window=1, kind=None, attrs=None)[source]#

Write a time series under observables/.

Parameters:

name (str) – The name, a standard one (see standard_observables) or, for the user’s extra properties, “user/<name>”.
steps (array_like of int) – The MD steps.
times (array_like of float) – The times, in time_unit.
values (array_like) – The values, with the first index running over the samples.
unit, time_unit (str) – The units, in H5MD syntax, as the code produced them.
average_window (int = 1) – The number of MD steps each value averages over; 1 for instantaneous values.
kind (str = None) – “correlation” for a correlation function accumulated by the code.
attrs (dict = None) – Extra attributes for the value dataset.

parameters(values, group='parameters')[source]#

Store the run’s parameters as attributes, nesting dicts as groups.

Parameters:: values (dict) – E.g. {“ensemble”: “NVT”, “temperature”: 298.15, “temperature_unit”: “K”, “timestep”: 1.0, …}.

particle_appender(name, n_particles, unit, time_unit='fs', group='all', attrs=None)[source]#

An appender for a per-particle property, frame by frame.

Positions, velocities and forces are written as float32, chunked by frame, so a long trajectory is streamed rather than held in memory.

Parameters:

name (str) – “position”, “velocity” or “force”.
n_particles (int) – The number of particles in the group.
unit, time_unit (str) – The units, in H5MD syntax.
group (str = “all”) – The particle group.
attrs (dict = None) – Extra attributes for the value dataset, e.g. {“seamm_unwrapped”: 1}.

Returns:

Call append(step, time, frame) with an n_particles × 3 array.

Return type:

object

particles(group='all', ids=None, species=None, masses=None, charges=None, cell=None, periodic=True, length_unit='Angstrom', mass_unit='amu', charge_unit='e')[source]#

Describe a group of particles and its box.

Parameters:

group (str = “all”) – The particle group: “all” for the atoms, “molecules” for the centres of mass of the molecules.
ids (array_like of int) – The SEAMM ids of the atoms, or the indices of the molecules.
species (array_like of int = None) – The atomic numbers.
masses (array_like of float = None) – The masses.
charges (array_like of float = None) – The charges, omitted when the model has none (MLFFs).
cell (3×3 array_like = None) – The cell vectors as rows, for a box that does not change. It is written for each frame of the positions when the file is closed, since MDAnalysis reads the box frame by frame. Leave it out for a non-periodic system, or when the box changes (see box_appender()).
periodic (bool = True) – Whether the system is periodic. A non-periodic box has no edges, because MDAnalysis treats any box with edges as periodic.
length_unit, mass_unit, charge_unit (str) – The units, in H5MD syntax.

provenance(values, segments=None)[source]#

Store SEAMM’s provenance under parameters/seamm.

Parameters:

values (dict) – The job and step ids, configuration, model chemistry, code and its version, package versions, hashes of model files, …
segments ([dict] = None) – One dict per segment of the run, with “first_step”, “last_step”, “wall_time” (s) and “node”.

seamm_md_analysis.h5md.standard_observables = {'density': (), 'hbond_count': (), 'hbond_energy': (), 'heat_flux': (3,), 'kinetic_energy': (), 'molecular_energy': (), 'pair_energy': (), 'potential_energy': (), 'pressure': (), 'pressure_tensor': (3, 3), 'species_mass_flux': (None, 3), 'temperature': (), 'total_dipole': (3,), 'total_energy': (), 'volume': ()}#: The standard observable names and their shapes; MD steps map their output onto these, and analysis steps ask for them.

seamm_md_analysis.io module#

Reading SEAMM’s H5MD files.

open_md() opens one run, given as one file or as the list of its segment files, and returns an MDData. Time series are read with h5py on their own time axis. MDAnalysis is used only for structural analysis, through MDData.universe(), because it attaches observables to frames by index and so misaligns any observable sampled more often than the frames.

class seamm_md_analysis.io.MDData(paths)[source]#

Bases: object

One MD run in SEAMM’s H5MD layout.

Use it as a context manager, or call close(). Several segment files are read as one run: time series and frames are joined in order, and a sample repeated at the boundary between segments is used once.

box(group='all', units='Å')[source]#: The box at every recorded time, as a TimeSeries of 3×3 matrices.

cell(group='all', units='Å')[source]#

The box of the first frame as a 3×3 array of cell vectors (rows).

Returns None for a non-periodic system.

charges(group='all')[source]#

close()[source]#

property creator#: The program that wrote the file and its version.

property first#: The first (or only) h5py file.

has_observable(name)[source]#

has_particles(name, group='all')[source]#: Whether a per-particle property, e.g. “position”, was sampled.

ids(group='all')[source]#

property ion_charge_scale#: The factor the charges of the ions were scaled by, 1 if they were not.

iter_frames(group='all', fields=('position',), chunk=100, units=None, time_units=None)[source]#

Stream frames of a trajectory in chunks.

Parameters:

group (str = “all”) – The particle group.
fields ([str]) – The properties to read, e.g. (“position”, “velocity”). They must have been sampled at the same steps.
chunk (int = 100) – The number of frames per chunk.
units ({str: str} = None) – Units to convert each field to, e.g. {“position”: “Å”}.
time_units (str = None) – Units for the times.

Yields:

(numpy.ndarray, numpy.ndarray, {str (numpy.ndarray})) – The steps, the times and the values of each field, with the shape (n, n_particles, 3), for each chunk of frames.

masses(group='all')[source]#

n_frames(name='position', group='all')[source]#: The number of frames of a per-particle property in all segments.

observable(name, units=None, time_units=None)[source]#

Read an observable on its own time axis.

Parameters:

name (str) – E.g. “temperature”, “pressure_tensor” or “user/<name>”.
units, time_units (str or pint.Unit = None) – Convert to these units; by default they are as the code wrote them.

Return type:

TimeSeries

observables()[source]#: The names of the observables, including “user/<name>”.

property parameters#: The run’s parameters, without SEAMM’s provenance.

particle_groups()[source]#: The particle groups, e.g. [“all”, “molecules”].

particle_series(name, group='molecules', units=None, time_units=None)[source]#

Read a per-particle property whole, e.g. the molecular centres of mass.

Use this only when it is small; iter_frames() streams large trajectories.

Returns:: values has the shape (n_frames, n_particles, 3).
Return type:: TimeSeries

property periodic#: Whether the system is periodic, from the box of the particles.

property provenance#: SEAMM’s provenance (parameters/seamm).

species(group='all')[source]#

universe(configuration=None, names=None)[source]#

An MDAnalysis Universe for the atoms, for structural analysis.

The topology comes from the file (species, masses and charges) and, if given, from the molsystem configuration (atom names, bonds, and molecules as residues). MDAnalysis is imported here only, so it is needed only by the steps that use it (pip install seamm-md-analysis[structure]).

Do not use its per-frame observables (ts.data): MDAnalysis attaches them to frames by index, ignoring their own steps.

Parameters:

configuration (molsystem._Configuration = None) – The configuration that was simulated, in the order of the atoms in the file.
names ([str] = None) – Atom names; by default the element symbols.

Return type:

MDAnalysis.Universe

class seamm_md_analysis.io.TimeSeries(name: str, values: ndarray, times: ndarray, steps: ndarray, units: object, time_units: object, average_window: int = 1, kind: str = None, attrs: dict = <factory>)[source]#

Bases: object

A time series or a per-particle trajectory read from an H5MD file.

Variables:

name (str) – The observable or property, e.g. “temperature” or “molecules/position”.
values (numpy.ndarray) – The values, with the first index running over the samples.
times (numpy.ndarray) – The times, in time_units.
steps (numpy.ndarray) – The MD steps.
units (pint.Unit) – The units of the values.
time_units (pint.Unit) – The units of the times.
average_window (int) – The number of MD steps each value averages over; 1 for instantaneous values.
kind (str or None) – “correlation” for a correlation function accumulated by the code.

attrs: dict#

average_window: int = 1#

property dt#: The sampling interval as a pint quantity, or None for fewer than 2.

kind: str = None#

name: str#

steps: ndarray#

time_units: object#

times: ndarray#

to(units=None, time_units=None)[source]#

A copy in other units.

Parameters:

units (str or pint.Unit = None) – The units of the values; None leaves them as they are.
time_units (str or pint.Unit = None) – The units of the times; None leaves them as they are.

Return type:

TimeSeries

units: object#

values: ndarray#

seamm_md_analysis.io.md_run_sources = ('the last MD run', 'the last MD run of each loop iteration')#: The choices for which earlier MD runs an analysis reads

seamm_md_analysis.io.open_md(paths)[source]#

Open an MD run’s H5MD file, or the list of its segment files.

Parameters:: paths (str, pathlib.Path or list) – One file, or the segments of one run in order.
Return type:: MDData

seamm_md_analysis.io.select_md_runs(records, source='the last MD run', name=None)[source]#

Choose the MD runs an analysis reads from the records in _md_runs.

Parameters:

records ([dict]) – The records MD steps added to _md_runs, oldest first.
source (str) – “the last MD run”, or “the last MD run of each loop iteration”: every run by the same MD step as the last one, which in a loop are the runs of the iterations, treated as independent replicas.
name (str = None) – Only consider runs that sampled this observable or per-particle property, e.g. “heat_flux” or “molecules/position”.

Returns:

The chosen records, oldest first; empty if there are none.

Return type:

[dict]

seamm_md_analysis.legacy module#

Readers for the text files the LAMMPS step has written so far.

These are the files that SEAMM’s H5MD output will replace: the !MolSSI trajectory 2.0 time series (*.trj), the !MolSSI vector_trajectory 2.0 per-molecule series (com_*.trj) and the LAMMPS dump files (*.dump_trj). The readers remain so that older jobs can still be analysed.

seamm_md_analysis.legacy.iter_dump_frames(path)[source]#

Stream the frames of a LAMMPS custom dump file, gzipped or not.

Parameters:: path (str or pathlib.Path) – The file; a name ending in “.gz” is read with gzip.
Yields:: dict – For each frame: “step”; “time” and “units” if the dump has them; “cell”, the 3×3 cell vectors as rows; “periodic”, from the boundary flags; and “columns”, {name: numpy.ndarray} for each per-atom column, sorted by the atom id. Columns that are not numbers (e.g. “element”) are strings.

seamm_md_analysis.legacy.iter_vector_frames(path)[source]#

Stream a !MolSSI vector_trajectory 2.0 file, such as molecular COMs.

Parameters:: path (str or pathlib.Path) – The file.
Yields:: dict, then (int, numpy.ndarray(n_rows, n_values)) for each frame – First the metadata (“code”, “type”, “dt”, “tunits”, “nsteps”, “fields”), then the step and the values of each frame.

seamm_md_analysis.legacy.lammps_box_to_cell(bounds, tilts=None)[source]#

The cell vectors (rows) from a LAMMPS dump’s box bounds.

Parameters:

bounds (3×2 array_like) – The (lo, hi) pairs from the ITEM: BOX BOUNDS lines. For a triclinic box these are LAMMPS’s bounding-box values.
tilts ((xy, xz, yz) = None) – The tilt factors of a triclinic box.

Return type:

numpy.ndarray(3, 3)

seamm_md_analysis.legacy.read_dump_trajectory(path)[source]#

Read a LAMMPS dump file with an id and three values per atom.

Parameters:

path (str or pathlib.Path) – The file.

Returns:

dict – “timesteps”, and when the dump has times, “times”, “dt” and “tunits” (“ps” for LAMMPS metal units, otherwise “fs”).
numpy.ndarray(n_frames, n_atoms, 3) – The values.

Note

Each frame looks like this:

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
256
ITEM: BOX BOUNDS pp pp pp
0.0000000000000000e+00 1.5900375866964442e+01
0.0000000000000000e+00 1.5900375866964442e+01
0.0000000000000000e+00 1.5900375866964442e+01
ITEM: ATOMS id xu yu zu
1 20.0993 14.7149 -5.20836
2 8.78739 12.8425 -25.9197
...

optionally preceded by ITEM: UNITS and ITEM: TIME.

seamm_md_analysis.legacy.read_trajectory(path)[source]#

Read a !MolSSI trajectory 2.0 time series, such as the heat flux.

The file starts with a header line holding JSON metadata, then a line of column names, then one line per sample. When the names do not include the first column, that column is the LAMMPS timestep and is returned as “tstep”.

Parameters:

path (str or pathlib.Path) – The file.

Returns:

dict – The metadata: “code”, “type”, “dt”, “tunits”, “nsteps”, and “fields”, the names of the columns.
{str (numpy.ndarray}) – The columns, by name.

seamm_md_analysis.legacy.read_vector_trajectory(path)[source]#

Read a !MolSSI vector_trajectory 2.0 file, such as molecular COMs.

Parameters:

path (str or pathlib.Path) – The file.

Returns:

dict – The metadata: “code”, “type”, “dt”, “tunits”, “nsteps”, and “fields”, the names of the values per row.
numpy.ndarray(n_frames, n_rows, n_values) – The values.

seamm_md_analysis.plots module#

Plots of correlation functions, Green–Kubo integrals, Helfand moments and mean square displacements, for seamm_util Figures (written as .graph files).

seamm_md_analysis.plots.add_component_traces(plot, x_axis, y_axis, rows, ts, name, yunits, err=None, fit=None, hover=None, labels=[('xx', 'red', 'rgba(255,0,0,0.1)'), ('yy', 'green', 'rgba(0,255,0,0.1)'), ('zz', 'blue', 'rgba(0,0,255,0.1)'), ('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)')], prefix='', fit_yunits=None)[source]#

Add the traces for each component of a curve: fit, error band and data.

Parameters:

plot (seamm_util.Plot) – The plot to add the traces to.
x_axis, y_axis – The axes of the plot.
rows (numpy.ndarray(n_components, m)) – The curves, one row per component.
ts ([float]) – The times, in ps, as a list.
name (str) – The name of the data traces, followed by the component label.
yunits (str) – The units of the data.
err (numpy.ndarray(n_components, m)) – The errors, drawn as a band hidden by default (optional).
fit ([{str: any}]) – For each component, a dict with the fitted curve “xs” and “ys” (optional).
hover (callable(label, fit) -> str) – The hover text of a fit trace.
labels ([(str, str, str)]) – The label, line color and band color of each component.
prefix (str) – Put in front of each trace name, for example the species.
fit_yunits (str) – The units of the fit traces, if different from yunits.

seamm_md_analysis.plots.add_helfand_trace(plot, x_axis, y_axis, species, M, ts, err=None, fit=None, labels=[('x', 'red', 'rgba(255,0,0,0.1)'), ('y', 'green', 'rgba(0,255,0,0.1)'), ('z', 'blue', 'rgba(0,0,255,0.1)'), ('', 'black', 'rgba(0,0,0,0.1)')])[source]#

Add the single-molecule Helfand moments of one species to a plot.

Parameters:

plot (seamm_util.Plot) – The plot to add the traces to.
x_axis, y_axis – The axes of the plot.
species (str) – The label of the species.
M (numpy.ndarray(m, 3 or 4)) – The Helfand moments, in Å^2.
ts ([float]) – The times, in ps.
err (numpy.ndarray(m, 3 or 4)) – The errors of the moments (optional).
fit ([{str: any}]) – The fits of the curves, with “D_s”, “scale”, “xs” and “ys” (optional).

seamm_md_analysis.plots.add_msd_trace(plot, x_axis, y_axis, species, msd, ts, err=None, fit=None, labels=[('x', 'red', 'rgba(255,0,0,0.1)'), ('y', 'green', 'rgba(0,255,0,0.1)'), ('z', 'blue', 'rgba(0,0,255,0.1)'), ('', 'black', 'rgba(0,0,0,0.1)')])[source]#

Add the mean square displacement of one species to a plot.

Parameters:

plot (seamm_util.Plot) – The plot to add the traces to.
x_axis, y_axis – The axes of the plot.
species (str) – The label of the species.
msd (numpy.ndarray(n, 3 or 4)) – The MSD in x, y and z, and optionally the total, in Å^2.
ts ([float]) – The times, in ps.
err (numpy.ndarray(n, 3 or 4)) – The errors of the MSD (optional).
fit ([{str: any}]) – The fits of the curves, with “D_s”, “scale”, “xs” and “ys” (optional).

seamm_md_analysis.plots.plot_GK_integrals(figure, x, ts, err=None, fit=None, labels=[('xx', 'red', 'rgba(255,0,0,0.1)'), ('yy', 'green', 'rgba(0,255,0,0.1)'), ('zz', 'blue', 'rgba(0,0,255,0.1)'), ('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)')], _range=None)[source]#

Plot the Green–Kubo running integrals of the thermal conductivity.

Parameters:

figure (seamm_util.Figure) – The figure that contains the plots.
x (numpy.ndarray(6, m)) – The running integrals, in W/m/K.
ts ([float]) – The times, in ps.
err (numpy.ndarray(6, m)) – Their standard errors (optional).
fit ([{str: any}]) – Fits to the curves (optional).
_range ([float]) – The range of the time axis.

seamm_md_analysis.plots.plot_components(figure, name, rows, ts, trace_name, ylabel, yunits, err=None, fit=None, labels=[('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)'), ('(xx-yy)/2', 'red', 'rgba(255,0,0,0.1)'), ('(yy-zz)/2', 'blue', 'rgba(0,0,255,0.1)'), ('', 'black', 'rgba(0,0,0,0.1)')], symbol='', fit_units=None, fit_yunits=None, _range=None)[source]#

Plot the components of a curve, such as running Green–Kubo integrals.

A general form of the heat-flux plots below, for any quantity and units.

Parameters:

figure (seamm_util.Figure) – The figure that contains the plots.
name (str) – The name of the plot in the figure.
rows (numpy.ndarray(n_components, m)) – The curves.
ts ([float]) – The times, in ps.
trace_name (str) – The name of the data traces, followed by the component label.
ylabel (str) – The label of the y axis.
yunits (str) – The units of the curves.
err (numpy.ndarray(n_components, m)) – Their standard errors (optional).
fit ([{str: any}]) – For each component, a dict with the fitted curve “xs” and “ys”, and the fitted “value” and its “stderr” for the hover text (optional).
labels ([(str, str, str)]) – The label, line color and band color of each component.
symbol (str) – The symbol of the fitted value in the hover text, e.g. “η”.
fit_units (str) – The units of the fitted value in the hover text; by default yunits.
fit_yunits (str) – The units of the fitted curves, if different from yunits.
_range ([float]) – The range of the time axis.

seamm_md_analysis.plots.plot_correlation_functions(figure, CCF, ts, err=None, fit=None, labels=[('xx', 'red', 'rgba(255,0,0,0.1)'), ('yy', 'green', 'rgba(0,255,0,0.1)'), ('zz', 'blue', 'rgba(0,0,255,0.1)'), ('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)')])[source]#

Plot the heat-flux correlation functions.

Parameters:

figure (seamm_util.Figure) – The figure that contains the plots.
CCF (numpy.ndarray(6, m)) – The correlation functions, in W^2/m^4.
ts ([float]) – The times, in ps.
err (numpy.ndarray(6, m)) – Their standard errors (optional).
fit ([{str: any}]) – Fits to the curves (optional).

seamm_md_analysis.plots.plot_helfand_moments(figure, M, ts, err=None, fit=None, labels=[('xx', 'red', 'rgba(255,0,0,0.1)'), ('yy', 'green', 'rgba(0,255,0,0.1)'), ('zz', 'blue', 'rgba(0,0,255,0.1)'), ('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)')])[source]#

Plot the Helfand moments of the thermal conductivity.

Parameters:

figure (seamm_util.Figure) – The figure that contains the plots.
M (numpy.ndarray(6, m)) – The Helfand moments, in W/m/K*ps.
ts ([float]) – The times, in ps.
err (numpy.ndarray(6, m)) – Their standard errors (optional).
fit ([{str: any}]) – Fits to the curves (optional).

seamm_md_analysis.plots.plot_helfand_slopes(figure, x, ts, err=None, fit=None, labels=[('xx', 'red', 'rgba(255,0,0,0.1)'), ('yy', 'green', 'rgba(0,255,0,0.1)'), ('zz', 'blue', 'rgba(0,0,255,0.1)'), ('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)')], _range=None)[source]#

Plot the local slopes of the Helfand moments of the thermal conductivity.

Parameters:

figure (seamm_util.Figure) – The figure that contains the plots.
x (numpy.ndarray(6, m)) – The slopes, in W/m/K.
ts ([float]) – The times, in ps.
err (numpy.ndarray(6, m)) – Their standard errors (optional).
fit ([{str: any}]) – Fits to the curves (optional).
_range ([float]) – The range of the time axis.

seamm_md_analysis.plots.shear_labels = [('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)'), ('(xx-yy)/2', 'red', 'rgba(255,0,0,0.1)'), ('(yy-zz)/2', 'blue', 'rgba(0,0,255,0.1)'), ('', 'black', 'rgba(0,0,0,0.1)')]#: Labels and colors of the five shear components used for the viscosity, in the order of seamm_md_analysis.transport.SHEAR_COMPONENTS, and their average.

seamm_md_analysis.plots.tensor_labels = [('xx', 'red', 'rgba(255,0,0,0.1)'), ('yy', 'green', 'rgba(0,255,0,0.1)'), ('zz', 'blue', 'rgba(0,0,255,0.1)'), ('xy', 'rgb(127,127,0)', 'rgba(127,127,0,0.1)'), ('xz', 'rgb(255,0,255)', 'rgba(255,0,255,0.1)'), ('yz', 'rgb(0,255,255)', 'rgba(0,255,255,0.1)')]#: Labels and colors (line, error band) of the six tensor components.

seamm_md_analysis.plots.vector_labels = [('x', 'red', 'rgba(255,0,0,0.1)'), ('y', 'green', 'rgba(0,255,0,0.1)'), ('z', 'blue', 'rgba(0,0,255,0.1)'), ('', 'black', 'rgba(0,0,0,0.1)')]#: Labels and colors of the x, y and z components and their total.

seamm_md_analysis.requirements module#

What an analysis needs from the MD, and what an MD step will provide.

An MD step describes what a run will write to its H5MD file with Provides, from its parameters, before anything runs. An analysis step describes what it needs with Requirements. check() compares the two, so that a missing or too coarse sampling is caught before the MD runs rather than after it.

The names are those of the H5MD file: an observable such as “heat_flux” or “pressure_tensor”, or a per-particle property as “<group>/<name>”, e.g. “molecules/position” or “all/velocity”.

Intervals are given in time, in MD steps, or both, since an MD step may know only one of them before it runs (the timestep may depend on the structure).

class seamm_md_analysis.requirements.Interval(time: object = None, steps: int = None)[source]#

Bases: object

A sampling interval, in time and/or in MD steps.

Variables:

time (pint.Quantity or None) – The interval in time, if known.
steps (int or None) – The interval in MD steps, if known.

steps: int = None#

time: object = None#

class seamm_md_analysis.requirements.Need(names: list, max_time: object = None, max_steps: int = None, why: str = '', sample_time: object = None)[source]#

Bases: object

One thing an analysis needs.

Variables:

names ([str]) – The acceptable names, in order of preference, e.g. [“molecules/position”, “all/position”].
max_time (pint.Quantity or str or None) – The longest acceptable sampling interval in time.
max_steps (int or None) – The longest acceptable sampling interval in MD steps.
why (str) – What it is needed for, for messages.
sample_time (pint.Quantity or str or None) – The interval to use when an analysis turns the sampling on, if finer than the longest acceptable one; by default max_time.

max_steps: int = None#

max_time: object = None#

names: list#

sample_time: object = None#

why: str = ''#

class seamm_md_analysis.requirements.Provides(data: dict = <factory>, timestep: object = None, step: str = '')[source]#

Bases: object

What an MD run will write to its H5MD file.

Variables:

data ({str: Interval}) – Each observable and per-particle property, with its sampling interval.
timestep (pint.Quantity or None) – The timestep, if known before the run.
step (str) – The MD step, for messages, e.g. “LAMMPS / NVT dynamics”.

data: dict#

interval(name)[source]#: The interval of a name, completed from the timestep if possible.

step: str = ''#

timestep: object = None#

class seamm_md_analysis.requirements.Requirements(needs: list = <factory>, periodic: bool = False, dynamics: bool = False)[source]#

Bases: object

What an analysis needs from the MD.

Variables:

needs ([Need]) – The data, each with its acceptable sampling.
periodic (bool) – Whether the system must be periodic.
dynamics (bool) – Whether the analysis uses the dynamics – transport properties such as the viscosity, diffusion and conductivities – so that a thermostat that perturbs the motion (Langevin) or none (NVE) should not be used.

dynamics: bool = False#

needs: list#

periodic: bool = False#

seamm_md_analysis.requirements.check(requirements, provides)[source]#

What the MD will not provide that the analysis needs.

A sampling interval that cannot be compared before the run (e.g. one given in steps when the timestep is not yet known) is not reported.

Parameters:

requirements (Requirements) – What the analysis needs.
provides (Provides) – What the MD run will write.

Returns:

A description of each gap; empty if there are none.

Return type:

[str]

seamm_md_analysis.state module#

The state the MD runs sampled, and whether they can be treated as isotropic.

Every analysis step reports the state its results belong to – temperature, pressure, volume and density, each with its error over the runs – because a property measured at the wrong pressure or density belongs to that state, not the one intended:

an NPT run’s mean density is the equilibrium density at its pressure;
an NVT run’s pressure is what it is; where the bulk modulus K is known (from the volume fluctuations of an NPT run, such as the equilibration before), the equilibrium density at the target pressure is estimated as ρ_eq ≈ ρ (1 + (P_target − P) / K).

A run is treated as isotropic – one value per property from the average of its components – only when its cell is cubic and no fix makes it anisotropic (walls, fields, deformation, flow, added forces).

seamm_md_analysis.state.AMU_PER_A3_IN_G_PER_ML = 1.6605390666#: amu/Å³ in g/mL

seamm_md_analysis.state.ANISOTROPIC_FIXES = ('wall', 'efield', 'deform', 'flow', 'addforce', 'aveforce', 'setforce', 'spring', 'smd', 'indent', 'gravity', 'move', 'oneway')#: The fix styles (or their beginnings) that make a run anisotropic

seamm_md_analysis.state.KB_ATM_A3 = 136.25946212681964#: Boltzmann’s constant in atm·Å³/K

class seamm_md_analysis.state.RunState(ensemble: str = '', temperature: float = None, pressure: float = None, volume: float = None, density: float = None, target_pressure: float = None)[source]#

Bases: object

The mean state of one run.

Variables:

ensemble (str) – “NVT”, “NPT”, “NVE”, …
density (temperature, pressure, volume,) – K, atm, Å³ and g/mL; None if not known.
target_pressure (float) – The barostat’s pressure for NPT, atm; None otherwise.

density: float = None#

ensemble: str = ''#

pressure: float = None#

target_pressure: float = None#

temperature: float = None#

volume: float = None#

class seamm_md_analysis.state.State(ensemble: str = '', temperature: Estimate = None, pressure: Estimate = None, volume: Estimate = None, density: Estimate = None, target_pressure: float = None, bulk_modulus: float = None, equilibrium_density: Estimate = None, notes: list = <factory>)[source]#

Bases: object

The state of a set of runs, each quantity with its error over the runs.

Variables:

ensemble (str)
density (temperature, pressure, volume,) – K, atm, Å³, g/mL.
target_pressure (float) – The pressure aimed at, atm, if known.
bulk_modulus (float) – K, atm, if known.
equilibrium_density (Estimate or None) – For runs at constant volume, the estimated density at the target pressure.
notes ([str]) – What the output should say about the state.

bulk_modulus: float = None#

density: Estimate = None#

ensemble: str = ''#

equilibrium_density: Estimate = None#

notes: list#

pressure: Estimate = None#

target_pressure: float = None#

temperature: Estimate = None#

volume: Estimate = None#

seamm_md_analysis.state.bulk_modulus(md)[source]#

The bulk modulus from the volume fluctuations of an NPT run, in atm.

K = k_B T ⟨V⟩ / (⟨V²⟩ − ⟨V⟩²), over the run after its equilibration. Returns None for a run without a volume series.

seamm_md_analysis.state.is_cubic(cell, tolerance=1e-06)[source]#: Whether a cell (3×3, vectors as rows) is cubic.

seamm_md_analysis.state.isotropy(md)[source]#

Whether a run can be treated as isotropic, and why.

Parameters:: md (seamm_md_analysis.MDData) – The run.
Returns:: Whether it is isotropic, and the reason, for the output.
Return type:: (bool, str)

seamm_md_analysis.state.run_state(md)[source]#

The mean temperature, pressure, volume and density of a run.

The means are over the run after its equilibration. The volume and density of a run at constant volume come from its cell and masses.

seamm_md_analysis.state.state_table(state)[source]#

The state as a table and notes, for the output.

Returns:: The table, with a line saying what the errors are, and the notes.
Return type:: str

seamm_md_analysis.state.summarize_state(runs, target_pressure=None, bulk_modulus=None)[source]#

The state over runs, with the equilibrium density where it can be estimated.

Parameters:

runs ([RunState]) – The state of each run, from run_state().
target_pressure (float) – The pressure aimed at, atm: by default that of NPT runs, else 1 atm when a bulk modulus is given.
bulk_modulus (float) – K, atm, for the equilibrium density of runs at constant volume, e.g. from bulk_modulus() of an NPT equilibration.

Return type:

State

seamm_md_analysis.stats module#

Statistics of time series and of replicas.

class seamm_md_analysis.stats.ComponentCheck(chi2: float = None, dof: int = 0, p: float = None, agree: bool = None, text: str = '')[source]#

Bases: object

Whether the components of a property agree within their errors.

Variables:

chi2 (float or None) – χ² of the component means about their average, weighted by their standard errors; None without errors.
dof (int) – The degrees of freedom, the number of components less one.
p (float or None) – The probability of a χ² at least this large if the components agree.
agree (bool or None) – Whether they agree, at p ≥ 0.01.
text (str) – One line for the output.

agree: bool = None#

chi2: float = None#

dof: int = 0#

p: float = None#

text: str = ''#

class seamm_md_analysis.stats.Estimate(value: float, stderr: float = None, n: int = 1, units: str = '', source: str = 'runs')[source]#

Bases: object

A value with its standard error, from n independent samples.

The samples are independent runs or blocks (source). Results stored for later steps carry the standard error; tables show the 95% confidence interval, from Student’s t, and say so.

Variables:

value (float)
stderr (float or None) – The standard error of the value; None if it cannot be estimated (one sample).
n (int) – The number of independent samples.
units (str)
source (str) – What the samples are: “runs” or “blocks”.

property ci95#: The half-width of the 95% confidence interval, or None.

describe_error()[source]#: What the ± in a table means, for a note under it.

formatted(precision=2)[source]#

The value and its 95% interval as text, e.g. (“0.72”, “0.19”).

A value without an error, or one whose error is negligible – the same in every sample, like the volume of runs at constant volume or the dipole of a rigid molecule – is given to six significant figures, without ±.

classmethod from_samples(values, units='', source='runs')[source]#: The mean of independent samples and its standard error.

n: int = 1#

source: str = 'runs'#

stderr: float = None#

units: str = ''#

value: float#

seamm_md_analysis.stats.block_average(series, n_blocks=10)[source]#

The mean of a series and its standard error from block averages.

The series is cut into n_blocks equal, consecutive blocks (dropping the remainder at the start). The blocks should be longer than the correlation time of the series for the error to be reliable.

Parameters:

series (numpy.ndarray(n) or numpy.ndarray(n, …)) – The series, with time along the first axis.
n_blocks (int = 10) – The number of blocks.

Returns:

mean, stderr – The mean over the blocks used and the standard error of the mean.
numpy.ndarray(n_blocks, …) – The block means.

seamm_md_analysis.stats.component_check(per_sample, names=None)[source]#

Compare the components of a property, for a one-line diagnostic.

Parameters:

per_sample (array_like(n_samples, n_components)) – The value of each component in each run or block.
names ([str]) – The components’ names, for the text.

Return type:

ComponentCheck

seamm_md_analysis.stats.detect_equilibration(series, nskip=1)[source]#

Find where a series has equilibrated, with pymbar.

Parameters:

series (numpy.ndarray(n)) – The series.
nskip (int = 1) – Test every nskip-th point as the start of the equilibrated region.

Returns:

t0 (int) – The index where the equilibrated region starts.
g (float) – The statistical inefficiency of the equilibrated region.
n_eff (float) – The effective number of uncorrelated samples in it.

seamm_md_analysis.stats.fmt_err(value, err, precision=2, as_float=False, fallback_decimals=3)[source]#

Format a value and its error with the error to precision digits.

Parameters:

value (float) – The value.
err (float) – Its error, for example twice the standard error for a 95% interval.
precision (int = 2) – The number of significant digits of the error.
as_float (bool = False) – Return rounded floats instead of strings.
fallback_decimals (int = 3) – The decimals of the value when the error cannot be used (zero, negative or not a number). The error is then “–“.

Returns:

The formatted value and error.

Return type:

(str, str) or (float, float)

seamm_md_analysis.stats.isotropic_average(per_sample)[source]#

One value for an isotropic property from its components.

The components are averaged within each run (or block), and the error is that of these averages over the runs: the components are correlated (e.g. the normal-stress differences of the viscosity), so their own spread would overstate the precision.

Parameters:: per_sample (array_like(n_samples, n_components)) – The value of each component in each run or block.
Return type:: Estimate

seamm_md_analysis.stats.mean_and_stderr(values, axis=0)[source]#

The mean of replicas and the standard error of that mean.

Parameters:

values (array-like) – The replica values along axis.
axis (int = 0) – The axis that runs over the replicas.

Returns:

numpy.ndarray or float – The mean.
numpy.ndarray or float – The standard error of the mean (sample standard deviation over sqrt(n)); NaN for a single replica.

seamm_md_analysis.stats.mean_and_stderr_correlated(series)[source]#

The mean of a correlated series and its standard error.

The error is the naive standard error inflated by the statistical inefficiency, sqrt(g * var / n).

Parameters:: series (numpy.ndarray(n)) – The series.
Returns:: The mean and its standard error.
Return type:: float, float

seamm_md_analysis.stats.statistical_inefficiency(series)[source]#

The statistical inefficiency g of a series, with pymbar.

The number of uncorrelated samples is about n / g.

Parameters:: series (numpy.ndarray(n)) – The series.
Returns:: The statistical inefficiency, at least 1.
Return type:: float

seamm_md_analysis.stats.t_factor(n, level=0.95)[source]#: The factor turning a standard error from n samples into a confidence interval: Student’s t with n - 1 degrees of freedom (2.78 for 5 samples at 95%, 1.96 for many). NaN for fewer than two samples.

seamm_md_analysis.step_support module#

Support for SEAMM analysis steps that read MD runs.

An analysis step either runs its own MD, in a sub-flowchart, as replicas (“mode A”), or analyses MD that ran earlier in the flowchart (“mode B”), found through the workspace variable _md_runs. These functions hold what both modes need and that every analysis step shares: finding the MD sub-steps in a flowchart, checking that they will sample what the analysis needs, turning sampling on, and finding the files to read.

They use only the methods of SEAMM’s nodes and flowcharts, so this library does not depend on SEAMM. An MD sub-step is any node with an md_provides() method.

seamm_md_analysis.step_support.NOT_CANONICAL = ('Berendsen', 'velocity rescaling')#: Thermostats that do not give the canonical fluctuations

seamm_md_analysis.step_support.STATE_PARAMETERS = ('T', 'T0', 'T1', 'P', 'P0', 'P1', 'Pinitial', 'Pfinal', 'temperature', 'pressure')#: The parameters of MD sub-steps that set the temperature or pressure

seamm_md_analysis.step_support.STATE_REPORTED = '_md_state_reported'#: The workspace variable recording the runs whose state was last printed

seamm_md_analysis.step_support.STOCHASTIC = ('Langevin', 'canonical sampling, langevin dynamics (csld)')#: Thermostats that perturb the motion of the particles

seamm_md_analysis.step_support.cell_widths(cell)[source]#

The perpendicular widths of a periodic cell, Å.

Parameters:: cell (array_like(3, 3)) – The cell vectors as rows.

seamm_md_analysis.step_support.check_md_flowchart(node, preceding, requirements, source)[source]#

Check, before anything runs, that the MD will provide what an analysis needs.

Besides the sampling, it checks where the analysis is relative to the loop around the MD, and the thermostat. Where it can, a problem comes with a fix.

Parameters:

node (seamm.Node) – The analysis step.
preceding ([seamm.Node]) – The steps that run before it, in order.
requirements (seamm_md_analysis.Requirements) – What the analysis needs.
source (str) – The value of its “md source” parameter.

Returns:

A description of each problem (a seamm.Problem where seamm has them); empty if none.

Return type:

[str]

seamm_md_analysis.step_support.each_iteration = 'the last MD run of each loop iteration'#

every run of the MD before the step, which in a Loop are its iterations’ runs, and otherwise the one run

Type:: The default source

seamm_md_analysis.step_support.md_inputs(node, P, name=None)[source]#

The H5MD files an analysis in mode B reads, one per independent run.

Parameters:

node (seamm.Node) – The analysis step, for the workspace variables and paths.
P (dict) – Its parameters, as values.
name (str = None) – For runs from _md_runs, only those that sampled this observable or per-particle property.

Return type:

[pathlib.Path]

seamm_md_analysis.step_support.md_substeps(nodes)[source]#

The MD sub-steps among nodes and inside their sub-flowcharts, in order.

Parameters:: nodes ([seamm.Node] or seamm.Flowchart) – The nodes to search, or a flowchart to search in execution order.
Returns:: The nodes that have md_provides.
Return type:: [seamm.Node]

seamm_md_analysis.step_support.molecule_extents(configuration)[source]#: The length of each molecule – the largest distance between two of its atoms, made whole across the cell by following its bonds – in Å, in the order of configuration.find_molecules(). The configuration is not changed.

seamm_md_analysis.step_support.npt_reference(node, paths)[source]#

The bulk modulus and target pressure from an NPT run before the analysed runs.

The last NPT run recorded in _md_runs before the first of the analysed runs – typically the equilibration before the production runs – gives the bulk modulus from its volume fluctuations, and the pressure it was run at.

Parameters:

node (seamm.Node) – The analysis step, for the workspace variables.
paths ([str or pathlib.Path]) – The md.h5md files of the analysed runs.

Returns:

The bulk modulus and the target pressure, in atm.

Return type:

(float, float) or (None, None)

seamm_md_analysis.step_support.preceding_viscosity(node)[source]#

The viscosity published by the last Viscosity step run, or None.

Parameters:: node (seamm.Node) – The step asking, for the workspace variables.
Returns:: The record written by publish_viscosity().
Return type:: dict or None

seamm_md_analysis.step_support.print_details(printer, text)[source]#

Write detailed output – the tables of each run, block or component – to the step’s own output (step.out) only, not to the job’s.

The step’s printer writes to step.out and, through the job’s printer, to job.out; this sends the text to the step.out handler alone.

Parameters:

printer (logging.Logger) – The step’s printer.
text (str) – The text, already formatted.

Returns:

Whether it was written (False if the step has no step.out handler).

Return type:

bool

seamm_md_analysis.step_support.publish_viscosity(node, value, stderr=None, T=None, method='Green-Kubo')[source]#

Record a viscosity for the steps after this one, such as Diffusivity.

Parameters:

node (seamm.Node) – The step that calculated it.
value (float) – The shear viscosity, in Pa s.
stderr (float = None) – Its standard error, in Pa s.
T (float = None) – The temperature, in K.
method (str = “Green-Kubo”) – How it was calculated.

Returns:

The record, also stored in the workspace variable _viscosity.

Return type:

dict

seamm_md_analysis.step_support.read_files = 'files'#: The choice that reads the files given by the user

seamm_md_analysis.step_support.report_state(node, printer, paths, indent=' ')[source]#

Print the state of the analysed runs and return it.

Parameters:

node (seamm.Node) – The analysis step.
printer (logging.Logger) – The step’s printer.
paths ([str or pathlib.Path]) – The md.h5md files of the runs.
indent (str) – The indentation of the table.

Returns:

None if the state could not be determined (it is noted, not raised).

Return type:

seamm_md_analysis.State or None

seamm_md_analysis.step_support.request_sampling(flowchart, requirements)[source]#

Turn on, in the last MD sub-step of a flowchart, what an analysis needs.

Parameters:

flowchart (seamm.Flowchart) – The analysis step’s sub-flowchart.
requirements (seamm_md_analysis.Requirements) – What the analysis needs.

Returns:

What was turned on, and what is still missing.

Return type:

([str], [str])

seamm_md_analysis.step_support.run_own_md = 'run the MD in the sub-flowchart'#: The choice by which an analysis step used to run its own MD, in a sub-flowchart; no longer offered (phase 8): replicas are a Loop around the MD

seamm_md_analysis.step_support.size_warnings(configuration, min_molecules=200, widths_per_molecule=3.0)[source]#

Whether a system is large enough for liquid properties.

Parameters:

configuration (molsystem._Configuration) – The structure about to be simulated.
min_molecules (int) – The fewest molecules for good statistics (~200-300).
widths_per_molecule (float) – The narrowest width of the cell should be at least this many times the length of the longest molecule.

Returns:

Warnings for the output; empty if the system is large enough.

Return type:

[str]

seamm_md_analysis.step_support.source_choices = ('the last MD run', 'the last MD run of each loop iteration', 'files')#: All choices of the source parameter

seamm_md_analysis.step_support.source_parameters(analysis='analysis')[source]#

The parameter definitions an analysis step adds for its source of MD.

Parameters:: analysis (str) – What the step calculates, for the help text, e.g. “the diffusivity”.
Returns:: {str – “md source”, “md files” and “number of blocks”.
Return type:: dict}

seamm_md_analysis.step_support.state_data(state)[source]#: The state as results for store_results, keyed as in state_results.

seamm_md_analysis.step_support.state_of_runs(node, paths)[source]#

The state of the analysed runs: T, P, V and ρ with their errors, and the equilibrium density where an NPT run before them gives the bulk modulus.

Parameters:

node (seamm.Node) – The analysis step.
paths ([str or pathlib.Path]) – The md.h5md files of the analysed runs.

Return type:

seamm_md_analysis.State

seamm_md_analysis.step_support.state_results = {'P': {'description': 'The mean pressure', 'dimensionality': 'scalar', 'type': 'float', 'units': 'atm'}, 'P,stderr': {'description': 'The standard error of the pressure', 'dimensionality': 'scalar', 'type': 'float', 'units': 'atm'}, 'T': {'description': 'The mean temperature', 'dimensionality': 'scalar', 'type': 'float', 'units': 'K'}, 'T,stderr': {'description': 'The standard error of the temperature', 'dimensionality': 'scalar', 'type': 'float', 'units': 'K'}, 'V': {'description': 'The mean volume', 'dimensionality': 'scalar', 'type': 'float', 'units': 'Å^3'}, 'bulk modulus': {'description': 'The bulk modulus used for the equilibrium density, from an NPT run', 'dimensionality': 'scalar', 'type': 'float', 'units': 'atm'}, 'density': {'description': 'The mean density', 'dimensionality': 'scalar', 'type': 'float', 'units': 'g/mL'}, 'density,stderr': {'description': 'The standard error of the density', 'dimensionality': 'scalar', 'type': 'float', 'units': 'g/mL'}, 'equilibrium density': {'description': 'The estimated density at the target pressure of runs at constant volume', 'dimensionality': 'scalar', 'type': 'float', 'units': 'g/mL'}, 'equilibrium density,stderr': {'description': 'The standard error of the estimated equilibrium density', 'dimensionality': 'scalar', 'type': 'float', 'units': 'g/mL'}}#: The results describing the state of the runs, for each analysis step’s metadata

seamm_md_analysis.step_support.viscosity_variable = '_viscosity'#: The workspace variable holding the result of the last Viscosity step

seamm_md_analysis.transport module#

Transport coefficients: Helfand moments, mean square displacements, and fits of Green–Kubo integrals and Helfand moments.

class seamm_md_analysis.transport.GreenKuboFit(value, stderr, a, a_err, tau, tau_err, xs, ys)#

Bases: tuple

A fit of a Green–Kubo integral by one or two rising exponentials.

a#: Alias for field number 2

a_err#: Alias for field number 3

stderr#: Alias for field number 1

tau#: Alias for field number 4

tau_err#: Alias for field number 5

value#: Alias for field number 0

xs#: Alias for field number 6

ys#: Alias for field number 7

class seamm_md_analysis.transport.LinearFit(slope, stderr, xs, ys)#

Bases: tuple

A straight-line fit: the slope, its standard error, and the fitted curve.

slope#: Alias for field number 0

stderr#: Alias for field number 1

xs#: Alias for field number 2

ys#: Alias for field number 3

seamm_md_analysis.transport.SHEAR_COMPONENTS = ('xy', 'xz', 'yz', '(xx-yy)/2', '(yy-zz)/2')#: The five independent components of the shear stress of an isotropic fluid.

seamm_md_analysis.transport.YEH_HUMMER_XI = 2.837297#: The Yeh–Hummer constant for a cubic cell (Yeh and Hummer, 2004).

seamm_md_analysis.transport.axb(x, a, b)[source]#: A straight line, a * x + b.

seamm_md_analysis.transport.collective_helfand_moments(J, m=None)[source]#

The Helfand moments of a collective vector flux, such as the heat flux.

\[M_{ab}(k) = \left\langle \int_0^{k} J_a \, dt \int_0^{k} J_b \, dt \right\rangle\]

averaged over n - m time origins, in units of the sampling interval squared, for the six components in seamm_md_analysis.correlate.TENSOR_COMPONENTS.

Parameters:

J (numpy.ndarray(3, n)) – The x, y and z components of the flux.
m (int) – The number of lags; by default min(n // 20, 10000).

Returns:

The Helfand moments.

Return type:

numpy.ndarray(6, m)

seamm_md_analysis.transport.component_helfand_moments(X, m=None, demean=True)[source]#

The Helfand moment of each of several independent series.

\[M_i(k) = \left\langle \left( \int_0^{k} X_i \, dt \right)^2 \right\rangle\]

averaged over n - m time origins, in units of the sampling interval squared. For the shear stress (shear_stress_components()) this gives the Helfand moments of the viscosity.

Parameters:

X (numpy.ndarray(k, n)) – The k series, each of n samples.
m (int) – The number of lags; by default min(n // 20, 10000).
demean (bool = True) – Subtract the mean of each series first. A non-zero mean adds a term growing as the square of the time, which is not part of the transport coefficient. When the true mean is known to be zero, as for the shear stress of a fluid, do not: subtracting the sample mean of a finite series biases the moments low.

Returns:

The Helfand moments.

Return type:

numpy.ndarray(k, m)

seamm_md_analysis.transport.effective_correlation_functions(C, L)[source]#

The derivative of the Schur complement of the integrals: the correlation function of the heat flux less the part carried by the species, with the coefficients L_JJ^-1 L_Je taken as fixed at each lag.

Parameters:: C, L (numpy.ndarray(…, m, 1 + k, 1 + k)) – The correlation functions and their integrals, from flux_correlation_matrices() or averages of them.
Return type:: numpy.ndarray(…, m)

seamm_md_analysis.transport.exp_1(t, a, tau)[source]#: A single rising exponential, a * (1 - exp(-t / tau)).

seamm_md_analysis.transport.exp_2(t, a1, tau1, a2, tau2)[source]#: The sum of two rising exponentials.

seamm_md_analysis.transport.fit_green_kubo_integral(y, xs, sigma=None)[source]#

Fit a Green–Kubo running integral by one or two rising exponentials.

The plateau of the integral, a = sum of the amplitudes, is the transport coefficient. A single exponential a (1 - exp(-t / tau)) is fitted first over the initial rise; a second exponential is then added over tau / 2 to 20 tau. The second exponential is kept only if both amplitudes are positive, the two time constants differ, the slower one is short enough for the plateau to be reached within the fitted data (3 tau), and the value is no more than twice that of the single exponential. Otherwise the single exponential is used.

The standard error of the value includes the covariance of the two amplitudes, which are strongly correlated.

Parameters:

y (numpy.ndarray(n)) – The running integral.
xs (numpy.ndarray(n) or [float]) – The times.
sigma (numpy.ndarray(n)) – The standard error of y (optional). Without it, or if any value is zero or not finite (as for the spread of a single run), the fit is unweighted and the errors come from the residuals.

Returns:

The plateau value and its standard error, the amplitudes and time constants with their errors, and the fitted curve.

Return type:

GreenKuboFit

seamm_md_analysis.transport.fit_helfand_moment(y, xs, sigma=None, start=1, logger=<Logger seamm_md_analysis.transport (WARNING)>)[source]#

A straight-line fit to a Helfand moment, from time start to the end.

Parameters:

y (numpy.ndarray(n)) – The Helfand moment.
xs (numpy.ndarray(n) or [float]) – The times.
sigma (numpy.ndarray(n)) – The standard error of y, to weight the fit (optional). Without it, or if any value is zero or not finite, the fit is unweighted.
start (float = 1) – The time, in the units of xs, before which the curve is ignored.

Returns:

The slope, its standard error, and the fitted line.

Return type:

LinearFit

seamm_md_analysis.transport.fit_slope(y, xs, sigma=None, start=0.1, end=0.9)[source]#

A straight-line fit to the part of a curve between two fractions.

Parameters:

y (numpy.ndarray(n)) – The curve, for example an MSD or a Helfand moment.
xs (numpy.ndarray(n) or [float]) – The times.
sigma (numpy.ndarray(n)) – The standard error of y, to weight the fit (optional). Without it, or if any value is zero or not finite, the fit is unweighted.
start (float = 0.1) – The fraction of the curve to skip at the beginning.
end (float = 0.9) – The fraction of the curve at which to stop.

Returns:

The slope, its standard error, and the fitted line over the range.

Return type:

LinearFit

seamm_md_analysis.transport.flux_correlation_matrices(J, fluxes, m=None)[source]#

The correlation functions of the heat flux and the species’ mass fluxes, and their running integrals, for each direction.

For direction d the series are the heat flux J_d and the mass fluxes of the species in the same direction. The correlation functions are symmetrized, (<a(t) b(0)> + <b(t) a(0)>) / 2, as the Onsager coefficients of two fluxes odd in time are. Being linear in the data, the matrices of several runs can be averaged before schur_complement() is applied, which is much less noisy than the complement of each run.

Parameters:

J (numpy.ndarray(3, n)) – The x, y and z components of the heat flux.
fluxes (numpy.ndarray(k, 3, n)) – The mass fluxes of k species, with k one less than the number of species.
m (int) – The number of lags to keep; all of them by default.

Returns:

numpy.ndarray(3, m, 1 + k, 1 + k) – The correlation functions, the heat flux first.
numpy.ndarray(3, m, 1 + k, 1 + k) – Their running integrals, in units of the sampling interval.

seamm_md_analysis.transport.flux_helfand_matrices(J, fluxes, m=None)[source]#

The Helfand moments of the heat flux and the species’ mass fluxes, for each direction.

The matrix form of collective_helfand_moments(): for direction d, <G_a(k) G_b(k)> for the integrals G of the heat flux J_d and of the mass fluxes of the species in the same direction. Like flux_correlation_matrices(), the matrices of several runs can be averaged before schur_complement() is applied.

Parameters:

J (numpy.ndarray(3, n)) – The x, y and z components of the heat flux.
fluxes (numpy.ndarray(k, 3, n)) – The mass fluxes of k species, with k one less than the number of species.
m (int) – The number of lags; by default min(n // 20, 10000).

Returns:

The Helfand moments, in units of the sampling interval squared, the heat flux first.

Return type:

numpy.ndarray(3, m, 1 + k, 1 + k)

seamm_md_analysis.transport.helfand_slope(y, xs, sigma)[source]#

The local slope of a Helfand moment, by central differences.

The slope of a Helfand moment is the running Green–Kubo integral, so it is zero at the origin.

Parameters:

y (numpy.ndarray(n)) – The Helfand moment.
xs (numpy.ndarray(n) or [float]) – The times.
sigma (numpy.ndarray(n)) – The standard error of y.

Returns:

numpy.ndarray(n) – The slope: zero at the first point, a central difference in between and a backward difference at the last point.
xs – The times, unchanged.
numpy.ndarray(n) – The error of the slope, from the sum of the errors of the two points.

seamm_md_analysis.transport.msd_by_species(xyz, species, average=True)[source]#

The mean square displacement of each species, averaged over molecules.

Parameters:

xyz (numpy.ndarray(n, n_molecules, 3)) – The unwrapped positions at each step.
species ([[int]]) – The indices of the molecules of each species.
average (bool = True) – Add a fourth column with the total MSD (the sum over x, y and z).

Returns:

[numpy.ndarray(n, 3 or 4)] – The MSD of each species.
[numpy.ndarray(n, 3 or 4)] – The standard deviation over the molecules.

seamm_md_analysis.transport.multicomponent_correlation_functions(J, fluxes, m=None)[source]#

The Green–Kubo input for the heat flux of a mixture.

In a mixture the heat flux is correlated with the interdiffusion of the species, and the Green–Kubo integral of its autocorrelation function gives too large a thermal conductivity. The thermal conductivity comes instead from the Schur complement of the matrix of Onsager coefficients (Bertossa, Grasselli, Ercole and Baroni, Phys. Rev. Lett. 122, 255901 (2019)):

\[\kappa \propto L_{ee} - L_{eJ} L_{JJ}^{-1} L_{Je}, \qquad L_{ab}(t) = \int_0^t \langle a(\tau) b(0) \rangle \, d\tau\]

with e the heat flux and J the mass fluxes of all the species but one (the last is fixed by the total momentum, which is zero). Their scale cancels, so any units, or the momenta of the species, will do. The complement is taken for each lag and each direction, coupling the fluxes in the same direction, as for an isotropic fluid.

This is for one run. The complement of a short run is noisy, since L_JJ is small and poorly known; for several runs or blocks, average their flux_correlation_matrices() and take the complement of the average.

Parameters:

J (numpy.ndarray(3, n)) – The x, y and z components of the heat flux.
fluxes (numpy.ndarray(k, 3, n)) – The mass fluxes of k species, with k one less than the number of species.
m (int) – The number of lags to keep; all of them by default.

Returns:

numpy.ndarray(6, m) – The correlation functions, in the order of seamm_md_analysis.correlate.TENSOR_COMPONENTS. The diagonal ones are the derivatives of the corrected integrals (effective_correlation_functions()); the off-diagonal ones are those of the heat flux alone, which are zero either way for an isotropic system.
numpy.ndarray(6, m) – Their running integrals, in units of the sampling interval.

seamm_md_analysis.transport.multicomponent_helfand_moments(J, fluxes, m=None)[source]#

The Helfand moments of the heat flux of a mixture.

The Helfand-moment form of multicomponent_correlation_functions(): the Schur complement, for each lag and direction, of the matrix of Helfand moments of the heat flux and the mass fluxes of all the species but one. At long times the moments grow as twice the time times the Onsager coefficients, so the slope of the complement gives the thermal conductivity of the mixture. As there, for several runs average their flux_helfand_matrices() and take the complement of the average.

Parameters:

J (numpy.ndarray(3, n)) – The x, y and z components of the heat flux.
fluxes (numpy.ndarray(k, 3, n)) – The mass fluxes of k species, with k one less than the number of species.
m (int) – The number of lags; by default min(n // 20, 10000).

Returns:

The Helfand moments, in units of the sampling interval squared, in the order of seamm_md_analysis.correlate.TENSOR_COMPONENTS. The off-diagonal ones are those of the heat flux alone, as from collective_helfand_moments().

Return type:

numpy.ndarray(6, m)

seamm_md_analysis.transport.schur_complement(L, coefficients=False)[source]#

L_ee - L_eJ L_JJ^-1 L_Je for each lag: the thermal conductivity of a mixture.

The matrices are of the Green-Kubo integrals or Helfand moments of the heat flux e, first, and the mass fluxes J of all the species but one (see multicomponent_correlation_functions()).

Parameters:

L (numpy.ndarray(…, m, 1 + k, 1 + k)) – The matrices for each lag, with the heat flux first.
coefficients (bool = False) – Also return the coefficients L_JJ^-1 L_Je.

Returns:

numpy.ndarray(…, m) – The Schur complement; L_ee at lags where L_JJ is singular, such as lag 0.
numpy.ndarray(…, m, k) – The coefficients, if asked for; zero where L_JJ is singular.

seamm_md_analysis.transport.shear_stress_components(pressure_tensor)[source]#

The five independent shear components of a pressure tensor.

For an isotropic fluid, Pxy, Pxz, Pyz, (Pxx - Pyy)/2 and (Pyy - Pzz)/2 all have the same autocorrelation function, so each gives the shear viscosity by Green–Kubo, and averaging the five improves the statistics over using only the three off-diagonal components (Daivis and Evans, 1994). The order is that of SHEAR_COMPONENTS.

Parameters:: pressure_tensor (numpy.ndarray(n, 3, 3)) – The pressure tensor at each sample.
Returns:: The five components.
Return type:: numpy.ndarray(5, n)

seamm_md_analysis.transport.species_helfand_moments(v, species, m=None, max_bytes=268435456)[source]#

Single-molecule Helfand moments from velocities, for each species.

For each molecule and direction, the square of the integral of its velocity over a window of k steps, averaged over n - m time origins, then averaged over the molecules of the species. A fourth column holds the total, the sum over x, y and z.

The molecules are taken a batch at a time, so the memory used is bounded by max_bytes whatever the length of the run: the integrals and the FFTs need several times the velocities of the molecules they work on, which for every molecule of a long run at once was ~8 times the whole array. The batches’ means and spreads are combined exactly.

Parameters:

v (numpy.ndarray(n, n_molecules, 3)) – The velocities at each step.
species ([[int]]) – The indices of the molecules of each species.
m (int) – The number of lags; by default min(n // 2, 5000).
max_bytes (int) – About how much memory a batch of molecules may use.

Returns:

[numpy.ndarray(m, 4)] – The Helfand moments of each species, in units of velocity times the sampling interval, squared.
[numpy.ndarray(m, 4)] – The standard deviation over the molecules.

seamm_md_analysis.transport.squared_displacement_sums(x)[source]#

Sum over all time origins of the squared displacement at each lag.

\[S(k) = \sum_{i=0}^{n-1-k} (x[i+k] - x[i])^2\]

computed with FFTs for every column.

Parameters:: x (numpy.ndarray(n, …)) – The positions, with time along the first axis.
Returns:: S for lags 0 .. n-1; divide by n - k for the mean square displacement.
Return type:: numpy.ndarray(n, …)

seamm_md_analysis.transport.windowed_displacement_products(Ga, Gb, m, n_origins)[source]#

Sum over time origins of the products of displacements.

\[S(k) = \sum_{i=0}^{N_0 - 1} (G_a[i+k] - G_a[i]) (G_b[i+k] - G_b[i])\]

for lags k = 0 .. m-1 and a fixed set of N_0 origins, computed with FFTs. With G the running integral of a flux, S/N_0 is the Helfand moment.

Each series is split into a straight line and the rest before the FFTs, and the line’s contribution is added back exactly. A flux with a non-zero mean has a running integral that drifts linearly, and without this the rounding error would grow as the square of the length of the series.

Parameters:

Ga, Gb (numpy.ndarray(n, …)) – The series, with time along the first axis. They may be the same array.
m (int) – The number of lags.
n_origins (int) – The number of time origins, N_0; n_origins + m - 1 must not exceed n.

Returns:

S for each lag and column.

Return type:

numpy.ndarray(m, …)

seamm_md_analysis.transport.yeh_hummer_correction(T, L, viscosity)[source]#

The Yeh–Hummer finite-size correction to a diffusion constant.

\[\Delta D = \frac{\xi k_B T}{6 \pi \eta L}\]

with xi = 2.837297 for a cubic cell. Add it to the diffusion constant computed in the periodic cell.

Parameters:

T (float) – The temperature, in K.
L (float) – The length of the cubic cell, in m.
viscosity (float) – The shear viscosity, in Pa s.

Returns:

The correction, in m^2/s.

Return type:

float

seamm_md_analysis.units module#

Unit strings in H5MD files and their translation to and from pint.

H5MD writes a unit as space-separated factors, each a unit symbol with an optional integer power, with negative powers for the denominator: kcal mol-1, Angstrom ps-1, g cm-3. MDAnalysis accepts only its own spellings for lengths, times, velocities and forces (Angstrom, Angstrom ps-1, kcal mol-1 Angstrom-1), so SEAMM writes those, and this module translates between them and pint.

pint’s application registry is used unless a registry is given, so that the quantities mix with SEAMM’s, which sets its registry as the application registry.

seamm_md_analysis.units.to_h5md(unit, ureg=None)[source]#

The H5MD unit string for a pint unit or unit expression.

Factors with positive powers come first, then the negative ones, with “mol” first among those, so that forces are written as “kcal mol-1 Angstrom-1” as MDAnalysis expects.

Parameters:

unit (str, pint.Unit or pint.Quantity) – The unit, e.g. “Å/fs” or ureg.kcal / ureg.mol.
ureg (pint.UnitRegistry = None) – The registry, by default pint’s application registry.

Returns:

The unit in H5MD syntax, or “” for a dimensionless unit.

Return type:

str

seamm_md_analysis.units.to_pint(unit, ureg=None)[source]#

The pint unit for an H5MD unit string.

Parameters:

unit (str or None) – The unit in H5MD syntax, e.g. “kcal mol-1”. None, “” and “1” are dimensionless.
ureg (pint.UnitRegistry = None) – The registry, by default pint’s application registry.

Return type:

pint.Unit

Module contents#

Analysis of molecular dynamics output for SEAMM.

Time-series statistics, correlation functions, Green–Kubo and Helfand transport coefficients, plots, SEAMM’s H5MD files, and readers for older MD output.