Module contents

mdevaluate.open(directory='', topology='*.tpr', trajectory='*.xtc', index_file=None, cached=False, reindex=True, verbose=True, nojump=False, ignore_index_timestamps=False)

Open a simulation from a directory.

Parameters:

• directory – Directory of the simulation.
• topology (opt.) – Descriptor of the topology file (tpr or gro). By default a tpr file is used, if there is exactly one in the directoy.
• trajectory (opt.) – Descriptor of the trajectory (xtc file).
• index_file (opt.) – An index file that will be loaded alongside with the topology.
• cached (opt.) – If the trajectory reader should be cached. Can be True, an integer or None. If this is True maxsize is 128, otherwise this is used as maxsize for the cache, None means infinite cache (this is a potential memory leak!).
• reindex (opt.) – Regenerate the xtc-index if necessary.
• verbose (opt.) – Be verbose about the opened files.
• nojump (opt.) – If nojump matrixes should be generated. They will alwyas be loaded if present
• ignore_index_timestamps (opt.) – Ignore timestamps in xtc index file. If True the index file will be loaded regardless of the timestamp

Returns:

A Coordinate object of the simulation.

Example

Open a simulation located in ‘/path/to/sim’, where the trajectory is located in a sub-directory ‘/path/to/sim/out’ and named for Example ‘nojump_traj.xtc’. All read frames will be cached in memory.

>>> open('/path/to/sim', trajectory='out/nojump*.xtc', cached=None)

The file descriptors can use unix style pathname expansion to define the filenames. The default patterns use the recursive placeholder ** which matches the base or any subdirctory, thus files in subdirectories with matching file type will be found too. For example: ‘out/nojump*.xtc’ would match xtc files in a subdirectory out that start with nojump and end with .xtc.

For more details see: https://docs.python.org/3/library/glob.html

mdevaluate.open_energy(energyfile)

Open an energy file with EnergyReader.

mdevaluate.autosave

mdevaluate.autosave.autosave_data(nargs, kwargs_keys=None, version=None)

Enable autosaving of results for a function.

Parameters:

• nargs – Number of args which are relevant for the calculation.
• kwargs_keys (opt.) – List of keyword arguments which are relevant for the calculation.
• version (opt.) – An optional version number of the decorated function, which replaces the checksum of the function code, hence the checksum does not depend on the function code.

mdevaluate.autosave.disable()

Disable autosave.

class mdevaluate.autosave.disabled

Bases: object

A context manager that disbales the autosave module within its context.

Example

import mdevaluate as md md.autosave.enable(‘data’) with md.autosave.disabled():

# Autosave functionality is disabled within this context. md.correlation.shifted_correlation(

…

)

# After the context is exited, autosave will work as before.

mdevaluate.autosave.enable(dir, load_data=True, verbose=True)

Enable auto saving results of functions decorated with autosave_data().

Parameters:

• dir – Directory where the data should be saved.
• load_data (opt., bool) – If data should also be loaded.

mdevaluate.autosave.get_directory(reader)

Get the autosave directory for a trajectory reader.

mdevaluate.autosave.get_filename(function, checksum, description, *args)

Get the autosave filename for a specific function call.

mdevaluate.autosave.load_data(filename)

Load data from a npz file.

mdevaluate.autosave.notify(msg)

mdevaluate.autosave.save_data(filename, checksum, data)

Save data and checksum to a file.

mdevaluate.autosave.verify_file(filename, checksum)

Verify if the file matches the function call.

mdevaluate.atoms

exception mdevaluate.atoms.AtomMismatch

Bases: Exception

class mdevaluate.atoms.AtomSubset(atoms, selection=None, description='')

Bases: object

atom_names

indices

residue_ids

residue_names

subset(atom_name=None, residue_name=None, residue_id=None, indices=None)

Return a subset of the system. The selection is specified by one or more of the keyworss below. Names are matched as a regular expression with re.match.

Parameters:

• atom_name – Specification of the atom name
• residue_name – Specification of the resiude name
• residue_id – Residue ID or list of IDs
• indices – List of atom indices

summary

class mdevaluate.atoms.Atoms(atoms, indices=None, masses=None, charges=None)

Bases: object

Basic container class for atom information.

Parameters:

• atoms – N tuples of residue id, residue name and atom name.
• indices (optional) – Dictionary of named atom index groups.

residue_ids

Indices of the atoms residues

residue_names

Names of the atoms residues

atom_names

Names of the atoms

indices

Dictionary of named atom index groups, if specified

subset(*args, **kwargs)

Return a subset of these atoms with all atoms selected.

All arguments are passed to the AtomSubset.subset() method directly.

mdevaluate.atoms.center_of_mass(position, mass=None)

mdevaluate.atoms.compare_regex(list, exp)

Compare a list of strings with a regular expression.

mdevaluate.atoms.distance_to_atoms(ref, atoms, box=None)

Get the minimal distance from atoms to ref. The result is an array of with length == len(atoms)

mdevaluate.atoms.distance_to_atoms_cKDtree(ref, atoms, box=None, thickness=None)

Get the minimal distance from atoms to ref. The result is an array of with length == len(atoms) Can be faster than distance_to_atoms. Thickness defaults to box/5. If this is too small results may be wrong. If box is not given then periodic boundary conditions are not applied!

mdevaluate.atoms.from_grofile(grofile, index_file=None)

Load atoms from a GROMACS coordinate file.

Parameters:

• grofile (str) – Filename of the grofile
• index_file (str, optional) – Filename of the gromacs indexfile

Returns:

All atoms in grofile

Return type:

AtomSubset

mdevaluate.atoms.from_tprfile(tprfile, index_file=None)

Load atoms from a compiled tpr file.

Parameters:

• tprile (str) – Filename of the tpr file
• index_file (opt.) – Index file that should be loaded alongside.

Returns:

All atoms in tprfile

Return type:

AtomSubset

mdevaluate.atoms.gyration_radius(position)

Calculates a list of all radii of gyration of all molecules given in the coordinate frame, weighted with the masses of the individual atoms.

Parameters:position – Coordinate frame object

..math::

R_G = left(frac{sum_{i=1}^{n} m_i |vec{r_i} - vec{r_{COM}}|^2 }{sum_{i=1}^{n} m_i } rigth)^{frac{1}{2}}

mdevaluate.atoms.layer_of_atoms(atoms, thickness, plane_offset=array([0, 0, 0]), plane_normal=array([1, 0, 0]))

mdevaluate.atoms.next_neighbors(atoms, query_atoms=None, number_of_neighbors=1, distance_upper_bound=inf, distinct=False)

Find the N next neighbors of a set of atoms.

Parameters:

• atoms – The reference atoms and also the atoms which are queried if query_atoms is net provided
• query_atoms (opt.) – If this is not None, these atoms will be queried
• number_of_neighbors (int, opt.) – Number of neighboring atoms to find
• distance_upper_bound (float, opt.) – Upper bound of the distance between neighbors
• distinct (bool, opt.) – If this is true, the atoms and query atoms are taken as distinct sets of atoms

mdevaluate.coordinates

class mdevaluate.coordinates.CoordinateFrame

Bases: numpy.ndarray

atom_names

box

charges

indices

masses

nojump

pbc

residue_ids

residue_names

selection

time

volume

whole

class mdevaluate.coordinates.Coordinates(frames, atom_filter=None, atom_subset: mdevaluate.atoms.AtomSubset = None, mode=None)

Bases: object

Coordinates represent trajectory data, which is used for evaluation functions.

Atoms may be selected by specifing a atom_subset or a atom_filter.

clear_cache()

Clears the frame cache, if it is enabled.

description

get_frame(fnr)

Returns the fnr-th frame.

get_mode(mode)

mode

nojump

pbc

subset(atom_name=None, residue_name=None, residue_id=None, indices=None)

Return a subset of the system. The selection is specified by one or more of the keyworss below. Names are matched as a regular expression with re.match.

Parameters:

• atom_name – Specification of the atom name
• residue_name – Specification of the resiude name
• residue_id – Residue ID or list of IDs
• indices – List of atom indices

whole

class mdevaluate.coordinates.CoordinatesFilter(coordinates, atom_filter)

Bases: object

atom_subset

class mdevaluate.coordinates.CoordinatesKDTree(frames, selector=None, boxsize=None, maxcache=128, ckdtree=True)

Bases: object

A KDTree of coordinates frames. The KDtrees are cached by a functools.lru_cache(). Uses scipy.spatial.cKDTree by default, since it’s significantly faster. Make sure to use scipy 0.17 or later or switch to the normal KDTree, since cKDTree has a memory leak in earlier versions.

cache_info

Return info about the state of the cache.

clear_cache()

Clear the LRU cache.

class mdevaluate.coordinates.CoordinatesMap(coordinates, function)

Bases: object

description

nojump

pbc

subset(atom_name=None, residue_name=None, residue_id=None, indices=None)

Return a subset of the system. The selection is specified by one or more of the keyworss below. Names are matched as a regular expression with re.match.

Parameters:

• atom_name – Specification of the atom name
• residue_name – Specification of the resiude name
• residue_id – Residue ID or list of IDs
• indices – List of atom indices

whole

class mdevaluate.coordinates.MeanCoordinates(frames, atom_filter=None, mean=1)

Bases: mdevaluate.coordinates.Coordinates

len()

exception mdevaluate.coordinates.UnknownCoordinatesMode

Bases: Exception

mdevaluate.coordinates.centers_of_mass(c, *, masses=None)

A- 1 B- 2 A- 1 C 3 A- B- A- C A- B- A- C

Example: rd = XTCReader(‘t.xtc’) coordinates = Coordinates(rd) com = centers_of_mass(coordinates, (1.0, 2.0, 1.0, 3.0))

mdevaluate.coordinates.map_coordinates(func)

mdevaluate.coordinates.polar_coordinates(x, y)

Convert cartesian to polar coordinates.

mdevaluate.coordinates.pore_coordinates(coordinates, origin, sym_axis='z')

Map coordinates of a pore simulation so the pore has cylindrical symmetry.

Parameters:

• coordinates – Coordinates of the simulation
• origin – Origin of the pore which will be the coordinates origin after mapping
• sym_axis (opt.) – Symmtery axis of the pore, may be a literal direction ‘x’, ‘y’ or ‘z’ or an array of shape (3,)

mdevaluate.coordinates.radial_selector(frame, coordinates, rmin, rmax)

Return a selection of all atoms with radius in the interval [rmin, rmax].

mdevaluate.coordinates.rotate_axis(coords, axis)

Rotate a set of coordinates to a given axis.

mdevaluate.coordinates.spatial_selector(frame, transform, rmin, rmax)

Select a subset of atoms which have a radius between rmin and rmax. Coordinates are filtered by the condition:

rmin <= transform(frame) <= rmax

Parameters:

• frame – The coordinates of the actual trajectory
• transform – A function that transforms the coordinates of the frames into the one-dimensional spatial coordinate (e.g. radius).
• rmin – Minimum value of the radius
• rmax – Maximum value of the radius

mdevaluate.coordinates.spherical_coordinates(x, y, z)

Convert cartesian to spherical coordinates.

mdevaluate.coordinates.spherical_radius(frame, origin=None)

Transform a frame of cartesian coordinates into the sperical radius. If origin=None the center of the box is taken as the coordinates origin.

mdevaluate.coordinates.vectors(coordinates, atoms_a, atoms_b, normed=False, box=None)

Compute the vectors between the atoms of two subsets.

Parameters:

• coordinates – The Coordinates object the atoms will be taken from
• atoms_a – Mask or indices of the first atom subset
• atoms_b – Mask or indices of the second atom subset
• normed (opt.) – If the vectors should be normed
• box (opt.) – If not None, the vectors are calcualte with PBC

The defintion of atoms_a/b can be any possible subript of a numpy array. They can, for example, be given as a masking array of bool values with the same length as the frames of the coordinates. Or they can be a list of indices selecting the atoms of these indices from each frame.

It is possible to compute the mean of several atoms before calculating the vectors, by using a two-dimensional list of indices. The following code computes the vectors between atoms 0, 3, 6 and the mean coordinate of atoms 1, 4, 7 and 2, 5, 8:

>>> inds_a = [0, 3, 6]
>>> inds_b = [[1, 4, 7], [2, 5, 8]]
>>> vectors(coords, inds_a, inds_b)
array([
    coords[0] - (coords[1] + coords[2])/2,
    coords[3] - (coords[4] + coords[5])/2,
    coords[6] - (coords[7] + coords[8])/2,
])

mdevaluate.correlation

mdevaluate.correlation.coherent_scattering_function(onset, frame, q)

Calculate the coherent scattering function.

mdevaluate.correlation.correlation(function, frames)

mdevaluate.correlation.isf(start, frame, q, box=None)

Incoherent intermediate scattering function. To specify q, use water_isf = functools.partial(isf, q=22.77) # q has the value 22.77 nm^-1

Parameters:q – length of scattering vector

mdevaluate.correlation.log_indices(first, last, num=100)

mdevaluate.correlation.msd(start, frame)

Mean square displacement

mdevaluate.correlation.multi_subensemble_correlation(selector_function)

selector_function has to expect a frame and to return either valid indices (as with subensemble_correlation) or a multidimensional array whose entries are valid indices

e.g. slice(10,100,2)

e.g. [1,2,3,4,5]

e.g. [[[0,1],[2],[3]],[[4],[5],[6]] -> shape: 2,3 with list of indices of varying length

e.g. [slice(1653),slice(1653,None,3)]

e.g. [np.ones(len_of_frames, bool)]

in general using slices is the most efficient. if the selections are small subsets of a frame or when many subsets are empty using indices will be more efficient than using masks.

mdevaluate.correlation.non_gaussian(onset, frame)

Calculate the Non-Gaussian parameter : ..math:

lpha_2 (t) =

rac{3}{5} rac{langle r_i^4(t) angle}{langle r_i^2(t) angle^2} - 1

mdevaluate.correlation.overlap(onset, frame, crds_tree, radius)

Compute the overlap with a reference configuration defined in a CoordinatesTree.

Parameters:

• onset – Initial frame, this is only used to get the frame index
• frame – The current configuration
• crds_tree – A CoordinatesTree of the reference configurations
• radius – The cutoff radius for the overlap

This function is intended to be used with shifted_correlation(). As usual the first two arguments are used internally and the remaining ones should be defined with functools.partial().

If the overlap of a subset of the system should be calculated, this has to be defined through a selection of the reference configurations in the CoordinatesTree.

Example

>>> shifted_correlation(
...     partial(overlap, crds_tree=CoordinatesTree(traj), radius=0.11),
...     traj
... )

mdevaluate.correlation.rotational_autocorrelation(onset, frame, order=2)

Compute the rotaional autocorrelation of the legendre polynamial for the given vectors.

Parameters:

• frame (onset,) – CoordinateFrames of vectors
• order (opt.) – Order of the legendre polynomial.

Returns:

Skalar value of the correltaion function.

mdevaluate.correlation.set_has_counter(func)

mdevaluate.correlation.shifted_correlation(function, frames, index_distribution=<function log_indices>, correlation=<function correlation>, segments=10, window=0.5, skip=None, average=False)

Calculate the time series for a correlation function.

The times at which the correlation is calculated are determined automatically by the function given as index_distribution. The default is a logarithmic distribution.

Parameters:

• function – The function that should be correlated
• frames – The coordinates of the simulation data
• index_distribution (opt.) – A function that returns the indices for which the timeseries will be calculated
• correlation (function, opt.) – The correlation function
• segments (int, opt.) – The number of segments the time window will be shifted
• window (float, opt.) – The fraction of the simulation the time series will cover
• skip (float, opt.) – The fraction of the trajectory that will be skipped at the beginning, if this is None the start index of the frames slice will be used, which defaults to 0.
• counter (bool, opt.) – If True, returns length of frames (in general number of particles specified)
• average (bool, opt.) – If True, returns averaged correlation function

Returns:

A list of length N that contains the indices of the frames at which the time series was calculated and a numpy array of shape (segments, N) that holds the (non-avaraged) correlation data

if has_counter == True: adds number of counts to output tupel.

if average is returned it will be weighted.

Return type:

tuple

Example

Calculating the mean square displacement of a coordinates object named coords:

>>> indices, data = shifted_correlation(msd, coords)

mdevaluate.correlation.subensemble_correlation(selector_function, correlation_function=<function correlation>)

mdevaluate.correlation.susceptibility(time, correlation, **kwargs)

Calculate the susceptibility of a correlation function.

Parameters:

• time – Timesteps of the correlation data
• correlation – Value of the correlation function
• **kwargs (opt.) – Additional keyword arguments will be passed to filon_fourier_transformation().

mdevaluate.correlation.van_hove_distinct(onset, frame, bins, box=None, use_dask=True, comp=False, bincount=True)

Compute the distinct part of the Van Hove autocorrelation function.

..math::

G(r, t) = sum_{i, j} delta(|\vec r_i(0) - \vec r_j(t)| - r)

mdevaluate.correlation.van_hove_self(start, end, bins)

Compute the self part of the Van Hove autocorrelation function.

..math::

G(r, t) = sum_i delta(|\vec r_i(0) - \vec r_i(t)| - r)

mdevaluate.distribution

mdevaluate.distribution.distance_distribution(atoms, bins)

mdevaluate.distribution.fast_averaged_rdf(from_coords, bins, to_coords=None, times=10, exclude=0, **kwargs)

mdevaluate.distribution.mixing_ratio_distribution(atoms_a, atoms_b, bins_ratio, bins_density, weights_a=None, weights_b=None, weights_ratio=None)

Compute the distribution of the mixing ratio of two sets of atoms.

mdevaluate.distribution.next_neighbor_distribution(atoms, reference=None, number_of_neighbors=4, bins=None, normed=True)

Compute the distribution of next neighbors with the same residue name.

mdevaluate.distribution.pbc_spm_rdf(atoms_a, atoms_b=None, bins=None, box=None, exclude=0, returnx=False, **kwargs)

mdevaluate.distribution.pbc_tree_rdf(atoms_a, atoms_b=None, bins=None, box=None, exclude=0, returnx=False, **kwargs)

mdevaluate.distribution.radial_density(atoms, bins, symmetry_axis=(0, 0, 1), origin=(0, 0, 0), height=1, returnx=False)

Calculate the radial density distribution.

This function is meant to be used with time_average.

Parameters:

• atoms – Set of coordinates.
• bins – Bin specification that is passed to numpy.histogram. This needs to be a list of bin edges if the function is used within time_average.
• symmetry_axis (opt.) – Vector of the symmetry axis, around which the radial density is calculated, default is z-axis.
• origin (opt.) – Origin of the rotational symmetry, e.g. center of the pore.
• height (opt.) – Height of the pore, necessary for correct normalization of the density.
• returnx (opt.) – If True, the x ordinate of the distribution is returned.

mdevaluate.distribution.rdf(atoms_a, atoms_b=None, bins=None, box=None, kind=None, chunksize=50000, returnx=False, **kwargs)

Compute the radial pair distribution of one or two sets of atoms.

\[g_{AB}(r) = \frac{1}{\langle \rho_B\rangle N_A}\sum\limits_{i\in A}^{N_A} \sum\limits_{j\in B}^{N_B}\frac{\delta(r_{ij} -r)}{4\pi r^2}\]

For use with time_average(), define bins through the use of partial(), the atom sets are passed to time_average(), if a second set of atoms should be used specify it as coordinates_b and it will be passed to this function.

Parameters:

• atoms_a – First set of atoms, used internally
• atoms_b (opt.) – Second set of atoms, used internally
• bins – Bins of the radial distribution function
• box (opt.) – Simulations box, if not specified this is taken from atoms_a.box
• kind (opt.) – Can be ‘inter’, ‘intra’ or None (default).
• chunksize (opt.) – For large systems (N > 1000) the distaces have to be computed in chunks so the arrays fit into memory, this parameter controlls the size of these chunks. It should be as large as possible, depending on the available memory.
• returnx (opt.) – If True the x ordinate of the histogram is returned.

mdevaluate.distribution.shell_density(atoms, shell_radius, bins, shell_thickness=0.5, symmetry_axis=(0, 0, 1), origin=(0, 0, 0))

Compute the density distribution on a cylindrical shell.

Parameters:

• atoms – The coordinates of the atoms
• shell_radius – Inner radius of the shell
• bins – Histogramm bins, this has to be a two-dimensional list of bins: [angle, z]
• shell_thickness (opt.) – Thicknes of the shell, default is 0.5
• symmetry_axis (opt.) – The symmtery axis of the pore, the coordinates will be rotated such that this axis is the z-axis
• origin (opt.) – Origin of the pore, the coordinates will be moved such that this is the new origin.

Returns:

Two-dimensional density distribution of the atoms in the defined shell.

mdevaluate.distribution.spatial_density(atoms, bins, weights=None)

Compute the spatial density distribution.

mdevaluate.distribution.tetrahedral_order(atoms, reference_atoms=None)

mdevaluate.distribution.tetrahedral_order_distribution(atoms, reference_atoms=None, bins=None)

mdevaluate.distribution.time_average(function, coordinates, coordinates_b=None, pool=None)

Compute the time average of a function.

Parameters:

• function – The function that will be averaged, it has to accept exactly one argument which is the current atom set
• coordinates – The coordinates object of the simulation
• pool (multiprocessing.Pool, opt.) – A multiprocessing pool which will be used for cocurrent calculation of the averaged function

mdevaluate.distribution.time_histogram(function, coordinates, bins, hist_range, pool=None)

mdevaluate.evaluation

mdevaluate.functions

mdevaluate.functions.cole_cole(w, A, b, t0)

mdevaluate.functions.cole_davidson(w, A, b, t0)

mdevaluate.functions.colen(d, X, t8, A)

mdevaluate.functions.colenQ(d, X, Qb, g)

mdevaluate.functions.havriliak_negami(ω, A, β, α, τ)

Imaginary part of the Havriliak-Negami function.

\[\chi_{HN}(\omega) = \Im\left(\frac{A}{(1 + (i\omega\tau)^\alpha)^\beta}\right)\]

mdevaluate.functions.kww(t, A, τ, β)

mdevaluate.functions.kww_1e(A, τ, β)

mdevaluate.pbc

mdevaluate.pbc.nojump(frame, usecache=True)

Return the nojump coordinates of a frame, based on a jump matrix.

mdevaluate.pbc.pbc_backfold_compact(act_frame, box_matrix)

mimics “trjconv … -pbc atom -ur compact”

folds coords of act_frame in wigner-seitz-cell (e.g. dodecahedron)

mdevaluate.pbc.pbc_backfold_rect(act_frame, box_matrix)

mimics “trjconv … -pbc atom -ur rect”

folds coords of act_frame in cuboid

mdevaluate.pbc.pbc_diff(v1, v2=None, box=None)

mdevaluate.pbc.pbc_diff_numba

mdevaluate.pbc.pbc_diff_old(v1, v2, box)

Calculate the difference of two vestors, considering optional boundary conditions.

mdevaluate.pbc.pbc_diff_rect(v1, v2, box)

Calculate the difference of two vectors, considering periodic boundary conditions.

mdevaluate.pbc.pbc_diff_tric(v1, v2=None, box=None)

difference vector for arbitrary pbc

Args: box_matrix: CoordinateFrame.box

mdevaluate.pbc.pbc_dist(a1, a2, box=None)

mdevaluate.pbc.pbc_extend(c, box)

in: c is frame, box is frame.box out: all atoms in frame and their perio. image (shape => array(len(c)*27,3))

mdevaluate.pbc.pbc_kdtree(v1, box, leafsize=32, compact_nodes=False, balanced_tree=False)

kd_tree with periodic images box - whole matrix rest: optional optimization

mdevaluate.pbc.pbc_kdtree_query(v1, v2, box, n=1)

kd_tree query with periodic images

mdevaluate.pbc.pbc_points(coordinates, box, thickness=0, index=False, inclusive=True, center=None)

Returns the points their first periodic images. Does not fold them back into the box. Thickness 0 means all 27 boxes. Positive means the box+thickness. Negative values mean that less than the box is returned. index=True also returns the indices with indices of images being their originals values. inclusive=False returns only images, does not work with thickness <= 0

mdevaluate.pbc.whole(frame)

Apply -pbc whole to a CoordinateFrame.

mdevaluate.reader

Module that provides different readers for trajectory files.

It also provides a common interface layer between the file IO packages, namely pygmx and mdanalysis, and mdevaluate.

class mdevaluate.reader.BaseReader(rd)

Bases: object

Base class for trajectory readers.

filename

nojump_matrixes

class mdevaluate.reader.CachedReader(rd, maxsize)

Bases: mdevaluate.reader.BaseReader

A reader that has a least-recently-used cache for frames.

cache_info

Get Information about the lru cache.

clear_cache()

Clear the cache of the frames.

class mdevaluate.reader.DelayedReader(filename, reindex=False, ignore_index_timestamps=False)

Bases: mdevaluate.reader.BaseReader

filename

class mdevaluate.reader.EnergyReader(edrfile)

Bases: object

A reader for Gromacs energy files.

exception mdevaluate.reader.NoReaderAvailabelError

Bases: Exception

exception mdevaluate.reader.NojumpError

Bases: Exception

mdevaluate.reader.generate_nojump_matrixes(trajectory)

Create the matrixes with pbc jumps for a trajectory.

mdevaluate.reader.is_writeable(fname)

Test if a directory is actually writeable, by writing a temporary file.

mdevaluate.reader.load_nojump_matrixes(reader)

mdevaluate.reader.nojump_filename(reader)

mdevaluate.reader.open(topology, trajectory, cached=False, index_file=None, reindex=False, ignore_index_timestamps=False)

Open a trajectory file with the apropiate reader.

Parameters:

• filename (str) – Trajectory file to open, the reader will be chosen according to the file extension.
• cached (opt.) – If Reader should be cached with lru_cache. If this is True, maxsize for the cache is 128, otherwise the argument is passed as maxsize. Use cached=None to get an unbound cache.
• reindex (opt.) – Regenerate the index of the xtc-file
• nojump (opt.) – If nojump matrixes should be generated.

mdevaluate.reader.open_with_mdanalysis(topology, trajectory, cached=False)

Open a the topology and trajectory with mdanalysis.

mdevaluate.reader.open_with_pygmx(topology, trajectory, cached=False, reindex=False, ignore_index_timestamps=False, index_file=None)

Open a topology and trajectory with pygmx.

mdevaluate.reader.parse_jumps(trajectory)

mdevaluate.reader.save_nojump_matrixes(reader, matrixes=None)

mdevaluate.utils

Collection of utility functions.

mdevaluate.utils.Fqt_from_Grt(data, q)

Calculate the ISF from the van Hove function for a given q value by fourier transform.

\[F_q(t) = \int\limits_0^\infty dr \; G(r, t) \frac{\sin(qr)}{qr}\]

Parameters:

• data – Input data can be a pandas dataframe with columns ‘r’, ‘time’ and ‘G’ or an array of shape (N, 3), of tuples (r, t, G).
• q – Value of q.

Returns:

If input data was a dataframe the result will be returned as one too, else two arrays will be returned, which will contain times and values of Fq(t) respectively.

mdevaluate.utils.Sq_from_gr(r, gr, q, ρ)

Compute the static structure factor as fourier transform of the pair correlation function. [Yarnell]

\[\begin{split}S(q) - 1 = \\frac{4\\pi \\rho}{q}\\int\\limits_0^\\infty (g(r) - 1)\\,r \\sin(qr) dr\end{split}\]

Parameters:

• r – Radii of the pair correlation function
• gr – Values of the pair correlation function
• q – List of q values
• ρ – Average number density

[Yarnell]

Yarnell, J. L., Katz, M. J., Wenzel, R. G., & Koenig, S. H. (1973). Physical Review A, 7(6), 2130–2144. http://doi.org/10.1017/CBO9781107415324.004

mdevaluate.utils.coherent_histogram(func, coord_a, coord_b, bins, distinct=False)

Compute a coherent histogram over two arrays, equivalent to coherent_sum. For numpy arrays ofthis is equal to:

N, d = x.shape
M, d = y.shape
bins = np.arange(1, 5, 0.1)
coherent_histogram(f, x, y, bins) == histogram(f(x.reshape(N, 1, d), x.reshape(1, M, d)), bins=bins)

Parameters:

• func – The function is called for each two items in both arrays, this should return a scalar value.
• coord_b (coord_a,) – The two arrays.
• bins – The bins used for the histogram must be distributed regular on a linear scale.

mdevaluate.utils.coherent_sum(func, coord_a, coord_b)

Perform a coherent sum over two arrays \(A, B\).

\[\frac{1}{N_A N_B}\sum_i\sum_j f(A_i, B_j)\]

For numpy arrays this is equal to:

N, d = x.shape
M, d = y.shape
coherent_sum(f, x, y) == f(x.reshape(N, 1, d), x.reshape(1, M, d)).sum()

Parameters:

• func – The function is called for each two items in both arrays, this should return a scalar value.
• coord_b (coord_a,) – The two arrays.

mdevaluate.utils.filon_fourier_transformation(time, correlation, frequencies=None, derivative='linear', imag=True)

Fourier-transformation for slow varrying functions. The filon algorithmus is described in detail in ref [Blochowicz], ch. 3.2.3.

Parameters:

• time – List of times where the correlation function was sampled.
• correlation – Values of the correlation function.
• frequencies (opt.) – List of frequencies where the fourier transformation will be calculated. If None the frequencies will be choosen based on the input times.
• derivative (opt.) – Approximation algorithmus for the derivative of the correlation function. Possible values are: ‘linear’, ‘stencil’ or a list of derivatives.
• imag (opt.) – If imaginary part of the integral should be calculated.

If frequencies are not explicitly given they will be evenly placed on a log scale in the interval [1/tmax, 0.1/tmin] where tmin and tmax are the smallest respectively the biggest time (greater than 0) of the provided times. The frequencies are cut off at high values by one decade, since the fourier transformation deviates quite strongly in this regime.

[Blochowicz]

T. Blochowicz, Broadband dielectric spectroscopy in neat and binary molecular glass formers, Ph.D. thesis, Universität Bayreuth (2003)

mdevaluate.utils.five_point_stencil(xdata, ydata)

Calculate the derivative dy/dx with a five point stencil. This algorith is only valid for equally distributed x values.

Parameters:

• xdata – x values of the data points
• ydata – y values of the data points

Returns:

Values where the derivative was estimated and the value of the derivative at these points.

This algorithm is only valid for values on a regular grid, for unevenly distributed data it is only an approximation, albeit a quite good one.

See: https://en.wikipedia.org/wiki/Five-point_stencil

mdevaluate.utils.histogram(data, bins)

Compute the histogram of the given data. Uses numpy.bincount function, if possible.

mdevaluate.utils.mask2indices(mask)

Return the selected indices of an array mask. If the mask is two-dimensional, the indices will be calculated for the second axis.

Example

>>> mask2indices([True, False, True, False])
array([0, 2])
>>> mask2indices([[True, True, False], [True, False, True]])
array([[0, 1], [0, 2]])

mdevaluate.utils.moving_average(A, n=3)

Compute the running mean of an array. Uses the second axis if it is of higher dimensionality.

Parameters:

• data – Input data of shape (N, )
• n – Number of points over which the data will be averaged

Returns:

Array of shape (N-(n-1), )

Supports 2D-Arrays. Slower than runningmean for small n but faster for large n.

mdevaluate.utils.norm

mdevaluate.utils.quick1etau(t, C, n=7)

Estimate the time for a correlation function that goes from 1 to 0 to decay to 1/e.

If successful, returns tau as fine interpolation with a kww fit. The data is reduce to points around 1/e to remove short and long times from the kww fit! t is the time C is C(t) the correlation function n is the minimum number of points around 1/e required

mdevaluate.utils.runningmean(data, nav)

Compute the running mean of a 1-dimenional array.

Parameters:

• data – Input data of shape (N, )
• nav – Number of points over which the data will be averaged

Returns:

Array of shape (N-(nav-1), )

mdevaluate.utils.singledispatchmethod(func)

A decorator to define a genric instance method, analogue to functools.singledispatch.

mdevaluate.utils.superpose(x1, y1, x2, y2, N=100, damping=1.0)