MDMC.trajectory_analysis.observables package

Submodules

MDMC.trajectory_analysis.observables.concurrency_tools module

Tools for concurrency in observable calculation.

MDMC.trajectory_analysis.observables.concurrency_tools.core_batch(generator: Iterable[T]) → Iterable[tuple[T, ...]][source]

Batch generator according to the number of available cores, OMP_NUM_THREADS.

Parameters:: generator (Iterable[T]) – The generator to batch.
Yields:: tuple[T] – Batches of size OMP_NUM_THREADS.

See also

itertools.batched: Standard implementation from 3.12.

Examples

>>> core_batch(range(10))
on 1 core produces [0], [1], [2], [3], [4], [5], [6], [7], [8], [9]
on 4 cores produces [0, 1, 2, 3], [4, 5, 6, 7], [8, 9].

MDMC.trajectory_analysis.observables.concurrency_tools.create_executor() → ThreadPoolExecutor[source]

Create a ThreadPoolExecutor with the relevant number of workers.

Returns:: A thread pool executor with max_workers=`OMP_NUM_THREADS` or 1 if not set.
Return type:: ThreadPoolExecutor

MDMC.trajectory_analysis.observables.fqt module

Module for Intermediate Scattering Function class.

class MDMC.trajectory_analysis.observables.fqt.AbstractFQt[source]

Bases: SQwMixins, Observable

Abstract class for total, coherent, and incoherent intermediate scattering functions.

Equations used for calculating this are based on Kneller et al. [Kneller].

See also

SQwMixins: Properties for MD frames & Q are found in the SQwMixins class.

References

[Kneller]

Kneller et al. Comput. Phys. Commun. 91 (1995) 191-214.

property FQt

apply_resolution(resolution: Resolution) → FQt[source]

Apply instrument resolution to an FQt object.

Parameters:: resolution (Resolution) – The Resolution object to apply to FQt.
Returns:: The FQt object with resolution applied.
Return type:: FQt

calculate_SQw(energy: list[float], resolution: Resolution = None) → ndarray[source]

Calculate S(Q, w) from F(Q, t), accounting for instrument resolution.

In order to obtain len(energy) values in energy, we reflect the intermediate scattering function in time to give it dimensions of (len(self.Q), 2 * (len(self.t)) - 2). This uses the fact it is even in time, and the number of time points is chosen to be 1 greater than the number of energy points [Rapaport, The Art of Molecular Dynamics Simulation (2nd Edition), 2004, page 142].

The numpy implementation of the FFT gives frequencies arranged so that the first len(energy) points in the energy dimension correspond to positive frequencies, and the remaining points have negative frequency.

Parameters:

energy (list of floats) – The list of energy (E) points at which S(Q, w) will be calculated.
resolution (Resolution (default None)) – The instrument resolution object which will be applied to FQt.

Returns:

The S(Q, w) calculated from F(Q, t).

Return type:

numpy.ndarray

calculate_from_MD(MD_input: CompactTrajectory, verbose: int = 0, **settings: dict) → None[source]

Calculate the intermediate scattering function from a trajectory.

independent_variables can either be set previously or defined within **settings.

Parameters:

MD_input (CompactTrajectory) – CompactTrajectory object to compute F(Q, t) from.
verbose (int, optional) –
The level of verbosity:
- Verbose level 0 gives no information.
- Verbose level 1 gives final time for the whole method.
- Verbose level 2 gives final time and also a progress bar.
- Verbose level 3 gives final time, a progress bar, and time per step.
**settings –

n_Q_vectors (int)
The maximum number of Q_vectors for any Q value. The greater the number of Q_vectors, the more accurate the calculation, but the longer it will take.

dimensions (list, tuple, numpy.ndarray)
A 3 element tuple or array of float specifying the dimensions of the Universe in units of Ang

calculate_rho_config(config: ndarray, single_Q_vectors: list) → ndarray[source]

Calculate density over an entire configuration.

Parameters:

config (np.ndarray) – An array of atom positions.
single_Q_vectors (list) – A list of Q-vectors for a single value of Q.

Returns:

An array of rho values for each timestep summed over the atoms in the system, corresponding to each Q.

Return type:

np.ndarray

property dependent_variables: dict

Get the dependent variables.

This is FQt, the intermediate scattering function (in arb).

Returns:: The dependent variables.
Return type:: dict

property dependent_variables_structure: dict[str, list]

The order of indexing of ‘FQt’ dependent variables in terms of ‘Q’ and ‘t’.

The purpose of this method is to ensure consistency between different readers/methods which create FQt objects.

Returns:: The shape of the SQw dependent variable.
Return type:: dict[str, list]

Notes

Explicitly: we have that self.FQt[Q_index, t_index] is the data point for given indices of self.Q and self.t.

It also means that: np.shape(self.FQt)=(np.size(self.Q), np.size(self.t)).

property errors: dict

Get or set the errors on the dependent variables.

This is the intermediate scattering function error (in arb).

Returns:: The errors on the dependent_variables.
Return type:: dict

property independent_variables: dict

Get or set the independent variables.

These are the frequency Q (in Ang^-1) and time t (in fs).

Returns:: The independent variables.
Return type:: dict

property t: ndarray

Get or set the times of the intermediate scattering function (in fs).

Returns:: 1D array of times in fs.
Return type:: numpy.array

property uniformity_requirements: dict[str, dict[str, bool]]

Get restrictions required for computing E & Q.

Capture the current limitations on the time ‘t’ and reciprocal lattice points ‘Q’ within the intermediate scattering function Observables.

If using FFT, then ‘t’ must be uniform and start at zero, otherwise it has no restrictions. ‘Q’ must be uniform but does not need to start at zero.

Returns:: Dictionary of uniformity restrictions for ‘t’ and ‘Q’.
Return type:: dict[str, dict[str, bool]]

class MDMC.trajectory_analysis.observables.fqt.FQt[source]

Bases: AbstractFQt

Class to process the intermediate scattering function for the total dynamic structure factor.

name = 'IntermediateScatteringFunction'

MDMC.trajectory_analysis.observables.fqt.calc_incoherent_scatt_length(element)[source]

Calculate incoherent scattering length from incoherent scattering cross section.

Parameters:: element (str) – A string representing the chemical symbol of the element, e.g ‘He’ for Helium.
Returns:: Incoherent scattering length of chemical symbol passed in.
Return type:: float

MDMC.trajectory_analysis.observables.fqt.calculate_rho(positions: ndarray, Q_vectors: list) → Generator[ndarray, None, None][source]

Calculate t dependent number density in reciprocal space for all Q vectors.

As rho is the sum of the contributions for all of the specified Q vectors, these Q vectors should have the same Q value.

Parameters:

positions (numpy.ndarray) – An array of atomic positions for which the reciprocal space number density should be calculated.
Q_vectors (list) – A list of one or more Q vectors with the same Q value.

Yields:

np.ndarray – Reciprocal space number density for each Q-vector.

MDMC.trajectory_analysis.observables.fqt.get_point_group(dimensions: ndarray) → str[source]

Get the Hermann-Mauguin point group for the universe.

Note

Currently, MDMC can only create universes with mutually orthogonal sides, so this method can only produce cubic, tetragonal, or orthorhombic universe point groups.

Parameters:: dimensions (numpy.ndarray) – An array of length 3, containing the dimensions of the universe.
Returns:: The point group symbol for the universe.
Return type:: str

Notes

For tetragonal universes, an additional identifier (x), (y), or (z) is added to identify which of the sides is unique.

MDMC.trajectory_analysis.observables.fqt.wyckoff_symmetries(point: tuple[float, float, float], point_group: str) → set[tuple[float, float, float]][source]

Get the Wyckoff symmetries for a point based on its point group.

Parameters:

point (tuple[float, float, float]) – A tuple of length 3 which corresponds to the coordinates (x, y, z) of a point.
point_group (str) –
The point group of the universe in Hermann-Mauguin notation. Currently accepted groups are:
- ’m-3m’ (cubic)
- ’4/mmm’ (tetragonal)
- ’mmm’ (orthorhombic)

Returns:

A calculated set of the symmetries for the point.

Return type:

set[tuple[float, float, float]]

MDMC.trajectory_analysis.observables.fqt_coh module

Module for coherent FQt class.

class MDMC.trajectory_analysis.observables.fqt_coh.FQtCoherent[source]

Bases: AbstractFQt

Class for processing the intermediate scattering function for coherent dynamic structure factor.

name = 'FQt_coh'

MDMC.trajectory_analysis.observables.fqt_incoh module

Module for incoherent FQt class.

class MDMC.trajectory_analysis.observables.fqt_incoh.FQtIncoherent[source]

Bases: AbstractFQt

Class for processing intermediate scattering function for incoherent dynamic structure factor.

name = 'IncoherentIntermediateScatteringFunction'

MDMC.trajectory_analysis.observables.obs module

Module defining a class for processing observables from MD trajectories.

class MDMC.trajectory_analysis.observables.obs.Observable[source]

Bases: ABC

Abstract class that defines methods common to all observable data containers.

Observable data can either be from a file or calculated from MD and stored in the data property, along with the associated uncertainty. The bool property from_MD states the source of the information.

reader

The file reader for reading experimental data.

Type:: ObservableReader

abstract calculate_from_MD(MD_input: CompactTrajectory | list[CompactTrajectory], verbose: int = 0, **parameters: dict) → None[source]

Calculate the observable using input from an MD simulation.

Parameters:

MD_input (CompactTrajectory | list[CompactTrajectory]) – Some input from an MD simulation, commonly a CompactTrajectory.
verbose (int) – Enables verbose printing of the calculation.
**parameters – Additional parameters required for calculation specific Observable objects.

property data: dict

Get the independent, dependent and error data.

Returns:: The independent, dependent and error data.
Return type:: dict

abstract property dependent_variables: dict

The dependent variables.

Returns:: The dependent variables.
Return type:: dict

abstract property dependent_variables_structure: dict

Structure of the dependent variables with respect to the independent variables.

Specifically, the order in which the dependent variables are indexed with regards to the independent variables.

The purpose of this method is to ensure that all Observable s of a particular type are created with dependent_variables that are consistent regardless of how they were created (e.g. by different Reader s).

Returns:: The np.shape of the dependent variables.
Return type:: dict

Examples

If dep_var1[indep_var1_index, indep_var2_index, ...] == data point for values of the independent_variables with the stated indices, then the relevant entry in the returned dict should be: {'dependent_variable1': [independent_variable1, independent_variable2, ...]}

Note

This would also correspond to numpy.shape of the dependent variable being:

np.shape(dependent_variable1)=(np.size(independent_variable1),
                               np.size(independent_variable2), ...)

abstract property errors: dict

The errors on the dependent variables.

Returns:: The errors on the dependent_variables.
Return type:: dict

abstract property independent_variables: dict

The independent variables.

Returns:: The independent variables.
Return type:: dict

abstract maximum_frames() → int[source]

The max no. of CompactTrajectory frames able to calculate the dependent_variables.

Returns:: The maximum number of frames.
Return type:: int

abstract minimum_frames(dt: float = None) → int[source]

The no. of CompactTrajectory frames needed to calculate the dependent_variables.

Parameters:: dt (float, optional) – The time separation of frames in fs, default is None.
Returns:: The minimum number of frames.
Return type:: int

property name: str

Get or set the module name that used for factory instantiation.

Returns:: The name of the module in which the Observable is located.
Return type:: str

property origin: Literal['experiment', 'MD']

Get or set the origin of the observable.

Returns:: The origin of the Observable, either 'experiment' or 'MD'.
Return type:: {‘experiment’, ‘MD’}

read_from_file(reader: str, file_name: str) → None[source]

Read in experimental data from a file using a specified reader.

Parameters:

reader (str) – The name of the required file reader.
file_name (str) – The name of the file.

abstract property uniformity_requirements: dict[str, dict[str, bool]]

Get the current limitations on independent_variables of the Observable.

It captures if the independent_variables are required to be uniform or to start at zero The keys of the returned dictionary should be the variables that have such a restriction, with the associated values being a dictionary with booleans if the variables are ‘uniform’ or ‘zeroed’.

Variables without any requirements do not need to be included, but can be included.

Note

If there are no uniformity requirements it is okay to return ‘None’.

Returns:: Dictionary of independent variables with their uniformity restrictions represented as booleans.
Return type:: dict[str, dict[str, bool]]

property use_FFT: bool

Get or set whether to use FFT when calculating from MD.

Returns:: Whether to use FFT.
Return type:: bool

MDMC.trajectory_analysis.observables.obs_factory module

Factory class for generating observables.

class MDMC.trajectory_analysis.observables.obs_factory.ObservableFactory[source]

Bases: object

Provide a factory for creating an Observable.

Any module within the observables submodule can be created with a string of the class name, as long as it is a subclass of Observable.

classmethod create_observable(name: str) → Observable[source]

Create an Observable object from a module name.

The Observable object must be registered with the ObservableFactory.

Parameters:: name (str) – The name of the module with which the Observable is registered.
Returns:: An Observable object.
Return type:: Observable

classmethod get_observable(name: str) → type[Observable][source]

Get an Observable class from a registry name.

Parameters:: name (str) – The name of the module with which the Observable is registered.
Returns:: A subclass of Observable.
Return type:: cls

classmethod register(names: str | Iterable[str]) → Callable[source]

A class level decorator for registering Observable classes.

The names of the modules with which the Observable is registered should be the parameter passed to the decorator.

Parameters:: names (str) – The names of the modules with which the Observable is registered.
Returns:: Wrapped class having been added to registry.
Return type:: Callable

Examples

To register the SQw class with ObservableFactory:

@ObservableFactory.register('SQw')
class SQw(Observable):

registry: dict[str, Observable] = {'CoherentDynamicStructureFactor': <class 'MDMC.trajectory_analysis.observables.sqw.SQwCoherent'>, 'CoherentIntermediateScatteringFunction': <class 'MDMC.trajectory_analysis.observables.fqt_coh.FQtCoherent'>, 'DynamicStructureFactor': <class 'MDMC.trajectory_analysis.observables.sqw.SQw'>, 'FQt': <class 'MDMC.trajectory_analysis.observables.fqt.FQt'>, 'FQtCoh': <class 'MDMC.trajectory_analysis.observables.fqt_coh.FQtCoherent'>, 'FQtCoherent': <class 'MDMC.trajectory_analysis.observables.fqt_coh.FQtCoherent'>, 'FQtIncoherentFQtIncoh': <class 'MDMC.trajectory_analysis.observables.fqt_incoh.FQtIncoherent'>, 'FQt_coh': <class 'MDMC.trajectory_analysis.observables.fqt_coh.FQtCoherent'>, 'FQt_incoh': <class 'MDMC.trajectory_analysis.observables.fqt_incoh.FQtIncoherent'>, 'IncoherentDynamicStructureFactor': <class 'MDMC.trajectory_analysis.observables.sqw.SQwIncoherent'>, 'IncoherentIntermediateScatteringFunction': <class 'MDMC.trajectory_analysis.observables.fqt_incoh.FQtIncoherent'>, 'IntermediateScatteringFunction': <class 'MDMC.trajectory_analysis.observables.fqt.FQt'>, 'PDF': <class 'MDMC.trajectory_analysis.observables.pdf.PairDistributionFunction'>, 'PairDistributionFunction': <class 'MDMC.trajectory_analysis.observables.pdf.PairDistributionFunction'>, 'SQw': <class 'MDMC.trajectory_analysis.observables.sqw.SQw'>, 'SQwCoh': <class 'MDMC.trajectory_analysis.observables.sqw.SQwCoherent'>, 'SQwCoherent': <class 'MDMC.trajectory_analysis.observables.sqw.SQwCoherent'>, 'SQwIncoherentSQwIncoh': <class 'MDMC.trajectory_analysis.observables.sqw.SQwIncoherent'>, 'SQw_coh': <class 'MDMC.trajectory_analysis.observables.sqw.SQwCoherent'>, 'SQw_incoh': <class 'MDMC.trajectory_analysis.observables.sqw.SQwIncoherent'>}

MDMC.trajectory_analysis.observables.pdf module

Module for calculating the total pair distribution function (PDF).

class MDMC.trajectory_analysis.observables.pdf.PairDistributionFunction[source]

Bases: Observable

Class for processing a pair distribution function (PDF).

We derive our definitions for this from [Keen]:

We employ the following mathematical form for the total pair distribution function (PDF):

\[G(r) = \sum_{i,j}^{N_{elements}} c_ic_jb_ib_j(g_{ij}(r) - 1)\]

where \(c_i\) is the number concentration of element \(i\), \(b_i\) is the (coherent) scattering length of element \(i\). (This corresponds to equation 8 in the above publication)

The partial pair distribution, \(g_{ij}\), is:

\[g_{ij}(r) = \frac{h_{ij}(r)}{4 \pi r^2 \rho_{j} \Delta{r}}\]

where \(h_{ij}`\) is the histogram of distances of \(j\) element atoms around atoms of element \(i\), with bins of size \(\Delta{r}\), and \(\rho_{j}\) is the number density of atoms of element \(j\). As \(g_{ij}(0) = 0\), it is evident that \(G(0) = -\sum_{i,j}^{N_{elements}} c_ic_jb_ib_j\). (This corresponds to equation 10 in the above publication)

The total PDF is contained in PDF and the partial pair PDFs (if calculated or imported) are contained in partial_pdfs.

References

[Keen]

Keen, D. A., “A comparison of various commonly used correlation functions for describing total scattering” J. Appl. Cryst. 34 (2001) 172-177.

DOI: https://doi.org/10.1107/S0021889800019993

property PDF

property PDF_err

calculate_from_MD(MD_input: CompactTrajectory, verbose: int = 0, **settings: dict)[source]

Calculate the pair distribution function, \(G(r)\) from a CompactTrajectory.

The partial pair distribution (from __) for a pair i-j, \(g_{ij}\), is:

\[g_{ij}(r) = \frac{h_{ij}(r)}{4 \pi r^2 \rho_{j} \Delta{r}}\]

where \(h_{ij}`\) is the histogram of distances of \(j\) element atoms around atoms of element \(i\), with bins of size \(\Delta{r}\), and \(\rho_{j}\) is the number density of atoms of element \(j\). As \(g_{ij}(0) = 0\), it is evident that \(G(0) = -\sum\limits_{i,j}^{N_{elements}} c_ic_jb_ib_j\).

The total pair distribution function (pdf.PDF) has the form:

\[G(r) = \sum_{i,j}^{N_{elements}} c_ic_jb_ib_j(g_{ij}(r) - 1)\]

where \(c_i\) is the proportion of element \(i\) in the material, \(b_i\) is the (coherent) scattering length of element \(i\)

This corresponds to the equation (10) in the above paper.

Independent variables can either be set previously or defined within settings.

A number of frames can be specified, from which the PDF and its error are calculated. If the number of frames is too large relative to the run length, the samples will be correlated, which will result in an underestimate of the error.

Parameters:

MD_input (CompactTrajectory) – A single CompactTrajectory object.
verbose (int) – Verbose print settings. Not currently implemented for PDF.
**settings – Extra options.
n_frames (int) –
The number of frames from which the PDF and its error are calculated. These frames are selected uniformly, and the step is taken to be n_frames / total number of frames rounded to the nearest positive integer.

If this is not passed, 1% of the total number of frames are used (rounded up).
use_average (bool) –
If set to False (default), only the last frame of the trajectory is used and n_frames will be ignored.

If set to True then the mean value for PDF is calculated across selected frames from the trajectory. Also, the errors are set to the standard deviation calculated over the multiple frames.
subset (list of tuples) –
The subset of element pairs from which the PDF is calculated.

This can be used to calculate the partial PDFs of a multicomponent system. If this is not passed, all combinations of elements are used i.e. the PDF is the total PDF.
b_coh (dict) –
A dictionary containing coherent scattering length values for one or more elements.

This can be used to calculate the PDF of a system where one or more elements has a coherent scattering length different from the coherent scattering lengths from periodictable.
r_min (float) – The minimum r (atomic separation) in Angstrom for which the PDF will be calculated.
r_max (float) – The maximum r (atomic separation) in Angstrom for which the PDF will be calculated.
r_step (float) – The step size of r (atomic separation) for which the PDF will be calculated.
r (numpy.ndarray) – The uniform r values in Angstrom for which the PDF will be calculated.
dimensions (array-like) – A 3 element array-like (list, tuple) with the dimensions of the Universe in Angstrom.

Notes

If this, r_min, r_max, and r_step are passed then these will create a range for the independent variable r, which will overwrite any r which has previously been defined. These cannot be passed if r is passed.

Examples

To calculate the O-O partial PDF from a simulation of water, use the subset keyword:

pdf.calculate_from_MD(trajectory, subset=[(O, O)])

To calculate the sum of the H-O and O-O partial PDFs:

pdf.calculate_from_MD(trajectory, subset=[(O, O), (H, O)])

To calculate the total PDF for sodium chloride with 37Cl:

pdf.calculate_from_MD(trajectory, b_coh={'Cl':3.08})

To calculate the total PDF for r values of [1., 2., 3., 4.]:

pdf.calculate_from_MD(trajectory, r=[1., 2., 3., 4.])

To calculate the total PDF over an average of 5 frames:

pdf.calculate_from_MD(trajectory, use_average=True, n_frames=5)

property dependent_variables: dict

Get the dependent variables.

This is the pair distribution function (in barn).

Returns:: The dependent variables.
Return type:: dict

property dependent_variables_structure: dict[str, list]

The shape of the ‘PDF’ dependent variable in terms of ‘r’.

Where np.shape(self.PDF) == (np.size(self.r)).

Returns:: The shape of the PDF dependent variable.
Return type:: dict[str, list[str]]

property errors: dict

Get or set the errors on the dependent variables.

This is the pair distribution function (in barn)

Returns:: The errors on the dependent_variables.
Return type:: dict

property independent_variables: dict

Get or set the independent variables.

This is the atomic separation distance r (in Ang).

Returns:: The independent variables.
Return type:: dict

maximum_frames() → None[source]

The maximum number of frames needed to calculate the dependent_variables.

Returns:: There is no hard limit on the number of frames that can be used, so return None.
Return type:: None

minimum_frames(dt: float = None) → int[source]

The minimum number of frames needed to calculate the dependent_variables.

For PDF, this is 1.

Parameters:: dt (float, optional) – The time separation of frames in fs, default is None, not used.
Returns:: The minimum number of frames.
Return type:: int

name = 'PairDistributionFunction'

property r: float | None

Get or set the value of the atomic separation distance (in Ang).

Returns:: The atomic separation distance.
Return type:: float

property uniformity_requirements: dict[str, dict[str, bool]]

Define the current limitations on the atomic separation distance ‘r’.

The requirement is that ‘r’ must be uniform, but it does not have to start at zero.

Returns:: Dictionary of uniformity restrictions for ‘r’.
Return type:: dict[str, dict[str, bool]]

MDMC.trajectory_analysis.observables.sqw module

Module for AbstractSQw and total SQw class.

class MDMC.trajectory_analysis.observables.sqw.AbstractSQw[source]

Bases: SQwMixins, Observable

An abstract class for total, coherent and incoherent dynamic structure factors.

The equations used for calculating this are based on Kneller et al. [Kneller].

Note

Properties for MD frames & Q are found in the SQwMixins class.

References

[Kneller]

Kneller et al. Comput. Phys. Commun. 91 (1995) 191-214.

property E

property SQw

property SQw_err

static calculate_E(nE: int, dt: float) → ndarray[source]

Compute energy from MD.

Calculate an array of nE uniformly spaced energy values from the time separation of the CompactTrajectory frames, dt.

The frequencies are determined by the fast Fourier transform, as implemented by numpy, for 2 * nE points in time which are truncated to only include nE positive frequencies.

As we are dealing with frequency rather than angular frequency here, the relation to between energy is given by:

\[E = h \nu\]

Parameters:

nE (int) – The number of energy values to be calculated.
dt (float) – The step size between frames in fs.

Returns:

An array of float specifying the energy in units of meV.

Return type:

numpy.ndarray

See also

numpy.fft.fftfreq: Frequency calculation.

calculate_dt() → float[source]

Calculate timestep required by experimental dataset assuming uniform spacing.

Note

This may be different from the time separation that the user has given as an input, as it only depends on the current values for self.E.

The relationship between time and energy comes from the numpy implementation of the FFT for 2 * nE points where:

\[\begin{split}\nu_{max} &=& \frac{n_E - 1}{2 n_E \Delta t} \\\\ \therefore \Delta t &=& \frac{h (n_E - 1)}{2 n_E E_{max}}\end{split}\]

Returns:: The time separation required by the current values of self.E.
Return type:: float

See also

numpy.fft.fft: Numpy FFT implementation.

calculate_from_MD(MD_input: CompactTrajectory, verbose: int = 0, **settings: dict)[source]

Calculate the dynamic structure factor, S(Q, w) from a CompactTrajectory.

If the CompactTrajectory has more frames than the self.maximum_frames() that can be used to recreate the grid of energy points, it can slice the CompactTrajectory into sub-trajectories of length self.maximum_frames(), with the slicing specified through the settings use_average and cont_slicing.

The independent_variable Q can either be set previously or defined within **settings.

Parameters:

MD_input (CompactTrajectory) – An MDMC CompactTrajectory from which to calculate SQw.
verbose (int, optional) –
The level of verbosity:
- Verbose level 0 gives no information.
- Verbose level 1 gives final time for the whole method.
- Verbose level 2 gives final time and also a progress bar.
- Verbose level 3 gives final time, a progress bar, and time per step.
**settings (dict) – Extra options.
n_Q_vectors (int) – The maximum number of Q_vectors for any Q value. The greater the number of Q_vectors, the more accurate the calculation, but the longer it will take.
dimensions (list, tuple or numpy.ndarray) – A 3 element tuple or array of float specifying the dimensions of the Universe in units of Ang.
energy_resolution (dict) –
Specify energy resolution and function in units of ueV (micro eV), in the format of the one-line dict {‘function’: value}, where function is the resolution function and value is the desired FWHM.

For example, to pass a Gaussian resolution of 80ueV we use {‘gaussian’: 80}.

Currently accepted functions are ‘gaussian’ and ‘lorentzian’ Can also be ‘lazily’ given as float, in which case it is assumed to be Gaussian.
Q_values (numpy.ndarray) – 1D array of Q float (in Ang^-1).
use_average (bool) –
Optional parameter if a list of more than one Trajectory is used.

If set to True (default) then the mean value for S(Q, w) is calculated. Also, the errors are set to the standard deviation calculated over the list of CompactTrajectory objects.
cont_slicing (bool) –
Flag to decide between two possible behaviours when the number of MD_steps is larger than the minimum required to calculate the observables.

If False (default) then the CompactTrajectory is sliced into non-overlapping sub-CompactTrajectory blocks for each of which the observable is calculated.

If True, then the CompactTrajectory is sliced into as many non-identical sub-CompactTrajectory blocks as possible (with overlap allowed).

calculate_resolution_functions(dt: float) → dict[source]

Generate resolution function in momentum and time.

Parameters:

dt (float) –

The time spacing for the inverse Fourier transform (in fs).

Ideally this should be the same as the frame separation expected when applying this function.

Returns:

A dictionary with the key ‘SQw’ corresponding to function which accepts arrays of time and momentum (respectively) and returns a 2D array of values for the instrument resolution.

Return type:

dict

Notes

There are several caveats to this function a user should be aware of:

This uses the SQw values of the Observable it is called from, and so should only be called for an observable which has been created from relevant resolution data, i.e. a vanadium sample.
If this resolution function is used on data outside its original range, then it will use nearest neighbour extrapolation.
The input will be reflected in the time/energy domain as symmetry about 0 is assumed.

If for whatever reason this is not appropriate for the data in question, this function should not be used.

property dependent_variables: dict

Get the dependent variables.

This is:

SQw, the dynamic structure factor (in arb)

Returns:: The dependent variables.
Return type:: dict

property dependent_variables_structure: dict[str, list]

The order in which the ‘SQw’ dependent variable is indexed in terms of ‘Q’ and ‘E’.

The purpose of this method is to ensure consistency between different readers/methods which create SQw objects.

Returns:: The shape of the SQw dependent variable.
Return type:: dict[str, list]

Notes

Explicitly: we have that self.SQw[Q_index, E_index] is the data point for given indices of self.Q and self.E.

It also means that: np.shape(self.SQw)=(np.size(self.Q), np.size(self.E)).

property errors: dict

Get or set the errors on the dependent variables.

The dynamic structure factor (in arb)

Returns:: The errors on the dependent_variables.
Return type:: dict

property independent_variables: dict

Get or set the independent variables.

These are:

the frequency \(Q\) (in Ang^-1)
energy \(E\) (in``meV``)

Returns:: The independent variables.
Return type:: dict

property recreated_Q: ndarray

Get the indices of the recreated Q_values.

Returns:: Array of recreated Q values.
Return type:: numpy.ndarray

property uniformity_requirements: dict[str, dict[str, bool]]

Get restrictions required for computing E & Q.

Captures the current limitations on the energy ‘E’ and reciprocal lattice points ‘Q’ within the dynamic structure factor Observables. If using FFT, then ‘E’ must be uniform and start at zero, otherwise it has no restrictions. ‘Q’ must be uniform but does not need to start at zero.

Returns:: Dictionary of uniformity restrictions for ‘E’ and ‘Q’.
Return type:: dict[str, dict[str, bool]]

validate_energy(time_step: float = None)[source]

Ensure dt is valid for computing energies.

Asserts that the user set frame separation dt leads to energy separation that matches that of the experiment. If not, it changes the time step and trajectory step to fix this. The time step value is prioritised here.

Parameters:

time_step (float, optional) – User specified length of time for each update of the atoms trajectories in the simulation, default is None.

Returns:

valid (bool) – Whether dt is valid.
traj_step (Optional[int]) – Ideal traj_step if provided.
time_step (Optional[float]) – Ideal time_step if provided.
dt (float) – Required dt.

property w

class MDMC.trajectory_analysis.observables.sqw.SQw[source]

Bases: AbstractSQw

A class for the total dynamic structure factor.

Calculation is done in the respective FQt object, and this is just a reference to get the correct FQt object.

name = 'SQw'

class MDMC.trajectory_analysis.observables.sqw.SQwCoherent[source]

Bases: AbstractSQw

A class for the coherent dynamic structure factor.

name = 'SQw_coh'

class MDMC.trajectory_analysis.observables.sqw.SQwIncoherent[source]

Bases: AbstractSQw

A class for the incoherent dynamic structure factor.

name = 'SQw_incoh'

class MDMC.trajectory_analysis.observables.sqw.SQwMixins[source]

Bases: object

A mixin class for properties used by both SQw and FQt objects.

property Q

maximum_frames() → int | None[source]

Compute maximum number of CompactTrajectory frames to calculate dependent_variables.

Returns:: The maximum number of frames.
Return type:: int

Notes

Depends on self.use_FFT.

If True, it is the number of energy steps + 1, in order to allow for a reflection in time which only counts the end points once.

Otherwise, there is no limit and all frames will contribute to the calculation.

minimum_frames(dt: float = None) → int[source]

Compute minimum number of CompactTrajectory frames to calculate dependent_variables.

Parameters:: dt (float, optional) – The time separation of frames in fs, default is None.
Returns:: The minimum number of frames.
Return type:: int

Notes

Depends on self.use_FFT.

If self.use_FFT == True, it is the number of energy steps + 1, in order to allow for a reflection in time which only counts the end points once.

If self.use_FFT == False, there is not a hard minimum on number of frames. However, to distinguish our smallest differences in energy \(F(Q,t)\) needs to cover at least a time period \(T_{min}\) such that:

\[T_{min} \sim \frac{h}{\Delta E_{min}}\]

Due to the aforementioned reflection in the time domain, to cover a period of \(T_{min}\) we only need \(N\) frames:

\[N = \frac{T_{min}}{2 \Delta t} + 1 = \frac{h}{2 \Delta t \Delta E_{min}} + 1\]

Module contents

Modules calculating observables from molecular dynamics trajectories.

Observables

DynamicStructureFactor / SQw CoherentDynamicStructureFactor / SQwCoherent / SQwCoh / SQw_coh IncoherentDynamicStructureFactor / SQwIncoherentSQwIncoh / SQw_incoh IntermediateScatteringFunction / FQt IncoherentIntermediateScatteringFunction / FQtIncoherentFQtIncoh / FQt_incoh CoherentIntermediateScatteringFunction / FQtCoherent / FQtCoh / FQt_coh PDF / PairDistributionFunction

Examples

Observables can be instantiated using the names above.

For instance an SQw observable can be instantiated using either of the aliases:

from MDMC.trajectory_analysis import observables
sqw = observables.SQw()                     # This line...
sqw = observables.DynamicStructureFactor()  # ...is equivalent to this line

MDMC.trajectory_analysis.observables.CoherentDynamicStructureFactor: alias of SQwCoherent

MDMC.trajectory_analysis.observables.CoherentIntermediateScatteringFunction: alias of FQtCoherent

MDMC.trajectory_analysis.observables.DynamicStructureFactor: alias of SQw

class MDMC.trajectory_analysis.observables.FQt[source]

Bases: AbstractFQt

Class to process the intermediate scattering function for the total dynamic structure factor.

name = 'IntermediateScatteringFunction'

MDMC.trajectory_analysis.observables.FQtCoh: alias of FQtCoherent

class MDMC.trajectory_analysis.observables.FQtCoherent[source]

Bases: AbstractFQt

Class for processing the intermediate scattering function for coherent dynamic structure factor.

name = 'FQt_coh'

MDMC.trajectory_analysis.observables.FQtIncoherentFQtIncoh: alias of FQtIncoherent

MDMC.trajectory_analysis.observables.FQt_coh: alias of FQtCoherent

MDMC.trajectory_analysis.observables.FQt_incoh: alias of FQtIncoherent

MDMC.trajectory_analysis.observables.IncoherentDynamicStructureFactor: alias of SQwIncoherent

MDMC.trajectory_analysis.observables.IncoherentIntermediateScatteringFunction: alias of FQtIncoherent

MDMC.trajectory_analysis.observables.IntermediateScatteringFunction: alias of FQt

MDMC.trajectory_analysis.observables.PDF: alias of PairDistributionFunction

class MDMC.trajectory_analysis.observables.PairDistributionFunction[source]

Bases: Observable

Class for processing a pair distribution function (PDF).

We derive our definitions for this from [Keen]:

We employ the following mathematical form for the total pair distribution function (PDF):

\[G(r) = \sum_{i,j}^{N_{elements}} c_ic_jb_ib_j(g_{ij}(r) - 1)\]

where \(c_i\) is the number concentration of element \(i\), \(b_i\) is the (coherent) scattering length of element \(i\). (This corresponds to equation 8 in the above publication)

The partial pair distribution, \(g_{ij}\), is:

\[g_{ij}(r) = \frac{h_{ij}(r)}{4 \pi r^2 \rho_{j} \Delta{r}}\]

where \(h_{ij}`\) is the histogram of distances of \(j\) element atoms around atoms of element \(i\), with bins of size \(\Delta{r}\), and \(\rho_{j}\) is the number density of atoms of element \(j\). As \(g_{ij}(0) = 0\), it is evident that \(G(0) = -\sum_{i,j}^{N_{elements}} c_ic_jb_ib_j\). (This corresponds to equation 10 in the above publication)

The total PDF is contained in PDF and the partial pair PDFs (if calculated or imported) are contained in partial_pdfs.

References

[Keen]

Keen, D. A., “A comparison of various commonly used correlation functions for describing total scattering” J. Appl. Cryst. 34 (2001) 172-177.

DOI: https://doi.org/10.1107/S0021889800019993

property PDF

property PDF_err

calculate_from_MD(MD_input: CompactTrajectory, verbose: int = 0, **settings: dict)[source]

Calculate the pair distribution function, \(G(r)\) from a CompactTrajectory.

The partial pair distribution (from __) for a pair i-j, \(g_{ij}\), is:

\[g_{ij}(r) = \frac{h_{ij}(r)}{4 \pi r^2 \rho_{j} \Delta{r}}\]

where \(h_{ij}`\) is the histogram of distances of \(j\) element atoms around atoms of element \(i\), with bins of size \(\Delta{r}\), and \(\rho_{j}\) is the number density of atoms of element \(j\). As \(g_{ij}(0) = 0\), it is evident that \(G(0) = -\sum\limits_{i,j}^{N_{elements}} c_ic_jb_ib_j\).

The total pair distribution function (pdf.PDF) has the form:

\[G(r) = \sum_{i,j}^{N_{elements}} c_ic_jb_ib_j(g_{ij}(r) - 1)\]

where \(c_i\) is the proportion of element \(i\) in the material, \(b_i\) is the (coherent) scattering length of element \(i\)

This corresponds to the equation (10) in the above paper.

Independent variables can either be set previously or defined within settings.

A number of frames can be specified, from which the PDF and its error are calculated. If the number of frames is too large relative to the run length, the samples will be correlated, which will result in an underestimate of the error.

Parameters:

MD_input (CompactTrajectory) – A single CompactTrajectory object.
verbose (int) – Verbose print settings. Not currently implemented for PDF.
**settings – Extra options.
n_frames (int) –
The number of frames from which the PDF and its error are calculated. These frames are selected uniformly, and the step is taken to be n_frames / total number of frames rounded to the nearest positive integer.

If this is not passed, 1% of the total number of frames are used (rounded up).
use_average (bool) –
If set to False (default), only the last frame of the trajectory is used and n_frames will be ignored.

If set to True then the mean value for PDF is calculated across selected frames from the trajectory. Also, the errors are set to the standard deviation calculated over the multiple frames.
subset (list of tuples) –
The subset of element pairs from which the PDF is calculated.

This can be used to calculate the partial PDFs of a multicomponent system. If this is not passed, all combinations of elements are used i.e. the PDF is the total PDF.
b_coh (dict) –
A dictionary containing coherent scattering length values for one or more elements.

This can be used to calculate the PDF of a system where one or more elements has a coherent scattering length different from the coherent scattering lengths from periodictable.
r_min (float) – The minimum r (atomic separation) in Angstrom for which the PDF will be calculated.
r_max (float) – The maximum r (atomic separation) in Angstrom for which the PDF will be calculated.
r_step (float) – The step size of r (atomic separation) for which the PDF will be calculated.
r (numpy.ndarray) – The uniform r values in Angstrom for which the PDF will be calculated.
dimensions (array-like) – A 3 element array-like (list, tuple) with the dimensions of the Universe in Angstrom.

Notes

If this, r_min, r_max, and r_step are passed then these will create a range for the independent variable r, which will overwrite any r which has previously been defined. These cannot be passed if r is passed.

Examples

To calculate the O-O partial PDF from a simulation of water, use the subset keyword:

pdf.calculate_from_MD(trajectory, subset=[(O, O)])

To calculate the sum of the H-O and O-O partial PDFs:

pdf.calculate_from_MD(trajectory, subset=[(O, O), (H, O)])

To calculate the total PDF for sodium chloride with 37Cl:

pdf.calculate_from_MD(trajectory, b_coh={'Cl':3.08})

To calculate the total PDF for r values of [1., 2., 3., 4.]:

pdf.calculate_from_MD(trajectory, r=[1., 2., 3., 4.])

To calculate the total PDF over an average of 5 frames:

pdf.calculate_from_MD(trajectory, use_average=True, n_frames=5)

property dependent_variables: dict

Get the dependent variables.

This is the pair distribution function (in barn).

Returns:: The dependent variables.
Return type:: dict

property dependent_variables_structure: dict[str, list]

The shape of the ‘PDF’ dependent variable in terms of ‘r’.

Where np.shape(self.PDF) == (np.size(self.r)).

Returns:: The shape of the PDF dependent variable.
Return type:: dict[str, list[str]]

property errors: dict

Get or set the errors on the dependent variables.

This is the pair distribution function (in barn)

Returns:: The errors on the dependent_variables.
Return type:: dict

property independent_variables: dict

Get or set the independent variables.

This is the atomic separation distance r (in Ang).

Returns:: The independent variables.
Return type:: dict

maximum_frames() → None[source]

The maximum number of frames needed to calculate the dependent_variables.

Returns:: There is no hard limit on the number of frames that can be used, so return None.
Return type:: None

minimum_frames(dt: float = None) → int[source]

The minimum number of frames needed to calculate the dependent_variables.

For PDF, this is 1.

Parameters:: dt (float, optional) – The time separation of frames in fs, default is None, not used.
Returns:: The minimum number of frames.
Return type:: int

name = 'PairDistributionFunction'

property r: float | None

Get or set the value of the atomic separation distance (in Ang).

Returns:: The atomic separation distance.
Return type:: float

property uniformity_requirements: dict[str, dict[str, bool]]

Define the current limitations on the atomic separation distance ‘r’.

The requirement is that ‘r’ must be uniform, but it does not have to start at zero.

Returns:: Dictionary of uniformity restrictions for ‘r’.
Return type:: dict[str, dict[str, bool]]

class MDMC.trajectory_analysis.observables.SQw[source]

Bases: AbstractSQw

A class for the total dynamic structure factor.

Calculation is done in the respective FQt object, and this is just a reference to get the correct FQt object.

name = 'SQw'

MDMC.trajectory_analysis.observables.SQwCoh: alias of SQwCoherent

class MDMC.trajectory_analysis.observables.SQwCoherent[source]

Bases: AbstractSQw

A class for the coherent dynamic structure factor.

name = 'SQw_coh'

MDMC.trajectory_analysis.observables.SQwIncoherentSQwIncoh: alias of SQwIncoherent

MDMC.trajectory_analysis.observables.SQw_coh: alias of SQwCoherent

MDMC.trajectory_analysis.observables.SQw_incoh: alias of SQwIncoherent