StructureFactor

class mdcraft.analysis.structure.StructureFactor(groups: AtomGroup | tuple[AtomGroup], groupings: str | tuple[str] = 'atoms', *, mode: str = None, form: str = 'exp', dimensions: ndarray[float] | Quantity | Quantity = None, n_points: int = 32, n_surfaces: int = None, n_surface_points: int = 8, q_max: float | Quantity | Quantity = None, wavevectors: ndarray[float] = None, sort: bool = True, unique: bool = True, parallel: bool = False, verbose: bool = True, **kwargs)

Bases: NumbaAnalysisBase

Serial and parallel implementations to calculate the static structure factor \(S(q)\) or partial structure factor \(S_{\alpha\beta}(q)\) for species \(\alpha\) and \(\beta\).

The static structure factor is a measure of how a material scatters incident radiation, and can be computed directly from a molecular dynamics simulation trajectory using

\[\begin{split}S(q)&=\frac{1}{N}\left\langle\sum_{j=1}^N\sum_{k=1}^N \exp{[-i\mathbf{q}\cdot(\mathbf{r}_j-\mathbf{r}_k)]} \right\rangle\\&=\frac{1}{N}\left\langle\left( \sum_{j=1}^N\sin{(\mathbf{q}\cdot\mathbf{r}_j)}\right)^2 +\left(\sum_{j=1}^N\cos{(\mathbf{q}\cdot\mathbf{r}_j)} \right)^2\right\rangle\end{split}\]

where \(N\) is the total number of particles or centers of mass, \(\mathbf{q}\) is the scattering wavevector, and \(\mathbf{r}_i\) is the position of the \(i\)-th particle.

For multicomponent systems, the equation above can be generalized to get the partial structure factors

\[\begin{split}S_{\alpha\beta}(q)&=\frac{2-\delta_{\alpha\beta}}{N} \left\langle\sum_{j=1}^{N_\alpha}\sum_{k=1}^{N_\beta} \exp{[-i\mathbf{q}\cdot(\mathbf{r}_j-\mathbf{r}_k)]} \right\rangle\\ &=\frac{2-\delta_{\alpha\beta}}{N}\left\langle \sum_{j=1}^{N_\alpha}\cos{(\mathbf{q}\cdot\mathbf{r}_j)} \sum_{k=1}^{N_\beta}\cos{(\mathbf{q}\cdot\mathbf{r}_k)} +\sum_{j=1}^{N_\alpha}\sin{(\mathbf{q}\cdot\mathbf{r}_j)} \sum_{k=1}^{N_\beta}\sin{(\mathbf{q}\cdot\mathbf{r}_k)} \right\rangle\end{split}\]

where \(\delta_{ij}\) is the Kronecker delta, \(N_\alpha\) and \(N_\beta\) are the numbers of particles or centers of mass for species \(\alpha\) and \(\beta\).

The partial structure factors \(S_{\alpha\beta}(q)\) and the static structure factor \(S(q)\) are related via

\[S(q)=\sum_{\alpha=1}^{N_\mathrm{g}} \sum_{\beta=\alpha}^{N_\mathrm{g}}S_{\alpha\beta}(q)\]

Parameters:

groupsMDAnalysis.AtomGroup or array-like

Groups of atoms that share the same grouping type. If mode=None, all atoms in the universe must be assigned to a group. If mode="pair", there must be exactly one or two groups.

groupingsstr or array-like, default: "atoms"

Determines whether the centers of mass are used in lieu of individual atom positions. If groupings is a str, the same value is used for all groups.

Note

If the desired grouping is not "atoms",

• the trajectory file should have segments (or chains) containing residues (or molecules) and residues containing atoms, and

• residues and segments should be locally unwrapped at the simulation box edges, if not already, using MDAnalysis.transformations.wrap.unwrap, MDAnalysis.core.groups.AtomGroup.unwrap(), or MDAnalysis.lib.mdamath.make_whole().

Valid values:

• "atoms": Atom positions (generally or for coarse-grained simulations).

• "residues": Residues’ centers of mass (for atomistic simulations).

• "segments": Segments’ centers of mass (for atomistic polymer simulations).

modestr, keyword-only, optional

Evaluation mode.

Valid values:

• None: The static structure factor is computed.

• "pair": The partial structure factor is computed between the groups in groups.

• "partial": The partial structure factors for all unique pairs in groups is computed.

formstr, keyword-only, default: "exp"

Expression used to evaluate the structure factors.

Valid values:

• "exp": Exponential form. Slightly faster due to fewer mathematical operations.

• "trig": Trigonometric form. Slightly slower but doesn’t have overflow issues.

dimensionsarray-like, openmm.unit.Quantity, or pint.Quantity, keyword-only, optional

System dimensions \((L_x,\,L_y,\,L_z)\). Only used to determine the scattering wavevectors when wavevectors is not specified.

Shape: \((3,)\).

Reference unit: \(\mathrm{Å}\).

n_pointsint, keyword-only, default: 32

Number of points \(n_\mathrm{points}\) in the scattering wavevector grid to generate

\[\mathbf{q}=2\pi\left(\frac{a}{L_x},\,\frac{b}{L_y},\, \frac{c}{L_z}\right)\]

where \(a\), \(b\), and \(c\) are integers from \(0\) up to \(n_\mathrm{points}-1\).

Additional wavevectors can be introduced via n_surfaces and n_surface_points for more accurate structure factors at small wavenumbers. Alternatively, the desired wavevectors can be specified directly in wavevectors.

n_surfacesint, keyword-only, optional

Number of spherical surfaces in the first octant that intersect with the grid wavevectors along the three coordinate axes for which to introduce extra wavevectors for more accurate structure factor values. Only available if the system is perfectly cubic.

n_surface_pointsint, keyword-only, default: 8

Number of extra wavevectors to introduce per spherical surface. Has no effect if n_surfaces is not specified.

q_maxfloat, openmm.unit.Quantity, or pint.Quantity, keyword-only, optional

Maximum wavenumber \(q_\mathrm{max}\).

Reference unit: \(\mathrm{Å}^{-1}\).

wavevectorsarray-like, openmm.unit.Quantity, or pint.Quantity, keyword-only, optional

Scattering wavevectors for which to compute structure factors. Has precedence over n_points, n_surfaces, and n_surface_points if specified.

Shape: \((N_q,\,3)\).

Reference unit: \(\mathrm{Å}^{-1}\).

sortbool, keyword-only, default: True

Determines whether the results are sorted by the wavenumbers.

uniquebool, keyword-only, default: True

Determines whether structure factors for the same wavenumber are grouped and averaged.

parallelbool, keyword-only, default: False

Determines whether the analysis is performed in parallel.

verbosebool, keyword-only, default: True

Determines whether detailed progress is shown.

**kwargs

Additional keyword arguments to pass to MDAnalysis.analysis.base.AnalysisBase.

Attributes:

universeMDAnalysis.Universe

MDAnalysis.core.universe.Universe object containing all information describing the system.

results.unitsdict

Reference units for the results. For example, to get the reference units for results.wavenumbers, call results.units["wavenumbers"].

results.pairstuple

All unique pairs of indices of the groups of atoms in groups. The ordering coincides with the column indices in results.ssf.

results.wavenumbersnumpy.ndarray

Wavenumbers \(q\).

Shape: \((N_q,)\).

Reference unit: \(\mathrm{Å}^{-1}\).

results.ssfnumpy.ndarray

Static structure factor \(S(q)\) or partial structure factors \(S_{\alpha\beta}(q)\).

Shape: \((1,\,N_q)\) or \((C(N_\mathrm{g}+1,\,2),\,N_q)\).

Methods

get_supported_backends	Tuple with backends supported by the core library for a given class.
run	Performs the calculation.
save	Saves results to a binary or archive file in NumPy format.

classmethod get_supported_backends()

Tuple with backends supported by the core library for a given class. User can pass either one of these values as backend=... to run() method, or a custom object that has apply method (see documentation for run()):

• ‘serial’: no parallelization

• ‘multiprocessing’: parallelization using multiprocessing.Pool

• ‘dask’: parallelization using dask.delayed.compute(). Requires installation of mdanalysis[dask]

If you want to add your own backend to an existing class, pass a backends.BackendBase subclass (see its documentation to learn how to implement it properly), and specify unsupported_backend=True.

Returns:

tuple

names of built-in backends that can be used in run(backend=...)()

Added in version 2.8.0: ..

property parallelizable

Boolean mark showing that a given class can be parallelizable with split-apply-combine procedure. Namely, if we can safely distribute _single_frame() to multiple workers and then combine them with a proper _conclude() call. If set to False, no backends except for serial are supported.

Note

If you want to check parallelizability of the whole class, without explicitly creating an instance of the class, see _analysis_algorithm_is_parallelizable. Note that you setting it to other value will break things if the algorithm behind the analysis is not trivially parallelizable.

Returns:

bool

if a given AnalysisBase subclass instance is parallelizable with split-apply-combine, or not

Added in version 2.8.0: ..

run(start: int = None, stop: int = None, step: int = None, frames: slice | ndarray[int] = None, verbose: bool = None, *, n_threads: int = None, **kwargs) → NumbaAnalysisBase

Performs the calculation.

Parameters:

startint, optional

Starting frame for analysis.

stopint, optional

Ending frame for analysis.

stepint, optional

Number of frames to skip between each analyzed frame.

framesslice or array-like, optional

Index or logical array of the desired trajectory frames.

verbosebool, optional

Determines whether detailed progress is shown.

n_threadsint, keyword-only, optional

Number of threads to use for analysis.

**kwargs

Additional keyword arguments to pass to MDAnalysis.lib.log.ProgressBar.

Returns:

selfNumbaAnalysisBase

Analysis object with results.

save(file: str | TextIO, archive: bool = True, compress: bool = True, **kwargs) → None

Saves results to a binary or archive file in NumPy format.

Parameters:

filestr or file

Filename or file-like object where the data will be saved. If file is a str, the .npy or .npz extension will be appended automatically if not already present.

archivebool, default: True

Determines whether the results are saved to a single archive file. If True, the data is stored in a .npz file. Otherwise, the data is saved to multiple .npy files.

compressbool, default: True

Determines whether the .npz file is compressed. Has no effect when archive=False.

**kwargs

Additional keyword arguments to pass to numpy.save(), numpy.savez(), or numpy.savez_compressed(), depending on the values of archive and compress.