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: Generator) → Generator[List, None, None][source]

Turn a generator into a new generator that yields in batches according to the number of available cores, OMP_NUM_THREADS.

Example: >>> core_batch((i for i in 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].

Parameters:: generator (Generator) – The generator to batch.
Returns:: A generator which iterates in batches of size OMP_NUM_THREADS.
Return type:: Generator

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

Creates a ThreadPoolExecutor with the relevant number of workers (according to the number of cores specified).

Returns:: A thread pool executor with max_workers=`OMP_NUM_THREADS`
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

An abstract class for total, coherent, and incoherent intermediate scattering functions. Equations used for calculating this are based on Kneller et al. Comput. Phys. Commun. 91 (1995) 191-214.

Note that properties for MD frames & Q are found in the SQwMixins class.

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.
Return type:: The FQt object with resolution applied.

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

Calculates 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]

Calculates the intermediate scattering function from a trajectory.

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

Parameters:

MD_input (CompactTrajectory) – a single CompactTrajectory object.
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

this is FQt, the intermediate scattering function (in arb)

Returns:: The dependent variables
Return type:: dict
Type:: Get the dependent variables

property dependent_variables_structure: dict[str, list]

The order in which the ‘FQt’ dependent variable is indexed in terms of ‘Q’ and ‘t’. 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))

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]

property errors: dict

Get or set the errors on the dependent variables, the intermediate scattering function (in arb)

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

property independent_variables: dict

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

Returns:: The independent variables
Return type:: dict
Type:: Get or set the independent variables

property t: ndarray

Get or set the times of the intermediate scattering function in units of fs

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

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

Captures 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

A class for containing, calculating and reading the intermediate scattering function for the total dynamic structure factor

name = 'IntermediateScatteringFunction'

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

Takes the incoherent scattering cross section of element and calculates the incoherent scattering length to be returned.

Parameters:: element (string) – 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]

Calculates 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

Returns:

A generator for the reciprocal space number density for each Q-vector.

Return type:

Generator[np.ndarray]

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

Gets the Hermann-Mauguin point group for the universe. Currently, MDMC can only create universes with mutually orthogonal sides, so this method can only produce cubic, tetragonal, or orthorhombic universe point groups.

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

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

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

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

Parameters:

point (tuple) – 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]

MDMC.trajectory_analysis.observables.fqt_coh module

Module for coherent FQt class

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

Bases: AbstractFQt

A class for containing, calculating and reading the intermediate scattering function for the 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

A class for containing, calculating and reading the intermediate scattering function for the incoherent dynamic structure factor

name = 'IncoherentIntermediateScatteringFunction'

MDMC.trajectory_analysis.observables.obs module

Module defining a class for storing, calculating and reading in observables from molecular dynamics 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]

Calculates the observable using input from an MD simulation

Parameters:

MD_input (Object) – 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

The 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. Example: 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 that this would also correspond to numpy.shape of the dependent variable being: np.shape(dependent_variable1)=(np.size(independent_variable1), np.size(independent_variable2), …)

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

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 maximum number of CompactTrajectory frames that can be used to calculate the dependent_variables

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

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

The minimum number 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: str

Get or set the origin of the observable

Returns:: The origin of the Observable, either 'experiment' or 'MD'
Return type:: str

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

Reads 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]]

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

Provides 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]

Creates 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]

Gets 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) → 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

Example

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

A class for containing, calculating and reading a pair distribution function (PDF).

We derive our definitions for this from the following publication: “A comparison of various commonly used correlation functions for describing total scattering” Keen, D. A. (2001). J. Appl. Cryst. 34, 172-177. DOI: https://doi.org/10.1107/S0021889800019993

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.

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 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_{i,j}^{N_{elements}} c_ic_jb_ib_j\).

This corresponds to the equation (8) in the following paper: “A comparison of various commonly used correlation functions for

describing total scattering”

Keen, D. A. (2001). J. Appl. Cryst. 34, 172-177. DOI: https://doi.org/10.1107/S0021889800019993

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 –

n_framesint
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 to the nearest positive integer).

use_averagebool
Optional parameter. 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. If set to False (default), only the last frame of the trajectory is used and n_frames will be ignored.

subsetlist 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_cohdict
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_minfloat
The minimum r (atomic separation) in Angstrom for which the PDF will be calculated. If this, 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. This cannot be passed if r is passed.

r_maxfloat
The maximum r (atomic separation) in Angstrom for which the PDF will be calculated. If this, r_min, 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. This cannot be passed if r is passed.

r_stepfloat
The step size of r (atomic separation) for which the PDF will be calculated. If this, r_min, and r_max are passed then these will create a range for the independent variable r, which will overwrite any r which has previously been defined. This cannot be passed if r is passed.

rnumpy.ndarray
The uniform r values in Angstrom for which the PDF will be calculated. This cannot be passed if r_min, r_max, and r_step are passed.

dimensionsarray-like
A 3 element array-like (list, tuple) with the dimensions of the Universe in Angstrom.

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

these are PDF, the pair distribution function (in barn)

Returns:: The dependent variables
Return type:: dict
Type:: Get the dependent variables

property dependent_variables_structure: dict[str, list]

The shape of the ‘PDF’ dependent variable in terms of ‘r’’: np.shape(self.PDF)=(np.size(self.r))

Returns:: The shape of the PDF dependent variable
Return type:: dict

property errors: dict

Get or set the errors on the dependent variables, the pair distribution function (in arb)

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

property independent_variables: dict

this is the atomic separation distance r (in Ang)

Returns:: The independent variables
Return type:: dict
Type:: Get or set the independent variable

maximum_frames() → None[source]

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 CompactTrajectory frames needed to calculate the dependent_variables 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)

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

Defines the current limitations on the atomic separation distance ‘r’ of the PairDistributionFunction ``Observable. 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. Comput. Phys. Commun. 91 (1995) 191-214.

Note that properties for MD frames & Q are found in the SQwMixins class.

property E

property SQw

property SQw_err

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

Calculates an array of nE uniformly spaced energy values from the time separation of the CompactTrajectory frames, dt. The frequencie are determined by the Fast Fourier Transform, as implemented by numpy, for 2 * nE points in time which we then crop 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

calculate_dt() → float[source]

Calculates the time separation of frames required by the experimental dataset, assuming uniform spacing. Note that 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

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 –

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

``energy_resolution` (dict)
Optionally 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. e.g. 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 (array)
1D array of Q float (in Ang^-1). (optional)

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]

Generates a resolution function in momentum and time that can be used in the calculation of SQw. Note that 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.

Note that if this resolution function is used on data outside its original range, then it will use nearest neighbour extrapolation. Additionally, 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.

Parameters:: dt (float) – The time spacing to use when performing the inverse Fourier transform in units of 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

property dependent_variables: dict

this is SQw, the dynamic structure factor (in arb)

Returns:: The dependent variables
Return type:: dict
Type:: Get the dependent variables

property dependent_variables_structure: dict[str, list]

The order in which the ‘SQw’ dependent variable is indexed in terms of ‘Q’ and ‘E’. 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))

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]

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

these are the frequency Q (in Ang^-1) and energy E (in``meV``)

Returns:: The independent variables
Return type:: dict
Type:: Get or set the independent variables

property recreated_Q: Get the indices of the recreated Q_values.

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

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]

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:: contains a boolean, and two floats (if traj_step and time_step have values) or two NoneType (if traj_step and time_step were passed in as None).
Return type:: tuple
Raises:: AssertionError –

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]

The maximum number of CompactTrajectory frames that can be used to calculate the dependent_variables 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.

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

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

The minimum number of CompactTrajectory frames needed to calculate the dependent_variables 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\]

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

Module contents

Modules calculating observables from molecular dynamics trajectories

Observables (/ indicates alias)

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

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

A class for containing, calculating and reading 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

A class for containing, calculating and reading the intermediate scattering function for the 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

A class for containing, calculating and reading a pair distribution function (PDF).

We derive our definitions for this from the following publication: “A comparison of various commonly used correlation functions for describing total scattering” Keen, D. A. (2001). J. Appl. Cryst. 34, 172-177. DOI: https://doi.org/10.1107/S0021889800019993

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.

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 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_{i,j}^{N_{elements}} c_ic_jb_ib_j\).

This corresponds to the equation (8) in the following paper: “A comparison of various commonly used correlation functions for

describing total scattering”

Keen, D. A. (2001). J. Appl. Cryst. 34, 172-177. DOI: https://doi.org/10.1107/S0021889800019993

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 –

n_framesint
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 to the nearest positive integer).

use_averagebool
Optional parameter. 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. If set to False (default), only the last frame of the trajectory is used and n_frames will be ignored.

subsetlist 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_cohdict
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_minfloat
The minimum r (atomic separation) in Angstrom for which the PDF will be calculated. If this, 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. This cannot be passed if r is passed.

r_maxfloat
The maximum r (atomic separation) in Angstrom for which the PDF will be calculated. If this, r_min, 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. This cannot be passed if r is passed.

r_stepfloat
The step size of r (atomic separation) for which the PDF will be calculated. If this, r_min, and r_max are passed then these will create a range for the independent variable r, which will overwrite any r which has previously been defined. This cannot be passed if r is passed.

rnumpy.ndarray
The uniform r values in Angstrom for which the PDF will be calculated. This cannot be passed if r_min, r_max, and r_step are passed.

dimensionsarray-like
A 3 element array-like (list, tuple) with the dimensions of the Universe in Angstrom.

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

these are PDF, the pair distribution function (in barn)

Returns:: The dependent variables
Return type:: dict
Type:: Get the dependent variables

property dependent_variables_structure: dict[str, list]

The shape of the ‘PDF’ dependent variable in terms of ‘r’’: np.shape(self.PDF)=(np.size(self.r))

Returns:: The shape of the PDF dependent variable
Return type:: dict

property errors: dict

Get or set the errors on the dependent variables, the pair distribution function (in arb)

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

property independent_variables: dict

this is the atomic separation distance r (in Ang)

Returns:: The independent variables
Return type:: dict
Type:: Get or set the independent variable

maximum_frames() → None[source]

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 CompactTrajectory frames needed to calculate the dependent_variables 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)

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

Defines the current limitations on the atomic separation distance ‘r’ of the PairDistributionFunction ``Observable. 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