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)[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.

Returns:

neutral, ions and other: lists of molecule indices; charge: the total charge of each molecule; z: the nearest integer charge of each.

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), '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]#

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)[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.

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.stats module#

Statistics of time series and of replicas.

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.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.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.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.check_md_flowchart(node, preceding, requirements, source)[source]#

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

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; empty if none.

Return type:

[str]

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.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.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.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 of the analysis step’s source parameter that runs its own MD

seamm_md_analysis.step_support.source_choices = ('run the MD in the sub-flowchart', '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.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.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.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.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)[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.

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).

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.