Index¶

mdapy.ackland_jones_analysis module¶

class mdapy.ackland_jones_analysis.AcklandJonesAnalysis(data: DataFrame, box: Box, verlet_list: ndarray, distance_list: ndarray)[source]¶

Bases: object

Perform Ackland-Jones Analysis (AJA) to identify crystal structures.

The Ackland-Jones Analysis is a structure identification method that classifies local atomic environments based on coordination patterns and bond angles. It can distinguish between FCC, BCC, HCP, and other common crystal structures.

Parameters:

data (pl.DataFrame) – Atomic data containing at least the columns ‘x’, ‘y’, ‘z’ for positions.
box (Box) – Simulation box object.
verlet_list (np.ndarray) – Neighbor list array of shape (N, max_neigh) containing neighbor indices.
distance_list (np.ndarray) – Distance list array of shape (N, max_neigh) containing neighbor distances.

data¶

Input atomic data.

Type:: pl.DataFrame

box¶

Simulation box.

Type:: Box

verlet_list¶

Neighbor indices.

Type:: np.ndarray

distance_list¶

Neighbor distances.

Type:: np.ndarray

aja¶

Structure type array of shape (N,) assigned after calling compute(). Possible values:

0: Unknown/other structure
1: FCC
2: HCP
3: BCC
4: ICO (icosahedral)

Type:: np.ndarray

Notes

This analysis requires a pre-computed neighbor list. Typically, at least 14 nearest neighbors are needed for accurate classification.

References

Ackland, G. J., & Jones, A. P. (2006). Applications of local crystal structure measures in experiment and simulation. Physical Review B, 73(5), 054104.

compute() → None[source]¶

Perform the Ackland-Jones structure analysis.

This method computes the structure type for each atom based on the neighbor list and assigns the result to the aja attribute.

Notes

After calling this method, the aja attribute will contain integer structure identifiers for each atom.

mdapy.angular_distribution_function module¶

class mdapy.angular_distribution_function.AngularDistributionFunction(data: DataFrame, box: Box, rc_dict: Dict[str, List[float]], nbin: int, verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray)[source]¶

Bases: object

Calculate Angular Distribution Function (ADF) for bond angle analysis.

The angular distribution function quantifies the distribution of bond angles between atomic triplets (i-j-k), where i is the central atom. This is useful for characterizing local atomic structure and coordination environments.

Parameters:

data (pl.DataFrame) – Atomic data containing positions and element information. Must have columns: ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Simulation box object with lattice parameters and boundaries.
rc_dict (dict of str to list of float) – Dictionary mapping triplet patterns to cutoff radii. Format: {‘A-B-C’: [rc_AB_min, rc_AB_max, rc_AC_min, rc_AC_max]} where A, B, C are element symbols and A is the central atom. Example: {‘O-H-H’: [0, 2.0, 0, 2.0]} for water molecules.
nbin (int) – Number of bins for angle discretization (0-180 degrees).
verlet_list (np.ndarray) – Verlet neighbor list, shape (N, max_neigh).
distance_list (np.ndarray) – Neighbor distances, shape (N, max_neigh).
neighbor_number (np.ndarray) – Number of neighbors for each atom, shape (N,).

bond_angle_distribution¶

Bond angle distribution for each triplet type, shape (Npair, nbin). Available after calling compute().

Type:: np.ndarray

r_angle¶

Angle bin centers in degrees, shape (nbin,). Available after calling compute().

Type:: np.ndarray

ele_unique¶

Unique element symbols in the system, sorted alphabetically.

Type:: list of str

pair_list¶

Type indices for each triplet pattern, shape (Npair, 3).

Type:: np.ndarray

rc_list¶

Cutoff radii for each triplet pattern, shape (Npair, 4).

Type:: np.ndarray

Notes

The angular distribution function is calculated by:

For each central atom i of type A
Find neighbors j of type B within [rc_AB_min, rc_AB_max]
Find neighbors k of type C within [rc_AC_min, rc_AC_max]
Calculate angle θ between j-i-k bonds
Histogram the angles into bins

Examples

Calculate ADF for water molecules (H-O-H angle):

>>> from mdapy import System
>>> system = System("water.xyz")
>>> system.build_neighbor(rc=3.0, max_neigh=50)
>>> from mdapy.angular_distribution_function import AngularDistributionFunction
>>> adf = AngularDistributionFunction(
...     system.data,
...     system.box,
...     rc_dict={"O-H-H": [0, 2.0, 0, 2.0]},  # O is central atom
...     nbin=40,
...     verlet_list=system.verlet_list,
...     distance_list=system.distance_list,
...     neighbor_number=system.neighbor_number,
... )
>>> adf.compute()
>>> fig, ax = adf.plot_bond_angle_distribution()

compute()[source]¶

Compute the angular distribution function.

This method calculates the bond angle distribution for all specified triplet patterns. Results are stored in bond_angle_distribution and r_angle.

plot_bond_angle_distribution(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶

Plot bond angle distribution.

Parameters:

fig (Figure, optional) – Existing matplotlib figure. If None, a new figure is created.
ax (Axes, optional) – Existing matplotlib axes. If None, new axes are created.

Returns:

fig (Figure) – Matplotlib figure object.
ax (Axes) – Matplotlib axes object.

Notes

The plot shows bond angle distribution as a function of angle (0-180°) for all triplet patterns specified in rc_dict. Each pattern is shown as a separate line with markers.

mdapy.atomic_strain module¶

class mdapy.atomic_strain.AtomicStrain(rc: float, ref: System, affine: bool = False, max_neigh: int | None = None)[source]¶

Bases: object

Calculate atomic-level shear and volumetric strain tensors.

This class computes local strain at each atom by comparing neighbor distances between a reference (undeformed) configuration and a current (deformed) configuration. The method calculates the deformation gradient tensor F and derives strain measures from it.

Parameters:

rc (float) – Cutoff radius for neighbor search in Angstroms.
ref (System) – Reference (undeformed) atomic system.
affine (bool, default=False) – Whether to apply affine transformation to map current box to reference box. If True, removes global box deformation to isolate local atomic strain.
max_neigh (int, optional) – Maximum number of neighbors to consider. If None, automatically determined.

ref¶

Reference system with neighbor list built.

Type:: System

rc¶

Cutoff radius used for neighbor calculations.

Type:: float

max_neigh¶

Maximum neighbors parameter.

Type:: int or None

affine¶

Affine transformation flag.

Type:: bool

repeat¶

Replication factors (nx, ny, nz) if reference box is too small.

Type:: tuple of int

Notes

The atomic strain is calculated using the following algorithm:

Build matrices V and W for each atom i:

\[ \begin{align}\begin{aligned}V_{mn} = \sum_j \Delta r^{ref}_{jn} \Delta r^{ref}_{jm}\\W_{mn} = \sum_j \Delta r^{ref}_{jn} \Delta r^{cur}_{jm}\end{aligned}\end{align} \]

where the sum is over all neighbors j within cutoff rc, and \(\Delta r^{ref}\) and \(\Delta r^{cur}\) are neighbor distance vectors in reference and current configurations.
Calculate deformation gradient tensor:

\[F = (W \cdot V^{-1})^T\]
Compute strain tensor (Green-Lagrange strain):

\[\varepsilon = \frac{1}{2}(F^T F - I)\]
Extract strain measures:
- Shear strain (von Mises equivalent strain):
  
  \[\eta_s = \sqrt{\varepsilon_{xy}^2 + \varepsilon_{xz}^2 + \varepsilon_{yz}^2 + \frac{1}{6}[(\varepsilon_{xx}-\varepsilon_{yy})^2 + (\varepsilon_{xx}-\varepsilon_{zz})^2 + (\varepsilon_{yy}-\varepsilon_{zz})^2]}\]
- Volumetric strain (hydrostatic strain):
  
  \[\eta_v = \frac{1}{3}(\varepsilon_{xx} + \varepsilon_{yy} + \varepsilon_{zz})\]

The affine transformation option is useful when the simulation box itself undergoes deformation (e.g., under stress). Setting affine=True removes this global deformation to isolate local atomic-level strain.

References

Shimizu, F., Ogata, S., & Li, J. (2007). Theory of shear banding in metallic glasses and molecular dynamics calculations. Materials Transactions, 48(11), 2923-2927.

Examples

Calculate atomic strain between two configurations:

>>> from mdapy import System
>>> # Load reference (undeformed) configuration
>>> ref_system = System("reference.dump")
>>> # Load current (deformed) configuration
>>> cur_system = System("deformed.dump")
>>> # Initialize strain calculator
>>> from mdapy.atomic_strain import AtomicStrain
>>> strain = AtomicStrain(rc=5.0, ref=ref_system)
>>> # Compute strain
>>> strain.compute(cur_system)
>>> # Results are stored in cur_system.data
>>> print(cur_system.data["shear_strain"])
>>> print(cur_system.data["volumetric_strain"])

With affine transformation (removes global box deformation):

>>> strain = AtomicStrain(rc=5.0, ref=ref_system, affine=True)
>>> strain.compute(cur_system)

compute(current: System)[source]¶

Compute atomic strain for the current configuration.

This method calculates shear and volumetric strain for each atom by comparing neighbor distances with the reference configuration. Results are added as new columns ‘shear_strain’ and ‘volumetric_strain’ to the current system’s data.

Parameters:: current (System) – Current (deformed) atomic system. Must have the same number of atoms as the reference system.
Raises:: AssertionError – If current system has different number of atoms than reference.

Notes

The computed strain values are added directly to current.data:

shear_strain: Von Mises equivalent shear strain at each atom
volumetric_strain: Hydrostatic (volumetric) strain at each atom

Positive volumetric strain indicates expansion, negative indicates compression. Shear strain is always non-negative and indicates local distortion.

If affine=True was set during initialization, the current configuration is first transformed to remove global box deformation before calculating local atomic strain.

Examples

>>> strain = AtomicStrain(rc=5.0, ref=ref_system)
>>> strain.compute(cur_system)
>>> # Access results
>>> shear = cur_system.data["shear_strain"].to_numpy()
>>> volumetric = cur_system.data["volumetric_strain"].to_numpy()
>>> print(f"Max shear strain: {shear.max():.4f}")
>>> print(f"Mean volumetric strain: {volumetric.mean():.4f}")

mdapy.atomic_temperature module¶

class mdapy.atomic_temperature.AtomicTemperature(data: DataFrame, verlet_list: ndarray, distance_list: ndarray, rc: float, factor: float = 1.0)[source]¶

Bases: object

Calculate atomic temperature from velocity fluctuations, in the units of A/fs.

This class computes the local atomic temperature for each atom by analyzing the velocity fluctuations of the atom and its neighbors. The temperature is calculated from the kinetic energy distribution.

Parameters:

data (pl.DataFrame) – Atomic data containing at least ‘vx’, ‘vy’, ‘vz’ velocity columns and either ‘amass’ (atomic mass) or ‘element’ columns.
verlet_list (np.ndarray) – Neighbor list array of shape (N, max_neigh).
distance_list (np.ndarray) – Distance list array of shape (N, max_neigh).
rc (float) – Cutoff radius for neighbor consideration.
factor (float, default=1.0) – Scaling factor for velocities (e.g., for unit conversion).

data¶

Input atomic data.

Type:: pl.DataFrame

verlet_list¶

Neighbor indices.

Type:: np.ndarray

distance_list¶

Neighbor distances.

Type:: np.ndarray

rc¶

Cutoff radius.

Type:: float

factor¶

Velocity scaling factor.

Type:: float

T¶

Computed atomic temperatures (K) after calling compute().

Type:: np.ndarray

Notes

The atomic temperature represents the local kinetic temperature based on velocity fluctuations. It differs from the thermodynamic temperature and is useful for identifying hot spots, shock fronts, or temperature gradients.

compute() → None[source]¶

Compute atomic temperatures.

This method calculates the atomic temperature for each atom based on velocity fluctuations and stores the result in the T attribute.

Raises:

AssertionError – If velocity columns (‘vx’, ‘vy’, ‘vz’) are not present in data.
ValueError – If neither ‘amass’ nor ‘element’ columns are present, or if an unknown element symbol is encountered.

mdapy.bond_analysis module¶

class mdapy.bond_analysis.BondAnalysis(data: DataFrame, box: Box, rc: float, nbin: int, verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray)[source]¶

Bases: object

Analyze bond length and bond angle distributions in atomic systems.

This class computes histograms of bond lengths (distances between bonded atoms) and bond angles (angles formed by atomic triplets i-j-k where i is the central atom). These distributions are useful for characterizing local structure and identifying coordination patterns.

Parameters:

data (pl.DataFrame) – Atomic data containing positions. Must have columns: ‘x’, ‘y’, ‘z’.
box (Box) – Simulation box object with lattice parameters and boundaries.
rc (float) – Cutoff radius for bond analysis in Angstroms. Pairs within this distance are considered bonded.
nbin (int) – Number of histogram bins for both bond length (0 to rc) and bond angle (0 to 180 degrees) distributions.
verlet_list (np.ndarray) – Verlet neighbor list, shape (N, max_neigh).
distance_list (np.ndarray) – Neighbor distances, shape (N, max_neigh).
neighbor_number (np.ndarray) – Number of neighbors for each atom, shape (N,).

bond_length_distribution¶

Histogram of bond lengths, shape (nbin,). Available after compute().

Type:: np.ndarray

bond_angle_distribution¶

Histogram of bond angles, shape (nbin,). Available after compute().

Type:: np.ndarray

r_length¶

Bin centers for bond length histogram in Angstroms, shape (nbin,). Available after compute().

Type:: np.ndarray

r_angle¶

Bin centers for bond angle histogram in degrees, shape (nbin,). Available after compute().

Type:: np.ndarray

rc¶

Cutoff radius used.

Type:: float

nbin¶

Number of bins used.

Type:: int

Notes

Bond Length Distribution:

For each atom i, counts all neighbors j within cutoff rc. Each distance r_ij is binned into a histogram with bins from 0 to rc.

Bond Angle Distribution:

For each atom i (central atom), considers all pairs of neighbors (j, k) within cutoff rc. The angle θ is calculated as:

\[\theta = \arccos\left(\frac{\vec{r}_{ij} \cdot \vec{r}_{ik}} {|\vec{r}_{ij}| |\vec{r}_{ik}|}\right)\]

where \(\vec{r}_{ij} = \vec{r}_j - \vec{r}_i\) is the distance vector from atom i to atom j. Angles are binned from 0° to 180°.

Examples

Analyze bond structure in a system:

>>> from mdapy import System
>>> system = System("structure.xyz")
>>> system.build_neighbor(rc=5.0, max_neigh=50)
>>> from mdapy.bond_analysis import BondAnalysis
>>> ba = BondAnalysis(
...     system.data,
...     system.box,
...     rc=5.0,
...     nbin=50,
...     verlet_list=system.verlet_list,
...     distance_list=system.distance_list,
...     neighbor_number=system.neighbor_number,
... )
>>> ba.compute()
>>> # Plot distributions
>>> fig1, ax1 = ba.plot_bond_length_distribution()
>>> fig2, ax2 = ba.plot_bond_angle_distribution()

compute()[source]¶

Compute bond length and bond angle distributions.

This method calculates histograms of all bond lengths and bond angles in the system. Results are stored in bond_length_distribution and bond_angle_distribution, with corresponding bin centers in r_length and r_angle.

Both distributions use equal bin spacing: bond lengths from 0 to rc, bond angles from 0° to 180°.

plot_bond_length_distribution(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶

Plot bond length distribution.

Parameters:

fig (Figure, optional) – Existing matplotlib figure. If None, a new figure is created.
ax (Axes, optional) – Existing matplotlib axes. If None, new axes are created.

Returns:

fig (Figure) – Matplotlib figure object.
ax (Axes) – Matplotlib axes object.

plot_bond_angle_distribution(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶

Plot bond angle distribution.

Parameters:

fig (Figure, optional) – Existing matplotlib figure. If None, a new figure is created.
ax (Axes, optional) – Existing matplotlib axes. If None, new axes are created.

Returns:

fig (Figure) – Matplotlib figure object.
ax (Axes) – Matplotlib axes object.

mdapy.bond_stiffness module¶

Bond stiffness vs bond length, modeled after ATAT’s fitsvsl.

Each chemical bond (i,j) between species (A,B) is treated as a 3D harmonic spring whose force-constant matrix in the bond-aligned frame is diag(k_l, k_t, k_t) — one longitudinal (stretch) constant and two identical transverse (bending) constants. Both are assumed to be universal functions of bond length only for a given element pair (the “bond stiffness vs bond length” approximation used by ATAT svsl / fitsvsl; see van de Walle & Ceder, Rev. Mod. Phys. 74, 11, 2002, and Pinsook & Ceder, Phys. Rev. B 60, 11997, 1999):

\[ \begin{align}\begin{aligned}k_l(r) = c^l_0 + c^l_1 r + c^l_2 r^2 + \dots\\k_t(r) = c^t_0 + c^t_1 r + c^t_2 r^2 + \dots\end{aligned}\end{align} \]

Workflow¶

Build the equilibrium neighbor list from the input mdapy.System. Bonds within rc_bond are grouped into distance shells (NN1, NN2, …) using shell_tol; each shell gets its own (k_l, k_t) polynomial. For BCC and similar lattices where the second-neighbor shell carries non-negligible stiffness, set rc_bond large enough to include 2NN.
Generate \(6N\) perturbations: each atom displaced by \(\pm\delta\) along x, y, z. Repeat at n_lattice isotropic strains symmetrically distributed across [-max_strain, +max_strain] to scan a range of bond lengths (needed for poly_order >= 1).
Compute forces with the user-supplied calculator (any mdapy.calculator.CalculatorMP subclass — typically mdapy.NEP).
Solve a single ordinary-least-squares system whose unknowns are the polynomial coefficients (c^l_q, c^t_q) for every (element_pair, shell, q) triple, and whose rows are the per-atom force components for every perturbation. The harmonic-pair model gives the force on atom \(k\) from a perturbation \(d_p\) of atom \(p\) as

\[F_k = \sum_{(i,j)\ \text{containing}\ p,k} \bigl[k_l\, (d_{rel}\!\cdot\!\hat u)\hat u + k_t\, (d_{rel} - (d_{rel}\!\cdot\!\hat u)\hat u)\bigr]\]

with the appropriate sign so that Newton’s third law holds. This is the same formulation used by ATAT fitsvsl and gives a well-conditioned, statistically efficient fit.
Visualise: BondStiffness.plot() reproduces the ATAT-style stiffness vs bond length figure — one panel per shell, with one colour per element pair, scatter for raw data, polynomial overlay.
Write the result in (multi-shell-extended) ATAT slspring.out format for cross-validation or downstream use with svsl.

Citation¶

Using this class should cite the ATAT fitsvsl method papers:

A. van de Walle, M. Asta. Self-driven lattice-model Monte Carlo simulations of alloy thermodynamic properties and phase diagrams. Modelling Simul. Mater. Sci. Eng. 10, 521 (2002).

E. J. Wu, G. Ceder, A. van de Walle. Using bond-length-dependent transferable force constants to predict vibrational entropies in Au-Cu, Au-Pd and Cu-Pd alloys. Phys. Rev. B 67, 134103 (2003).

class mdapy.bond_stiffness.BondStiffness(system: System, calculator: CalculatorMP, rc_bond: float | None = None, shell_tol: float = 0.1, delta: float = 0.05, poly_order: int = 1, n_lattice: int = 3, max_strain: float = 0.02, central_diff: bool = True, rcond: float = 1e-06)[source]¶

Bases: object

Fit bond longitudinal and transverse stiffness vs bond length.

Parameters:

system (mdapy.System) – Reference structure (already relaxed; atoms in equilibrium positions). Must have an element column.
calculator (CalculatorMP) – Force calculator (e.g. mdapy.NEP). Must implement get_forces(data, box) -> (N, 3).
rc_bond (float, optional) – Bond cutoff in Å. All atom pairs with distance below this threshold are considered bonds. Pairs are then automatically grouped into distance shells (NN1, NN2, …) with tolerance shell_tol; each shell is fit independently. If omitted, defaults to 1.05 × the smallest non-zero pairwise distance in system — a 1NN-only cutoff. For BCC and similar lattices where 2NN contributes meaningfully, pass an rc_bond large enough to include both shells (e.g. rc_bond = 1.2 * a for BCC catches NN1 ≈ a√3/2 and NN2 = a).
shell_tol (float, default=0.1) – Distance tolerance (Å) for grouping bonds into the same shell. Two bonds within shell_tol of each other share a shell index; bonds that move within their shell across strains provide the (r, k) data points for a length-dependent fit.
delta (float, default=0.05) – Magnitude of atomic displacement used to probe forces, in Å. Use central differences (+delta and -delta per axis) so the linear-response coefficient is extracted exactly even with a small anharmonic part.
poly_order (int, default=1) – Polynomial order of the k(r) fit. 0 gives a single constant per (element pair, shell). For poly_order >= 1 set n_lattice >= 2 to actually scan bond lengths.
n_lattice (int, default=3) – Number of isotropic strains at which the perturbation set is repeated. Strains are spaced linearly and symmetrically across [-max_strain, +max_strain] (so n_lattice=1 gives [0], n_lattice=2 gives [-max_strain, +max_strain], n_lattice=3 gives [-max_strain, 0, +max_strain], etc.). This contrasts with ATAT fitsvsl’s -ns which only stretches; we scan symmetrically since extrapolating to compressed states is just as physical as expanded ones.
max_strain (float, default=0.02) – Maximum isotropic strain magnitude. With the default n_lattice=3 this scans bond lengths over ±2 %.
central_diff (bool, default=True) – If True, displace each atom by both \(+\delta\) and \(-\delta\) per axis and use central differences. Doubles the force calculations but makes the linear response exact in the limit of small \(\delta\). If False, displace only \(+\delta\) (forward differences; cheaper, less accurate).
rcond (float, default=1e-6) – Relative cutoff for singular values in the OLS pseudo-inverse. Singular values below rcond * max(sv) are treated as zero, regularising the fit when polynomial-feature columns are nearly collinear (e.g. poly_order >= 1 with a tiny strain range).

bond_table¶

One row per bond per perturbation strain: columns element_a, element_b, shell, r, k_long, k_trans, strain. Set by compute().

Type:: polars.DataFrame

shells¶

Mean diameter (in Å) of each distance shell, sorted from nearest to farthest. len(shells) is the number of shells retained.

Type:: list[float]

k_long, k_trans

Fitted polynomial coefficients (lowest order first) keyed by (element_a, element_b, shell). Element pair is sorted alphabetically so ("Al","Cr",0) and ("Cr","Al",0) always map to the same key. Set by compute().

Type:: dict[tuple[str, str, int], numpy.ndarray]

compute() → BondStiffness[source]¶

Run the full pipeline and populate k_long, k_trans, bond_table. Returns self.

The fit is set up exactly as in ATAT fitsvsl: a single ordinary least-squares system whose unknowns are the polynomial coefficients (c^l_0,...,c^l_p, c^t_0,...,c^t_p) for each element pair, NOT per-bond stiffnesses. This is the well- conditioned formulation that gives a numerically stable answer even when many bonds of the same chemical pair share identical lengths.

write_slspring(path: str) → None[source]¶

Write fitted (k_l, k_t) polynomials to path in ATAT slspring.out format (extended for multi-shell output).

With one shell, the file is byte-for-byte compatible with what ATAT fitsvsl writes. With multiple shells, each (element pair, shell) becomes a separate block annotated with # shell <id> d=<diameter> so you can tell them apart while keeping the section structure ATAT svsl parses.

File layout:

<element_a> <element_b>     # shell 0 d=<NN1 distance>
<n_coeffs>
c_l0
...
<n_coeffs>
c_t0
...
<element_a> <element_b>     # shell 1 d=<NN2 distance>
...

plot(which: str = 'both', ax=None, ncol: int | None = None)[source]¶

Plot fitted bond stiffness vs bond length.

Produces one panel per element pair, with all distance shells (NN1, NN2, …) drawn on the same axes. Different shells use different scatter markers; each shell’s polynomial fit is overlaid as a continuous curve within that shell’s r range. This reproduces the canonical “stiffness vs bond length” figure of Wu, Ceder, van de Walle (MSMSE 10, 521, 2002), where all shells of a chemical pair share a single panel and the x-axis spans the full r range covered by the chosen rc_bond.

To get a wider x-range like the published figures (~2.5 to 5.0 Å), pass an rc_bond covering NN1..NN4 (e.g. ~1.4 × the lattice parameter for fcc).

Parameters:

which ({'both', 'long', 'trans'}, default='both') – Which stiffness mode to plot. 'both' shows longitudinal and transverse on the same panel (different line styles).
ax (matplotlib Axes or list of Axes, optional) – Pre-existing axes to plot into. If omitted, a fresh figure is created via mdapy.set_figure().
ncol (int, optional) – Number of subplot columns when multiple element pairs are plotted and ax is None. Defaults to min(3, n_pairs).

Returns:

fig, axes (matplotlib Figure, list of Axes) – One axes per element pair.

generate_perturbed_structures(output_dir: str | None = None) → List['System'][source]¶

Generate the list of perturbed mdapy.System objects used internally by compute(). If output_dir is given, also dump each as an ATAT-format str.out file (and a single str_unpert.out) ready for an external DFT or potential run, plus a perturb_index.txt describing the perturbation applied to each.

Returns the list of System objects in the same order they are used by compute(). Note this enumerates only the n_lattice=1 (unstrained) case; for multi-strain workflows call this method with different self.max_strain / n_lattice settings.

mdapy.box module¶

Bases: object

The Box class represents a parallelepiped simulation cell that can be either orthogonal or triclinic (non-orthogonal). It supports periodic boundary conditions and provides methods for coordinate transformations and box manipulations.

Parameters:

box (int, float, Iterable[float], np.ndarray, or Box) –
Defines the box vectors or size:
- int/float: Creates a cubic box with given edge length
- Iterable[float] of length 3: Creates an orthogonal box with given dimensions [Lx, Ly, Lz]
- np.ndarray of shape (3,3): Full box matrix with vectors as rows
- np.ndarray of shape (4,3): Box matrix with origin as 4th row
- Box: Creates a deep copy of an existing Box object
boundary (Iterable[int] or np.ndarray, optional) –
Boundary condition flags for each dimension (x, y, z):
- 1: Periodic boundary condition
- 0: Fixed (non-periodic) boundary
Defaults to [1, 1, 1] (all periodic).
origin (Iterable[float] or np.ndarray, optional) – The origin position of the box in 3D space. Defaults to [0.0, 0.0, 0.0].

Examples

Creating different types of boxes:

import numpy as np
from mdapy.box import Box

# Cubic box with edge length 10
cubic_box = Box(10)

# Orthogonal box with different dimensions
ortho_box = Box([10, 20, 30])

# Triclinic box from full matrix
matrix = np.array([[10, 0, 0],
                  [5, 10, 0],
                  [0, 0, 10]])
triclinic_box = Box(matrix)

# Box with custom origin and boundary conditions
custom_box = Box([20, 20, 20],
                boundary=[1, 1, 0],  # z-direction non-periodic
                origin=[5, 5, 5])

# Copy an existing box
box_copy = Box(cubic_box)

Applying periodic boundary conditions:

# Create a box
box = Box(10)

# Apply PBC to a displacement vector
rij = np.array([12, -8, 5])  # Vector that may cross boundaries
rij_pbc = box.pbc(rij)
print(f"Original: {rij}")
print(f"After PBC: {rij_pbc}")

Notes

The box matrix convention follows the standard where box vectors are rows:

\[\begin{split}\mathbf{B} = \begin{pmatrix} a_x & a_y & a_z \\ b_x & b_y & b_z \\ c_x & c_y & c_z \end{pmatrix}\end{split}\]

set_box(box: int | float | Iterable[float] | ndarray) → None[source]¶

Set box matrix.

Parameters:: box (int, float, Iterable[float], or 3x3 np.ndarray) – Box specification in various formats.

property volume: float¶

no-index:

Get the volume of the simulation box.

Returns:: float – Box volume in cubic units.

property box: ndarray¶

no-index:

Get the box matrix.

Returns:: np.ndarray – 3x3 matrix with box vectors as rows.

property inverse_box: ndarray¶

no-index:

Get the inverse box matrix.

Returns:: np.ndarray – 3x3 inverse box matrix.

property origin: ndarray¶

no-index:

Get the origin position of the box.

Returns:: np.ndarray – 3D origin position.

set_origin(origin: Iterable[float] | ndarray) → None[source]¶

Set origin coordinates.

Parameters:: origin (Iterable[float] or np.ndarray) – Origin position.

property boundary: ndarray¶

no-index:

Get the boundary condition flags.

Returns:: np.ndarray – Array of boundary flags (1=periodic, 0=fixed).

set_boundary(boundary: Iterable[int] | ndarray) → None[source]¶

Set boundary conditions.

Parameters:: boundary (Iterable[int] or np.ndarray) – Boundary condition flags (1=periodic, 0=fixed).

property triclinic: bool¶

no-index:

Check if the box is triclinic.

Returns:: bool – True if triclinic, False if orthogonal.

is_general_box(tol: float = 1e-06) → bool[source]¶

Determine whether the simulation box is a general (non-lower-triangular) box.

A box is considered general if it does not satisfy the lower-triangular form required by LAMMPS, i.e.,

[[ax, 0, 0], [bx, by, 0], [cx, cy, cz]]

within a given numerical tolerance.

Parameters:: tol (float, optional) – Numerical tolerance used to judge whether a value is effectively zero. Default is 1e-6.
Returns:: bool – True if the box is a general box; False if it is already in lower-triangular form.

align_to_lammps_box() → Tuple[Box, ndarray][source]¶

Transform the box to LAMMPS-compatible format.

Returns:

Tuple[Box, np.ndarray] –

Aligned Box object in LAMMPS format
3x3 rotation matrix used for the transformation

pbc(rij: ndarray) → ndarray[source]¶

Apply periodic boundary conditions to a displacement vector.

Parameters:: rij (np.ndarray) – 3D displacement vector in Cartesian coordinates.
Returns:: np.ndarray – Wrapped displacement vector following minimum image convention.

get_thickness() → ndarray[source]¶

Calculate the perpendicular thickness along each box dimension.

Returns:: np.ndarray – Array of thicknesses [tx, ty, tz] for each dimension.

check_small_box(rc: float) → ndarray[source]¶

Check if the box satisfies the minimum image convention for a given cutoff.

Parameters:: rc (float) – Cutoff radius for interactions.
Returns:: np.ndarray – Integer array [nx, ny, nz] indicating required replications.

mdapy.build_bond module¶

mdapy.build_bond.build_bond(verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray, type_list: ndarray, cutoff_matrix: ndarray) → ndarray[source]¶: Build bond pairs from precomputed neighbor information.

mdapy.build_lattice module¶

mdapy.build_lattice.build_crystal(name, structure: str, a: float, miller1: Tuple[int, int, int] | None = None, miller2: Tuple[int, int, int] | None = None, miller3: Tuple[int, int, int] | None = None, nx: int = 1, ny: int = 1, nz: int = 1, c: float | None = None) → System[source]¶

Build a crystal supercell. Output matches atomsk’s --create <structure> byte-for-byte (after sorting on Cartesian position).

Parameters:

name (str or sequence of str) – Element symbols. Pass a single string to broadcast it to every basis atom (single-species crystal); pass a tuple/list whose length equals one of the structure’s allowed species counts for an ordered multi-species structure (NaCl, GaN, SrTiO3, …). Each element must be a recognised symbol from mdapy.data.atomic_numbers.
structure (str) –
Crystal type. Case-insensitive; common synonyms accepted ("nacl" → "rocksalt", "zb" → "zincblende", "b2" → "cscl", "hex_diamond" → "lonsdaleite", …). Supported structures and their allowed species counts:

Cubic (Miller-rotatable, 3-index [h, k, l]):
- sc (1): simple cubic, e.g. α-Po
- fcc (1 or 2): Cu, Al, Ni; or ordered FCC
- bcc (1 or 2): Fe, W; B2 alloy with two species
- diamond (1 or 2): C, Si; zincblende for two species
- cscl / B2 (2): CsCl, NiAl
- rocksalt / B1 (2): NaCl, MgO
- zincblende / B3 (2): GaAs, ZnS
- fluorite (2): CaF2, UO2
- l1_2 (2): Ni3Al, Cu3Au
- perovskite (3, B-A-O ordering): SrTiO3 → ("Ti","Sr","O")
Hexagonal (Miller-rotatable, 3-index [u, v, w] or 4-index [h, k, i, l]):
- hcp (1 or 2): Mg, Ti, Zr
- wurtzite / B4 (1 or 2): GaN, ZnO
- lonsdaleite / hex-diamond (1): hexagonal carbon polytype
- graphite (1 or 2): graphite, hex-BN — c required
- graphene (1 or 2): single-layer honeycomb — c required
a (float) – Primary lattice constant in Å (cubic edge, hexagonal in-plane spacing).
miller1 (sequence of int, optional) –
Crystallographic axes for the three sides of the supercell. For cubic structures pass three integers [h, k, l]; for hexagonal structures pass either Miller-Bravais [h, k, i, l] (with h + k + i == 0) or 3-index [u, v, w] directly. When all three are None, the canonical orientation is used.

Cubic Miller indices must be mutually orthogonal in Cartesian space; left-handed frames are accepted and yield a mirror-image crystal (a lattice symmetry for every supported cubic structure). Hexagonal Miller indices only need to be linearly independent — non-orthogonal frames are returned as triclinic boxes in lower-triangular form. This matches atomsk’s --create ... orient semantics.
miller2 (sequence of int, optional) –
Crystallographic axes for the three sides of the supercell. For cubic structures pass three integers [h, k, l]; for hexagonal structures pass either Miller-Bravais [h, k, i, l] (with h + k + i == 0) or 3-index [u, v, w] directly. When all three are None, the canonical orientation is used.

Cubic Miller indices must be mutually orthogonal in Cartesian space; left-handed frames are accepted and yield a mirror-image crystal (a lattice symmetry for every supported cubic structure). Hexagonal Miller indices only need to be linearly independent — non-orthogonal frames are returned as triclinic boxes in lower-triangular form. This matches atomsk’s --create ... orient semantics.
miller3 (sequence of int, optional) –
Crystallographic axes for the three sides of the supercell. For cubic structures pass three integers [h, k, l]; for hexagonal structures pass either Miller-Bravais [h, k, i, l] (with h + k + i == 0) or 3-index [u, v, w] directly. When all three are None, the canonical orientation is used.

Cubic Miller indices must be mutually orthogonal in Cartesian space; left-handed frames are accepted and yield a mirror-image crystal (a lattice symmetry for every supported cubic structure). Hexagonal Miller indices only need to be linearly independent — non-orthogonal frames are returned as triclinic boxes in lower-triangular form. This matches atomsk’s --create ... orient semantics.
nx (int, default=1) – Replication counts along the (possibly Miller-rotated) axes. Total atom count = len(basis) * nx * ny * nz.
ny (int, default=1) – Replication counts along the (possibly Miller-rotated) axes. Total atom count = len(basis) * nx * ny * nz.
nz (int, default=1) – Replication counts along the (possibly Miller-rotated) axes. Total atom count = len(basis) * nx * ny * nz.
c (float, optional) – Out-of-plane lattice constant in Å for hexagonal structures. For hcp, wurtzite and lonsdaleite defaults to a * sqrt(8/3) (ideal close-packing). For graphite and graphene it must be supplied — no sensible default exists (interlayer / vacuum spacing).

Returns:

System – A System with the per-atom DataFrame populated with x, y, z, element and an orthogonal-or-triclinic simulation box matching atomsk’s convention for the chosen structure.

Examples

Single-species cubic — copper FCC:

>>> cu = build_crystal("Cu", "fcc", a=3.615, nx=5, ny=5, nz=5)

Two-species cubic — NaCl:

>>> nacl = build_crystal(("Na", "Cl"), "rocksalt", a=5.64, nx=4, ny=4, nz=4)

Three-species perovskite (B, A, O ordering):

>>> sto = build_crystal(("Ti", "Sr", "O"), "perovskite", a=3.905)

Hexagonal — wurtzite GaN with explicit c:

>>> gan = build_crystal(("Ga", "N"), "wurtzite", a=3.19, c=5.18)

Hexagonal diamond (lonsdaleite) — single-species wurtzite:

>>> lons = build_crystal("C", "lonsdaleite", a=2.51, c=4.12)

Single-layer graphene (vacuum-padded along c):

>>> g = build_crystal("C", "graphene", a=2.46, c=20.0)

Cubic Miller-rotated — FCC [111] slab orientation:

>>> cu_111 = build_crystal(
...     "Cu", "fcc", a=3.615,
...     miller1=(1, -1, 0), miller2=(1, 1, -2), miller3=(1, 1, 1),
... )

Hexagonal Miller-Bravais — Mg prismatic (11-20) plane:

>>> mg = build_crystal(
...     "Mg", "hcp", a=3.21, c=5.21,
...     miller1=(1, -1, 0, 0), miller2=(1, 1, -2, 0), miller3=(0, 0, 0, 1),
... )

mdapy.build_lattice.build_hea(element_list: Tuple[str], element_ratio: Tuple[float], structure: str, a: float, miller1: Tuple[int, int, int] | None = None, miller2: Tuple[int, int, int] | None = None, miller3: Tuple[int, int, int] | None = None, nx: int = 1, ny: int = 1, nz: int = 1, c: float | None = None, random_seed: int | None = None) → System[source]¶

Build a random high-entropy alloy (HEA) on a single sublattice.

Parameters:

element_list (tuple of str) – Element symbols, e.g. ('Cu', 'Al', 'Mg'). Must be unique.
element_ratio (tuple of float) – Corresponding atomic ratios; must sum to 1.
structure (str) – Single-species crystal type — typically 'fcc', 'bcc' or 'hcp'. Any structure from build_crystal() whose allowed species count includes 1 will work; only one sublattice is randomly substituted.
a (float) – Primary lattice constant in Å.
miller1 (sequence of int, optional) – Crystallographic axes (3-index for cubic, 3- or 4-index for hexagonal). See build_crystal() for the conventions.
miller2 (sequence of int, optional) – Crystallographic axes (3-index for cubic, 3- or 4-index for hexagonal). See build_crystal() for the conventions.
miller3 (sequence of int, optional) – Crystallographic axes (3-index for cubic, 3- or 4-index for hexagonal). See build_crystal() for the conventions.
nx (int, default=1) – Replication counts along each box axis.
ny (int, default=1) – Replication counts along each box axis.
nz (int, default=1) – Replication counts along each box axis.
c (float, optional) – Out-of-plane lattice constant for hexagonal structures.
random_seed (int, optional) – Seed for the per-atom element assignment.

Returns:

System – Crystal with the element column randomly populated by element_list according to element_ratio.

Examples

>>> sys_ = build_hea(
...     ("Cr", "Co", "Ni"), (0.1, 0.2, 0.7),
...     "fcc", 3.526, nx=5, ny=5, nz=5, random_seed=1,
... )

mdapy.build_lattice.build_hea_fromsystem(system: System, element_list: Tuple[str], element_ratio: Tuple[float], random_seed: int | None = None) → System[source]¶

Generate element per atom for a system,

Parameters:

system (System) – System need to be modified.
element_list (Tuple[str]) – List of element symbols, e.g., ('Cu', 'Al', 'Mg').
element_ratio (Tuple[float]) – Corresponding atomic ratios. The sum must equal 1.
random_seed (int, optional) – Random seed for reproducible element assignment.

Returns:

System – A System object containing the generated element information.

mdapy.build_lattice.build_partial_dislocation_fcc(name: str, a: float, nx: int, ny: int, nz: int, element_list: Tuple[str] | None = None, element_ratio: Tuple[float] | None = None, random_seed: int | None = None) → System[source]¶

Generate a FCC strcuture with a pair partial dislocations along y axis. The crystalline orientation is: x->[1-10], y->[11-2], z->[111]. If element_list and element_ratio are given, it will generate a HEA. Similar to Construct a dislocation by superimposing two crystals in atomsk: https://atomsk.univ-lille.fr/tutorial_Al_edge.php

Parameters:

name (str) – Element symbol (e.g., ‘Cu’, ‘Al’, ‘Mg’, ‘C’).
a (float) – Lattice constant in Angstroms.
nx (int, default=1) – Number of repetitions along x-axis direction.
ny (int, default=1) – Number of repetitions along y-axis direction.
nz (int, default=1) – Number of repetitions along z-axis direction.
element_list (Tuple[str], optional) – List of element symbols, e.g., ('Cu', 'Al', 'Mg').
element_ratio (Tuple[float], optional) – Corresponding atomic ratios. The sum must equal 1.
random_seed (int, optional) – Random seed for reproducible element assignment.

Returns:

System – A System object containing the dislocations.

mdapy.calculator module¶

Abstract Base Class for Static Calculators¶

This module defines the abstract base class for implementing static calculators in the mdapy project. All calculator implementations should inherit from this class and implement the required abstract methods.

Classes¶

CalculatorMP : Abstract base class for MD calculators

class mdapy.calculator.CalculatorMP[source]¶

Bases: ABC

Abstract base class for static property calculators.

This class defines the interface that all calculators must implement. Calculators are responsible for computing physical properties such as energies, forces, stresses, and virials for a given atomic configuration.

Notes

All subclasses must implement the five abstract methods defined here. The data parameter should contain atomic positions and elements, while the box parameter describes the simulation cell geometry.

Examples

>>> class MyCalculator(CalculatorMP):
...     def get_energies(self, data, box):
...         # Implementation here
...         pass
...
...     # ... implement other methods

results = {}¶

abstractmethod get_energies(data: DataFrame, box: Box) → ndarray[source]¶

Calculate per-atom potential energies.

Parameters:

data (pl.DataFrame) –
Polars DataFrame containing atomic information. Must include:
- ’x’, ‘y’, ‘z’: atomic coordinates (float)
- ’element’: atomic species information
box (Box) – Simulation box object containing cell dimensions and boundary conditions

Returns:

np.ndarray – 1D array of shape (N,) containing potential energy for each atom, where N is the number of atoms

abstractmethod get_energy(data: DataFrame, box: Box) → float[source]¶

Calculate total potential energy of the system.

Parameters:

data (pl.DataFrame) –
Polars DataFrame containing atomic information. Must include:
- ’x’, ‘y’, ‘z’: atomic coordinates (float)
- ’element’: atomic species information
box (Box) – Simulation box object containing cell dimensions and boundary conditions

Returns:

float – Total potential energy of the system (sum of all per-atom energies)

Notes

This is typically implemented as return self.get_energies(data, box).sum()

abstractmethod get_forces(data: DataFrame, box: Box) → ndarray[source]¶

Calculate forces acting on each atom.

Parameters:

data (pl.DataFrame) –
Polars DataFrame containing atomic information. Must include:
- ’x’, ‘y’, ‘z’: atomic coordinates (float)
- ’element’: atomic species information
box (Box) – Simulation box object containing cell dimensions and boundary conditions

Returns:

np.ndarray – 2D array of shape (N, 3) containing force components [fx, fy, fz] for each atom, where N is the number of atoms

abstractmethod get_stress(data: DataFrame, box: Box) → ndarray[source]¶

Calculate stress tensor of the system.

Parameters:

data (pl.DataFrame) –
Polars DataFrame containing atomic information. Must include:
- ’x’, ‘y’, ‘z’: atomic coordinates (float)
- ’element’: atomic species information
box (Box) – Simulation box object containing cell dimensions and boundary conditions

Returns:

np.ndarray – 1D array of shape (6,) containing stress tensor components in Voigt notation: [σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy]

Notes

The stress tensor is symmetric, so only 6 independent components are returned. Stress is typically computed from virial: σ = -(virial + virial^T) / (2 * volume)

abstractmethod get_virials(data: DataFrame, box: Box) → ndarray[source]¶

Calculate per-atom virial tensors.

Parameters:

data (pl.DataFrame) –
Polars DataFrame containing atomic information. Must include:
- ’x’, ‘y’, ‘z’: atomic coordinates (float)
- ’element’: atomic species information
box (Box) – Simulation box object containing cell dimensions and boundary conditions

Returns:

np.ndarray – 2D array of shape (N, 9) containing virial tensor components for each atom. Components ordered as: [v_xx, v_xy, v_xz, v_yx, v_yy, v_yz, v_zx, v_zy, v_zz] where N is the number of atoms

mdapy.centro_symmetry_parameter module¶

class mdapy.centro_symmetry_parameter.CentroSymmetryParameter(data: DataFrame, box: Box, N: int, verlet_list: ndarray)[source]¶

Bases: object

Calculate the Centro-Symmetry Parameter (CSP) for detecting crystal defects.

The Centro-Symmetry Parameter is a structural metric that measures the local symmetry of an atom’s environment. It is particularly useful for identifying defects in centro-symmetric structures like FCC and BCC crystals.

Parameters:

data (pl.DataFrame) – Atomic data containing at least ‘x’, ‘y’, ‘z’ position columns.
box (Box) – Simulation box object.
N (int) – Number of nearest neighbors to consider. Must be a positive even number. Typical values are 12 for FCC and 8 for BCC structures.
verlet_list (np.ndarray) – Neighbor list array of shape (n_atoms, max_neigh) containing neighbor indices.

data¶

Input atomic data.

Type:: pl.DataFrame

box¶

Simulation box.

Type:: Box

N¶

Number of neighbors used in calculation.

Type:: int

verlet_list¶

Neighbor indices.

Type:: np.ndarray

csp¶

Centro-symmetry parameters after calling compute(). Shape (n_atoms,). Lower values indicate higher symmetry (perfect crystal), while higher values indicate defects, surfaces, or disordered regions.

Type:: np.ndarray

Notes

The CSP is calculated as:

\[CSP = \sum_{i=1}^{N/2} |\mathbf{r}_i + \mathbf{r}_{i+N/2}|^2\]

where \(\mathbf{r}_i\) and \(\mathbf{r}_{i+N/2}\) are vectors to opposite nearest neighbors.

For perfect centro-symmetric crystals, CSP ≈ 0. Typical CSP values: - Perfect FCC/BCC bulk: 0-0.1 - Surfaces/interfaces: 1-10 - Dislocations: 10-50 - Amorphous regions: > 50

References

Kelchner, C. L., Plimpton, S. J., & Hamilton, J. C. (1998). Dislocation nucleation and defect structure during surface indentation. Physical Review B, 58(17), 11085.

compute() → None[source]¶

Compute the Centro-Symmetry Parameter for all atoms.

This method calculates the CSP value for each atom and stores the result in the csp attribute.

Notes

After calling this method, the csp attribute will contain a float array of CSP values, one for each atom.

mdapy.chill_plus module¶

class mdapy.chill_plus.ChillPlus(data: DataFrame, box: Box, cutoff: float = 3.5, verlet_list: ndarray[tuple[Any, ...], dtype[int32]] | None = None, distance_list: ndarray[tuple[Any, ...], dtype[float64]] | None = None, neighbor_number: ndarray[tuple[Any, ...], dtype[int32]] | None = None)[source]¶

Bases: object

CHILL+ structure identification for water phases.

Classifies each particle (oxygen / coarse-grained water bead) as one of six phases by examining the bond-orientation correlation \(c_{ij}\) between neighbouring molecules within cutoff:

0 – Other
1 – Hexagonal ice
2 – Cubic ice
3 – Interfacial ice
4 – Gas hydrate
5 – Interfacial gas hydrate

The cutoff defines the O-O distance below which two molecules are considered bonded (default 3.5 Å). The algorithm operates only on the particles supplied – pass oxygen / water-bead positions, not hydrogens.

Parameters:

data (pl.DataFrame) – Particle data with at least ‘x’, ‘y’, ‘z’ columns.
box (Box) – Simulation box.
cutoff (float, optional) – O-O cutoff radius in Å. Default 3.5.
verlet_list (np.ndarray, optional) – Pre-computed cutoff neighbour list. When omitted, an internal Neighbor build is performed at cutoff.
distance_list (np.ndarray, optional) – Pre-computed cutoff neighbour list. When omitted, an internal Neighbor build is performed at cutoff.
neighbor_number (np.ndarray, optional) – Pre-computed cutoff neighbour list. When omitted, an internal Neighbor build is performed at cutoff.

pattern¶

Structure code per particle (filled by compute()).

Type:: np.ndarray of int32, shape (N,)

References

Nguyen, A. H. & Molinero, V. Identification of Clathrate Hydrates, Hexagonal Ice, Cubic Ice, and Liquid Water in Simulations: the CHILL+ Algorithm. J. Phys. Chem. B 119 (2015) 9369-9376. DOI: 10.1021/jp510289t

compute() → None[source]¶: Run the CHILL+ classification, populating self.pattern.

mdapy.cluster_analysis module¶

class mdapy.cluster_analysis.ClusterAnalysis(rc: float | int | Dict[str, float], verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray, type_list: ndarray | None = None)[source]¶

Bases: object

Perform cluster analysis based on atomic connectivity within a cutoff distance.

This class identifies connected groups of atoms (“clusters”) by examining neighbor relationships. Two atoms belong to the same cluster if their interatomic distance is smaller than a given cutoff value.

The cutoff can be either:

a single scalar value (applied to all atom pairs), or
a dictionary defining different cutoffs between atom types, e.g. {"1-1": 1.2, "1-2": 1.5, "2-2": 1.8}.

Parameters:

rc (float, int, or dict of str to float) – The cutoff distance for bonding. If a dict is provided, the keys should be in the form "type1-type2" (e.g. "1-2"), and type_list must be given.
verlet_list (np.ndarray) – The neighbor list array, where each row contains neighbor indices for each atom.
distance_list (np.ndarray) – The corresponding neighbor distances for each atom.
neighbor_number (np.ndarray) – The number of neighbors for each atom.
type_list (np.ndarray, optional) – Atom type array. Required if rc is a dict specifying per-type cutoffs.

max_rc¶

The maximum cutoff distance among all pair types.

Type:: float

particleClusters¶

The cluster ID assigned to each atom after computation.

Type:: np.ndarray

cluster_number¶

The total number of clusters identified.

Type:: int

compute()[source]¶

Perform the actual cluster analysis.

This function assigns a unique cluster ID to each atom based on connectivity within the given cutoff distance.

Returns:

None – The results are stored in the attributes:

particleClusters: array of cluster IDs.
cluster_number: total number of clusters found.

mdapy.common_neighbor_analysis module¶

class mdapy.common_neighbor_analysis.CommonNeighborAnalysis(data: DataFrame, box: Box, verlet_list: ndarray | None = None, neighbor_number: ndarray | None = None, rc: float | None = None)[source]¶

Bases: object

Perform Common Neighbor Analysis (CNA) on atomic structures.

This class classifies each atom according to its local atomic environment using either an adaptive or fixed cutoff approach.

Two modes are supported:

Adaptive cutoff (default): determines neighbors based on the 14 nearest atoms.
Fixed cutoff: uses a user-specified cutoff distance rc to define neighbors.

After computation, each atom is assigned a coordination pattern stored in pattern:

0 = Other (unknown coordination)
1 = FCC (face-centered cubic)
2 = HCP (hexagonal close-packed)
3 = BCC (body-centered cubic)
4 = ICO (icosahedral coordination)

Parameters:

data (pl.DataFrame) – Atomic coordinates stored as a Polars DataFrame with columns [‘x’, ‘y’, ‘z’].
box (Box) – Simulation box object (supports periodic boundaries).
verlet_list (Optional[np.ndarray], default=None) – Precomputed neighbor indices (optional). If not provided, neighbors will be determined automatically.
neighbor_number (Optional[np.ndarray], default=None) – Number of neighbors for each atom (required if verlet_list is provided).
rc (Optional[float], default=None) – Cutoff radius for fixed-cutoff CNA. If None, adaptive cutoff is used.

pattern¶

Array of integer codes indicating the coordination type of each atom.

Type:: np.ndarray

Notes

Adaptive cutoff automatically finds neighbors based on the 14 nearest atoms.
Fixed cutoff requires the rc parameter and uses all neighbors within rc.

References

Faken D, Jónsson H. Systematic analysis of local atomic structure combined with 3D computer graphics[J]. Computational Materials Science, 1994, 2(2): 279-286.
Stukowski A. Structure identification methods for atomistic simulations of crystalline materials[J]. Modelling and Simulation in Materials Science and Engineering, 2012, 20(4): 045021.

compute()[source]¶

This method fills self.pattern with integer codes representing the local coordination environment of each atom.

The method automatically handles small simulation boxes by replicating atoms as needed to ensure sufficient neighbors.

mdapy.common_neighbor_parameter module¶

class mdapy.common_neighbor_parameter.CommonNeighborParameter(data: DataFrame, box: Box, rc: float, verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray)[source]¶

Bases: object

Calculate the Common Neighbor Parameter (CNP) for structural analysis.

The Common Neighbor Parameter quantifies the local structural order by analyzing the bonding topology between neighbors. It is effective for distinguishing between crystalline and disordered regions.

Parameters:

data (pl.DataFrame) – Atomic data containing at least ‘x’, ‘y’, ‘z’ position columns.
box (Box) – Simulation box object.
rc (float) – Cutoff radius for neighbor bonds. Must be positive.
verlet_list (np.ndarray) – Neighbor list array of shape (N, max_neigh).
distance_list (np.ndarray) – Distance list array of shape (N, max_neigh).
neighbor_number (np.ndarray) – Number of neighbors for each atom, shape (N,).

data¶

Input atomic data.

Type:: pl.DataFrame

box¶

Simulation box.

Type:: Box

rc¶

Cutoff radius.

Type:: float

verlet_list¶

Neighbor indices.

Type:: np.ndarray

distance_list¶

Neighbor distances.

Type:: np.ndarray

neighbor_number¶

Neighbor counts.

Type:: np.ndarray

cnp¶

Common Neighbor Parameters after calling compute(). Shape (N,). Values closer to 0 indicate more ordered/crystalline structures, while larger values indicate disordered regions.

Type:: np.ndarray

Notes

The CNP is calculated by analyzing the connectivity graph of nearest neighbors. For each atom, it examines which of its neighbors are also neighbors of each other, forming a local bonding topology signature.

Typical CNP interpretations: - Crystalline structures (FCC, BCC, HCP): Low CNP values (< 1) - Grain boundaries: Intermediate values (1-3) - Liquid/amorphous: Higher values (> 3)

References

Tsuzuki, H., Branicio, P. S., & Rino, J. P. (2007). Structural characterization of deformed crystals by analysis of common atomic neighborhood. Computer Physics Communications, 177(6), 518-523.

compute() → None[source]¶

Compute the Common Neighbor Parameter for all atoms.

This method calculates the CNP value for each atom based on its neighbor topology and stores the result in the cnp attribute.

Notes

After calling this method, the cnp attribute will contain a float array with CNP values for each atom.

mdapy.create_polycrystal module¶

class mdapy.create_polycrystal.CreatePolycrystal(unitcell: System, box: int | float | Iterable[float] | ndarray | Box, seed_number: int, seed_position: ndarray | None = None, theta_list: ndarray | None = None, randomseed: int | None = None, metal_overlap_dis: float | None = None, add_graphene: bool = False, metal_gra_overlap_dis: float = 3.0, face_threshold: float = 0.0, need_rotation: bool = True)[source]¶

Bases: object

Create polycrystalline structures using Voronoi tessellation method.

This class generates polycrystalline atomic structures by: 1. Creating Voronoi cells from seed positions 2. Filling each cell with rotated unit cell atoms 3. Optionally decorating grain boundaries with graphene layers 4. Removing overlapping atoms based on distance criteria

Parameters:

unitcell (System) – The unit cell structure to be replicated in each grain.
box (Union[int, float, Iterable[float], np.ndarray, Box]) – Simulation box dimensions. Can be a single value (cubic), array of 3 values, or Box object.
seed_number (int) – Number of grains (Voronoi seeds) to generate.
seed_position (Optional[np.ndarray], default=None) – Positions of Voronoi seeds with shape (seed_number, 3). If None, randomly generated within the box.
theta_list (Optional[np.ndarray], default=None) – Rotation angles (in degrees) for each grain with shape (seed_number, 3). Represents rotations around x, y, z axes. If None, randomly generated.
randomseed (Optional[int], default=None) – Random seed for reproducibility. If None, randomly chosen.
metal_overlap_dis (Optional[float], default=None) – Minimum allowed distance between metal atoms (Angstrom). Atoms closer than this will be removed.
add_graphene (bool, default=False) – Whether to add graphene layers at grain boundaries.
metal_gra_overlap_dis (float, default=3.0) – Minimum allowed distance between metal and carbon atoms (Angstrom).
face_threshold (float, default=0.0) – Minimum face area (Angstrom²) for graphene decoration. Faces smaller than this will be skipped.
need_rotation (bool, default=True) – Whether to apply random rotations to grains.

unitcell¶

The unit cell structure.

Type:: System

box¶

The simulation box.

Type:: Box

seed_number¶

Number of grains.

Type:: int

seed_position¶

Voronoi seed positions.

Type:: np.ndarray

theta_list¶

Rotation angles for each grain.

Type:: np.ndarray

con¶

Voronoi container with all cells.

Type:: Container

randomseed¶

Random seed used.

Type:: int

Examples

>>> from mdapy.build_lattice import build_crystal
>>> unit = build_crystal("Al", "fcc", 4.05)
>>> poly = CreatePolycrystal(unit, box=100, seed_number=10, metal_overlap_dis=2.0)
>>> system = poly.compute()
>>> system.write_xyz("polycrystal.xyz")

Notes

Triclinic boxes are not supported
Free boundary conditions are not supported

compute(verbose: bool = True) → System[source]¶

Compute and generate the polycrystalline structure.

Main workflow:

Generate Voronoi tessellation from seeds
Fill each cell with rotated unit cell atoms
Optionally add graphene at grain boundaries
Remove overlapping atoms
Wrap atoms into simulation box

Parameters:

verbose (bool, default=True) – If True, print detailed progress information. If False, suppress all output.

Returns:

System – The generated polycrystalline system with attributes:

data: polars DataFrame with columns [‘element’, ‘x’, ‘y’, ‘z’, ‘grain_id’, ‘type’]
box: simulation box

Examples

>>> poly = CreatePolycrystal(unitcell, box=100, seed_number=10)
>>> system = poly.compute(verbose=True)
>>> system.write_xyz("output.xyz")

mdapy.data module¶

This module provides fundamental atomic data used throughout MDAPY, including chemical symbols, atomic numbers, van der Waals radii, and atomic masses (copy from ase.data). All data are stored as immutable NumPy arrays or Python dictionaries for fast lookup and consistent use across simulations and analysis.

The data originate from reliable references such as IUPAC atomic weight tables and peer-reviewed compilations of van der Waals radii.

References

S. Alvarez, A cartography of the van der Waals territories, Dalton Trans., 2013, 42, 8617–8636. DOI: 10.1039/C3DT50599E
J. Meija et al., Atomic weights of the elements 2013 (IUPAC Technical Report), Pure and Applied Chemistry, 88(3), 265–291 (2016). DOI: 10.1515/pac-2015-0305

mdapy.data.chemical_symbols¶

List of all chemical element symbols in periodic table order. The first entry “X” is a placeholder for undefined species.

Type:: list[str]

mdapy.data.atomic_numbers¶

A mapping from element symbol to its atomic number, e.g. atomic_numbers["Fe"] == 26.

Type:: dict[str, int]

mdapy.data.vdw_radii¶

NumPy array of van der Waals radii (in Å). Van der Waals radii in [A] taken from: A cartography of the van der Waals territories S. Alvarez, Dalton Trans., 2013, 42, 8617-8636 DOI: 10.1039/C3DT50599E. Some elements have np.nan where data are unavailable. The array is read-only.

Type:: np.ndarray

mdapy.data.atomic_masses¶

NumPy array of standard atomic weights (in unified atomic mass units, u). Values follow the IUPAC 2013 report (Meija et al., 2016). For elements without stable isotopes, the mass of the most stable isotope is used. The array is read-only.

Type:: np.ndarray

Examples

>>> from mdapy import data
>>> data.atomic_numbers["C"]
6
>>> data.chemical_symbols[8]
'O'
>>> data.vdw_radii[8]
1.5
>>> data.atomic_masses[26]
55.845

mdapy.data.xray_form_factor: Dict[str, List[float]] = {'Ac': [35.6597, 0.589092, 23.1032, 3.65155, 12.5977, 18.599, 4.08655, 117.02, 13.5266], 'Ac3+': [35.1736, 0.579689, 22.1112, 3.41437, 8.19216, 12.9187, 7.05545, 25.9443, 13.4637], 'Ag': [19.2808, 0.6446, 16.6885, 7.4726, 4.8045, 24.6605, 1.0463, 99.8156, 5.179], 'Ag1+': [19.1812, 0.646179, 15.9719, 7.19123, 5.27475, 21.7326, 0.357534, 66.1147, 5.21572], 'Ag2+': [19.1643, 0.645643, 16.2456, 7.18544, 4.3709, 21.4072, 0.0, 0.0, 5.21404], 'Al': [6.4202, 3.0387, 1.9002, 0.7426, 1.5936, 31.5472, 1.9646, 85.0886, 1.1151], 'Al3+': [4.17448, 1.93816, 3.3876, 4.14553, 1.20296, 0.228753, 0.528137, 8.28524, 0.706786], 'Am': [36.6706, 0.483629, 24.0992, 3.20647, 17.3415, 14.3136, 3.49331, 102.273, 13.3592], 'Ar': [7.4845, 0.9072, 6.7723, 14.8407, 0.6539, 43.8983, 1.6442, 33.3929, 1.4445], 'As': [16.6723, 2.6345, 6.0701, 0.2647, 3.4313, 12.9479, 4.2779, 47.7972, 2.531], 'At': [35.3163, 0.68587, 19.0211, 3.97458, 9.49887, 11.3824, 7.42518, 45.4715, 13.7108], 'Au': [16.8819, 0.4611, 18.5913, 8.6216, 25.5582, 1.4826, 5.86, 36.3956, 12.0658], 'Au1+': [28.0109, 1.35321, 17.8204, 7.7395, 14.3359, 0.356752, 6.58077, 26.4043, 11.2299], 'Au3+': [30.6886, 1.2199, 16.9029, 6.82872, 12.7801, 0.212867, 6.52354, 18.659, 9.0968], 'B': [2.0545, 23.2185, 1.3326, 1.021, 1.0979, 60.3498, 0.7068, 0.1403, -0.1932], 'Ba': [20.3361, 3.216, 19.297, 0.2756, 10.888, 20.2073, 2.6959, 167.202, 2.7731], 'Ba2+': [20.1807, 3.21367, 19.1136, 0.28331, 10.9054, 20.0558, 0.77634, 51.746, 3.02902], 'Be': [1.5919, 43.6427, 1.1278, 1.8623, 0.5391, 103.483, 0.7029, 0.542, 0.0385], 'Be2+': [6.2603, 0.0027, 0.8849, 0.8313, 0.7993, 2.2758, 0.1647, 5.1146, -6.1092], 'Bi': [33.3689, 0.704, 12.951, 2.9238, 16.5877, 8.7937, 6.4692, 48.0093, 13.5782], 'Bi3+': [21.8053, 1.2356, 19.5026, 6.24149, 19.1053, 0.469999, 7.10295, 20.3185, 12.4711], 'Bi5+': [33.5364, 0.91654, 25.0946, 0.39042, 19.2497, 5.71414, 6.91555, 12.8285, -6.7994], 'Bk': [36.7881, 0.451018, 24.7736, 3.04619, 17.8919, 12.8946, 4.23284, 86.003, 13.2754], 'Br': [17.1789, 2.1723, 5.2358, 16.5796, 5.6377, 0.2609, 3.9851, 41.4328, 2.9557], 'Br1-': [17.1718, 2.2059, 6.3338, 19.3345, 5.5754, 0.2871, 3.7272, 58.1535, 3.1776], 'C': [2.31, 20.8439, 1.02, 10.2075, 1.5886, 0.5687, 0.865, 51.6512, 0.2156], 'Ca': [8.6266, 10.4421, 7.3873, 0.6599, 1.5899, 85.7484, 1.0211, 178.437, 1.3751], 'Ca2+': [15.6348, -0.0074, 7.9518, 0.6089, 8.4372, 10.3116, 0.8537, 25.9905, -14.875], 'Cd': [19.2214, 0.5946, 17.6444, 6.9089, 4.461, 24.7008, 1.6029, 87.4825, 5.0694], 'Cd2+': [19.1514, 0.597922, 17.2535, 6.80639, 4.47128, 20.2521, 0.0, 0.0, 5.11937], 'Ce': [21.1671, 2.81219, 19.7695, 0.226836, 11.8513, 17.6083, 3.33049, 127.113, 1.86264], 'Ce3+': [20.8036, 2.77691, 19.559, 0.23154, 11.9369, 16.5408, 0.612376, 43.1692, 2.09013], 'Ce4+': [20.3235, 2.65941, 19.8186, 0.21885, 12.1233, 15.7992, 0.144583, 62.2355, 1.5918], 'Cf': [36.9185, 0.437533, 25.1995, 3.00775, 18.3317, 12.4044, 4.24391, 83.7881, 13.2674], 'Cl': [11.4604, 0.0104, 7.1964, 1.1662, 6.2556, 18.5194, 1.6455, 47.7784, -9.5574], 'Cl1-': [18.2915, 0.0066, 7.2084, 1.1717, 6.5337, 19.5424, 2.3386, 60.4486, -16.378], 'Cm': [36.6488, 0.465154, 24.4096, 3.08997, 17.399, 13.4346, 4.21665, 88.4834, 13.2887], 'Co': [12.2841, 4.2791, 7.3409, 0.2784, 4.0034, 13.5359, 2.3488, 71.1692, 1.0118], 'Co2+': [11.2296, 4.1231, 7.3883, 0.2726, 4.7393, 10.2443, 0.7108, 25.6466, 0.9324], 'Co3+': [10.338, 3.90969, 7.88173, 0.238668, 4.76795, 8.35583, 0.725591, 18.3491, 0.286667], 'Cr': [10.6406, 6.1038, 7.3537, 0.392, 3.324, 20.2626, 1.4922, 98.7399, 1.1832], 'Cr2+': [9.54034, 5.66078, 7.7509, 0.344261, 3.58274, 13.3075, 0.509107, 32.4224, 0.616898], 'Cr3+': [9.6809, 5.59463, 7.81136, 0.334393, 2.87603, 12.8288, 0.113575, 32.8761, 0.518275], 'Cs': [20.3892, 3.569, 19.1062, 0.3107, 10.662, 24.3879, 1.4953, 213.904, 3.3352], 'Cs1+': [20.3524, 3.552, 19.1278, 0.3086, 10.2821, 23.7128, 0.9615, 59.4565, 3.2791], 'Cu': [13.338, 3.5828, 7.1676, 0.247, 5.6158, 11.3966, 1.6735, 64.8126, 1.191], 'Cu1+': [11.9475, 3.3669, 7.3573, 0.2274, 6.2455, 8.6625, 1.5578, 25.8487, 0.89], 'Cu2+': [11.8168, 3.37484, 7.11181, 0.244078, 5.78135, 7.9876, 1.14523, 19.897, 1.14431], 'Cval': [2.26069, 22.6907, 1.56165, 0.656665, 1.05075, 9.75618, 0.839259, 55.5949, 0.286977], 'Dy': [26.507, 2.1802, 17.6383, 0.202172, 14.5596, 12.1899, 2.96577, 111.874, 4.29728], 'Dy3+': [25.5395, 1.9804, 20.2861, 0.143384, 11.9812, 9.34972, 4.50073, 19.581, 0.68969], 'Er': [27.6563, 2.07356, 16.4285, 0.223545, 14.9779, 11.3604, 2.98233, 105.703, 5.92046], 'Er3+': [26.722, 1.84659, 19.7748, 0.13729, 12.1506, 8.36225, 5.17379, 17.8974, 1.17613], 'Eu': [24.6274, 2.3879, 19.0886, 0.1942, 13.7603, 13.7546, 2.9227, 123.174, 2.5745], 'Eu2+': [24.0063, 2.27783, 19.9504, 0.17353, 11.8034, 11.6096, 3.87243, 26.5156, 1.36389], 'Eu3+': [23.7497, 2.22258, 20.3745, 0.16394, 11.8509, 11.311, 3.26503, 22.9966, 0.759344], 'F': [3.5392, 10.2825, 2.6412, 4.2944, 1.517, 0.2615, 1.0243, 26.1476, 0.2776], 'F1-': [3.6322, 5.27756, 3.51057, 14.7353, 1.26064, 0.442258, 0.940706, 47.3437, 0.653396], 'Fe': [11.7695, 4.7611, 7.3573, 0.3072, 3.5222, 15.3535, 2.3045, 76.8805, 1.0369], 'Fe2+': [11.0424, 4.6538, 7.374, 0.3053, 4.1346, 12.0546, 0.4399, 31.2809, 1.0097], 'Fe3+': [11.1764, 4.6147, 7.3863, 0.3005, 3.3948, 11.6729, 0.0724, 38.5566, 0.9707], 'Fr': [35.9299, 0.646453, 23.0547, 4.17619, 12.1439, 23.1052, 2.11253, 150.645, 13.7247], 'Ga': [15.2354, 3.0669, 6.7006, 0.2412, 4.3591, 10.7805, 2.9623, 61.4135, 1.7189], 'Ga3+': [12.692, 2.81262, 6.69883, 0.22789, 6.06692, 6.36441, 1.0066, 14.4122, 1.53545], 'Gd': [25.0709, 2.25341, 19.0798, 0.181951, 13.8518, 12.9331, 3.54545, 101.398, 2.4196], 'Gd3+': [24.3466, 2.13553, 20.4208, 0.155525, 11.8708, 10.5782, 3.7149, 21.7029, 0.645089], 'Ge': [16.0816, 2.8509, 6.3747, 0.2516, 3.7068, 11.4468, 3.683, 54.7625, 2.1313], 'Ge4+': [12.9172, 2.53718, 6.70003, 0.205855, 6.06791, 5.47913, 0.859041, 11.603, 1.45572], 'H': [0.489918, 20.6593, 0.262003, 7.74039, 0.196767, 49.5519, 0.049879, 2.20159, 0.001305], 'H1-': [0.897661, 53.1368, 0.565616, 15.187, 0.415815, 186.576, 0.116973, 3.56709, 0.002389], 'He': [0.8734, 9.1037, 0.6309, 3.3568, 0.3112, 22.9276, 0.178, 0.9821, 0.0064], 'Hf': [29.144, 1.83262, 15.1726, 9.5999, 14.7586, 0.275116, 4.30013, 72.029, 8.58154], 'Hf4+': [28.8131, 1.59136, 18.4601, 0.128903, 12.7285, 6.76232, 5.59927, 14.0366, 2.39699], 'Hg': [20.6809, 0.545, 19.0417, 8.4484, 21.6575, 1.5729, 5.9676, 38.3246, 12.6089], 'Hg1+': [25.0853, 1.39507, 18.4973, 7.65105, 16.8883, 0.443378, 6.48216, 28.2262, 12.0205], 'Hg2+': [29.5641, 1.21152, 18.06, 7.05639, 12.8374, 0.284738, 6.89912, 20.7482, 10.6268], 'Ho': [26.9049, 2.07051, 17.294, 0.19794, 14.5583, 11.4407, 3.63837, 92.6566, 4.56796], 'Ho3+': [26.1296, 1.91072, 20.0994, 0.139358, 11.9788, 8.80018, 4.93676, 18.5908, 0.852795], 'I': [20.1472, 4.347, 18.9949, 0.3814, 7.5138, 27.766, 2.2735, 66.8776, 4.0712], 'I1-': [20.2332, 4.3579, 18.997, 0.3815, 7.8069, 29.5259, 2.8868, 84.9304, 4.0714], 'In': [19.1624, 0.5476, 18.5596, 6.3776, 4.2948, 25.8499, 2.0396, 92.8029, 4.9391], 'In3+': [19.1045, 0.551522, 18.1108, 6.3247, 3.78897, 17.3595, 0.0, 0.0, 4.99635], 'Ir': [27.3049, 1.59279, 16.7296, 8.86553, 15.6115, 0.417916, 5.83377, 45.0011, 11.4722], 'Ir3+': [30.4156, 1.34323, 15.862, 7.10909, 13.6145, 0.204633, 5.82008, 20.3254, 8.27903], 'Ir4+': [30.7058, 1.30923, 15.5512, 6.71983, 14.2326, 0.167252, 5.53672, 17.4911, 6.96824], 'K': [8.2186, 12.7949, 7.4398, 0.7748, 1.0519, 213.187, 0.8659, 41.6841, 1.4228], 'K1+': [7.9578, 12.6331, 7.4917, 0.7674, 6.359, -0.002, 1.1915, 31.9128, -4.9978], 'Kr': [17.3555, 1.9384, 6.7286, 16.5623, 5.5493, 0.2261, 3.5375, 39.3972, 2.825], 'La': [20.578, 2.94817, 19.599, 0.244475, 11.3727, 18.7726, 3.28719, 133.124, 2.14678], 'La3+': [20.2489, 2.9207, 19.3763, 0.250698, 11.6323, 17.8211, 0.336048, 54.9453, 2.4086], 'Li': [1.1282, 3.9546, 0.7508, 1.0524, 0.6175, 85.3905, 0.4653, 168.261, 0.0377], 'Li1+': [0.6968, 4.6237, 0.7888, 1.9557, 0.3414, 0.6316, 0.1563, 10.0953, 0.0167], 'Lu': [28.9476, 1.90182, 15.2208, 9.98519, 15.1, 0.261033, 3.71601, 84.3298, 7.97628], 'Lu3+': [28.4628, 1.68216, 18.121, 0.142292, 12.8429, 7.33727, 5.59415, 16.3535, 2.97573], 'Mg': [5.4204, 2.8275, 2.1735, 79.2611, 1.2269, 0.3808, 2.3073, 7.1937, 0.8584], 'Mg2+': [3.4988, 2.1676, 3.8378, 4.7542, 1.3284, 0.185, 0.8497, 10.1411, 0.4853], 'Mn': [11.2819, 5.3409, 7.3573, 0.3432, 3.0193, 17.8674, 2.2441, 83.7543, 1.0896], 'Mn2+': [10.8061, 5.2796, 7.362, 0.3435, 3.5268, 14.343, 0.2184, 41.3235, 1.0874], 'Mn3+': [9.84521, 4.91797, 7.87194, 0.294393, 3.56531, 10.8171, 0.323613, 24.1281, 0.393974], 'Mn4+': [9.96253, 4.8485, 7.97057, 0.283303, 2.76067, 10.4852, 0.054447, 27.573, 0.251877], 'Mo': [3.7025, 0.2772, 17.2356, 1.0958, 12.8876, 11.004, 3.7429, 61.6584, 4.3875], 'Mo3+': [21.1664, 0.014734, 18.2017, 1.03031, 11.7423, 9.53659, 2.30951, 26.6307, -14.421], 'Mo5+': [21.0149, 0.014345, 18.0992, 1.02238, 11.4632, 8.78809, 0.740625, 23.3452, -14.316], 'Mo6+': [17.8871, 1.03649, 11.175, 8.48061, 6.57891, 0.058881, 0.0, 0.0, 0.344941], 'N': [12.2126, 0.0057, 3.1322, 9.8933, 2.0125, 28.9975, 1.1663, 0.5826, -11.529], 'Na': [4.7626, 3.285, 3.1736, 8.8422, 1.2674, 0.3136, 1.1128, 129.424, 0.676], 'Na1+': [3.2565, 2.6671, 3.9362, 6.1153, 1.3998, 0.2001, 1.0032, 14.039, 0.404], 'Nb': [17.6142, 1.18865, 12.0144, 11.766, 4.04183, 0.204785, 3.53346, 69.7957, 3.75591], 'Nb3+': [19.8812, 0.019175, 18.0653, 1.13305, 11.0177, 10.1621, 1.94715, 28.3389, -12.912], 'Nb5+': [17.9163, 1.12446, 13.3417, 0.028781, 10.799, 9.28206, 0.337905, 25.7228, -6.3934], 'Nd': [22.6845, 2.66248, 19.6847, 0.210628, 12.774, 15.885, 2.85137, 137.903, 1.98486], 'Nd3+': [21.961, 2.52722, 19.9339, 0.199237, 12.12, 14.1783, 1.51031, 30.8717, 1.47588], 'Ne': [3.9553, 8.4042, 3.1125, 3.4262, 1.4546, 0.2306, 1.1251, 21.7184, 0.3515], 'Ni': [12.8376, 3.8785, 7.292, 0.2565, 4.4438, 12.1763, 2.38, 66.3421, 1.0341], 'Ni2+': [11.4166, 3.6766, 7.4005, 0.2449, 5.3442, 8.873, 0.9773, 22.1626, 0.8614], 'Ni3+': [10.7806, 3.5477, 7.75868, 0.22314, 5.22746, 7.64468, 0.847114, 16.9673, 0.386044], 'Np': [36.1874, 0.511929, 23.5964, 3.25396, 15.6402, 15.3622, 4.1855, 97.4908, 13.3573], 'Np3+': [35.7074, 0.502322, 22.613, 3.03807, 12.9898, 12.1449, 5.43227, 25.4928, 13.2544], 'Np4+': [35.5103, 0.498626, 22.5787, 2.96627, 12.7766, 11.9484, 4.92159, 22.7502, 13.2116], 'Np6+': [35.0136, 0.48981, 22.7286, 2.81099, 14.3884, 12.33, 1.75669, 22.6581, 13.113], 'O': [3.0485, 13.2771, 2.2868, 5.7011, 1.5463, 0.3239, 0.867, 32.9089, 0.2508], 'O1-': [4.1916, 12.8573, 1.63969, 4.17236, 1.52673, 47.0179, -20.307, -0.01404, 21.9412], 'Os': [28.1894, 1.62903, 16.155, 8.97948, 14.9305, 0.382661, 5.67589, 48.1647, 11.0005], 'Os4+': [30.419, 1.37113, 15.2637, 6.84706, 14.7458, 0.165191, 5.06795, 18.003, 6.49804], 'P': [6.4345, 1.9067, 4.1791, 27.157, 1.78, 0.526, 1.4908, 68.1645, 1.1149], 'Pa': [35.8847, 0.547751, 23.2948, 3.41519, 14.1891, 16.9235, 4.17287, 105.251, 13.4287], 'Pb': [31.0617, 0.6902, 13.0637, 2.3576, 18.442, 8.618, 5.9696, 47.2579, 13.4118], 'Pb2+': [21.7886, 1.3366, 19.5682, 0.488383, 19.1406, 6.7727, 7.01107, 23.8132, 12.4734], 'Pb4+': [32.1244, 1.00566, 18.8003, 6.10926, 12.0175, 0.147041, 6.96886, 14.714, 8.08428], 'Pd': [19.3319, 0.698655, 15.5017, 7.98929, 5.29537, 25.2052, 0.605844, 76.8986, 5.26593], 'Pd2+': [19.1701, 0.696219, 15.2096, 7.55573, 4.32234, 22.5057, 0.0, 0.0, 5.2916], 'Pd4+': [19.2493, 0.683839, 14.79, 7.14833, 2.89289, 17.9144, -7.9492, 0.005127, 13.0174], 'Pm': [23.3405, 2.5627, 19.6095, 0.202088, 13.1235, 15.1009, 2.87516, 132.721, 2.02876], 'Pm3+': [22.5527, 2.4174, 20.1108, 0.185769, 12.0671, 13.1275, 2.07492, 27.4491, 1.19499], 'Po': [34.6726, 0.700999, 15.4733, 3.55078, 13.1138, 9.55642, 7.02588, 47.0045, 13.677], 'Pr': [22.044, 2.77393, 19.6697, 0.222087, 12.3856, 16.7669, 2.82428, 143.644, 2.0583], 'Pr3+': [21.3727, 2.6452, 19.7491, 0.214299, 12.1329, 15.323, 0.97518, 36.4065, 1.77132], 'Pr4+': [20.9413, 2.54467, 20.0539, 0.202481, 12.4668, 14.8137, 0.296689, 45.4643, 1.24285], 'Pt': [27.0059, 1.51293, 17.7639, 8.81174, 15.7131, 0.424593, 5.7837, 38.6103, 11.6883], 'Pt2+': [29.8429, 1.32927, 16.7224, 7.38979, 13.2153, 0.263297, 6.35234, 22.9426, 9.85329], 'Pt4+': [30.9612, 1.24813, 15.9829, 6.60834, 13.7348, 0.16864, 5.92034, 16.9392, 7.39534], 'Pu': [36.5254, 0.499384, 23.8083, 3.26371, 16.7707, 14.9455, 3.47947, 105.98, 13.3812], 'Pu3+': [35.84, 0.484938, 22.7169, 2.96118, 13.5807, 11.5331, 5.66016, 24.3992, 13.1991], 'Pu4+': [35.6493, 0.481422, 22.646, 2.8902, 13.3595, 11.316, 5.18831, 21.8301, 13.1555], 'Pu6+': [35.1736, 0.473204, 22.7181, 2.73848, 14.7635, 11.553, 2.28678, 20.9303, 13.0582], 'Ra': [35.763, 0.616341, 22.9064, 3.87135, 12.4739, 19.9887, 3.21097, 142.325, 13.6211], 'Ra2+': [35.215, 0.604909, 21.67, 3.5767, 7.91342, 12.601, 7.65078, 29.8436, 13.5431], 'Rb': [17.1784, 1.7888, 9.6435, 17.3151, 5.1399, 0.2748, 1.5292, 164.934, 3.4873], 'Rb1+': [17.5816, 1.7139, 7.6598, 14.7957, 5.8981, 0.1603, 2.7817, 31.2087, 2.0782], 'Re': [28.7621, 1.67191, 15.7189, 9.09227, 14.5564, 0.3505, 5.44174, 52.0861, 10.472], 'Rh': [19.2957, 0.751536, 14.3501, 8.21758, 4.73425, 25.8749, 1.28918, 98.6062, 5.328], 'Rh3+': [18.8785, 0.764252, 14.1259, 7.84438, 3.32515, 21.2487, -6.1989, -0.01036, 11.8678], 'Rh4+': [18.8545, 0.760825, 13.9806, 7.62436, 2.53464, 19.3317, -5.6526, -0.0102, 11.2835], 'Rn': [35.5631, 0.6631, 21.2816, 4.0691, 8.0037, 14.0422, 7.4433, 44.2473, 13.6905], 'Ru': [19.2674, 0.80852, 12.9182, 8.43467, 4.86337, 24.7997, 1.56756, 94.2928, 5.37874], 'Ru3+': [18.5638, 0.847329, 13.2885, 8.37164, 9.32602, 0.017662, 3.00964, 22.887, -3.1892], 'Ru4+': [18.5003, 0.844582, 13.1787, 8.12534, 4.71304, 0.36495, 2.18535, 20.8504, 1.42357], 'S': [6.9053, 1.4679, 5.2034, 22.2151, 1.4379, 0.2536, 1.5863, 56.172, 0.8669], 'Sb': [19.6418, 5.3034, 19.0455, 0.4607, 5.0371, 27.9074, 2.6827, 75.2825, 4.5909], 'Sb3+': [18.9755, 0.467196, 18.933, 5.22126, 5.10789, 19.5902, 0.288753, 55.5113, 4.69626], 'Sb5+': [19.8685, 5.44853, 19.0302, 0.467973, 2.41253, 14.1259, 0.0, 0.0, 4.69263], 'Sc': [9.189, 9.0213, 7.3679, 0.5729, 1.6409, 136.108, 1.468, 51.3531, 1.3329], 'Sc3+': [13.4008, 0.29854, 8.0273, 7.9629, 1.65943, -0.28604, 1.57936, 16.0662, -6.6667], 'Se': [17.0006, 2.4098, 5.8196, 0.2726, 3.9731, 15.2372, 4.3543, 43.8163, 2.8409], 'Si': [6.2915, 2.4386, 3.0353, 32.3337, 1.9891, 0.6785, 1.541, 81.6937, 1.1407], 'Si4+': [4.43918, 1.64167, 3.20345, 3.43757, 1.19453, 0.2149, 0.41653, 6.65365, 0.746297], 'Sival': [5.66269, 2.6652, 3.07164, 38.6634, 2.62446, 0.916946, 1.3932, 93.5458, 1.24707], 'Sm': [24.0042, 2.47274, 19.4258, 0.196451, 13.4396, 14.3996, 2.89604, 128.007, 2.20963], 'Sm3+': [23.1504, 2.31641, 20.2599, 0.174081, 11.9202, 12.1571, 2.71488, 24.8242, 0.954586], 'Sn': [19.1889, 5.8303, 19.1005, 0.5031, 4.4585, 26.8909, 2.4663, 83.9571, 4.7821], 'Sn2+': [19.1094, 0.5036, 19.0548, 5.8378, 4.5648, 23.3752, 0.487, 62.2061, 4.7861], 'Sn4+': [18.9333, 5.764, 19.7131, 0.4655, 3.4182, 14.0049, 0.0193, -0.7583, 3.9182], 'Sr': [17.5663, 1.5564, 9.8184, 14.0988, 5.422, 0.1664, 2.6694, 132.376, 2.5064], 'Sr2+': [18.0874, 1.4907, 8.1373, 12.6963, 2.5654, 24.5651, -34.193, -0.0138, 41.4025], 'Ta': [29.2024, 1.77333, 15.2293, 9.37046, 14.5135, 0.295977, 4.76492, 63.3644, 9.24354], 'Ta5+': [29.1587, 1.50711, 18.8407, 0.116741, 12.8268, 6.31524, 5.38695, 12.4244, 1.78555], 'Tb': [25.8976, 2.24256, 18.2185, 0.196143, 14.3167, 12.6648, 2.95354, 115.362, 3.58324], 'Tb3+': [24.9559, 2.05601, 20.3271, 0.149525, 12.2471, 10.0499, 3.773, 21.2773, 0.691967], 'Tc': [19.1301, 0.864132, 11.0948, 8.14487, 4.64901, 21.5707, 2.71263, 86.8472, 5.40428], 'Te': [19.9644, 4.81742, 19.0138, 0.420885, 6.14487, 28.5284, 2.5239, 70.8403, 4.352], 'Th': [35.5645, 0.563359, 23.4219, 3.46204, 12.7473, 17.8309, 4.80703, 99.1722, 13.4314], 'Th4+': [35.1007, 0.555054, 22.4418, 3.24498, 9.78554, 13.4661, 5.29444, 23.9533, 13.376], 'Ti': [9.7595, 7.8508, 7.3558, 0.5, 1.6991, 35.6338, 1.9021, 116.105, 1.2807], 'Ti2+': [9.11423, 7.5243, 7.62174, 0.457585, 2.2793, 19.5361, 0.087899, 61.6558, 0.897155], 'Ti3+': [17.7344, 0.22061, 8.73816, 7.04716, 5.25691, -0.15762, 1.92134, 15.9768, -14.652], 'Ti4+': [19.5114, 0.178847, 8.23473, 6.67018, 2.01341, -0.29263, 1.5208, 12.9464, -13.28], 'Tl': [27.5446, 0.65515, 19.1584, 8.70751, 15.538, 1.96347, 5.52593, 45.8149, 13.1746], 'Tl1+': [21.3985, 1.4711, 20.4723, 0.517394, 18.7478, 7.43463, 6.82847, 28.8482, 12.5258], 'Tl3+': [30.8695, 1.1008, 18.3481, 6.53852, 11.9328, 0.219074, 7.00574, 17.2114, 9.8027], 'Tm': [28.1819, 2.02859, 15.8851, 0.238849, 15.1542, 10.9975, 2.98706, 102.961, 6.75621], 'Tm3+': [27.3083, 1.78711, 19.332, 0.136974, 12.3339, 7.96778, 5.38348, 17.2922, 1.63929], 'U': [36.0228, 0.5293, 23.4128, 3.3253, 14.9491, 16.0927, 4.188, 100.613, 13.3966], 'U3+': [35.5747, 0.52048, 22.5259, 3.12293, 12.2165, 12.7148, 5.37073, 26.3394, 13.3092], 'U4+': [35.3715, 0.516598, 22.5326, 3.05053, 12.0291, 12.5723, 4.7984, 23.4582, 13.2671], 'U6+': [34.8509, 0.507079, 22.7584, 2.8903, 14.0099, 13.1767, 1.21457, 25.2017, 13.1665], 'V': [10.2971, 6.8657, 7.3511, 0.4385, 2.0703, 26.8938, 2.0571, 102.478, 1.2199], 'V2+': [10.106, 6.8818, 7.3541, 0.4409, 2.2884, 20.3004, 0.0223, 115.122, 1.2298], 'V3+': [9.43141, 6.39535, 7.7419, 0.383349, 2.15343, 15.1908, 0.016865, 63.969, 0.656565], 'V5+': [15.6887, 0.679003, 8.14208, 5.40135, 2.03081, 9.97278, -9.576, 0.940464, 1.7143], 'W': [29.0818, 1.72029, 15.43, 9.2259, 14.4327, 0.321703, 5.11982, 57.056, 9.8875], 'W6+': [29.4936, 1.42755, 19.3763, 0.104621, 13.0544, 5.93667, 5.06412, 11.1972, 1.01074], 'Xe': [20.2933, 3.9282, 19.0298, 0.344, 8.9767, 26.4659, 1.99, 64.2658, 3.7118], 'Y': [17.776, 1.4029, 10.2946, 12.8006, 5.72629, 0.125599, 3.26588, 104.354, 1.91213], 'Y3+': [17.9268, 1.35417, 9.1531, 11.2145, 1.76795, 22.6599, -33.108, -0.01319, 40.2602], 'Yb': [28.6641, 1.9889, 15.4345, 0.257119, 15.3087, 10.6647, 2.98963, 100.417, 7.56672], 'Yb2+': [28.1209, 1.78503, 17.6817, 0.15997, 13.3335, 8.18304, 5.14657, 20.39, 3.70983], 'Yb3+': [27.8917, 1.73272, 18.7614, 0.13879, 12.6072, 7.64412, 5.47647, 16.8153, 2.26001], 'Zn': [14.0743, 3.2655, 7.0318, 0.2333, 5.1652, 10.3163, 2.41, 58.7097, 1.3041], 'Zn2+': [11.9719, 2.9946, 7.3862, 0.2031, 6.4668, 7.0826, 1.394, 18.0995, 0.7807], 'Zr': [17.8765, 1.27618, 10.948, 11.916, 5.41732, 0.117622, 3.65721, 87.6627, 2.06929], 'Zr4+': [18.1668, 1.2148, 10.0562, 10.1483, 1.01118, 21.6054, -2.6479, -0.10276, 9.41454]}¶: X-ray form factors (Cromer-Mann coefficients) from International Tables for Crystallography, Volume C (Brown et al., 2004), Table 6.1.1.4. Keys are element/ion symbols (str), values are lists of 9 float coefficients [a1, b1, a2, b2, a3, b3, a4, b4, c].

mdapy.data.neutron_form_factor: Dict[str, float | complex] = {'100Mo': 6.73, '102Pd': 7.77, '104Pd': 7.77, '105Pd': 5.5, '106Cd': 5.2, '106Pd': 6.4, '107Ag': 7.555, '108Cd': 5.4, '108Pd': 4.1, '109Ag': 4.165, '10B': (-0.1-1.066j), '110Cd': 5.9, '110Pd': 7.77, '111Cd': 6.5, '112Cd': 6.4, '112Sn': 6.1, '113Cd': (-8-5.73j), '113In': 5.39, '114Cd': 7.5, '114Sn': 6.2, '115In': (4.01-0.0562j), '115Sn': 6.1, '116Cd': 6.3, '116Sn': 5.93, '117Sn': 6.48, '118Sn': 6.07, '119Sn': 6.12, '11B': 6.65, '120Sn': 6.49, '120Te': 5.3, '121Sb': 5.71, '122Sn': 5.74, '122Te': 3.8, '123Sb': 5.38, '123Te': (-0.05-0.116j), '124Sn': 5.97, '124Te': 7.96, '125Te': 5.02, '126Te': 5.56, '128Te': 5.89, '12C': 6.6511, '130Ba': -3.6, '130Te': 6.02, '132Ba': 7.8, '134Ba': 5.7, '135Ba': 4.67, '136Ba': 4.91, '136Ce': 5.8, '137Ba': 6.83, '138Ba': 4.84, '138Ce': 6.7, '138La': 8.2, '139La': 8.24, '13C': 6.19, '140Ce': 4.84, '142Ce': 4.75, '142Nd': 7.7, '143Nd': 14.2, '144Nd': 2.8, '144Sm': -3.4, '145Nd': 14.2, '146Nd': 8.7, '147Sm': 14.3, '148Nd': 5.7, '148Sm': -3.4, '149Sm': (-19.2-11.7j), '14N': 9.37, '150Nd': 5.3, '150Sm': 14.3, '151Eu': (6.13-2.53j), '152Gd': 10.3, '152Sm': -5.0, '153Eu': 8.22, '154Gd': 10.3, '154Sm': 9.3, '155Gd': (6-17j), '156Dy': 6.1, '156Gd': 6.3, '157Gd': (-1.14-71.9j), '158Dy': 6.4, '158Gd': 9.2, '15N': 6.44, '160Dy': 6.7, '160Gd': 9.15, '161Dy': 10.3, '162Dy': -1.4, '162Er': 8.8, '163Dy': 5.0, '164Dy': (49.4-0.79j), '164Er': 8.2, '166Er': 10.6, '167Er': 3.0, '168Er': 7.4, '168Yb': (-4.07-0.62j), '16O': 5.803, '170Er': 9.6, '170Yb': 6.77, '171Yb': 9.66, '172Yb': 9.43, '173Yb': 9.56, '174Hf': 10.911, '174Yb': 19.3, '175Lu': 7.24, '176Hf': 6.61, '176Lu': (6.1-0.57j), '176Yb': 8.72, '177Hf': 0.81, '178Hf': 5.9, '179Hf': 7.46, '17O': 5.78, '180Hf': 13.2, '180Ta': 7.2, '180W': 5.3, '181Ta': 6.91, '182W': 6.97, '183W': 6.53, '184Os': 10.2, '184W': 7.48, '185Re': 9.0, '186Os': 11.617, '186W': -0.72, '187Os': 10.2, '187Re': 9.3, '188Os': 7.6, '189Os': 10.7, '18O': 5.84, '190Os': 11.0, '190Pt': 9.0, '192Os': 11.5, '192Pt': 9.9, '194Pt': 10.55, '195Pt': 8.83, '196Hg': 30.31, '196Pt': 9.89, '198Pt': 7.8, '199Hg': 16.9, '1H': -3.7406, '203Tl': 6.99, '204Pb': 9.9, '205Tl': 9.52, '206Pb': 9.22, '207Pb': 9.28, '208Pb': 9.5, '20Ne': 4.631, '21Ne': 6.66, '22Ne': 3.87, '233U': 10.1, '234U': 12.4, '235U': 10.47, '238Pu': 14.1, '238U': 8.402, '239Pu': 7.7, '240Pu': 3.5, '242Pu': 8.1, '244Cm': 9.5, '246Cm': 9.3, '248Cm': 7.7, '24Mg': 5.66, '25Mg': 3.62, '26Mg': 4.89, '28Si': 4.107, '29Si': 4.7, '2H': 6.671, '30Si': 4.58, '32S': 2.804, '33S': 4.74, '34S': 3.48, '35Cl': 11.65, '36Ar': 24.9, '36S': 3.1, '37Cl': 3.08, '38Ar': 3.5, '39K': 3.74, '3H': 4.792, '3He': (5.74-1.483j), '40Ar': 1.83, '40Ca': 4.8, '40K': 3.1, '41K': 2.69, '42Ca': 3.36, '43Ca': -1.56, '44Ca': 1.42, '46Ca': 3.6, '46Ti': 4.93, '47Ti': 3.63, '48Ca': 0.39, '48Ti': -6.08, '49Ti': 1.04, '4He': 3.26, '50Cr': -4.5, '50Ti': 6.18, '50V': 7.6, '51V': -0.402, '52Cr': 4.92, '53Cr': -4.2, '54Cr': 4.55, '54Fe': 4.2, '56Fe': 9.94, '57Fe': 2.3, '58Fe': 15.7, '58Ni': 14.4, '60Ni': 2.8, '61Ni': 7.6, '62Ni': -8.7, '63Cu': 6.43, '64Ni': -0.37, '64Zn': 5.22, '65Cu': 10.61, '66Zn': 5.97, '67Zn': 7.56, '68Zn': 6.03, '69Ga': 7.88, '6Li': (2-0.261j), '70Ge': 10.0, '70Zn': 6.1, '71Ga': 6.4, '72Ge': 8.51, '73Ge': 5.02, '74Ge': 7.58, '74Se': 0.8, '76Ge': 8.2, '76Se': 12.2, '77Se': 8.25, '78Se': 8.24, '79Br': 6.8, '7Li': -2.22, '80Se': 7.48, '81Br': 6.79, '82Se': 6.34, '84Sr': 7.1, '85Rb': 7.03, '86Kr': 8.1, '86Sr': 5.67, '87Rb': 7.23, '87Sr': 7.4, '88Sr': 7.15, '90Zr': 6.4, '91Zr': 8.7, '92Mo': 6.91, '92Zr': 7.4, '94Mo': 6.8, '94Zr': 8.2, '95Mo': 6.91, '96Mo': 6.2, '96Zr': 5.5, '97Mo': 7.24, '98Mo': 6.58, 'Ag': 5.922, 'Al': 3.449, 'Am': 8.3, 'Ar': 1.909, 'As': 6.58, 'Au': 7.63, 'B': (5.3-0.213j), 'Ba': 5.07, 'Be': 7.79, 'Bi': 8.532, 'Br': 6.795, 'C': 6.646, 'Ca': 4.7, 'Cd': (4.87-0.7j), 'Ce': 4.84, 'Cl': 9.577, 'Co': 2.49, 'Cr': 3.635, 'Cs': 5.42, 'Cu': 7.718, 'Dy': (16.9-0.276j), 'Er': 7.79, 'Eu': (7.22-1.26j), 'F': 5.654, 'Fe': 9.45, 'Ga': 7.288, 'Gd': (6.5-13.82j), 'Ge': 8.185, 'H': -3.739, 'He': 3.263, 'Hf': 7.7, 'Hg': 12.692, 'Ho': 8.01, 'I': 5.28, 'In': (4.065-0.0539j), 'Ir': 10.6, 'K': 3.67, 'Kr': 7.81, 'La': 8.24, 'Li': -1.9, 'Lu': 7.21, 'Mg': 5.375, 'Mn': -3.73, 'Mo': 6.715, 'N': 9.36, 'Na': 3.63, 'Nb': 7.054, 'Nd': 7.69, 'Ne': 4.566, 'Ni': 10.3, 'Np': 10.55, 'O': 5.803, 'Os': 10.7, 'P': 5.13, 'Pa': 9.1, 'Pb': 9.405, 'Pd': 5.91, 'Pm': 12.6, 'Pr': 4.58, 'Pt': 9.6, 'Ra': 10.01, 'Rb': 7.09, 'Re': 9.2, 'Rh': 5.88, 'Ru': 7.03, 'S': 2.847, 'Sb': 5.57, 'Sc': 12.29, 'Se': 7.97, 'Si': 4.1491, 'Sm': (0.8-1.65j), 'Sn': 6.225, 'Sr': 7.02, 'Ta': 6.91, 'Tb': 7.38, 'Tc': 6.8, 'Te': 5.8, 'Th': 10.31, 'Ti': -3.438, 'Tl': 8.776, 'Tm': 7.07, 'U': 8.417, 'V': -0.3824, 'W': 4.86, 'Xe': 4.92, 'Y': 7.75, 'Yb': 12.43, 'Zn': 5.68, 'Zr': 7.16}¶: Coherent neutron scattering lengths (in fm) from NIST, https://www.ncnr.nist.gov/resources/n-lengths/list.html. Keys are element / isotope symbols (str). Values are real numbers for elements without significant absorption and complex numbers (with imaginary part) for absorptive isotopes such as \({}^{10}\mathrm{B}\), Cd, Sm, Eu, Gd. The same value is used at every k since neutron form factors are k-independent.

mdapy.eam module¶

Embedded Atom Method (EAM) Potential Calculator¶

This module implements the Embedded Atom Method (EAM) for static simulations. The EAM potential is widely used for metallic systems and consists of two main components:

Embedding energy: \(F_i(\rho_i)\) - energy to embed atom i into the electron density \(\rho_i\)
Pair interaction: \(\phi_{ij}(r_{ij})\) - pair potential between atoms i and j

The total energy is given by:

\[E_{tot} = \sum_i F_i(\rho_i) + \frac{1}{2} \sum_{i,j \neq i} \phi_{ij}(r_{ij})\]

where the electron density at atom i is:

\[\rho_i = \sum_{j \neq i} \rho_j(r_{ij})\]

The forces on atom i are computed as:

\[\mathbf{F}_i = -\sum_{j \neq i} \left[ \left(\frac{dF_i}{d\rho_i} + \frac{dF_j}{d\rho_j}\right) \frac{d\rho_j}{dr_{ij}} + \frac{d\phi_{ij}}{dr_{ij}} \right] \frac{\mathbf{r}_{ij}}{r_{ij}}\]

The stress tensor is computed from the virial theorem:

\[\sigma_{\alpha\beta} = -\frac{1}{V} \sum_i \sum_{j \neq i} r_{ij,\alpha} F_{ij,\beta}\]

References

Daw, M. S., & Baskes, M. I. (1984). Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. Physical Review B, 29(12), 6443.

Foiles, S. M., Baskes, M. I., & Daw, M. S. (1986). Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys. Physical Review B, 33(12), 7983.

class mdapy.eam.EAM(filename: str, mass_list: List[float] | None = None)[source]¶

Bases: CalculatorMP

Embedded Atom Method (EAM) potential calculator.

This class implements the EAM potential for metallic systems, supporting both single-element and multi-element alloy systems. It reads EAM potential files in the LAMMPS alloy format and provides methods to calculate energies, forces, stresses, and virials.

Parameters:

filename (str) – Path to the EAM potential file in LAMMPS alloy format.
mass_list (List[float], optional) – Atomic masses in the same order as elements_list. When None (default) masses are looked up from the built-in periodic table via the element symbols read from the potential file.

Notes

After parsing, the instance exposes the potential as these attributes:

filename — path to the EAM file.
header — the first 3 lines of the file.
Nelements, elements_list — element count and symbols.
nrho / drho — size and spacing of the electron-density grid.
nr / dr — size and spacing of the distance grid.
rc — cutoff radius, in Å.
F_rho — embedding function F(ρ), shape (Nelements, nrho).
rho_r — electron density ρ(r), shape (Nelements, nr).
phi_r — pair potential φ(r), shape (Nelements, Nelements, nr).
r, rho — grid points, shape (nr,) and (nrho,).
results — dict storing computed energies, forces, stresses, virials.

write_eam_alloy(output_name: str | None = None)[source]¶

Write to an eam.alloy file.

Args:: output_name (Optional[str], optional): output filename, it will automatically generated if not given.

get_energies(data: DataFrame, box: Box) → ndarray[tuple[Any, ...], dtype[float64]][source]¶

Calculate per-atom potential energies.

Parameters:

data (pl.DataFrame) – Atomic configuration containing columns: ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Simulation box object containing boundary conditions and dimensions.

Returns:

NDArray[np.float64] – Per-atom potential energies in eV. Shape: (N,) where N is the number of atoms.

get_energy(data: DataFrame, box: Box) → float[source]¶

Calculate total system potential energy.

Parameters:

data (pl.DataFrame) – Atomic configuration containing columns: ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Simulation box object containing boundary conditions and dimensions.

Returns:

float – Total potential energy in eV.

get_forces(data: DataFrame, box: Box) → ndarray[tuple[Any, ...], dtype[float64]][source]¶

Calculate atomic forces.

Parameters:

data (pl.DataFrame) – Atomic configuration containing columns: ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Simulation box object containing boundary conditions and dimensions.

Returns:

NDArray[np.float64] – Atomic forces in eV/Å. Shape: (N, 3) where N is the number of atoms. Columns represent [Fx, Fy, Fz].

get_stress(data: DataFrame, box: Box) → ndarray[tuple[Any, ...], dtype[float64]][source]¶

Calculate stress tensor in Voigt notation.

Parameters:

data (pl.DataFrame) – Atomic configuration containing columns: ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Simulation box object containing boundary conditions and dimensions.

Returns:

NDArray[np.float64] – Stress tensor in Voigt notation [σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy]. Units: eV/Å³ or GPa depending on unit conversion. Shape: (6,).

get_virials(data: DataFrame, box: Box) → ndarray[tuple[Any, ...], dtype[float64]][source]¶

Calculate per-atom virial tensors.

Parameters:

data (pl.DataFrame) – Atomic configuration containing columns: ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Simulation box object containing boundary conditions and dimensions.

Returns:

NDArray[np.float64] – Per-atom virial tensors. Shape: (N, 9) where N is the number of atoms. Components: [W_xx, W_xy, W_xz, W_yx, W_yy, W_yz, W_zx, W_zy, W_zz]. Units: eV.

plot(fig: Figure | None = None, ax: List[Axes] | None = None) → Tuple[Figure, List[Axes]][source]¶

Plot EAM potential functions: embedding function, electron density, and pair potential.

Creates three subplots showing:

Embedding function F(ρ) vs electron density ρ
Electron density function ρ(r) vs distance r
Pair potential φ(r) vs distance r

Parameters:

fig (Figure, optional) – Matplotlib figure object. If None, creates a new figure.
ax (List[Axes], optional) – List of three matplotlib axes objects. If None, creates new axes.

Returns:

fig (Figure) – Matplotlib figure object containing the plots.
ax (List[Axes]) – List of three matplotlib axes objects [ax_F, ax_rho, ax_phi].

class mdapy.eam.EAMAverage(filename: str, concentration: List[float], output_name: str | None = None)[source]¶

Bases: EAM

Generate average EAM (A-atom) potential for alloy investigation.

This class creates an average EAM potential, which is useful in alloy investigation. The A-atom potential has a similar formula to the original EAM potential:

\[E_{i}=\sum_{i} F^{A}\left(\bar{\rho}_{i}\right)+\frac{1}{2} \sum_{i, j \neq i} \phi_{ij}^{A A},\]

\[F^{A}\left(\bar{\rho}_{i}\right)=\sum_{\alpha} c_{\alpha} F^{\alpha}\left(\bar{\rho}_{i}\right),\]

\[\phi_{ij}^{A A}=\sum_{\alpha, \beta } c_{\alpha} c_{\beta} \phi_{ij}^{\alpha \beta},\]

\[\quad \bar{\rho}_{i}=\sum_{j \neq i} \sum_{\alpha} c_{\alpha} \rho_{ij}^{\alpha},\]

where \(A\) denotes an average-atom.

Note

If you use this module in publication, you should also cite the original paper. Average-atom interatomic potential for random alloys

Parameters:

filename (str) – Filename of eam.alloy file.
concentration (List[float]) – Atomic ratio list, such as [0.5, 0.5]. The summation should equal 1.
output_name (str, optional) – Filename of generated average EAM potential. If None, auto-generates name.

Examples

>>> potential = EAMAverage("CuNiAl.eam.alloy", [0.25, 0.25, 0.5])

class mdapy.eam.EAMGenerator(elements_list: List[str], output_filename: str | None = None, nr: int = 2000, nrho: int = 2000, rst: float = 0.5)[source]¶

Bases: object

Generator for EAM (Embedded Atom Method) potential files.

This class creates EAM potential files in the eam.alloy format for single or multi-element systems. The implementation is based on the Zhou-Johnson-Wadley parameterization.

Parameters:

elements_list (List[str]) – List of element symbols (e.g., [‘Cu’, ‘Ni’, ‘Al’]).
output_filename (str, optional) – Name of the output file. If None, generates name from elements.
nr (int, optional) – Number of points in distance grid (default: 2000).
nrho (int, optional) – Number of points in electron density grid (default: 2000).
rst (float, optional) – Minimum distance for pair potential computation (default: 0.5 Å).

Raises:

ValueError – If any element is not in the supported elements list.

References

Zhou, X. W., Johnson, R. A., & Wadley, H. N. G. (2004). Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers. Physical Review B, 69(14), 144113. https://doi.org/10.1103/PhysRevB.69.144113

Examples

>>> gen = EAMGenerator(["Cu", "Ni"])
>>> # This automatically generates a 'CuNi.eam.alloy' file

SUPPORTED_ELEMENTS: Tuple[str, ...] = ('Cu', 'Ag', 'Au', 'Ni', 'Pd', 'Pt', 'Al', 'Pb', 'Fe', 'Mo', 'Ta', 'W', 'Mg', 'Co', 'Ti', 'Zr')¶

DEFAULT_NR: int = 2000¶

DEFAULT_NRHO: int = 2000¶

DEFAULT_RST: float = 0.5¶

mdapy.elastic module¶

mdapy.elastic.strain_from_deformation(deformation: ndarray) → ndarray[source]¶

Green-Lagrange strain from a deformation gradient F.

E = ½(Fᵀ F − I)

Matches pymatgen Strain.from_deformation.

Parameters:: deformation ((3,3) deformation matrix F)
Returns:: strain ((3,3) symmetric strain tensor in Voigt-ready form)

mdapy.elastic.strain_to_voigt(strain: ndarray) → ndarray[source]¶: 3×3 symmetric strain tensor → 6-vector Voigt notation. Convention (matches pymatgen): [e11, e22, e33, 2e23, 2e13, 2e12]

mdapy.elastic.stress_to_voigt(stress: ndarray) → ndarray[source]¶: 3×3 stress tensor → 6-vector Voigt notation. Convention: [s11, s22, s33, s23, s13, s12]

mdapy.elastic.apply_deformation_to_cell(cell: ndarray, deformation: ndarray) → ndarray[source]¶

Apply deformation matrix F to a cell (row-vector convention).

pymatgen applies F to lattice row vectors as: new_lattice = old_lattice @ F.T (because lattice rows are basis vectors, F acts on column vectors)

Parameters:

cell ((3,3) original cell, row vectors [a; b; c])
deformation ((3,3) deformation gradient F)

Returns:

new_cell ((3,3) deformed cell)

mdapy.elastic.apply_deformation_to_positions(positions: ndarray, old_cell: ndarray, new_cell: ndarray) → ndarray[source]¶: Map Cartesian positions into the deformed cell via fractional coordinates. Atoms stay at the same fractional coordinates — only the box changes.

class mdapy.elastic.DeformedStructureSet(system: System, norm_strains: Sequence[float] = (-0.01, -0.005, 0.005, 0.01), shear_strains: Sequence[float] = (-0.06, -0.03, 0.03, 0.06))[source]¶

Bases: object

Generate a set of deformed cells for elastic constant fitting.

Deformation strategy (identical to pymatgen):

Normal modes (0,0) (1,1) (2,2) : each strain in norm_strains

Shear modes (0,1) (0,2) (1,2) : each strain in shear_strains

Total configurations = 3xlen(norm_strains) + 3xlen(shear_strains)

Parameters:

system (optimized structure)
norm_strains (normal strain magnitudes (default same as pymatgen))
shear_strains (shear strain magnitudes (default same as pymatgen))

deformations¶

Type:: list of (3,3) deformation matrices F

deformed_systems¶

Type:: list of deformed systems

class mdapy.elastic.ElasticTensor(voigt: ndarray)[source]¶

Bases: object

6x6 elastic stiffness tensor in Voigt notation.

Build via the class method:: et = ElasticTensor.from_independent_strains(strain_list, stress_list, eq_stress)

voigt¶

Type:: (6,6) Cij matrix

classmethod from_independent_strains(strains: List[ndarray], stresses: List[ndarray], eq_stress: ndarray | None = None, tol: float = 1e-10) → ElasticTensor[source]¶

Least-squares fit of elastic constants from independent strain–stress pairs.

Exact port of pymatgen ElasticTensor.from_independent_strains.

Algorithm:

Convert all strains/stresses to Voigt 6-vectors.
Group by “strain state” — which Voigt component is active.
For each strain state i (0–5) and each response component j (0–5):
C_ij = slope of stress_j vs strain_i (linear polyfit, degree 1)
Include the zero-strain equilibrium point in every group’s fit.

Parameters:

strains (list of (3,3) strain tensors (output of strain_from_deformation))
stresses (list of (3,3) stress tensors in GPa)
eq_stress ((3,3) equilibrium stress tensor (at zero strain))
tol (zero-out entries smaller than this in the final tensor)

Returns:

ElasticTensor with .voigt attribute = (6,6) Cij in same units as stresses

print()[source]¶

vrh()[source]¶

Voigt-Reuss-Hill polycrystalline averages.

Converts the single-crystal elastic tensor into isotropic polycrystalline properties by averaging over all grain orientations.

Three averaging schemes:

Voigt : assumes uniform strain across grains → upper bound

Reuss : assumes uniform stress across grains → lower bound

Hill : arithmetic mean of Voigt and Reuss → best estimate

The Voigt-Reuss gap reflects single-crystal anisotropy: a larger gap means stronger directional dependence (e.g. HCP Mg has a wider gap than FCC Cu).

Returns:: dict with keys (all in GPa except nu which is dimensionless) – K_V, K_R, K_H : bulk modulus G_V, G_R, G_H : shear modulus E : Young’s modulus, Hill only nu : Poisson’s ratio, Hill only

print_vrh()[source]¶: Print Voigt-Reuss-Hill polycrystalline mechanical properties.

mdapy.elastic.get_elastic_constant(system: System, calc: CalculatorMP, norm_strains: Sequence[float] = (-0.01, -0.005, 0.005, 0.01), shear_strains: Sequence[float] = (-0.06, -0.03, 0.03, 0.06), fmax: float = 0.0001) → ElasticTensor[source]¶

Workflow to compute elastic constants from a system.

Parameters:

system (atomic structure)
calc (calculator for computing energy, force and stress)
norm_strains (normal strain magnitudes, defaults is (-0.01, -0.005, 0.005, 0.01))
shear_strains (shear strain magnitudes, defaults is (-0.06, -0.03, 0.03, 0.06))
fmax (converge limit for minimization, defaults is 1e-4)

Returns:

ElasticTensor with .voigt attribute = (6,6) Cij in same units as stresses

mdapy.identify_diamond_structure module¶

class mdapy.identify_diamond_structure.IdentifyDiamondStructure(data: DataFrame, box: Box, verlet_list: ndarray[tuple[Any, ...], dtype[int32]] | None = None)[source]¶

Bases: object

This modifier analyzes the coordination structure of particles to identify diamond lattice structures, including both cubic and hexagonal diamond polymorphs. The identification is performed using an extended Common Neighbor Analysis method that considers both first and second nearest neighbor shells.

The pattern array assigns one of the following structure type identifiers to each particle:

0 Other
1 Cubic diamond
2 Cubic diamond (1st neighbor)
3 Cubic diamond (2nd neighbor)
4 Hexagonal diamond
5 Hexagonal diamond (1st neighbor)
6 Hexagonal diamond (2nd neighbor)

Parameters:

data (pl.DataFrame) – Particle properties data frame containing columns: ‘x’, ‘y’, ‘z’ for positions.
box (Box) – Simulation cell object containing boundary conditions and dimensions.
verlet_list (NDArray[np.int32], optional) – Pre-computed neighbor list for each particle with at least 4 neighbors. If None, k-nearest neighbor search (k=4) is performed automatically. Shape: (N, k) where N is the number of particles and k ≥ 4.

data¶

Input particle configuration data.

Type:: pl.DataFrame

box¶

Simulation cell information.

Type:: Box

verlet_list¶

Neighbor list used for structure identification.

Type:: NDArray[np.int32] or None

pattern¶

Structure type classification for each particle. Available after calling compute(). Shape: (N,) where N is the number of particles.

Type:: NDArray[np.int32]

References

[1] Maras E, Trushin O, Stukowski A, et al. Global transition path search for dislocation formation in Ge on Si (001)[J]. Computer Physics Communications, 2016, 205: 13-21.

compute() → None[source]¶

Perform diamond structure identification on all particles. The method automatically handles different system configurations:

For systems with ≤4 particles and no periodic boundaries, all particles are classified as type 0 (Other)
For small systems (thickness < 15 angstroms) with periodic boundaries, the simulation cell is automatically replicated to ensure accurate neighbor identification
If no neighbor list is provided, a k-nearest neighbor search with k=4 is performed automatically

mdapy.identify_fcc_planar_faults module¶

class mdapy.identify_fcc_planar_faults.IdentifyFccPlanarFaults(structure_types: ndarray, ptm_indices: ndarray, cal_esf: bool = True)[source]¶

Bases: object

Identify planar faults in FCC structures including stacking faults, twin boundaries, and ESF.

This class identifies intrinsic and extrinsic stacking faults as well as coherent twin boundaries in FCC crystal structures. The algorithm is based on the method described in the OVITO documentation, with an additional capability to identify Extrinsic Stacking Faults (ESF).

Parameters:

structure_types (np.ndarray) – Array of structure type identifiers for each atom, shape (N,). Typically from PTM analysis where 1=FCC, 2=HCP, etc.
ptm_indices (np.ndarray) – PTM nearest neighbor indices for each atom, shape (N, 12). Contains the indices of the 12 nearest neighbors in FCC/HCP ordering.
cal_esf (bool, default=True) – Whether to identify Extrinsic Stacking Faults (ESF) in addition to intrinsic stacking faults and twin boundaries.

fault_types¶

Array of fault type identifiers for each atom, shape (N,). Available after calling compute().

Fault type codes:

0: Non-hcp atoms (e.g. perfect fcc or disordered)
1: Indeterminate hcp-like (isolated hcp-like atoms, not forming a planar defect)
2: Intrinsic stacking fault (ISF, two adjacent hcp-like layers)
3: Coherent twin boundary (TB, one hcp-like layer)
4: Multi-layer stacking fault (three or more adjacent hcp-like layers)
5: Extrinsic Stacking Fault (ESF, if cal_esf=True)

Type:: np.ndarray

structure_types¶

Input structure types array.

Type:: np.ndarray

ptm_indices¶

Input PTM neighbor indices array.

Type:: np.ndarray

cal_esf¶

Flag for ESF identification.

Type:: bool

References

OVITO Identify FCC Planar Faults Modifier. https://www.ovito.org/manual/reference/pipelines/modifiers/identify_fcc_planar_faults.html

compute()[source]¶

Compute fault type for each atom.

This method identifies planar faults by analyzing HCP-classified atoms and their neighbors. Results are stored in fault_types.

mdapy.knn module¶

class mdapy.knn.NearestNeighbor(data: DataFrame, box: Box, k: int)[source]¶

Bases: object

Perform a nearest-neighbor search for atoms within a periodic or non-periodic box.

This class computes the indices and distances of the k nearest neighbors for each atom, considering periodic boundary conditions (PBC) if applicable.

For small systems, where the number of atoms is less than k, the simulation box is automatically replicated according to its boundary conditions to ensure enough neighbors are available.

Parameters:

data (pl.DataFrame) – A Polars DataFrame containing atomic coordinates with columns "x", "y", and "z".
box (Box) – The simulation box, defined as an instance of mdapy.box.Box.
k (int) – The number of nearest neighbors to search for. Must be less than 25.

indices_py¶

A 2D integer array of shape (N, k), storing the indices of the nearest neighbors for each atom.

Type:: np.ndarray

distances_py¶

A 2D float array of shape (N, k), storing the distances to the corresponding neighbors.

Type:: np.ndarray

_enlarge_data¶

Internal replicated atomic data when periodic extension is required.

Type:: pl.DataFrame, optional

_enlarge_box¶

The enlarged simulation box corresponding to replicated atoms.

Type:: Box, optional

compute()[source]¶

Compute the nearest neighbors for all atoms in the system.

If the number of atoms is smaller than k and periodic boundaries are enabled, the system will be automatically replicated along the periodic directions to ensure sufficient neighbors are found.

Returns:

None – Results are stored in the following attributes:

indices_py: nearest neighbor indices.
distances_py: corresponding neighbor distances.

mdapy.lammps_potential module¶

mdapy.lammps_potential.silence(enabled: bool = True)[source]¶: Redirect stdout/stderr to /dev/null when enabled is True; otherwise no-op.

class mdapy.lammps_potential.LammpsPotential(pair_parameter: str, element_list: List[str], units: str = 'metal', centroid_stress: bool = False, cmdargs: List[str] | None = None, extra_commands: str | None = None, silence_lammps: bool = True)[source]¶

Bases: CalculatorMP

LAMMPS-based calculator that runs a single-point evaluation to obtain per-atom energies, forces, virials and the global stress.

Parameters:

pair_parameter (str) – The LAMMPS pair style / pair coeff commands as a single string. This string is passed directly to LAMMPS with commands_string.
element_list (List[str]) – List of element names supported by this potential. The index in this list defines the corresponding LAMMPS atom type (1-based).
units (str, optional) – Units for LAMMPS (default "metal"). Currently the code asserts that units == “metal”.
centroid_stress (bool, optional) – If True, uses compute centroid/stress/atom NULL in LAMMPS; otherwise uses compute stress/atom NULL.
cmdargs (list of str, optional) – Extra command-line arguments forwarded to lammps(cmdargs=...). Useful for accelerator packages, e.g. Kokkos: ["-k", "on", "g", "1", "-sf", "kk", "-pk", "kokkos", "newton", "on", "neigh", "half"]. These are appended after the default ["-echo", "none", "-log", "none", "-screen", "none"].
extra_commands (str, optional) – Extra LAMMPS commands executed right after the box/atoms are set up and before pair_parameter. Use this for things that must come before the pair style (e.g. "package kokkos newton on neigh half", "newton on", "atom_modify ..." overrides handled in-script).
silence_lammps (bool, optional) – If True (default), redirect LAMMPS stdout/stderr to /dev/null while the calculator runs. Set to False to see LAMMPS output for debugging.

calculate(data: DataFrame, box: Box) → None[source]¶

Run LAMMPS to calculate per-atom energies, forces and virials and compute global stress.

This function validates inputs, constructs a triclinic LAMMPS box, creates atoms, sets up computes, runs run 0, extracts LAMMPS computed quantities, converts units, reorders virials, and stores results in self.results.

Parameters:

data (polars.DataFrame) – Polars DataFrame with required columns: “x”, “y”, “z”, “element”.
box (Box) – Box object from mdapy.

Notes

The method relies on lammps Python bindings to exist and provide: - lammps(cmdargs=…), commands_string, create_atoms, - numpy.extract_atom(…), numpy.extract_compute(…), - numpy.extract_atom(“f”), and .close().
Virial unit conversion: virial = virial / 1e4 / 160.21766208 (converts LAMMPS reported units to eV).
The final global stress is computed as: stress = -(virial_tensor + virial_tensor.T) / (2 * box.volume) and returned in Voigt order [xx, yy, zz, yz, xz, xy].

get_energies(data: DataFrame, box: Box) → Any[source]¶

Return per-atom energies. If not already computed, triggers calculate().

Parameters:

data (polars.DataFrame)
box (Box)

Returns:

Any – Stored per-atom energies (as placed into self.results[“energies”]).

get_energy(data: DataFrame, box: Box) → Any[source]¶

Return total energy (sum of per-atom energies).

Parameters:

data (polars.DataFrame)
box (Box)

Returns:

Any – Sum of per-atom energies.

get_forces(data: DataFrame, box: Box) → Any[source]¶: Return per-atom forces; compute if necessary.

get_stress(data: DataFrame, box: Box) → Any[source]¶: Return global stress in Voigt order [xx, yy, zz, yz, xz, xy]; compute if necessary.

get_virials(data: DataFrame, box: Box) → Any[source]¶: Return per-atom virials (9 components) and compute if necessary.

mdapy.lammps_runner module¶

class mdapy.lammps_runner.LammpsRunner(system: System, pair_parameter: str, element_list: List[str], units: str = 'metal', cmdargs: List[str] | None = None, extra_commands: str | None = None, masses: dict | None = None, log_file: str | None = None, silence_lammps: bool = True)[source]¶

Bases: object

Persistent LAMMPS runner driving common simulation tasks (minimization with/without box relaxation, NVE/NVT/NPT MD) on top of an mdapy System.

Parameters:

system (System) – Initial mdapy System. Must contain x, y, z and element columns. If vx, vy, vz columns exist they are pushed into LAMMPS as initial velocities.
pair_parameter (str) – Multi-line LAMMPS string for pair_style + pair_coeff (and any related commands). Passed to commands_string.
element_list (list of str) – Element ordering. Defines the LAMMPS atom-type index (1-based). Every element in system.data must appear here.
units (str, optional) – LAMMPS units. Only "metal" is supported.
cmdargs (list of str, optional) – Extra args appended to lammps(cmdargs=...) after the default ["-echo", "none", "-log", "none", "-screen", "none"]. Use this for accelerator packages (Kokkos, OMP, GPU).
extra_commands (str, optional) – Extra LAMMPS commands executed once after the box/atoms are set up and before pair_parameter (e.g. "newton on" or "package omp 8").
masses (dict[str, float], optional) – Override masses (in g/mol). Keys are element symbols. Any element not in the dict falls back to data.atomic_masses.
log_file (str, optional) – If given, LAMMPS writes its log to this file (log <file>); otherwise no log file is created.
silence_lammps (bool, optional) – If True (default), redirect LAMMPS stdout/stderr to /dev/null across all calls. Set to False to see LAMMPS output for debugging.

minimize(etol: float = 1e-08, ftol: float = 1e-08, maxiter: int = 10000, maxeval: int = 100000, style: str = 'cg') → None[source]¶

Energy minimization at fixed cell.

Parameters:

etol (float) – Energy / force tolerances (LAMMPS minimize args).
ftol (float) – Energy / force tolerances (LAMMPS minimize args).
maxiter (int) – Max iterations and force evaluations.
maxeval (int) – Max iterations and force evaluations.
style (str) – Minimization style passed to min_style (cg, sd, fire, …).

minimize_box(etol: float = 1e-08, ftol: float = 1e-08, maxiter: int = 10000, maxeval: int = 100000, style: str = 'cg', pressure_bar: float = 0.0, relax_style: str = 'iso') → None[source]¶

Energy minimization with simultaneous box relaxation (fix box/relax).

Parameters:

pressure_bar (float) – Target external pressure (bar in metal units).
relax_style (str) – One of "iso", "aniso", "tri". For "tri" all six components are relaxed independently.

run_md(ensemble: str = 'nvt', nsteps: int = 1000, timestep: float = 0.001, temp: float | Tuple[float, float] = 300.0, tdamp: float = 0.1, pressure_bar: float | Tuple[float, float] = 0.0, pdamp: float = 1.0, press_style: str = 'iso', init_velocity_temp: float | None = None, velocity_seed: int = 12345, thermo_freq: int = 100, thermo_keywords: str | None = None, dump_file: str | None = None, dump_freq: int = 100, dump_keywords: str | None = None) → None[source]¶

Run molecular dynamics.

Parameters:

ensemble (str) – One of "nve", "nvt", "npt".
nsteps (int) – Number of integration steps.
timestep (float) – MD timestep (ps in metal units).
temp (float or (float, float)) – Target temperature for nvt/npt. A tuple (T_start, T_stop) gives a linear ramp; a scalar uses it for both.
tdamp (float) – Temperature damping (ps).
pressure_bar (float or (float, float)) – Target pressure (bar) for npt; tuple gives a ramp.
pdamp (float) – Pressure damping (ps).
press_style (str) – One of "iso", "aniso", "tri". Controls how the barostat couples cell components.
init_velocity_temp (float, optional) – If set, velocity all create <T> <seed> dist gaussian is issued before the run. If None and the system already had vx/vy/vz, those persist; otherwise atoms start at rest.
velocity_seed (int) – RNG seed for velocity create.
thermo_freq (int) – thermo output frequency. Set to 0 to keep the LAMMPS default behaviour.
thermo_keywords (str, optional) – Custom thermo_style custom <...> columns. Defaults to a common set including step/temp/pe/press/vol/cell.
dump_file (str, optional) – If set, a custom dump is created for the run.
dump_freq (int) – Dump frequency (every N steps). Ignored if dump_file is None.
dump_keywords (str, optional) – Columns for the dump (LAMMPS custom syntax). Defaults to "id type element x y z".

get_system() → System[source]¶

Build a fresh mdapy System from the current LAMMPS state.

Notes

Positions/velocities are returned in the current LAMMPS frame. For a non-orthogonal input box that required rotation at setup, the returned System lives in the rotated (LAMMPS) frame — call System.write_data and reload if you need the original frame.

close() → None[source]¶

mdapy.lindemann_parameter module¶

class mdapy.lindemann_parameter.LindemannParameter(pos_list: ndarray, only_global: bool = False)[source]¶

Bases: object

Calculate the Lindemann index for distinguishing melt processes.

This class computes the Lindemann index, which is useful for distinguishing the melting process and determining the melting points of nano-particles. The Lindemann index is defined as the root-mean-square bond-length fluctuation with the following mathematical expression:

\[\left\langle\sigma_{i}\right\rangle=\frac{1}{N_{p}(N_{p}-1)} \sum_{j \neq i} \frac{\sqrt{\left\langle r_{i j}^{2}\right\rangle_t-\left\langle r_{i j}\right\rangle_t^{2}}}{\left\langle r_{i j}\right\rangle_t}\]

where \(N_p\) is the number of particles, \(r_{ij}\) is the distance between atom \(i\) and \(j\), and the brackets \(\left\langle \right\rangle_t\) represent a time average.

Parameters:

pos_list (np.ndarray) –
Array of particle positions with shape (\(N_f\), \(N_p\), 3), where \(N_f\) is the number of frames and \(N_p\) is the number of particles.

Warning

The positions MUST be unwrapped coordinates! For systems with periodic boundary conditions, you must provide unwrapped positions to correctly calculate inter-particle distances across time. Wrapped coordinates will produce incorrect results due to discontinuities when particles cross periodic boundaries.
only_global (bool, optional) – If True, only compute the global Lindemann index (faster, parallel computation). If False, compute both local and global indices (slower, serial computation). Default is False.

lindemann_trj¶

Global Lindemann index for the entire trajectory.

Type:: float

lindemann_frame¶

Lindemann index per frame with shape (\(N_f\),). Only available when only_global=False.

Type:: np.ndarray

lindemann_atom¶

Local Lindemann index per atom with shape (\(N_f\), \(N_p\)). Only available when only_global=False.

Type:: np.ndarray

Notes

This implementation is partly based on the work at https://github.com/N720720/lindemann for calculating the Lindemann index.
Memory requirement: This calculation has high memory requirements. You can estimate the required memory using: \(2 \times 8 \times N_p^2 / 1024^3\) GB.
Parallelization: If only the global Lindemann index is needed, the calculation can be performed in parallel. The local Lindemann index runs serially due to dependencies between different frames.
Algorithm: The Welford method is used to update the variance and mean of \(r_{ij}\) for numerical stability.

Examples

>>> import mdapy as mp
>>> import numpy as np

>>> # Generate a random walk trajectory with 200 frames and 1000 particles
>>> pos_list = np.cumsum(
...     np.random.choice([-1.0, 0.0, 1.0], size=(200, 1000, 3)), axis=0
... )

>>> # Calculate only global Lindemann index (faster)
>>> LDMG = mp.LindemannParameter(pos_list, only_global=True)
>>> LDMG.compute()
>>> print(f"Global Lindemann index: {LDMG.lindemann_trj:.6f}")

>>> # Calculate both local and global Lindemann indices
>>> LDML = mp.LindemannParameter(pos_list)
>>> LDML.compute()
>>> print(f"Global Lindemann index: {LDML.lindemann_trj:.6f}")
>>> print(f"First 5 frame indices: {LDML.lindemann_frame[:5]}")

>>> # Verify consistency
>>> np.isclose(LDML.lindemann_trj, LDMG.lindemann_trj)
True

>>> # Plot evolution
>>> fig, ax = LDML.plot()

compute() → None[source]¶

Perform the Lindemann index calculation.

This method calls the C++ backend for efficient computation. If only_global is True, only the global Lindemann index is computed; otherwise, both local and global indices are calculated.

The computation results are stored in the following attributes:

lindemann_trj : Global Lindemann index
lindemann_frame : Lindemann index per frame (if only_global=False)
lindemann_atom : Local Lindemann index per atom (if only_global=False)

plot(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶

Plot the evolution of Lindemann index per frame.

Parameters:

fig (Figure, optional) – Matplotlib figure object. If None, a new figure will be created.
ax (Axes, optional) – Matplotlib axes object. If None, a new axes will be created.

Returns:

fig (Figure) – The matplotlib figure object.
ax (Axes) – The matplotlib axes object.

Raises:

RuntimeError – If compute() has not been called with only_global=False.

Examples

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> pos_list = np.random.randn(100, 500, 3).cumsum(axis=0)
>>> ldm = LindemannParameter(pos_list)
>>> ldm.compute()
>>> fig, ax = ldm.plot()
>>> plt.savefig("lindemann_evolution.png")
>>> plt.show()

Notes

This method requires that compute() has been called with only_global=False, as it visualizes the frame-by-frame evolution of the Lindemann index.

mdapy.load_save module¶

class mdapy.load_save.BuildSystem[source]¶

Bases: object

classmethod from_file(filename: str, format: str | None = None) → Tuple[DataFrame, Box, Dict[str, Any] | None][source]¶

Load system from a file (supports .gz compression).

Parameters:

filename (str) – Path to the input file (can be .gz compressed).
format (str, optional) – File format. If None, inferred from file extension. Supported formats: ‘data’, ‘lmp’, ‘dump’, ‘poscar’, ‘xyz’, self difined mp. Add ‘.gz’ suffix for compressed files (e.g., ‘data.gz’, ‘xyz.gz’).

Returns:

Tuple[pl.DataFrame, mdapy.box.Box, Optional[Dict[str, Any]]] – DataFrame with atom data, simulation box, and optional global info.

static from_ovito(atom: DataCollection) → Tuple[pl.DataFrame, Box, Dict[str, Any]][source]¶

Load system from OVITO DataCollection.

Parameters:: atom (ovito.data.DataCollection) – OVITO data collection object.
Returns:: Tuple[pl.DataFrame, mdapy.box.Box, Dict[str, Any]] – DataFrame with atom data, simulation box, and global info.
Raises:: ImportError – If OVITO is not installed.

static from_ase(atom: Atoms) → Tuple[pl.DataFrame, Box][source]¶

Load system from ASE Atoms.

Parameters:: atom (ase.Atoms) – ASE Atoms object.
Returns:: Tuple[pl.DataFrame, mdapy.box.Box] – DataFrame with atom data and simulation box.
Raises:: ImportError – If ASE is not installed.

static from_array(pos: ndarray, box: int | float | Iterable[float] | ndarray | Box) → Tuple[DataFrame, Box][source]¶

Create system from position array.

Parameters:

pos (np.ndarray) – Atom positions (N x 3).
box (Union[int, float, Iterable[float], np.ndarray, mdapy.box.Box]) – Simulation box definition.

Returns:

Tuple[pl.DataFrame, mdapy.box.Box] – DataFrame with positions and simulation box.

static from_data(data: DataFrame, box: int | float | Iterable[float] | ndarray | Box) → Tuple[DataFrame, Box][source]¶

Create system from DataFrame.

Parameters:

data (pl.DataFrame) – DataFrame with at least ‘x’, ‘y’, ‘z’ columns.
box (Union[int, float, Iterable[float], np.ndarray, mdapy.box.Box]) – Simulation box definition.

Returns:

Tuple[pl.DataFrame, mdapy.box.Box] – Input DataFrame and simulation box.

static read_mp(filename: str) → Tuple[DataFrame, Box, Dict[str, Any] | None][source]¶

Read system from mp file.

Parameters:: filename (str) – Path to mp file, such as model.mp.
Returns:: Tuple[pl.DataFrame, mdapy.box.Box, Optional[Dict[str, Any]]] – DataFrame with atom data, simulation box, and global info.

static read_xyz(filename: str) → Tuple[DataFrame, Box, Dict[str, Any] | None][source]¶

Read system from XYZ file (supports .gz compression).

Parameters:: filename (str) – Path to XYZ file (can be .xyz.gz).
Returns:: (pl.DataFrame, mdapy.box.Box, dict) – Per-atom data, simulation cell, and the trailing key=value metadata from the comment line as a dict (with lattice, properties, pbc, origin removed since they’re already represented in the box / column schema).

Notes

Both classical (legacy) and extended XYZ are supported. The extended-XYZ properties= string is parsed token by token, supporting the conventional aliases:

pos:R:3 → x y z
species:S:1 / element:S:1 → element
vel:R:3 / velo:R:3 → vx vy vz
forces:R:3 / force:R:3 → fx fy fz

Multi-column user-defined fields (e.g. custom:R:5) become custom_0 … custom_4.

Performance: when the per-atom rows are uniformly single-space separated, parsing dispatches to pl.read_csv (fast, vectorised); otherwise it falls back to a Python str.split() pass that tolerates tabs, runs of spaces and trailing whitespace.

static read_poscar(filename: str) → Tuple[DataFrame, Box, Dict[str, Any] | None][source]¶

Read a VASP POSCAR file (supports .gz compression).

Format (line indices below are 0-based):

line 0: free-form comment.
line 1: universal scale factor (multiplies both the lattice and the Cartesian atom positions).
lines 2–4: three lattice vectors a, b, c (one per row).
line 5: element-symbol line (e.g. "Al Cu") or the per-species count line (e.g. "1 1") — distinguished by whether the first non-whitespace char is alphabetic vs. numeric.
line 6: per-species count line (only when line 5 is symbols).
optional Selective dynamics header (first non-blank char S / s).
coordinate-type line: first non-blank char in {D, d} → Direct, {C, c, K, k} → Cartesian / K-point.
N atom rows (3 floats each, plus 3 T/F flags when SD is on).
optional lattice-velocity / ion-velocity blocks.

Returns:: Tuple[pl.DataFrame, mdapy.box.Box, dict]

static read_data(filename: str) → Tuple[DataFrame, Box, Dict[str, Any] | None][source]¶

Read a single-frame LAMMPS data file (supports .gz compression).

Only atom_style atomic and atom_style charge are supported. Other styles (bond, molecular, full, …) raise NotImplementedError.

The optional image-flag columns nx ny nz are preserved as ix iy iz int columns alongside the wrapped x y z.

Parameters:: filename (str) – Path to data file (can be .gz).
Returns:: (pl.DataFrame, mdapy.box.Box, dict) – Per-atom data, simulation cell, and an empty global-info dict.

static parse_dump_frame(lines: List[str], source: str = '<dump>') → Tuple[DataFrame, Box, Dict[str, Any]][source]¶: Parse one frame’s worth of LAMMPS dump lines (the 9 header lines + N atom rows). Used by both read_dump() (single-frame) and mdapy.Trajectory (multi-frame).

static read_dump(filename: str) → Tuple[DataFrame, Box, Dict[str, Any]][source]¶

Read a single-frame LAMMPS dump file (text style, supports .gz).

Multi-frame dumps are rejected with a clear error — use mdapy.Trajectory for those.

Coordinate column variants understood (see LAMMPS dump custom docs):

x y z wrapped Cartesian → kept as x y z
xs ys zs scaled (fractional) → unscaled to Cartesian
xu yu zu unwrapped Cartesian → kept as x y z
xsu ysu zsu scaled-unwrapped → unscaled to Cartesian

Image flags ix iy iz and the optional element column from dump_modify ... element ... are preserved verbatim.

BOX BOUNDS header forms recognised:

ITEM: BOX BOUNDS xx yy zz — orthogonal
ITEM: BOX BOUNDS xy xz yz xx yy zz — restricted triclinic
ITEM: BOX BOUNDS abc origin — general triclinic (newer LAMMPS)

class mdapy.load_save.SaveSystem[source]¶

Bases: object

static to_ase(data: pl.DataFrame, box: Box) → Atoms[source]¶

Convert system to ASE Atoms object.

Parameters:

data (pl.DataFrame) – Particle data with at least ‘x’, ‘y’, ‘z’, ‘element’ columns.
box (mdapy.box.Box) – Simulation box.

Returns:

ase.Atoms – ASE Atoms object.

Raises:

ImportError – If ASE is not installed.
ValueError – If required columns are missing.

static to_ovito(data: pl.DataFrame, box: Box, global_info: Dict[str, Any] | None = None) → DataCollection[source]¶

Save system to OVITO DataCollection.

Parameters:

data (pl.DataFrame) – Particle informations.
box (Box) – Simulation cell.
global_info (Optional[Dict[str, Any]], optional) – Global attributes to add.

Returns:

ovito.data.DataCollection – OVITO data collection object.

Raises:

ImportError – If OVITO is not installed.

write_mp(data: DataFrame, box: Box, global_info: Dict[str, str]) → None[source]¶

Write system to mp file.

Parameters:

output_name (str) – Output file path.
data (pl.DataFrame) – Atom data with ‘x’, ‘y’, ‘z’, ‘element’.
box (Box) – Box information.
global_info (Dict) – Global information.

static write_xyz(output_name: str, box: Box, data: DataFrame, classical: bool = False, compress: bool = False, **kargs)[source]¶

Write system to XYZ file (with optional compression).

Parameters:

output_name (str) – Output file path.
box (mdapy.box.Box) – Simulation box.
data (pl.DataFrame) – Atom data with ‘x’, ‘y’, ‘z’, ‘element’.
classical (bool, optional) – Use classical XYZ format.
compress (bool, optional) – Compress output to .gz format.
**kargs – Additional key-value pairs for extended XYZ.

static write_poscar(output_name: str, box: Box, data: DataFrame, reduced_pos: bool = False, selective_dynamics: bool = False, compress: bool = False)[source]¶

Write system to POSCAR file (with optional compression).

Parameters:

output_name (str) – Output file path.
box (mdapy.box.Box) – Simulation box.
data (pl.DataFrame) – Atom data with ‘element’, ‘x’, ‘y’, ‘z’.
reduced_pos (bool, optional) – Use reduced coordinates.
selective_dynamics (bool, optional) – Include selective dynamics.
compress (bool, optional) – Compress output to .gz format.

static write_data(output_name: str, box: Box, data: DataFrame, element_list: List | None = None, num_type: int | None = None, data_format: str = 'atomic', compress: bool = False)[source]¶

Write system to LAMMPS data file (with optional compression).

Parameters:

output_name (str) – Output file path.
box (mdapy.box.Box) – Simulation box.
data (pl.DataFrame) – Atom data with ‘type’, ‘x’, ‘y’, ‘z’.
element_list (Optional[List], optional) – Element names for atom types.
num_type (Optional[int], optional) – Number of atom types.
data_format (str, optional) – ‘atomic’ or ‘charge’.
compress (bool, optional) – Compress output to .gz format.

static write_dump(output_name: str, box: Box, data: DataFrame, timestep: float = 0.0, compress: bool = False)[source]¶

Write system to LAMMPS dump file (with optional compression).

Parameters:

output_name (str) – Output file path.
box (mdapy.box.Box) – Simulation box.
data (pl.DataFrame) – Atom data with ‘type’, ‘x’, ‘y’, ‘z’.
timestep (float, optional) – Timestep value.
compress (bool, optional) – Compress output to .gz format.

static write_dump_frame_to(op, box: Box, data: DataFrame, timestep: float = 0.0) → None[source]¶: Write one LAMMPS dump frame to an already-open binary file handle op. Used by both write_dump() (single-frame) and mdapy.Trajectory (multi-frame).

mdapy.md_elastic module¶

Finite-temperature elastic constants from MD stress-strain¶

Implements the Aidan Thompson finite-T recipe (LAMMPS examples/ELASTIC/T):

(Optional) NPT pre-relax to equilibrium volume V(T) at target pressure.
Equilibrate the cell at T (NVT, langevin); average reference stress sigma_0.
For each Voigt direction d in {1..6} and each sign s in {-1,+1}: - Reset the cell from a saved restart. - change_box delta = s*delta*L0 in direction d (with proper xy/xz/yz coupling). - NVT (isothermal) or NVE (adiabatic) for n_equil + n_run steps. - Time-average the pressure tensor -> sigma_pm.
C_id = -(sigma_pos[i] - sigma_neg[i]) / (2*delta) for i,d in 1..6.
Symmetrize: C_ij = (C_ij + C_ji) / 2.

Each LAMMPS instance is a single subprocess; the 12 deformations within one temperature run on a ProcessPoolExecutor so a serial-Kokkos LAMMPS build parallelises to N cores. MDElastic.scan_parallel further runs all temperatures concurrently (1 NPT + 12 deformations per T, 13*N tasks total).

Driver uses the LAMMPS Python module (from lammps import lammps). pair_style / pair_coeff are passed by the caller — for serial CPU runs use -k on -sf kk only (no g 1 and no -pk kokkos GPU args).

class mdapy.md_elastic.MDElasticResult(cij_voigt: ndarray, V_eq: float, T_actual: float, stress_ref: ndarray, ensemble: str, temperature: float)[source]¶

Bases: object

Output of one MDElastic.run() call.

cij_voigt¶

Full 6x6 elastic stiffness tensor in GPa (already symmetrized).

Type:: (6,6) ndarray

V_eq¶

Equilibrium cell volume at T (Angstrom^3).

Type:: float

T_actual¶

Time-averaged temperature during reference run (K).

Type:: float

stress_ref¶

Time-averaged reference stress in GPa (should be ~0 if relax_volume).

Type:: (6,) ndarray

ensemble¶

“isothermal” (NVT) or “adiabatic” (NVE after equilibration).

Type:: str

temperature¶

Target temperature (K).

Type:: float

cubic_average() → Tuple[float, float, float][source]¶

Cubic-symmetry-averaged (C11, C12, C44) in GPa.

SQS / defected cells are not strictly cubic; this is the standard averaging convention (used by the LAMMPS reference script).

vrh() → dict[source]¶: Voigt-Reuss-Hill polycrystalline averages computed on the full 6x6.

born_stable_cubic() → bool[source]¶: Born stability for cubic (uses cubic_average).

print() → None[source]¶

class mdapy.md_elastic.MDElastic(system: System, pair_style: str, pair_coeff: str, elements: Sequence[str], temperature: float, pressure: float = 0.0, pressure_coupling: str = 'iso', delta: float = 0.01, timestep: float = 0.001, n_equil: int = 10000, n_run: int = 5000, nevery: int = 10, nrepeat: int = 10, ensemble: str = 'isothermal', thermostat: str = 'nose-hoover', thermostat_damp: float | None = None, seed: int = 87287, relax_volume: bool = True, n_relax: int = 50000, pdamp: float | None = None, accelerator: str = 'kokkos', n_workers: int | None = None, n_threads_ref: int = 1, n_threads_def: int = 1, n_gpus_ref: int = 0, n_gpus_def: int = 0, log_dir: str | PathLike | None = None, lammps_cmdargs: List[str] | None = None, work_dir: str | None = None, quiet: bool = False)[source]¶

Bases: object

Finite-T elastic constants via stress-strain MD (LAMMPS).

Parameters:

system (System) – mdapy System with element labels; will be written to a LAMMPS data file internally.
pair_style (str) – LAMMPS pair_style argument, e.g. "eam/alloy", "nep", "nep/kk", "meam".
pair_coeff (str) – Full pair_coeff line after pair_coeff, e.g. "* * NiCoCr.lammps.eam Ni Co Cr". Element order here must match elements below.
elements (list of str) – Element symbol order used to map system.data["element"] to LAMMPS atom types (1..N). Must match the species list in pair_coeff.
temperature (float) – Target temperature (K).
pressure (float, optional) – Target pressure for the optional NPT relax stage (bar). Default 0.
pressure_coupling ({"iso", "aniso"}, optional) – fix npt coupling. iso keeps the cell shape isotropic (rescales all three lattice vectors together — appropriate for nominally cubic SQS / chemically-disordered cells). aniso lets each axis relax independently (use for non-cubic structures). Default iso.
delta (float, optional) – Voigt strain magnitude. LAMMPS reference uses 2e-2; 5e-3..2e-2 is the useful range. Default 0.01.
timestep (float, optional) – MD timestep (ps). Default 0.001.
n_equil (int, optional) – Equilibration steps per stage. Default 10000.
n_run (int, optional) – Averaging steps per stage. Default 5000.
nevery (int, optional) – Parameters of LAMMPS fix ave/time. The averaging window is nevery * nrepeat. Defaults 10 and 10.
nrepeat (int, optional) – Parameters of LAMMPS fix ave/time. The averaging window is nevery * nrepeat. Defaults 10 and 10.
ensemble ({"isothermal", "adiabatic"}, optional) – Integration during the averaging run. Default “isothermal”.
thermostat ({"nose-hoover", "langevin"}, optional) – Default “nose-hoover” (Nose-Hoover NVT, lower stress noise).
thermostat_damp (float, optional) – Thermostat tau in ps. Default 100 * timestep.
seed (int, optional) – RNG seed. Default 87287.
relax_volume (bool, optional) – Run NPT pre-relax. Default True.
n_relax (int, optional) – NPT pre-relax steps. The first half equilibrates the barostat; the second half time-averages box dimensions (lx, ly, lz). The averaged values are reported as V_eq and used as len0 for the deformation strains, giving a noise-free V(T) without rescaling the actual cell. Default 50000.
pdamp (float, optional) – Barostat tau in ps. Default 1000 * timestep.
accelerator ({"kokkos", "omp", "gpu", "none"}, optional) – LAMMPS accelerator backend. kokkos (default) injects -k on [t N] [g M] -sf kk and uses pair_style/kk variants. omp uses the USER-OMP package: -pk omp N -sf omp. gpu uses the GPU package: -pk gpu M -sf gpu. none runs serial / MPI LAMMPS without any per-process accelerator. Pair-style suffix is added automatically by -sf X so the user’s pair_style string need not be modified.
n_workers (int, optional) – Number of deformation worker processes. The 12 deformations are scheduled across this pool. Default min(12, os.cpu_count()). Set to 1 to run sequentially in the main process.
n_threads_ref (int, optional) – OpenMP threads per LAMMPS process for the reference / deformation phases. Used by accelerator="kokkos" (as -k on t N) and accelerator="omp" (as -pk omp N). Default 1.
n_threads_def (int, optional) – OpenMP threads per LAMMPS process for the reference / deformation phases. Used by accelerator="kokkos" (as -k on t N) and accelerator="omp" (as -pk omp N). Default 1.
n_gpus_ref (int, optional) – GPUs per LAMMPS process. Used by accelerator="kokkos" (as -k on g M) and accelerator="gpu" (as -pk gpu M). Default 0.
n_gpus_def (int, optional) – GPUs per LAMMPS process. Used by accelerator="kokkos" (as -k on g M) and accelerator="gpu" (as -pk gpu M). Default 0.
log_dir (str or Path, optional) – Per-stage log files written here. The reference (NPT + ref NVT) log is log_dir/ref_T{T}.log (one per temperature); each deformation gets log_dir/def_T{T}_d{dir}_{sign}.log. Each log carries thermo every 100 steps with columns step temp press lx ly lz vol pe ke etotal pxx pyy pzz pxy pxz pyz.
lammps_cmdargs (list of str, optional) – Extra args appended to lammps(cmdargs=...). The accelerator prefix is injected automatically based on accelerator / n_threads_* / n_gpus_*.
work_dir (str, optional) – Directory for data + restart files. Default a tmpdir cleaned up on completion.
quiet (bool) – Suppress mdapy progress prints.

run() → MDElasticResult[source]¶

Execute the protocol for a single temperature.

Reference NPT+NVT runs in the main process; the 12 deformations run on a ProcessPoolExecutor of size min(n_workers, 12).

Run MDElastic at multiple temperatures with two-phase parallelism.

Phase 1: each temperature’s reference (NPT + NVT) runs in its own process; up to n_workers_ref simultaneously. Each process can use n_threads_ref OpenMP threads (passed via kwargs). Phase 2: each (T, direction, sign) deformation runs in its own process; up to n_workers_def simultaneously, each typically single-threaded (n_threads_def via kwargs).

Typical 16-core CPU + 4 temperatures:

scan_parallel(
    ..., n_workers_ref=4, n_threads_ref=4,    # 4*4 = 16 cores
         n_workers_def=12, n_threads_def=1,   # 12 procs single-thread
)

Parameters:

system (as named.)
temperatures (as named.)
work_dir (str or Path) – Persistent directory for LAMMPS data + per-T restart files.
n_workers_ref (int, optional) – Number of reference processes. Default len(temperatures).
n_workers_def (int, optional) – Number of deformation processes. Default min(12 * len(temperatures), os.cpu_count()).
n_workers (int, optional) – Convenience alias: sets both n_workers_ref and n_workers_def. Ignored if either is set explicitly.
log_dir (str or Path, optional) – Directory for per-stage LAMMPS log files. Each T’s stage 1 (NPT + ref NVT + averaging) goes to log_dir/ref_T{T}.log; each deformation goes to log_dir/def_T{T}_d{dir}_{sign}.log. Useful for checking convergence + debugging.
**kwargs (forwarded to MDElastic.__init__ for each T (including) – n_threads_ref / n_threads_def / thermostat / etc.).

Returns:

pl.DataFrame with one row per T, columns – T, V_eq, T_actual, C11, C12, C44 (cubic averages), K_VRH, G_VRH, E_VRH, nu_VRH, stable, stress_ref_max, c{ij} (full 6x6 entries).

static thermal_expansion(df: DataFrame, n_unit_cells: int) → DataFrame[source]¶

Compute lattice constant and thermal expansion from a scan_parallel result DataFrame.

Parameters:

df (pl.DataFrame) – Output of MDElastic.scan_parallel. Must contain T and V_eq columns.
n_unit_cells (int) – Number of conventional unit cells in the simulation cell. For a 3x3x3-of-conventional-FCC SQS (108 atoms with 4 atoms/unit cell), pass n_unit_cells=27.

Returns:

pl.DataFrame – Original columns plus a (lattice constant in Å, equal to (V / n_unit_cells)**(1/3)), alpha_V (volumetric thermal expansion in 1/K, centred finite-difference of V(T)) and alpha_L (linear thermal expansion, alpha_V / 3 under the cubic-symmetry assumption).

static parse_log(log_path: str | PathLike, skip_fraction: float = 0.5) → dict[source]¶

Parse a per-stage LAMMPS log written by MDElastic and return time-averages of the thermo columns over the equilibrated portion.

Reads the thermo_style custom step temp press lx ly lz vol pe ke etotal pxx pyy pzz pxy pxz pyz block emitted by stage 1 / stage 2 runs. The first skip_fraction of rows is discarded as equilibration; the rest is averaged.

Parameters:

log_path (path to ref_T{T}.log or similar.)
skip_fraction (float in [0, 1), default 0.5.)

Returns:

dict with keys T_avg, press_avg, a_avg (Å, from lx ly lz),
V_avg (Å^3), n_samples, plus stress_voigt 6-vector
[pxx, pyy, pzz, pyz, pxz, pxy] in bar.

mdapy.mean_squared_displacement module¶

class mdapy.mean_squared_displacement.MeanSquaredDisplacement(pos_list: ndarray, mode: Literal['window', 'direct'] = 'window')[source]¶

Bases: object

Mean Squared Displacement (MSD) calculator with optional FFT acceleration.

The MSD quantifies particle motion over time:

\[MSD(t) = \langle \lvert \mathbf{r}(t_0 + t) - \mathbf{r}(t_0) \rvert^2 \rangle\]

where \(\mathbf{r}(t)\) is the particle position at time \(t\) and \(\langle \cdot \rangle\) denotes ensemble averaging.

Two computation modes are available:

window (FFT-accelerated, O(N log N)):
Efficient for long trajectories, based on the Wiener-Khinchin theorem:

\[MSD(m) = S_1(m) - 2 S_2(m)\]

with

\[S_1(m) = \frac{1}{N-m} \sum_{t=0}^{N-m-1} \Big[ \mathbf{r}^2(t) + \mathbf{r}^2(t+m) \Big]\]

\[S_2(m) = \frac{1}{N-m} \sum_{t=0}^{N-m-1} \mathbf{r}(t) \cdot \mathbf{r}(t+m)\]

FFT is used to efficiently compute the autocorrelation \(S_2(m)\).
direct (O(N)):
Computes MSD relative to the initial frame:

\begin{eqnarray*} MSD(t) =& \dfrac{1}{N} \sum_{i=1}^{N} (r_i(t) - r_i(0))^2 \\ \end{eqnarray*}

Note

Input positions must be unwrapped. Wrapped coordinates across periodic boundaries will produce incorrect MSD.
For diffusive motion in 3D: \(MSD(t) \approx 6 D t\) where \(D\) is the diffusion coefficient.
Using pyFFTW (if installed) can accelerate FFT computation significantly for large trajectories.

pos_list¶

Particle positions, shape [Nframe, N, 3] (unwrapped).

Type:: np.ndarray

mode¶

Calculation mode (‘window’ or ‘direct’).

Type:: str

particle_msd¶

Per-particle MSD, shape [Nframe, N].

Type:: Optional[np.ndarray]

msd¶

Ensemble-averaged MSD, shape [Nframe].

Type:: Optional[np.ndarray]

References

Vania Calandrini, Eric Pellegrini, Paolo Calligari, Konrad Hinsen, Gerald R. Kneller. “Nmoldyn – interfacing spectroscopic experiments, molecular dynamics simulations and models for time correlation functions.” École thématique de la Société Française de la Neutronique, 12:201–232, 2011.

pyFFTW: https://github.com/pyFFTW/pyFFTW

compute() → None[source]¶: Compute MSD for all particles and time frames.

plot(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶

Plot MSD vs. frame number.

Parameters:

fig (Optional[Figure]) – Existing matplotlib figure.
ax (Optional[Axes]) – Existing matplotlib axes.

Returns:

Tuple[Figure, Axes] – Figure and axes for further customization.

mdapy.minimizer module¶

This module implements the Fast Inertial Relaxation Engine (FIRE2) algorithm for atomic structure optimization. It is a direct translation and adaptation of the FIRE2 method from ASE (Atomic Simulation Environment).

References

J. Guénolé, W.G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (FIRE) for energy minimization in atomistic simulations and its implementation in LAMMPS, Computational Materials Science 175 (2020) 109584. https://doi.org/10.1016/j.commatsci.2020.109584
S. Echeverri Restrepo, P. Andric, ABC-FIRE: Accelerated Bias-Corrected Fast Inertial Relaxation Engine, Computational Materials Science 218 (2023) 111978. https://doi.org/10.1016/j.commatsci.2022.111978

When optimize_cell=True, this implementation behaves equivalently to ase.optimize.FIRE2 combined with ase.constraints.UnitCellFilter, allowing simultaneous optimization of both atomic positions and cell shape.

class mdapy.minimizer.FIRE(system: System, dt: float = 0.1, maxstep: float = 0.2, dtmax: float = 1.0, dtmin: float = 0.002, Nmin: int = 20, finc: float = 1.1, fdec: float = 0.5, astart: float = 0.25, fa: float = 0.99, use_abc: bool = False, optimize_cell: bool = False, mask=None, cell_factor=None, hydrostatic_strain: bool = False, constant_volume: bool = False, scalar_pressure: float = 0.0)[source]¶

Bases: object

Implementation of the FIRE2 (and optional ABC-FIRE) energy minimization algorithm.

This optimizer relaxes atomic structures using a velocity-Verlet–like integration with adaptive time-stepping, as proposed in the FIRE method. It can optionally optimize both the atomic positions and simulation cell (similar to ASE’s UnitCellFilter).

Parameters:

system (mdapy.system.System) – The system to be relaxed. It must provide get_force(), get_energy(), and get_virials() methods, and an associated Box object.
dt (float, optional) – Initial time step for the integration (default: 0.1).
maxstep (float, optional) – Maximum allowed atomic displacement per step (default: 0.2 Å).
dtmax (float, optional) – Maximum time step allowed during optimization (default: 1.0).
dtmin (float, optional) – Minimum time step (default: 2e-3).
Nmin (int, optional) – Minimum number of consecutive “good” steps (dot(v, f) > 0) before increasing the time step (default: 20).
finc (float, optional) – Factor by which the time step is increased after successful steps (default: 1.1).
fdec (float, optional) – Factor by which the time step is decreased after bad steps (default: 0.5).
astart (float, optional) – Initial value of the velocity-force mixing coefficient α (default: 0.25).
fa (float, optional) – Multiplicative decay factor for α after successful steps (default: 0.99).
use_abc (bool, optional) – Whether to use the ABC-FIRE algorithm (default: False).
optimize_cell (bool, optional) – Whether to optimize the simulation cell along with atomic positions (default: False). Equivalent to wrapping ASE’s UnitCellFilter.
mask (array_like of shape (3,3) or (6,), optional) – Tensor mask controlling which components of strain are relaxed. 1 means free, 0 means fixed. Default is fully relaxed.
cell_factor (float, optional) – Scale factor applied to cell degrees of freedom. Equivalent to ASE’s cell_factor. Default is the number of atoms.
hydrostatic_strain (bool, optional) – If True, only isotropic (hydrostatic) deformation is allowed (default: False).
constant_volume (bool, optional) – If True, constrain relaxation to approximately constant volume (default: False).
scalar_pressure (float, optional) – External scalar pressure (in eV/Å³) added to the enthalpy term (default: 0.0).

Notes

When optimize_cell=True, the deformation gradient is treated as additional three “virtual atoms”, following the formalism from Tadmor et al., Phys. Rev. B 59, 235 (1999).
This class modifies both atomic positions and the simulation box to minimize the potential energy (or enthalpy if pressure is applied).

Examples

>>> from mdapy.minimizer import FIRE
>>> from mdapy import build_crystal
>>> from mdapy.nep import NEP
>>> system = build_crystal("Al", "fcc", 4.05)
>>> system.calc = NEP("Al.txt")  # A NEP model
>>> fire = FIRE(system, optimize_cell=True)
>>> fire.run(steps=1000, fmax=1e-4)

get_forces() → ndarray[source]¶

Compute the combined atomic and (optionally) cell forces.

Returns:: np.ndarray – Force array of shape (N, 3) when optimize_cell=False, or (N+3, 3) when optimize_cell=True. The last three rows correspond to the cell forces derived from the virial stress tensor.

update_data_box(extended_dr: ndarray)[source]¶

Update atomic positions and simulation cell based on the provided displacement.

Parameters:: extended_dr (np.ndarray) – Displacement array of shape (N, 3) if not optimizing the cell, or (N+3, 3) when optimize_cell=True.

run(steps: int, fmax=0.0001, show_process=True) → bool[source]¶

Run the FIRE relaxation process.

Parameters:

steps (int) – Maximum number of optimization steps.
fmax (float, optional) – Convergence criterion for maximum force (default: 1e-4 eV/Å).
show_process (bool, optional) – Whether to print per-step progress (default: True).

Returns:

bool – Return whether the minimization process is converged.

Notes

This function adaptively updates velocity, time step, and the mixing parameter α according to the FIRE2 algorithm. If use_abc=True, it applies the ABC-FIRE bias correction and per-direction displacement capping.

mdapy.neighbor module¶

class mdapy.neighbor.Neighbor(rc: float, box: Box, data: DataFrame, max_neigh: int | None = None)[source]¶

Bases: object

Construct the neighbor list for all atoms within a cutoff distance rc.

The Neighbor class builds the Verlet neighbor list and corresponding distance list for each atom in the system.

Parameters:

rc (float) – Cutoff radius for neighbor searching. Must be positive.
box (Box) – Simulation box represented by a mdapy.box.Box instance. Supports both orthogonal and triclinic boxes.
data (pl.DataFrame) – Atomic data containing at least the columns "x", "y", and "z", representing atomic coordinates.
max_neigh (int, optional) – Pre-allocated per-atom slot count. The fast path (used when max_neigh is given) writes directly into a (N, max_neigh) buffer and is significantly faster than the dynamic-sizer fallback. The C++ kernel guards each write against the slot bound, so an over-tight max_neigh cannot corrupt memory; if the true coordination exceeds max_neigh for any atom, compute() raises ValueError reporting the required size. When omitted, the dynamic-sizer kernel runs (slower, no size hint required).

rc¶

Neighbor search cutoff radius.

Type:: float

box¶

Simulation box used for the neighbor search.

Type:: Box

data¶

Input atomic data.

Type:: pl.DataFrame

max_neigh¶

Maximum number of neighbors allowed per atom (if specified).

Type:: Optional[int]

N¶

Number of atoms in the input data.

Type:: int

verlet_list¶

Integer array of shape (N, max_neigh) or dynamically sized, storing the neighbor indices for each atom.

Type:: np.ndarray

distance_list¶

Float array of the same shape as verlet_list, storing the corresponding neighbor distances.

Type:: np.ndarray

neighbor_number¶

Integer array of length N, storing the number of neighbors for each atom.

Type:: np.ndarray

_enlarge_data¶

Internal replicated atomic data when periodic extension is required.

Type:: pl.DataFrame, optional

_enlarge_box¶

The enlarged simulation box corresponding to replicated atoms.

Type:: Box, optional

compute()[source]¶

Build the neighbor list and compute interatomic distances.

Modifies self in place; self.verlet_list, self.distance_list and self.neighbor_number are populated. For systems whose box is too small for rc along periodic directions, self._enlarge_data and self._enlarge_box are also populated and the neighbor arrays index into them.

mdapy.nep module¶

class mdapy.nep.NEP(filename: str)[source]¶

Bases: CalculatorMP

NEP calculator for mdapy static calculation framework.

This class provides an interface to compute atomic properties using a pre-trained NEP (Neuroevolution Potential) machine learning model. It can calculate energies, forces, stresses, virials, descriptors, and latent space representations.

Parameters:: filename (str) – Path to the NEP model file (typically ‘nep.txt’ or similar)

calc¶

Underlying C++ NEP calculator object (NEP_CPU: https://github.com/brucefan1983/NEP_CPU)

Type:: nepcal.NEPCalculator

results¶

Dictionary storing computed results with keys:

‘energies’: per-atom potential energies
‘forces’: per-atom forces (N×3 array)
‘virials’: per-atom virials (N×9 array)
‘stress’: system stress tensor (6-component Voigt notation)

Type:: dict

Raises:: FileNotFoundError – If the specified model file does not exist

Examples

>>> from mdapy import build_crystal
>>> from mdapy.nep import NEP
>>> # Load NEP model
>>> calc = NEP("nep.txt")
>>> system = build_crystal("Cu", "fcc", 3.615)
>>> system.calc = calc
>>> # Calculate properties
>>> energy = system.get_energies()  # per-atom potential energy
>>> forces = system.get_forces()  # per-atom force
>>> virial = system.get_virials()  # per-atom virials
>>> stress = system.get_stress()  # system stress
>>> total_energy = system.get_energy()  # system potential energy
>>> descriptor = calc.get_descriptor(system.data, system.box)  # descriptor
>>> laten_space = calc.get_latentspace(system.data, system.box)  # latentspace

setAtoms(data: DataFrame, box: Box) → Tuple[ndarray, ndarray, ndarray, ndarray, ndarray][source]¶

Prepare atomic configuration data for NEP calculator.

This method converts the mdapy data format to the format required by the underlying NEP calculator, including type mapping and coordinate extraction.

Parameters:

data (pl.DataFrame) –
DataFrame containing atomic information with columns:
- ’x’, ‘y’, ‘z’: atomic coordinates
- ’element’: chemical symbols (e.g., ‘Cu’, ‘Au’)
box (Box) – Simulation box object

Returns:

type_list (np.ndarray) – Integer array mapping each atom to its type index
x (np.ndarray) – X-coordinates of all atoms
y (np.ndarray) – Y-coordinates of all atoms
z (np.ndarray) – Z-coordinates of all atoms
box_array (np.ndarray) – Box matrix (3×3 array)

Raises:

AssertionError – If required columns are missing or if elements are not in the model

Notes

The NEP model must have been trained with all elements present in the data.

calculate(data: DataFrame, box: Box)[source]¶

Calculate energies, forces, virials, and stress for the system.

This method performs the main NEP calculation and stores results in the self.results dictionary. It’s automatically called by other getter methods if results are not already cached.

Parameters:

data (pl.DataFrame) – DataFrame containing atomic positions and elements
box (Box) – Simulation box object

Notes

Results are stored in self.results with keys:

‘energies’: per-atom potential energies (N,)
‘forces’: per-atom forces (N, 3)
‘virials’: per-atom virials (N, 9)
‘stress’: system stress tensor (6,) in Voigt notation
‘charge’ : per-atom charge (N,), only available for qNEP
‘bec’ : per-atom bec (N, 9), only available for qNEP

The stress tensor is computed from virials as: σ = -(W + W^T) / (2V) where W is the total virial and V is volume.

get_energies(data: DataFrame, box: Box) → ndarray[source]¶

Get per-atom potential energies.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N,) containing energy for each atom

get_energy(data: DataFrame, box: Box) → float[source]¶

Get total potential energy of the system.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

float – Total potential energy (sum of per-atom energies)

get_virials(data: DataFrame, box: Box) → ndarray[source]¶

Get per-atom virial tensors.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N, 9) with virial components for each atom Ordered as: [v_xx, v_xy, v_xz, v_yx, v_yy, v_yz, v_zx, v_zy, v_zz]

get_forces(data: DataFrame, box: Box) → ndarray[source]¶

Get forces acting on each atom.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N, 3) containing force components [fx, fy, fz]

get_stress(data: DataFrame, box: Box) → ndarray[source]¶

Get stress tensor of the system.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Stress tensor in Voigt notation: [σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy]

get_charges(data: DataFrame, box: Box) → ndarray[source]¶

Get per-atom charges for qNEP.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N,) containing charge for each atom

get_bec(data: DataFrame, box: Box) → ndarray[source]¶

Get per-atom bec for qNEP.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N, 9) with bec components for each atom Ordered as: [bec_xx, bec_xy, bec_xz, bec_yx, bec_yy, bec_yz, bec_zx, bec_zy, bec_zz]

get_descriptor(data: DataFrame, box: Box) → ndarray[source]¶

Get atomic descriptors from the NEP model.

Descriptors are the learned feature representations in the hidden layers of the neural network, before the final energy prediction.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N, num_ndim) containing descriptor vectors for each atom, where num_ndim is the descriptor dimension

Notes

Descriptors can be useful for analyzing structural similarities.

get_latentspace(data: DataFrame, box: Box) → ndarray[source]¶

Get latent space representations from the NEP model.

The latent space is the representation in an intermediate layer of the neural network, capturing compressed structural information.

Parameters:

data (pl.DataFrame) – Atomic configuration data
box (Box) – Simulation box

Returns:

np.ndarray – Array of shape (N, num_nlatent) containing latent vectors for each atom, where num_nlatent is the latent space dimension

mdapy.nep4ase module¶

class mdapy.nep4ase.NEP4ASE(model_filename: str, atoms: Atoms | None = None)[source]¶

Bases: Calculator

NEP calculator compatible with ASE (Atomic Simulation Environment).

This class wraps the NEP calculator to work seamlessly with ASE’s calculator interface, allowing NEP models to be used in ASE workflows for geometry optimization, molecular dynamics, and other simulations.

Parameters:

model_filename (str) – Path to the NEP model file
atoms (Atoms, optional) – ASE Atoms object to attach the calculator to

calc¶

Underlying NEP calculator object

Type:: nepcal.NEPCalculator

results¶

Dictionary storing calculation results

Type:: dict

Examples

>>> from ase import Atoms
>>> from ase.optimize import BFGS
>>> from mdapy.nep import NEP4ASE
>>>
>>> # Create atoms object
>>> atoms = Atoms("Cu2", positions=[[0, 0, 0], [1.5, 0, 0]])
>>> atoms.set_cell([10, 10, 10])
>>> atoms.set_pbc(True)
>>>
>>> # Attach NEP calculator
>>> calc = NEP4ASE("nep.txt")
>>> atoms.calc = calc
>>>
>>> # Run geometry optimization
>>> opt = BFGS(atoms)
>>> opt.run(fmax=0.01)
>>>
>>> # Get energy and forces
>>> energy = atoms.get_potential_energy()
>>> forces = atoms.get_forces()

implemented_properties: List[str] = ['energy', 'energies', 'forces', 'stress', 'virials']¶: Properties calculator can handle (energy, forces, …)

set_nep(atoms: Atoms) → Tuple[ndarray, ndarray, ndarray, ndarray, ndarray][source]¶

Prepare ASE Atoms object for NEP calculation.

Converts ASE Atoms format to the format required by NEP calculator.

Parameters:

atoms (Atoms) – ASE Atoms object

Returns:

type_list (np.ndarray) – Integer type indices for each atom
x (np.ndarray) – X-coordinates
y (np.ndarray) – Y-coordinates
z (np.ndarray) – Z-coordinates
box (np.ndarray) – Cell matrix (3×3)

Raises:

AssertionError – If atoms contain elements not in the NEP model

calculate(atoms: Atoms = None, properties: List[str] = None, system_changes: List[str] = ['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms'])[source]¶

Perform calculation for requested properties.

This method is called by ASE when properties are requested. It calculates the specified properties and stores them in self.results.

Parameters:

atoms (Atoms, optional) – ASE Atoms object (uses self.atoms if None)
properties (list of str, optional) – List of properties to calculate (default: all implemented)
system_changes (list of str, optional) – List of changes since last calculation (for caching)

Notes

Special properties ‘descriptor’, ‘latentspace’, ‘charges’ and ‘bec’ are not in the standard ASE interface but are available through this calculator.

get_descriptor(atoms: Atoms = None, system_changes: List[str] = ['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms']) → ndarray[source]¶

Get atomic descriptors (not part of standard ASE interface).

Parameters:

atoms (Atoms, optional) – ASE Atoms object
system_changes (list of str, optional) – System changes since last calculation

Returns:

np.ndarray – Descriptor array of shape (N, num_ndim)

get_charges(atoms: Atoms = None, system_changes: List[str] = ['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms']) → ndarray[source]¶

Get atomic charges for qNEP model (not part of standard ASE interface).

Parameters:

atoms (Atoms, optional) – ASE Atoms object
system_changes (list of str, optional) – System changes since last calculation

Returns:

np.ndarray – Charge array of shape (N,)

get_bec(atoms: Atoms = None, system_changes: List[str] = ['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms']) → ndarray[source]¶

Get atomic bec for qNEP model (not part of standard ASE interface).

Parameters:

atoms (Atoms, optional) – ASE Atoms object
system_changes (list of str, optional) – System changes since last calculation

Returns:

np.ndarray – Bec array of shape (N, 9)

get_latentspace(atoms: Atoms = None, system_changes: List[str] = ['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms']) → ndarray[source]¶

Get latent space representations (not part of standard ASE interface).

Parameters:

atoms (Atoms, optional) – ASE Atoms object
system_changes (list of str, optional) – System changes since last calculation

Returns:

np.ndarray – Latent space array of shape (N, num_nlatent)

get_virials(atoms: Atoms = None, system_changes: List[str] = ['positions', 'numbers', 'cell', 'pbc', 'initial_charges', 'initial_magmoms']) → ndarray[source]¶

Get per-atom virials (not part of standard ASE interface).

Parameters:

atoms (Atoms, optional) – ASE Atoms object
system_changes (list of str, optional) – System changes since last calculation

Returns:

np.ndarray – Virial array of shape (N, 9)

mdapy.orthogonal_cell module¶

Convert a triclinic, fully-periodic mdapy.System to an equivalent orthogonal supercell.

This is the mdapy port of atomsk’s -orthogonal-cell option. The algorithm follows opt_orthocell.f90: along each Cartesian axis we search integer combinations of the original lattice vectors that produce a vector aligned with that axis, then keep the shortest such combination. Replicating the original cell over those integer ranges and filtering atoms inside the new box yields the orthogonal supercell.

A typical use case is converting an HCP / wurtzite / graphite primitive cell (120° hex) into its conventional orthogonal supercell — useful when downstream tooling does not handle triclinic boxes.

Usage¶

>>> import mdapy as mp
>>> hcp = mp.build_crystal("Mg", "hcp", a=3.21, c=5.21)
>>> ortho = mp.orthogonal_cell(hcp)              # 4-atom orthogonal cell
>>> ortho_min = mp.orthogonal_cell(hcp, find_minimal=True)

mdapy.orthogonal_cell.orthogonal_cell(system: System, find_minimal: bool = False, max_search: int = 20, tol: float = 1e-06) → System[source]¶

Convert a triclinic System to an equivalent orthogonal supercell.

Parameters:

system (System) – Input crystal. Must have all three boundaries periodic; an open boundary makes the orthogonal-supercell construction ill-defined and is rejected.
find_minimal (bool, default=False) – After building the orthogonal cell, also search for the smallest orthogonal sub-cell that reproduces the same crystal when replicated. Species-aware (won’t merge sites that carry different elements).
max_search (int, default=20) – Range of integer combinations searched along each lattice vector. Increase for very oblique cells where the canonical orthogonal supercell needs large coefficients; the cost grows as (2 max_search + 1)^3.
tol (float, default=1e-6) – Numerical tolerance for axis alignment, fractional-coordinate wrapping, and duplicate-atom detection.

Returns:

System – New System whose box is diagonal and whose atoms reproduce the input crystal exactly. The output is fully periodic.

Raises:

ValueError – If the input has any open boundary, if the input box is singular, or if no orthogonal supercell can be found within max_search.

Notes

Algorithm follows atomsk’s opt_orthocell.f90:

For each Cartesian axis, find the shortest non-zero integer combination of the input lattice vectors that produces a cartesian vector aligned with that axis.
Replicate the input cell by enough copies to cover the resulting orthogonal box.
Wrap atoms into the orthogonal box and deduplicate with tol precision; species-aware (an exact-position duplicate carrying a different element raises).
(Optional) reduce to the smallest periodic sub-cell.

Examples

HCP primitive (120° hex) → 4-atom orthogonal supercell:

>>> hcp = mp.build_crystal("Mg", "hcp", a=3.21, c=5.21)
>>> ortho = mp.orthogonal_cell(hcp)
>>> assert ortho.N == 4

Wurtzite GaN → orthogonal supercell with both species preserved:

>>> wz = mp.build_crystal(("Ga", "N"), "wurtzite", a=3.19, c=5.18)
>>> ortho = mp.orthogonal_cell(wz)
>>> sorted(set(ortho.data["element"].to_list()))
['Ga', 'N']

mdapy.parallel module¶

Unified parallelism control for mdapy.

mdapy uses several parallel backends:

OpenMP inside the C++ extension modules
Polars’ internal thread pool
Python multiprocessing (e.g. pigz)
Tachyon’s POSIX thread pool (rendering)

They are all driven by the environment variable MDAPY_NUM_THREADS.

Usage¶

Set the env var before importing mdapy:

import os
os.environ["MDAPY_NUM_THREADS"] = "1"  # serial
import mdapy

Rules¶

MDAPY_NUM_THREADS must be a positive integer (> 0) — otherwise mdapy raises ValueError at import time.
If the value exceeds os.cpu_count() a warning is emitted (oversubscribing usually hurts performance).
If unset, all parallel regions default to os.cpu_count() (use all cores).
Polars locks its thread pool at import time — that’s why mdapy mirrors the value to POLARS_MAX_THREADS in mdapy/__init__.py before polars is imported. This is the only env var mdapy writes.
OpenMP threads are passed explicitly into each C++ extension function via a num_threads(nt) clause, so mdapy never relies on OMP_NUM_THREADS. We deliberately do not mutate OMP_NUM_THREADS to avoid affecting other OpenMP libraries (torch, sklearn, scipy) running in the same process.

mdapy.parallel.get_num_threads() → int[source]¶

Resolve the thread count to use for any mdapy parallel region.

Reads MDAPY_NUM_THREADS (validated at mdapy import time). Falls back to os.cpu_count() when the env var is unset.

mdapy.phonon module¶

mdapy.pigz module¶

Parallel gzip compression using multiprocessing.

Simple, fast, and reliable: - Uses multiprocessing (true parallelism, no GIL) - Always uses maximum compression (level 9) - Focuses on: speed for large files, correctness for all files - Clean code, proper resource management

mdapy.pigz.compress_file(input_file, output_file=None)[source]¶

Compress a file to .gz format using multiprocessing.

This function provides parallel gzip compression for faster processing of large files. The number of worker processes is controlled by the MDAPY_NUM_THREADS environment variable (defaults to os.cpu_count()).

Parameters:

input_file (str) – Path to input file to compress.
output_file (str, optional) – Path to output file. If not specified, adds .gz to input filename.

Returns:

str – Path to the created compressed file.

Raises:

FileNotFoundError – If input file doesn’t exist.
ValueError – If input file already has .gz extension.
Exception – If compression fails for any reason.

Notes

Small files (<5 MB) use single-process compression.
Uses 512 KB chunks for optimal parallelism.

Examples

>>> compress_file("data.txt")
'data.txt.gz'

>>> compress_file("input.txt", "output.gz")
'output.gz'

mdapy.plotset module¶

Plotting Utilities for Scientific Publications¶

This module provides utilities for creating publication-quality figures with Matplotlib. It includes functions for configuring plot styles, creating figures with consistent formatting, and saving figures with uniform margins.

The module is designed to produce figures suitable for scientific papers and presentations, with sensible defaults for font sizes, line widths, tick marks, and color schemes.

Note: This module requires matplotlib to be installed. Install it with:: pip install matplotlib

Functions¶

set_figure : Create a figure with scientific style settings save_figure : Save a figure with uniform whitespace margins _pltset : Configure global Matplotlib style (internal) _cm2inch : Convert centimeters to inches (internal) _ensure_matplotlib : Check matplotlib availability (internal)

Examples

Basic usage for creating a simple plot:

>>> import numpy as np
>>> fig, ax = set_figure(figsize=(8.5, 7.0))
>>> x = np.linspace(0, 2 * np.pi, 100)
>>> ax.plot(x, np.sin(x), label="sin(x)")
>>> ax.set_xlabel("x")
>>> ax.set_ylabel("y")
>>> ax.legend()
>>> save_figure(fig, "output.png")

Creating a multi-panel figure:

>>> fig, axes = set_figure(figsize=(17, 7), nrow=1, ncol=2)
>>> for i, ax in enumerate(axes):
...     ax.plot(x, np.sin(x * (i + 1)))
...     ax.set_xlabel("x")
...     ax.set_ylabel("y")
>>> save_figure(fig, "multi_panel.pdf")

Create a Matplotlib figure and axes with scientific publication style.

Parameters:

figsize (tuple of (float or int, float or int), default=(8.5, 7.0)) – Figure dimensions in centimeters as (width, height).
figdpi (int, default=150) – Figure resolution in dots per inch.
nrow (int, default=1) – Number of subplot rows.
ncol (int, default=1) – Number of subplot columns.
color_cycler (list of str or tuple of str, optional) – Custom color palette for the figure.
**kwargs (Any) – Additional keyword arguments passed to _pltset.

Returns:

fig (Figure) – The created figure object.
ax (Axes or numpy.ndarray of Axes) – Single Axes or array of Axes objects.

mdapy.plotset.save_figure(fig: Figure, filename: str, dpi: int = 300, format: Literal['png', 'pdf', 'svg', 'eps', 'tiff'] = 'png', transparent: bool = True, pad_scale: float = 1.02) → None[source]¶

Save a figure with uniform whitespace margins on all sides.

Parameters:

fig (Figure) – The figure object to save.
filename (str) – Output filename.
dpi (int, default=300) – Resolution in dots per inch.
format ({'png', 'pdf', 'svg', 'eps', 'tiff'}, default='png') – Output file format.
transparent (bool, default=True) – Whether to use a transparent background.
pad_scale (float, default=1.02) – Scale factor to expand the bounding box.

mdapy.polyhedral_template_matching module¶

class mdapy.polyhedral_template_matching.PolyhedralTemplateMatching(structure: str, data: DataFrame, box: Box, rmsd_threshold: float = 0.1, verlet_list: ndarray | None = None)[source]¶

Bases: object

Polyhedral Template Matching (PTM) classifier for identifying local atomic structures.

This class implements the Polyhedral Template Matching algorithm to classify the local structural environment of atoms in a simulation dataset. It matches the neighborhood of each atom against predefined polyhedral templates for common crystal structures, providing robustness against thermal fluctuations and elastic strains compared to methods like Common Neighbor Analysis (CNA).

The supported structure types include:

FCC (Face-Centered Cubic)
HCP (Hexagonal Close-Packed)
BCC (Body-Centered Cubic)
ICO (Icosahedral)
SC (Simple Cubic)
DCUB (Cubic Diamond)
DHEX (Hexagonal Diamond)
Graphene

The algorithm computes a Root-Mean-Square Deviation (RMSD) for the best-matching template and assigns a structure type if below the specified threshold. Additional outputs include RMSD values, scale factors, orientations (as quaternions), and neighbor indices.

References

[1] Larsen P M, Schmidt S, Schiøtz J. Robust structural identification via polyhedral template matching[J]. Modelling and Simulation in Materials Science and Engineering, 2016, 24(5): 055007.

Parameters:

structure (str) – String specifying the structure types to consider, separated by hyphens (e.g., “fcc-hcp-bcc”). Supported values: “fcc”, “hcp”, “bcc”, “ico”, “sc”, “dcub”, “dhex”, “graphene”. Special values: “all” (all types), “default” (fcc, hcp, bcc).
data (pl.DataFrame) – DataFrame containing atomic positions with columns ‘x’, ‘y’, ‘z’. Optionally includes ‘type’ or ‘element’ for atomic types.
box (Box) – Simulation box object defining cell dimensions, origin, and boundary conditions.
rmsd_threshold (float, optional) – Maximum RMSD for a valid structure match. Particles exceeding this are classified as “Other”. Default is 0.1.
verlet_list (np.ndarray, optional) – Precomputed Verlet neighbor list (shape (N, 18)). If None, computed internally.

output¶

Array of shape (N, 8) containing per-atom results:

Column 0: Structure type (integer, 0=Other, 1=FCC, 2=HCP, 3=BCC, 4=ICO, 5=SC, 6=DCUB, 7=DHEX, 8=Graphene)
Column 1: Ordering type (interger, 0=Other, 1=L10, 2=L12 (A-site), 3=L12 (B-site), 4=B2, 5=zincblende / wurtzite)
Column 2: RMSD value
Column 3: Interatomic distance
Columns 4-7: Orientation quaternion (x, y, z, w)

Type:: np.ndarray

ptm_indices¶

Array of shape (N, 18) containing indices of neighboring atoms used in the template matching.

Type:: np.ndarray

compute() → None[source]¶

Perform the PTM computation and store results in self.output and self.ptm_indices.

Raises:: AssertionError – If invalid structure types are specified.

mdapy.potential_tool module¶

mdapy.potential_tool.run_gpumd(system: System, dirname: str, runin: str, nep_file: str, gpumd_path: str = 'gpumd') → None[source]¶

Warpper to run GPUMD.

Parameters:

system (System) – Use this to generate model.xyz.
dirname (str) – Run MD in this dirname. If it has been existing, raise error.
runin (str) – Use this to generate run.in file. Do not specify potential here, using nep_file parameter.
nep_file (str) – The path of nep.txt. Defaults to nep.txt.
gpumd_path (str) – The path of gpumd. Defaults to gpumd.

mdapy.potential_tool.rmse(predictions: ndarray, targets: ndarray) → float[source]¶

Compute Root-Mean-Square Error (RMSE).

Parameters:

predictions (np.ndarray) – Predicted values.
targets (np.ndarray) – Target reference values.

Returns:

float – RMSE value.

mdapy.potential_tool.read_thermo(path: str) → DataFrame[source]¶

Load GPUMD thermo.out file into a Polars DataFrame.

The file contains 18 columns: T, K, U, Pxx, Pyy, Pzz, Pyz, Pxz, Pxy, ax, ay, az, bx, by, bz, cx, cy, cz.

Parameters:: path (str) – Directory path containing thermo.out.
Returns:: pl.DataFrame – Thermo data with 18 columns.

mdapy.potential_tool.plot_nep_train(path: str, outname: str | None = None, figdpi: int | None = 300, **kargs) → Tuple[Figure, List[List[Axes]]][source]¶

Plot NEP training results, including energy/force/stress scatter plots and loss curves.

Parameters:

path (str) – Path containing NEP output files: loss.out, energy_train.out, force_train.out, stress_train.out.
outname (Optional[str], optional) – Filename to save figure, by default None.
figdpi (Optional[int], optional) – DPI of generated figure, by default 300.
**kargs – Extra arguments passed to set_figure().

Returns:

Tuple[Figure, List[List[Axes]]] – The figure and 2×2 axes list.

mdapy.potential_tool.get_sfe_fcc(name: str, a: float, calc: CalculatorMP) → float[source]¶

Compute stacking fault energy (SFE) for an FCC crystal.

Parameters:

name (str) – Element name.
a (float) – Lattice constant.
calc (CalculatorMP) – MDAPY calculator.

Returns:

float – Stacking fault energy in mJ/m².

mdapy.potential_tool.get_average_sfe_fcc_hea(N: int, element_list: List[str], element_ratio: List[float], a: float, calc: CalculatorMP) → ndarray[source]¶

Compute averaged SFE for random FCC HEA configurations.

Parameters:

N (int) – Number of random samples.
element_list (List[str]) – Element species.
element_ratio (List[float]) – Element ratios.
a (float) – Lattice constant.
calc (CalculatorMP) – MD calculator.

Returns:

np.ndarray – Array of: [i, mean_sfe up to sample i]

mdapy.potential_tool.get_eos(system: System, scale_start: float, scale_end: float, num: int) → ndarray[source]¶

Compute equation of state (EOS) by uniformly scaling volume.

Parameters:

system (System) – Input structure.
scale_start (float) – Initial scale factor (>0).
scale_end (float) – Final scale factor (> scale_start).
num (int) – Number of sampling points.

Returns:

np.ndarray – (num, 2) array of [volume_per_atom, energy_per_atom]

class mdapy.potential_tool.PCA(n_components: int)[source]¶

Bases: object

Simple PCA implementation (similar to sklearn PCA).

fit_transform(X: ndarray) → ndarray[source]¶

Fit PCA and return transformed coordinates.

Parameters:: X (np.ndarray) – Input array shape (n_samples, n_features).
Returns:: np.ndarray – Projected data with shape (n_samples, n_components).

mdapy.potential_tool.fps_sample(n_sample: int, descriptors: ndarray, start_idx: int = 0) → ndarray[source]¶

This function is used to sample the configurations using farthest point sampling method, based on the descriptors. It is helpful to select the structures during active learning process.

Parameters:

n_sample (int) – Number of structures one wants to select.
descriptors (np.ndarray) – Two dimensional ndarray, it can be any descriptors.
start_idx (int) – For deterministic results, fix the first sampled point index. Defaults to 0.

Returns:

sampled_indices (ndarray, shape (n_sample,))

mdapy.potential_tool.cfg2xyz(file_list: List[str] | str, type_dict: Dict[str, int], output_name: str = 'train.xyz', f_max: float = 25.0) → None[source]¶

Convert cfg file for MTP to xyz file for GPUMD, including energy, force and virial.

Parameters:

file_list (List[str] or str) – Single or multi cfg file.
type_dict (Dict[str, int]) – Map type from number to element, such as {0:’Al’, 1:’C’}.
output_name (str) – Output filename with append mode. Defaults to train.xyz.
f_max (float) – Force absolute maximum larger than this value will be filtered. Defaults to 25.0 eV/A.

mdapy.potential_tool.read_OUTCAR(filename: str) → Dict | bool[source]¶

This is only working for single-point calculation OUTCAR.

Args:: filename (str): outcar filename.
Returns:: Dict[str:Any]: system information, such as ‘Natom’, ‘lattice’, ‘energy’, ‘pos_force’, ‘symbols’, ‘virial’

mdapy.potential_tool.outcar2xyz(outcar_list: List[str] | str, output_path: str = 'train.xyz', mode: str = 'w', print_no_converge: bool = True)[source]¶

Convert OUTCAR file for VASP to xyz file for GPUMD, including energy, force and virial. It only works for single-point energy calculation generated OUTCAR.

Parameters:

outcar_list (List[str] or str) – Single or multi OUTCAR file.
output_path (str) – Output filename. Defaults to train.xyz.
mode (str) – Writting mode, such as ‘w’ and ‘a’. Defaults to ‘w’.
print_no_converge (bool) – Whether print non-converged outcars filename. Defaults to True

mdapy.potential_tool.outcars2xyz(outcar_list: List[str] | str, output_path: str = 'train.xyz', mode: str = 'w', print_no_converge: bool = True)[source]¶

Convert OUTCAR file for VASP to xyz file for GPUMD, including energy, force and virial.

Parameters:

outcar_list (List[str] or str) – Single or multi OUTCAR file.
output_path (str) – Output filename. Defaults to train.xyz.
mode (str) – Writting mode, such as ‘w’ and ‘a’. Defaults to ‘w’.
print_no_converge (bool) – Whether print non-converged outcars filename. Defaults to True

mdapy.qha_elastic module¶

mdapy.radial_distribution_function module¶

class mdapy.radial_distribution_function.RadialDistributionFunction(rc: float, nbin: int, box: Box, verlet_list: ndarray[tuple[Any, ...], dtype[int32]] | None = None, distance_list: ndarray[tuple[Any, ...], dtype[float64]] | None = None, neighbor_number: ndarray[tuple[Any, ...], dtype[int32]] | None = None, type_list: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None = None, streaming: bool = False, x: ndarray[tuple[Any, ...], dtype[float64]] | None = None, y: ndarray[tuple[Any, ...], dtype[float64]] | None = None, z: ndarray[tuple[Any, ...], dtype[float64]] | None = None)[source]¶

Bases: object

Radial distribution function \(g(r)\) of an atomic system.

For a multi-component system, the total RDF is the concentration-weighted sum of the partial RDFs:

\[g(r) = \sum_{\alpha \beta} c_{\alpha} c_{\beta} g_{\alpha \beta}(r)\]

where \(c_{\alpha}\) is the concentration of species \(\alpha\) and \(g_{\alpha \beta}(r) = g_{\beta \alpha}(r)\). The partial RDFs are normalised so that \(g_{\alpha \beta}(r) \to 1\) for an ideal homogeneous system.

Two execution paths are supported:

Verlet-list path (default): supply a pre-built neighbour list via verlet_list / distance_list / neighbor_number. Best when the cutoff is small relative to the box.
Streaming path (streaming=True): supply atomic positions via x / y / z. Pair distances are binned directly into the histogram with no neighbour list materialised, so memory stays \(\mathcal{O}(N)\) regardless of cutoff.

Parameters:

rc (float) – Cutoff distance for RDF calculation (in Angstroms).
nbin (int) – Number of histogram bins.
box (Box) – Simulation box.
verlet_list (NDArray[np.int32], optional) – Required when streaming=False.
distance_list (NDArray[np.float64], optional) – Required when streaming=False.
neighbor_number (NDArray[np.int32], optional) – Required when streaming=False.
type_list (NDArray, optional) – Per-atom species labels (int type ids or str element symbols). Used as keys in g_partial. If omitted, all atoms share label 0.
streaming (bool, default=False) – Use the streaming kernel (no Verlet list).
x (NDArray[np.float64], optional) – Atomic Cartesian coordinates. Required when streaming=True.
y (NDArray[np.float64], optional) – Atomic Cartesian coordinates. Required when streaming=True.
z (NDArray[np.float64], optional) – Atomic Cartesian coordinates. Required when streaming=True.

r¶

Bin centres, shape (nbin,).

Type:: NDArray[np.float64]

g_total¶

Total \(g(r)\), shape (nbin,).

Type:: NDArray[np.float64]

g_partial¶

Mapping (species_a, species_b) -> NDArray for \(a \leq b\) (upper triangle). Each value is the partial \(g_{\alpha\beta}(r)\).

Type:: dict

elements¶

Sorted list of unique species labels — the same labels used as keys in g_partial.

Type:: list

Ntype¶

Number of distinct species.

Type:: int

vol¶

Simulation box volume.

Type:: float

N¶

Total number of particles.

Type:: int

compute() → None[source]¶: Run the kernel and populate g_total, g_partial, and r.

plot(fig: 'Figure' | None = None, ax: 'Axes' | None = None) → Tuple[source]¶: Plot the total RDF \(g(r)\).

plot_partial(elements_list: List[str] | None = None, fig: 'Figure' | None = None, ax: 'Axes' | None = None) → Tuple[source]¶

Plot all partial RDFs \(g_{\alpha\beta}(r)\).

Parameters:: elements_list (list of str, optional) – Override the element labels used in the legend. Length must match Ntype. If None, uses the keys of g_partial directly.

mdapy.render module¶

mdapy.render.is_gpu_available() → bool[source]¶: Return True if this mdapy build supports GPU rendering (OptiX 7+).

mdapy.render.load_image(path: str) → ndarray[source]¶

Load an image file and return a (H, W, 4) uint8 RGBA numpy array.

Supports PNG, JPEG, BMP, TGA, GIF, PSD, HDR, PIC, PNM (via stb_image).

Parameters:: path (str) – Path to the image file.
Returns:: numpy.ndarray, shape (H, W, 4), dtype uint8, RGBA.

mdapy.render.save_image(path: str, img: ndarray) → None[source]¶

Save a (H, W, 4) uint8 RGBA image array to a file.

Format is determined by the file extension:

.png — RGBA, lossless (alpha preserved)
.jpg/.jpeg — RGB, lossy (alpha discarded, quality 95)
.bmp — RGB
.tga — RGBA
other — treated as .png

Parameters:

path (str) – Destination file path.
img (numpy.ndarray, shape (H, W, 4), dtype uint8) – RGBA image to save.

class mdapy.render.CameraParams(is_perspective: bool = True, field_of_view: float = 0.6981317007977318, position: Tuple[float, float, float] = (0.0, 0.0, 50.0), direction: Tuple[float, float, float] = (0.0, 0.0, -1.0), up: Tuple[float, float, float] = (0.0, 1.0, 0.0), znear: float = 0.0, dof_enabled: bool = False, dof_focal_len: float = 40.0, dof_aperture: float = 0.01)[source]¶

Bases: object

Camera parameters compatible with OVITO’s ViewProjectionParameters.

Perspective mode¶

field_of_view: vertical angle in radians. Tachyon zoom = 0.5 / tan(fov/2) — same formula as OVITO.

Orthographic mode¶

field_of_view: viewport half-height in world units. Tachyon zoom = 0.5 / field_of_view.

class mdapy.render.TachyonRender(backend: str = 'cpu', antialiasing: bool = True, aa_samples: int = 12, ao: bool = True, ao_samples: int = 12, ao_brightness: float = 0.8, ao_max_dist: float = 3.402823e+38, shadows: bool = True, direct_light_intensity: float = 0.9, background: tuple = (0.0, 0.0, 0.0))[source]¶

Bases: object

mdapy Tachyon ray-tracing renderer (Tachyon 0.99.5).

Supports two backends selectable via the backend parameter:

"cpu" (default): Classic multi-threaded CPU renderer. Always available. Uses Tachyon’s POSIX-thread parallelism.
"gpu": NVIDIA OptiX 7 GPU ray-tracing backend (RTX hardware acceleration). Only available when mdapy was compiled with CUDA + OptiX support. Raises RuntimeError at construction time if no CUDA GPU is found. Check availability with is_gpu_available().
"auto": Use GPU if available, fall back to CPU silently.

Parameters:

backend (str) – Rendering backend: "cpu", "gpu", or "auto". Default "cpu".
antialiasing (bool) – Enable anti-aliasing. Default True.
aa_samples (int) – Anti-aliasing samples per pixel. Default 12.
ao (bool) – Enable ambient occlusion. Default True.
ao_samples (int) – AO samples. Default 12.
ao_brightness (float) – Sky-light brightness for AO. Default 0.8.
ao_max_dist (float) – Maximum AO occlusion distance. Default unlimited.
shadows (bool) – Enable hard shadows (CPU only; GPU always uses shadows). Default True.
direct_light_intensity (float) – Directional light intensity. Default 0.9.
background (tuple) – Background colour (R, G, B) or (R, G, B, A) in [0, 1]. Default black.

Notes

Image size (width / height) is specified per-call in render() / render_system(), so the same renderer instance can produce images at different resolutions without re-initialisation.

Examples

>>> ren_cpu = TachyonRender()                     # CPU, default settings
>>> ren_gpu = TachyonRender(backend="gpu")         # GPU
>>> ren_auto = TachyonRender(backend="auto")       # auto-select
>>> # render at 1920×1080; save directly to PNG
>>> ren_cpu.render_system(sys, width=1920, height=1080,
...                       output_figure="out.png")
>>> # render with transparent background
>>> ren_cpu.render_system(sys, width=800, height=600,
...                       output_figure="out.png", transparent=True)
>>> print(ren_cpu.backend)
cpu

property backend: str¶

'cpu' or 'gpu'.

Type:: The active rendering backend

render(positions: ndarray, colors: ndarray, radii: ndarray, camera: CameraParams | None = None, bond_edges: ndarray | None = None, bond_colors: ndarray | None = None, bond_radius: float = 0.1, bond_color: tuple = (0.8, 0.8, 0.8, 1.0), box_edges: ndarray | None = None, box_edge_radius: float = 0.05, box_color: tuple = (1.0, 1.0, 1.0, 1.0), width: int = 800, height: int = 600, output_figure: str | None = None, transparent: bool = False) → ndarray | None[source]¶

Render spherical particles and optional bond / box cylinders.

Parameters:

positions ((N, 3) float64 Particle positions.)
colors ((N, 4) float32 Per-particle RGBA colour, values in [0, 1].)
radii ((N,) float32 Per-particle radius.)
camera (CameraParams or None. If None, a perspective camera) – is generated automatically from the bounding box.
bond_edges ((K, 2, 3) float64 or None. Pairs of endpoints for each) – bond cylinder. Pass None to skip drawing bonds.
bond_colors ((K, 4) float32 or None. Per-bond RGBA colours.) – If None, bond_color is used for all bond cylinders.
bond_radius (float. Cylinder radius for bonds. Default 0.1.)
bond_color (tuple (R, G, B) or (R, G, B, A), values in [0, 1].) – Default opaque grey.
box_edges ((M, 2, 3) float64 or None. Pairs of endpoints for each) – box edge cylinder. Pass None to skip drawing.
box_edge_radius (float. Cylinder radius for box edges. Default 0.05.)
box_color (tuple (R, G, B) or (R, G, B, A), values in [0, 1].) – Default opaque white (1, 1, 1, 1).
width (int. Output image size in pixels. Default 800 × 600.)
height (int. Output image size in pixels. Default 800 × 600.)
output_figure (str or None. If given, save the image to this path) – (format inferred from extension: png/jpg/bmp/tga) and return None. If None, return the RGBA numpy array.
transparent (bool. When True, pixels whose RGB matches the background) – colour are made fully transparent (alpha=0). Only meaningful for PNG output or when inspecting the array. Default False.

Returns:

numpy.ndarray, shape (H, W, 4), dtype uint8, RGBA image — or None if
output_figure is provided.

render_system(system: System, colors: ndarray | None = None, radii: ndarray | None = None, camera: CameraParams | None = None, draw_bond: bool = False, bond: ndarray | None = None, bond_radius: float = 0.1, bond_color: tuple = (0.8, 0.8, 0.8, 1.0), bond_color_mode: str = 'uniform', draw_box: bool = True, box_edge_radius: float = 0.05, box_color: tuple = (1.0, 1.0, 1.0, 1.0), default_radius: float = 1.0, width: int = 800, height: int = 600, output_figure: str | None = None, transparent: bool = False) → ndarray | None[source]¶

Render a mdapy.System object in one call.

Parameters:

system (mp.System The atomistic system to render.)
colors ((N, 4) float32 or None. If None, Jmol colours are) – assigned by element type.
radii ((N,) float32 or None. If None, default_radius is) – used for every particle.
camera (CameraParams or None. Auto-generated if not provided.)
draw_bond (bool. Whether to draw bond cylinders. Default False.)
bond ((Nbond, 2) int array or None. If None and) – draw_bond=True, system.bond is used.
bond_radius (float. Cylinder radius for bonds. Default 0.1.)
bond_color (tuple (R, G, B) or (R, G, B, A), values in [0, 1].) – Used when bond_color_mode='uniform'.
bond_color_mode (str. 'uniform' or 'atom'. In 'atom' mode) – each bond is split into two half-cylinders coloured by the connected atoms.
draw_box (bool. Whether to draw simulation-cell edges. Default True.)
box_edge_radius (float. Cylinder radius for cell edges. Default 0.05.)
box_color (tuple (R, G, B) or (R, G, B, A), values in [0, 1].) – Default opaque white (1, 1, 1, 1).
default_radius (float. Fallback radius when radii is None. Default 1.0.)
width (int. Output image size in pixels. Default 800 × 600.)
height (int. Output image size in pixels. Default 800 × 600.)
output_figure (str or None. If given, save image to this path and) – return None. Format inferred from extension.
transparent (bool. Transparent background (PNG recommended).)

mdapy.render.PRESET_VIEWS = ('perspective', 'orthographic', 'top', 'bottom', 'front', 'back', 'left', 'right')¶: All supported preset view names.

mdapy.render.preset_camera(view: str, positions: ndarray, fov_deg: float = 40.0, margin: float = 1.0, max_radius: float = 0.0) → CameraParams[source]¶

Build a camera that matches one of OVITO’s four default viewport orientations.

OVITO coordinate convention¶

OVITO uses a right-handed coordinate system where Z is the global up axis. The eight standard viewports are:

Name	Meaning	camera_dir (OVITO)	up
`"top"`	Look down the -Z axis → sees XY plane	(0, 0, -1)	(0,1,0)
`"front"`	Look along +Y axis → sees XZ plane	(0, 1, 0)	(0,0,1)
`"left"`	Look along +X axis → sees YZ plane	(1, 0, 0)	(0,0,1)
"perspective"``\| Tilted view from +X+Y+Z / ``"ortho" \| → 3D isometric look		(-1,-1,-1)/√3	(0,0,1)

All eight available view names¶

"perspective" – Perspective projection, tilted isometric camera
(matches OVITO Perspective viewport)
"orthographic" – Same direction as perspective but parallel projection
(matches OVITO Ortho viewport)
"top" – Orthographic, camera_dir=(0,0,-1), up=(0,1,0)
looks at the XY plane (matches OVITO Top viewport)
"bottom" – Orthographic, camera_dir=(0,0,+1), up=(0,1,0)
"front" – Orthographic, camera_dir=(0,+1,0), up=(0,0,1)
looks at the XZ plane, X points right, Z points up (matches OVITO Front viewport)
"back" – Orthographic, camera_dir=(0,-1,0), up=(0,0,1)
"left" – Orthographic, camera_dir=(+1,0,0), up=(0,0,1)
looks at the YZ plane, Y points right, Z points up (matches OVITO Left viewport)
"right" – Orthographic, camera_dir=(-1,0,0), up=(0,0,1)

param view:: View name. Must be one of PRESET_VIEWS.
type view:: str
param positions:: Particle positions in world coordinates.
type positions:: (N, 3) array_like
param fov_deg:: Vertical field of view in degrees for the perspective camera. Default is 40°.
type fov_deg:: float, optional
param margin:: Extra padding around the structure in world units (Å). Applied to all views. Default is 1.0 Å.
type margin:: float, optional
param max_radius:: Largest particle radius. Pass radii.max() so that edge atoms are not clipped. Default 0 (uses only atom centres).
type max_radius:: float, optional
returns:: CameraParams

Notes

About CameraParams.znear (orthographic cross-section clipping)

znear defaults to 0 and you normally do not need to change it. Setting it to a positive value shifts the camera plane forward along the view direction by that many world units, effectively clipping away the front portion of the structure. This is useful for cross-section renders:

cam = preset_camera("top", pos, max_radius=1.3)
cam.znear = 9.0  # discard atoms in front of the z=9 Å plane

mdapy.spatial_binning module¶

class mdapy.spatial_binning.SpatialBinning(data: DataFrame, box: Box, direction: str, bin_width: float, origin: List[float] | None = None)[source]¶

Bases: object

Spatial binning of atomic data in a simulation box.

This class divides a simulation box into spatial bins along specified directions and computes aggregated statistics for atomic properties within each bin.

Parameters:

data (pl.DataFrame) – DataFrame containing atomic data with at least ‘x’, ‘y’, ‘z’ columns.
box (Box) – Simulation box object containing box dimensions.
direction (str) – Binning direction(s). Supported values: ‘x’, ‘y’, ‘z’, ‘xy’, ‘xz’, ‘yz’, ‘xyz’.
bin_width (float) – Width of each bin in the same units as box dimensions.
origin (list of float, optional) – Origin point for binning [x, y, z]. If None, uses box.origin.

bin_volume¶

Volume of a single bin.

Type:: float

Notes

Currently only supports orthogonal (non-triclinic) boxes.

compute(name: str | List[str], operation: str | List[str], fill_empty: bool = False) → DataFrame[source]¶

Compute aggregated statistics for spatial bins.

Parameters:

name (str or list of str) – Column name(s) to aggregate.
operation (str or list of str) – Aggregation operation(s) to perform. Supported operations: ‘mean’, ‘sum’, ‘min’, ‘max’, ‘sum/binvol’, ‘count’.

Returns:

pl.DataFrame – Aggregated results with bin coordinates and statistics.

mdapy.spline module¶

Bases: object

Cubic spline interpolation on a strictly-increasing grid.

Constructs a piecewise-cubic \(s(x)\) that is \(C^2\) over the whole range, reproduces the sample points \((x_i, y_i)\) exactly, and satisfies the chosen boundary condition at the two endpoints. The grid need not be uniform — if it is, prefer the internal UniformCubicSpline (used by the EAM code, not exposed here) which is O(1) per lookup.

Parameters:

x (array_like) – 1-D array of x-coordinates. Must be strictly increasing and contain at least two points.
y (array_like) – 1-D array of y-coordinates, same length as x.
bc_type ({"not-a-knot", "natural", "clamped"}, default "not-a-knot") –
Boundary condition at the two endpoints:
- "not-a-knot" — the third derivative is continuous at x[1] and x[n-2] (equivalently, the first two and last two cubic pieces are each a single polynomial). Same default as scipy.interpolate.CubicSpline. Best for general data when the endpoint slopes are unknown.
- "natural" — \(s''(x_0) = s''(x_{n-1}) = 0\). Produces a minimum-curvature interpolant that flattens out at the ends.
- "clamped" — \(s'(x_0) = \texttt{dy0}\), \(s'(x_{n-1}) = \texttt{dyn}\). If dy0 and dyn are not given, they are estimated by fitting a quadratic through the first (last) three points and taking its analytic derivative at the endpoint.
dy0 (float, optional) – Endpoint first derivatives, only used when bc_type="clamped". Both must be provided together; if either is None the three-point estimates are used.
dyn (float, optional) – Endpoint first derivatives, only used when bc_type="clamped". Both must be provided together; if either is None the three-point estimates are used.

Notes

Evaluation is O(log n) per point via binary search. Batch evaluation is OpenMP-parallelised.

Out-of-range queries raise IndexError for scalar calls and return NaN element-wise for array calls. There is deliberately no silent extrapolation — cubic extrapolation past the last knot can swing wildly on smooth-looking data (see the EAM rho-clamping incident for a live example).

Examples

>>> import numpy as np
>>> x = np.linspace(0, 2 * np.pi, 13)
>>> y = np.sin(x)
>>> sp = Spline(x, y)                 # default: not-a-knot
>>> abs(sp.evaluate(np.pi / 4) - np.sin(np.pi / 4)) < 1e-4
True
>>> sp.derivative(0.0)                # should be ~cos(0) = 1
1.0000...

A clamped spline with user-supplied endpoint slopes — useful when you know the analytic derivative at the ends (here we know \(\cos(0) = 1\) and \(\cos(2\pi) = 1\)):

>>> sp_c = Spline(x, y, bc_type="clamped", dy0=1.0, dyn=1.0)

The natural spline, by contrast, forces \(s'' = 0\) at the ends, which is appropriate when you expect the data to flatten beyond the sample range.

Evaluate \(s(x)\) at scalar or array x.

Array inputs return an np.ndarray of the same length; entries outside the interpolation range become NaN. Scalar inputs raise IndexError if out of range.

Evaluate \(s'(x)\) at scalar or array x.

The derivative is computed analytically from the stored cubic coefficients, not by finite differencing.

Evaluate \(s''(x)\) at scalar or array x.

\(s''\) is piecewise-linear between the knots (a property of cubic splines), so this is exact up to floating-point rounding.

mdapy.sqs module¶

Special Quasirandom Structure (SQS) generation, modeled after ATAT mcsqs.

Take an mdapy.System containing a random alloy (e.g. produced by mdapy.build_hea()), enumerate pair / triplet / quadruplet clusters from the neighbor list, and run a Monte-Carlo swap engine to drive the cluster correlations to their fully-random target — i.e. an SQS. All cell shapes are supported (cubic, orthorhombic, triclinic) and small boxes that mdapy.Neighbor replicates internally are handled by folding ghost neighbor indices back to the primary cell.

Quick recipe¶

>>> import mdapy as mp
>>> sys_init = mp.build_hea(
...     ("Cr", "Co", "Ni"), (1/3,) * 3,
...     "fcc", 3.53, nx=3, ny=3, nz=3, random_seed=1,
... )
>>> sqs = mp.SQS(sys_init, cutoffs={2: 4.0, 3: 3.0}).compute()
>>> sqs.is_sqs()                      # prints a one-screen summary
>>> sqs.system.write_data("sqs.data") # write the resulting structure

How to choose `cutoffs`¶

cutoffs is a dict keyed by cluster body order (2 = pair, 3 = triplet, 4 = quadruplet). Only 2 is required; 3 and 4 are optional and add longer-range matching at the cost of more channels.

The right value is “just past the shell you want to constrain”. The first two coordination shells (NN1, NN2) for the common metallic lattices, with a the conventional lattice parameter:

fcc: NN1 = a / sqrt(2) ≈ 0.707 a, NN2 = a
bcc: NN1 = sqrt(3)/2 * a ≈ 0.866 a, NN2 = a
hcp (ideal c/a): NN1 = a, NN2 = sqrt(2) a

Practical choice — covering NN1+NN2 with pairs and NN1 alone with triplets is a good default:

Lattice	`cutoffs[2]`	`cutoffs[3]` (opt.)
fcc, `a=3.53`	`4.0` (NN1+NN2)	`3.0` (NN1 only)
bcc, `a=2.87`	`3.5` (NN1+NN2)	`2.7` (NN1 only)
hcp, `a=2.95`	`4.5` (NN1+NN2)	`3.1` (NN1 only)

If you don’t know a: build a quick neighbor list, plot the RDF, and pick cutoffs[2] halfway between the second and third RDF peaks.

What the parameters mean¶

cutoffs (required) — see above.
n_replicas — number of independent MC chains run in parallel; the one with the lowest objective wins. More replicas help when the landscape has multiple minima. Default 4.
max_steps — swap attempts per chain. Each chain remembers its best-ever state, so increasing max_steps is monotonically helpful (it cannot make the result worse). Set 0 to evaluate the input as-is without running MC — useful for diagnostics. Default 100000.
T — Metropolis temperature. Smaller is greedier; 0.01–0.1 work well for equiatomic HEAs. Default 0.05.
seed — RNG seed for reproducibility.

Mathematical background¶

We use the trigonometric per-site basis of van de Walle (CALPHAD 42, 13–18, 2013): for an m-component sublattice each species s maps to m-1 basis values phi_k(s). A channel groups cluster instances that share both a distance shell and the same function tuple (f_1, ..., f_n). The channel correlation is

pi = (1/N_inst) sum_inst prod_p phi_{f_p}(s_{a_p})

with target (under fully-random species placement) prod_p < phi_{f_p} >. For an equiatomic alloy every multi-body target is exactly zero.

The objective implemented is ATAT mcsqs’s d1-perfect-match formula: it combines a weighted residual sum with a reward for extending the longest range over which all correlations are matched within a small tolerance.

Citation¶

If you use this in a publication, please cite:

A. van de Walle, P. Tiwary, M. de Jong, D. L. Olmsted, M. Asta, A. Dick, D. Shin, Y. Wang, L.-Q. Chen, Z.-K. Liu. Efficient stochastic generation of special quasirandom structures. CALPHAD 42, 13–18 (2013).

class mdapy.sqs.SQS(system: System, cutoffs: Dict[int, float], n_replicas: int = 4, max_steps: int = 100000, T: float = 0.05, seed: int = 0)[source]¶

Bases: object

Generate a Special Quasirandom Structure from a random alloy template.

Parameters:

system (mdapy.System) – Input alloy. Must contain an element (or type) column with the same number of distinct species as the desired SQS. Atom positions and cell are taken as-is — only species labels are reshuffled.
cutoffs (dict[int, float]) – Cluster cutoffs in Å, keyed by cluster size: 2 (required) is the pair cutoff; 3 and 4 (optional) add triplet / quadruplet constraints. See the module docstring for how to choose them per crystal structure.
n_replicas (int, default=4) – Number of independent Monte-Carlo chains run in parallel via OpenMP.
max_steps (int, default=100000) – Swap attempts per chain. 0 evaluates the input only (no MC). The engine tracks the best-ever state per chain, so larger values are monotonically helpful.
T (float, default=0.05) – Metropolis temperature.
seed (int, default=0) – RNG seed; each replica gets seed plus a per-replica offset.

system¶

Output SQS structure (set by compute()). Same atoms / cell as the input but with the element / type columns reshuffled.

Type:: mdapy.System

objective¶

Final value of the Monte-Carlo objective (lower is better).

Type:: float

correlations¶

Per-channel correlations of the SQS, shape (n_channels,).

Type:: numpy.ndarray

channel_info¶

One entry per channel: n_pts, shell, diameter, funcs, n_instances, target, corr.

Type:: list[dict]

compute() → SQS[source]¶: Run the SQS optimization. Returns self for chaining.

is_sqs(tol: float = 0.03, verbose: bool = True)[source]¶

Verify whether self.system is a good SQS.

The verdict follows the formal SQS definition (Zunger 1990; van de Walle ATAT 2013): an SQS is a configuration whose cluster correlation residuals all sit below tolerance,

max_alpha |pi_alpha - pi_target_alpha| < tol

across every cluster channel (pair / triplet / quad) the optimiser was given. tol=0.03 is fine for ~100-atom equiatomic cells; tighten as N grows.

For convenience this method also computes Warren-Cowley max |alpha_ij| (full matrix, including the diagonal) at every pair shell within cutoffs[2] and lists them in the report. WCP is a derived physical quantity — it’s a linear transform of the pair correlations — and is not part of the verdict, mirroring what ATAT mcsqs does. It’s there so the user can inspect the SRO, especially in non-equiatomic / small-cell cases where the cell may not have enough degrees of freedom to drive every alpha to zero even when the correlation residuals all pass.

Parameters:

tol (float, default=0.03) – Threshold for the absolute-residual verdict.
verbose (bool, default=True) – Print a one-screen summary report.

Returns:

passed (bool) – True iff max |Delta pi| < tol.
info (dict) – Diagnostic numbers behind the verdict.

mdapy.steinhardt_bond_orientation module¶

class mdapy.steinhardt_bond_orientation.SteinhardtBondOrientation(box: Box, data: DataFrame, llist: ndarray, nnn: int, rc: float, average: bool, use_voronoi: bool, use_weight: bool, weight: ndarray, verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray, wl: bool, wlhat: bool, identify_liquid: bool, threshold: float, n_bond: int)[source]¶

Bases: object

Compute Steinhardt bond orientation order parameters. The Steinhardt bond orientation order parameters provide a way to characterize local structural order in molecular dynamics simulations. This class computes the spherical harmonics-based order parameters \(q_l\) and \(w_l\) for identifying crystal structures and distinguishing between solid and liquid phases.

Parameters:

box (Box) – Simulation box object containing boundary conditions and dimensions.
data (pl.DataFrame) – Polars DataFrame containing particle positions with columns ‘x’, ‘y’, ‘z’.
llist (np.ndarray) – Array of spherical harmonic degrees \(l\) to compute (e.g., [4, 6, 8]).
nnn (int) – Number of nearest neighbors to consider (0 to use cutoff distance instead).
rc (float) – Cutoff radius for neighbor search (ignored if nnn > 0 or use_voronoi=True).
average (bool) – Whether to compute averaged order parameters \(\bar{q}_l\) over neighbor shells.
use_voronoi (bool) – Use Voronoi tessellation to determine neighbors.
use_weight (bool) – Apply weights to neighbor contributions.
weight (np.ndarray) – Weight array for each neighbor pair (must match verlet_list shape if use_weight=True).
verlet_list (np.ndarray) – Precomputed neighbor list array.
distance_list (np.ndarray) – Precomputed distances for neighbor pairs.
neighbor_number (np.ndarray) – Number of neighbors for each particle.
wl (bool) – Compute third-order invariant \(w_l\) parameters.
wlhat (bool) – Compute normalized third-order invariant \(\hat{w}_l\) parameters.
identify_liquid (bool) – Enable solid-liquid classification (requires \(l=6\) in llist).
threshold (float) – Threshold for solid-liquid identification (default: 0.7 for normalized \(q_6\)).
n_bond (int) – Minimum number of “solid-like” bonds for solid classification.

qlm_r¶

Real part of spherical harmonics \(q_{lm}\), shape (\(N_{particles}\), \(N_l\), \(2 l_{max} + 1\)).

Type:: np.ndarray

qlm_i¶

Imaginary part of spherical harmonics \(q_{lm}\), shape (\(N_{particles}\), \(N_l\), \(2 l_{max} + 1\)).

Type:: np.ndarray

qnarray¶

Computed order parameters, shape (\(N_{particles}\) \(N_{columns}\)) where columns contain \(q_l\), and optionally \(w_l\) and \(\hat{w}_l\) values.

Type:: np.ndarray

solidliquid¶

Solid-liquid classification array (0=liquid, 1=solid), only computed if identify_liquid=True.

Type:: np.ndarray, optional

nbond¶

Number of solid-like bonds for each particle, only computed if identify_liquid=True.

Type:: np.ndarray, optional

Notes

For a particle \(i\), we calculate the quantity \(q_{lm}\) by summing the spherical harmonics between particle \(i\) and its neighbors \(j\) in a local region:

\[q_{lm}(i) = \frac{1}{N_b} \sum \limits_{j=1}^{N_b} Y_{lm}(\theta(\vec{r}_{ij}), \phi(\vec{r}_{ij}))\]

Then the \(q_l\) order parameter is computed by combining the \(q_{lm}\) in a rotationally invariant fashion to remove local orientational order:

\[q_l(i) = \sqrt{\frac{4\pi}{2l+1} \sum \limits_{m=-l}^{l} |q_{lm}(i)|^2 }\]

If the wl parameter is True, this class computes the quantity \(w_l\), defined as a weighted average over the \(q_{lm}(i)\) values using Wigner 3-j symbols (related to Clebsch-Gordan coefficients). The resulting combination is rotationally invariant:

\[\begin{split}w_l(i) = \sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix} l & l & l \\ m_1 & m_2 & m_3 \end{pmatrix} q_{lm_1}(i) q_{lm_2}(i) q_{lm_3}(i)\end{split}\]

If wlhat parameter to True will normalize the \(w_l\) order parameter as follows:

\[\begin{split}w_l(i) = \frac{ \sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix} l & l & l \\ m_1 & m_2 & m_3 \end{pmatrix} q_{lm_1}(i) q_{lm_2}(i) q_{lm_3}(i)} {\left(\sum \limits_{m=-l}^{l} |q_{lm}(i)|^2 \right)^{3/2}}\end{split}\]

If average is True, the class computes a variant of this order parameter that performs an average over the first and second shell combined . To compute this parameter, we perform a second averaging over the first neighbor shell of the particle to implicitly include information about the second neighbor shell. This averaging is performed by replacing the value \(q_{lm}(i)\) in the original definition by \(\overline{q}_{lm}(i)\), the average value of \(q_{lm}(k)\) over all the \(N_b\) neighbors \(k\) of particle \(i\), including particle \(i\) itself:

\[\overline{q}_{lm}(i) = \frac{1}{N_b} \sum \limits_{k=0}^{N_b} q_{lm}(k)\]

If use_weight is True, the contributions of each neighbor are weighted. Neighbor weights \(w_{ij}\) are defined for a Vornoi face area obtained from Voronoi neighbor or one with user-provided weights, and default to 1 if not otherwise provided. The formulas are modified as follows, replacing \(q_{lm}(i)\) with the weighted value \(q'_{lm}(i)\):

\[q'_{lm}(i) = \frac{1}{\sum_{j=1}^{N_b} w_{ij}} \sum \limits_{j=1}^{N_b} w_{ij} Y_{lm}(\theta(\vec{r}_{ij}), \phi(\vec{r}_{ij}))\]

For solid-liquid classification (only performed when \(l=6\) is included), particles are identified as solid if they have at least \(n_{\text{bond}}\) “solid-like” bonds. A bond between particles \(i\) and \(j\) is considered solid-like if their local bond orientations are correlated:

\[\mathrm{dot}_{6}(i,j) = \frac{\sum_{m=-6}^{6} q_{6m}(i) \cdot q_{6m}^*(j)}{\sqrt{\sum_{m=-6}^{6} |q_{6m}(i)|^2} \sqrt{\sum_{m=-6}^{6} |q_{6m}(j)|^2}} > \text{threshold}\]

where the default threshold is 0.7 for normalized \(q_6\) values. When solid bond is larger than n_bond, the atom is treated as solid atom. The default n_bond is 7, where 6-8 is generally good for FCC and BCC strutures.

References

Paul J. Steinhardt. Bond-orientational order in liquids and glasses. Physical Review B, 28(2):784-805, 1983. doi:10.1103/PhysRevB.28.784.

Wolfgang Lechner and Christoph Dellago. Accurate determination of crystal structures based on averaged local bond order parameters. The Journal of Chemical Physics, 129(11):114707, 2008. http://dx.doi.org/10.1063/1.2977970

Filion, M. Hermes, R. Ni, and M. Dijkstra. Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: a comparison of simulation techniques. The Journal of Chemical Physics, 133(24):244115, http://dx.doi.org/10.1063/1.3506838.

compute() → None[source]¶

Compute the Steinhardt bond orientation order parameters.

This method calculates the bond orientation order parameters and stores the results as instance attributes. If identify_liquid is True, it also performs solid-liquid classification.

Returns:

None – Results are stored in the following instance attributes:

qlm_r : Real part of spherical harmonics \(q_{lm}\)
qlm_i : Imaginary part of spherical harmonics \(q_{lm}\)
qnarray : Computed order parameters array
solidliquid : Solid-liquid classification (if identify_liquid=True)
nbond : Number of solid-like bonds per particle (if identify_liquid=True)

Raises:

AssertionError – If identify_liquid is True but l=6 is not in llist.
AssertionError – If identify_liquid is True but threshold or n_bond are not positive.
AssertionError – If rc is not positive when using cutoff-based neighbor search.
AssertionError – If use_weight is True but weight array shape doesn’t match verlet_list.

mdapy.structure_entropy module¶

This module implements the calculation of local structural entropy based on the method proposed by Piaggi and Parrinello (J. Chem. Phys. 147, 114112 (2017)), which quantifies the degree of local ordering in atomic configurations.

The local structural entropy measures how similar the local radial distribution around each atom is to the average distribution of its surroundings. It is a useful descriptor for identifying crystalline and disordered regions, phase transitions, or defects.

The implementation supports both global and density-weighted (normalized) entropy evaluation, and an optional neighbor-based averaging of the entropy field.

References

P.M. Piaggi and M. Parrinello, Entropy based fingerprint for local crystalline order, J. Chem. Phys. 147, 114112 (2017) https://doi.org/10.1063/1.4998408

class mdapy.structure_entropy.StructureEntropy(box: Box, verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray, rc: float, sigma: float, use_local_density: bool, average_rc: float = 0.0)[source]¶

Bases: object

Calculate the local structural entropy for each atom based on its neighbor environment within a specified cutoff radius.

The local structural entropy quantifies the degree of local disorder by comparing the local pair distribution to a reference distribution. Lower entropy values correspond to more ordered (crystalline) environments, while higher values indicate disordered or amorphous regions.

Parameters:

box (Box) – Simulation box object defining the system boundaries.
verlet_list (np.ndarray) – Neighbor list of atomic indices. Each row corresponds to one atom and contains indices of its neighbors.
distance_list (np.ndarray) – Distance list corresponding to the verlet_list. Each element stores the pairwise distances between atoms.
neighbor_number (np.ndarray) – Number of valid neighbors for each atom.
rc (float) – Cutoff radius (in Å) used for calculating the local pair distribution function.
sigma (float) – Gaussian smoothing width used in constructing the local radial distribution function.
use_local_density (bool) – If True, the local RDF is normalized by the local atomic density. If False, a global normalization is used.
average_rc (float, optional) – Cutoff radius for averaging entropy values over neighboring atoms. If zero or not specified, no averaging is performed.

entropy¶

Array of local structural entropy values for all atoms.

Type:: np.ndarray

entropy_ave¶

Averaged entropy values over neighbors within average_rc. Only computed if average_rc > 0.

Type:: np.ndarray, optional

Notes

The algorithm internally uses a precomputed neighbor list (verlet_list, distance_list, and neighbor_number).
The entropy field can be spatially averaged to reduce statistical noise or highlight larger-scale structural features.

compute()[source]¶

Perform the local structural entropy calculation.

This method computes the entropy value for each atom based on its local pair distribution function within rc. Optionally, if average_rc > 0, the resulting entropy field is averaged over neighbors within the given radius.

Returns:: None – Results are stored in the instance attributes: entropy and, if applicable, entropy_ave.

mdapy.structure_factor module¶

class mdapy.structure_factor.StructureFactor(data: DataFrame, box: Box, k_min: float, k_max: float, nbins: int, cal_partial: bool = False, atomic_form_factors: bool = False, mode: Literal['direct', 'debye', 'rdf'] = 'debye', rc: float | None = None, nbin_rdf: int = 200, window: bool = False)[source]¶

Bases: object

Static structure factor \(S(k)\) of an atomic system.

The class computes the partial Faber-Ziman structure factors \(A_{\alpha\beta}(k)\), the total \(S(k)\), and (optionally) their X-ray, neutron, and electron-diffraction equivalents along with the real-space derived quantities \(g(r)\), \(G(r)\), and \(R(r)\).

Two algorithms are available:

Debye method (mode='debye', default). The radial distribution function \(g_{\alpha\beta}(r)\) is computed once via the streaming RDF kernel up to rc (default \(L_{\max}/2\), where \(L_{\max}\) is the longest cell-vector length). Each partial is then a one-dimensional sin transform, integrated by numpy.trapezoid() on the bin-centre grid:

\[A_{\alpha\beta}(k) = 1 + \frac{4\pi\rho}{k} \int_0^{r_c} r\,[g_{\alpha\beta}(r) - 1]\,\sin(kr)\,w(r)\,dr,\]

with \(w(r) = \mathrm{sinc}(2\pi r/L_{\max})\) if window=True and \(w(r) = 1\) otherwise. The total \(S(k) = \sum_{\alpha\beta} c_\alpha c_\beta A_{\alpha\beta}(k)\). This implementation is bit-for-bit identical (max diff \(\lesssim 10^{-10}\)) to the reference code accompanying [Erhard2024].
Direct method (mode='direct'). Generates every reciprocal lattice vector \(\vec{k}\) in \((k_{\min}, k_{\max}]\) and sums

\[F_\alpha(\vec{k}) = \frac{1}{\sqrt{N}} \sum_{i\in\alpha} e^{i \vec{k}\cdot\vec{r}_i}, \qquad A_{\alpha\beta}(\vec{k}) = \mathrm{Re}\left[ F_\alpha^{*}(\vec{k})\,F_\beta(\vec{k})\right] / (c_\alpha c_\beta).\]

Spherically averaging \(A_{\alpha\beta}(\vec{k})\) in \(|\vec{k}|\)-bins yields \(A_{\alpha\beta}(k)\) of the same Faber-Ziman normalisation as the Debye method, modulo statistical noise from the discrete reciprocal-lattice grid. The direct method is much slower than the Debye method for large systems but resolves Bragg peaks in crystals where the Debye \(g(r)\)-based method is blurred by the finite radial cutoff.

The X-ray, neutron, and electron total structure factors are constructed from the partials following Erhard et al.:

\[S^{(X)}(k) = \frac{ \sum_{\alpha\beta} (2 - \delta_{\alpha\beta}) c_\alpha c_\beta f_\alpha(k) f_\beta(k)\,A_{\alpha\beta}(k) }{ \left[\sum_\alpha c_\alpha f_\alpha(k)\right]^2 },\]

with the appropriate form factors \(f(k)\):

X-ray: nine-coefficient Cromer-Mann fit \(f(k) = c + \sum_{i=1}^4 a_i\,\exp(-b_i\,(k/4\pi)^2)\) from [BrownITC].
Neutron: tabulated coherent scattering length from NIST, \(k\)-independent.
Electron: Mott-Bethe relation \(f_e(k) = (Z - f_X(k)) / (8 \pi^2 a_0 k^2)\), with \(a_0\) the Bohr radius in Angstrom.

Parameters:

data (pl.DataFrame) – Atomic data with at least x, y, z columns. For \(A_{\alpha\beta}\) and X-ray/neutron/electron derivations, element (preferred) or type is required.
box (Box) – Simulation box.
k_min (float) – Wavenumber range.
k_max (float) – Wavenumber range.
nbins (int) – Number of \(k\)-bins.
cal_partial (bool, default False) – Compute \(A_{\alpha\beta}(k)\). Required for X-ray / neutron / electron totals.
atomic_form_factors (bool, default False) – Compute X-ray total \(S^{(X)}(k)\). Implies cal_partial.
mode ({'debye', 'direct'}, default 'debye') – Algorithm selector. 'rdf' is accepted as an alias for 'debye'.
rc (float, optional) – Radial cutoff for the Debye-mode RDF. Defaults to \(L_{\max}/2\). Only used by mode='debye'.
nbin_rdf (int, default 200) – Number of radial bins for the Debye-mode RDF.
window (bool, default False) – Apply the Lorch window \(w(r) = \mathrm{sinc}(2\pi r/L_{\max})\) in the Debye-mode integral. Default off, matching the textbook / Erhard reference.

k¶

Wavenumber bin centres.

Type:: np.ndarray, shape (nbins,)

Sk¶

Total structure factor.

Type:: np.ndarray, shape (nbins,)

Sk_partial¶

Maps (species_a, species_b) -> A_alpha_beta(k) for the upper triangle (\(\alpha \leq \beta\)). Populated when cal_partial=True.

Type:: dict, optional

Sk_xray, Sk_neutron, Sk_electron

Total \(S^{(X)}\), \(S^{(N)}\), \(S^{(e)}\). Populated on demand by get_xray_structure_factor() etc., and automatically when atomic_form_factors=True.

Type:: np.ndarray, optional

References

[Erhard2024]

Erhard, L. C., Rohrer, J., Albe, K. & Deringer, V. L. “Modelling atomic and nanoscale structure in the silicon-oxygen system through active machine learning.” Nature Communications 15, 1927 (2024). https://doi.org/10.1038/s41467-024-45840-9. Reference Python implementation: https://github.com/LinusErhard/SiO_active_learning.

[BrownITC] (1,2)

Brown, P. J., Fox, A. G., Maslen, E. N., O’Keefe, M. A. & Willis, B. T. M. “Intensity of diffraction intensities,” in International Tables for Crystallography, Volume C: Mathematical, Physical, and Chemical Tables, Table 6.1.1.4, 4th ed. (2004).

Notes

The Faber-Ziman convention used here makes every \(A_{\alpha\beta}(k) \to 1\) at large \(k\). To convert to the OVITO structure-factor (Ashcroft-Langreth-style) partials use

\[S_{\alpha\beta}^{\mathrm{AL}}(k) = c_\alpha\,\delta_{\alpha\beta} + c_\alpha c_\beta\,[A_{\alpha\beta}(k) - 1].\]

compute() → None[source]¶

Run the configured algorithm and populate Sk and Sk_partial.

When atomic_form_factors=True this also populates Sk_xray. Neutron and electron totals stay None until the corresponding getters are called.

get_xray_structure_factor() → ndarray[source]¶

Return the X-ray total \(S^{(X)}(k)\).

Implements the Erhard convention

\[S^{(X)}(k) = \frac{ \sum_{\alpha\beta} (2 - \delta_{\alpha\beta}) c_\alpha c_\beta f_\alpha(k) f_\beta(k)\,A_{\alpha\beta}(k) }{\left[\sum_\alpha c_\alpha f_\alpha(k)\right]^2}\]

with \(f(k)\) the Cromer-Mann fit \(c + \sum_i a_i e^{-b_i (k/4\pi)^2}\) from [BrownITC].

get_neutron_structure_factor() → ndarray[source]¶

Return the neutron total \(S^{(N)}(k)\).

Form factors are the tabulated coherent scattering lengths from NIST (see mdapy.data.neutron_form_factor) and are \(k\)-independent. For absorptive isotopes (e.g. \(^{10}\mathrm{B}\), \(\mathrm{Cd}\), \(\mathrm{Sm}\)) the tabulated value is complex; the returned \(S^{(N)}(k)\) is then the modulus of the weighted sum.

get_electron_structure_factor() → ndarray[source]¶

Return the electron-diffraction total \(S^{(e)}(k)\).

Form factors follow the Mott-Bethe relation

\[f_e(k) = \frac{Z - f_X(k)}{8 \pi^2 a_0\,k^2}\]

with \(a_0\) the Bohr radius (in \(\AA\)) and \(f_X(k)\) the X-ray form factor.

get_xray_pair_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, g_X(r))\) reconstructed from \(S^{(X)}(k)\).

get_xray_reduced_pair_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, G_X(r) = 4\pi r \rho [g_X(r) - 1])\).

get_xray_radial_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, R_X(r) = 4\pi r^2 \rho g_X(r))\).

get_neutron_pair_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, g_N(r))\).

get_neutron_reduced_pair_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, G_N(r))\).

get_neutron_radial_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, R_N(r))\).

get_electron_pair_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, g_e(r))\).

get_electron_reduced_pair_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, G_e(r))\).

get_electron_radial_distribution_function(r: ndarray | None = None) → Tuple[ndarray, ndarray][source]¶: Return \((r, R_e(r))\).

plot(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶: Plot the total \(S(k)\) (or \(S^{(X)}(k)\) if atomic_form_factors=True).

plot_partial(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes][source]¶: Plot the partial Faber-Ziman \(A_{\alpha\beta}(k)\).

mdapy.system module¶

Core System Class for MDAPY¶

This module provides the central System class that serves as the main interface for atomic simulation data and analysis in mdapy. The System class integrates all analysis modules and provides a unified API for structural analysis, property calculations, and data manipulation.

Classes¶

System : Core class representing an atomic system

Bases: object

Core class representing an atomic system with integrated analysis capabilities.

The System class is the central data structure in mdapy, providing a unified interface for loading, manipulating, and analyzing atomic configurations. It integrates various structural analysis methods, property calculators, and visualization tools.

Parameters:

filename (str, optional) – Path to an atomic configuration file. Supports various formats including POSCAR, LAMMPS data/dump, XYZ, mp etc. Format is auto-detected or can be specified via the format parameter.
data (pl.DataFrame, optional) – Polars DataFrame containing atomic information. Must include columns ‘x’, ‘y’, ‘z’ for positions. Other common columns include ‘type’, ‘element’, ‘vx’, ‘vy’, ‘vz’, ‘fx’, ‘fy’, ‘fz’, etc.
pos (np.ndarray, optional) – Array of shape (N, 3) containing atomic positions. Used in conjunction with box parameter.
box (int, float, Iterable[float], np.ndarray, or Box, optional) – Simulation box specification. Required when using data or pos parameters. See Box for supported formats.
ase_atom (ase.Atoms, optional) – ASE Atoms object to initialize the system from.
ovito_atom (ovito.data.DataCollection, optional) – OVITO DataCollection object to initialize the system from.
format (str, optional) – File format specification when loading from filename. If None, format is auto-detected.
global_info (dict, optional) – Dictionary containing global system properties such as total energy, virial, stress tensor, timestep, etc. Defaults to empty dict.

verlet_list¶

Neighbor list array when neighbors are built.

Type:: np.ndarray, optional

distance_list¶

Distance list array when neighbors are built.

Type:: np.ndarray, optional

neighbor_number¶

Number of neighbors for each atom when neighbors are built.

Type:: np.ndarray, optional

rc¶

Cutoff radius used for neighbor list construction.

Type:: float, optional

Notes

The System class uses lazy evaluation for many analyses - neighbor lists and other expensive computations are only performed when needed. Results are cached and reused when possible.

At initialization, exactly one of the following must be provided: - filename - data and box - pos and box - ase_atom - ovito_atom

Examples

Creating a system from different sources:

from mdapy.system import System
import numpy as np

# Load from file (auto-detect format)
system = System("config.dump")

# Load from file with explicit format
system = System("POSCAR", format="poscar")

# Create from numpy array
pos = np.random.rand(100, 3) * 10
system = System(pos=pos, box=10)

# Create from polars DataFrame
import polars as pl

data = pl.DataFrame(
    {
        "x": np.random.rand(100),
        "y": np.random.rand(100),
        "z": np.random.rand(100),
        "type": np.ones(100, dtype=int),
    }
)
system = System(data=data, box=[10, 10, 10])

Performing structural analysis:

# Identify crystal structures
system.cal_polyhedral_template_matching()

# Calculate radial distribution function
rdf = system.cal_radial_distribution_function(rc=10.0, nbin=200)

# Perform cluster analysis
system.cal_cluster_analysis(rc=3.0)

Accessing and modifying data:

# Get positions
pos = system.data.select("x", "y", "z")

# Add a new column
new_data = system.data.with_columns(energy=np.random.rand(system.N))
system.update_data(new_data)

# Save to file
system.write("output.dump")

mdapy.tool_function module¶

mdapy.tool_function.average_by_neighbor(average_rc: float, data: DataFrame, property_name: str, verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray, include_self: bool = True, output_name: str | None = None) → DataFrame[source]¶

Compute the averaged property of each atom based on its neighbors within a given cutoff radius.

This function calculates the local average of a specified property (e.g., temperature, stress) for each atom by averaging over its neighboring atoms within the cutoff average_rc. The neighbor information is provided through the Verlet list and corresponding distances.

Parameters:

average_rc (float) – The cutoff radius used for averaging the property.
data (pl.DataFrame) – The input atomic data containing the property to be averaged.
property_name (str) – The name of the column in data to average.
verlet_list (np.ndarray) – A 2D array of neighbor indices. Each row corresponds to one atom and lists its neighbors.
distance_list (np.ndarray) – A 2D array of distances corresponding to the verlet_list entries.
neighbor_number (np.ndarray) – A 1D array indicating the actual number of neighbors for each atom.
include_self (bool, optional, default=True) – Whether to include the atom itself when computing the average.
output_name (str, optional) – The name of the output column to store the averaged property. If not provided, defaults to f"{property_name}_ave".

Returns:

pl.DataFrame – A new DataFrame containing the original data plus an additional column with the averaged property values.

mdapy.tool_function.sort_neighbor(verlet_list: ndarray, distance_list: ndarray, neighbor_number: ndarray, k: int)[source]¶

Sort the first k neighbors of each atom by ascending distance.

This function reorders the first k neighbor entries in the Verlet list and distance list so that the nearest neighbors appear first.

Parameters:

verlet_list (np.ndarray) – A 2D array of neighbor indices. Each row corresponds to one atom and lists its neighbors.
distance_list (np.ndarray) – A 2D array of distances corresponding to the verlet_list entries.
neighbor_number (np.ndarray) – A 1D array indicating the actual number of neighbors for each atom.
k (int) – The number of nearest neighbors to sort. Must be less than or equal to the minimum value in neighbor_number.

Raises:

AssertionError – If any atom has fewer than k neighbors.

Examples

>>> verlet = np.array([[3, 1, 2], [0, 2, 3]])
>>> distance = np.array([[0.3, 0.1, 0.2], [0.5, 0.2, 0.4]])
>>> neighbor_num = np.array([3, 3])
>>> sort_neighbor(verlet, distance, neighbor_num, k=2)
>>> verlet
array([[1, 2, 3],
       [2, 3, 0]])
>>> distance
array([[0.1, 0.2, 0.3],
       [0.2, 0.4, 0.5]])

mdapy.tool_function.wrap_pos(data: DataFrame, box: Box) → DataFrame[source]¶

Wrap position into box based on PBC.

Args:: data (pl.DataFrame): atom information. box (Box): box information
Returns:: pl.DataFrame: DataFrame with new wraped position.

mdapy.tool_function.replicate(data: DataFrame, box: Box, nx: int, ny: int, nz: int) → Tuple[DataFrame, Box][source]¶

Replicates atomic data along the x, y, and z directions.

This function creates a supercell by replicating the input atomic positions and associated data nx times along x, ny times along y, and nz times along z. The input DataFrame is expected to contain at least ‘x’, ‘y’, and ‘z’ columns for atomic coordinates.

Args:: data: Atomic information, including positions in ‘x’, ‘y’, ‘z’ columns. box: Simulation box information to be scaled. nx: Number of replications along the x direction. ny: Number of replications along the y direction. nz: Number of replications along the z direction.
Returns:: A tuple containing the replicated DataFrame (with updated positions and optionally renumbered ‘id’ column) and the scaled Box object.

mdapy.tool_function.split_xyz(input_file, output_dir='res', output_prefix=None, in_memory=True)[source]¶

Split a multi-frame XYZ file into individual frame files.

Parameters:

input_file (str) –
Path to the input XYZ file containing multiple frames. The file should follow standard XYZ format:
- Line 1: Number of atoms (integer)
- Line 2: Comment line
- Lines 3 to N+2: Atomic coordinates (element x y z …)
- Repeat for each frame
output_dir (str, optional) – Directory where individual frame files will be saved. The directory will be created if it doesn’t exist. Default is “res”.
output_prefix (str, optional) – Prefix for output filenames. If None, uses the input filename (without extension) as prefix. Output files are named as: {output_prefix}.{frame:0Nd}.xyz where N is the number of digits needed (e.g., 5 digits for 10000 frames). Default is None.
in_memory (bool, optional) – If input file is too big to load into memory, set this parameter to False. Default is True.

Returns:

None

Examples

>>> # Split trajectory.xyz into res/trajectory.00000.xyz, res/trajectory.00001.xyz, ...
>>> split_xyz("trajectory.xyz")

>>> # Specify custom output directory and prefix
>>> split_xyz("trajectory.xyz", output_dir="frames", output_prefix="md")
>>> # Output: frames/md.00000.xyz, frames/md.00001.xyz, ...

>>> # For a file with 100000 frames, output will use 6 digits
>>> split_xyz("large_traj.xyz")
>>> # Output: res/large_traj.000000.xyz, ..., res/large_traj.099999.xyz

Raises:

RuntimeError – If the input file cannot be opened or contains no valid frames.
OSError – If the output directory cannot be created.

mdapy.tool_function.generate_velocity(N, mass, temperature, remove_com=True, seed=None)[source]¶

Generate velocities from Maxwell-Boltzmann distribution at given temperature.

Parameters:

N (int) – Number of atoms.
mass (float or array_like) – Atomic mass in g/mol. Can be a scalar (all atoms same mass) or array of length N (different masses).
temperature (float) – Target temperature in Kelvin.
remove_com (bool, optional) – If True, remove center-of-mass velocity to ensure zero total momentum. Default is True.
seed (int, optional) – Random seed for reproducibility. Default is None.

Returns:

vel (ndarray, shape (N, 3)) – Velocities in Å/fs.

Notes

The velocity is generated from Maxwell-Boltzmann distribution:

P(v) ∝ exp(-m*v²/(2*kb*T))

The standard deviation of velocity component is:

σ_v = sqrt(kb*T/m)

Units:

Temperature: K
Mass: g/mol
Velocity: Å/fs
kb = 1.380649e-23 J/K

mdapy.trajectory module¶

class mdapy.trajectory.XYZTrajectory(filename: str | None = None, systems: List[System] | None = None, fast_mode: bool = False, verbose: bool = True)[source]¶

Bases: _TrajectoryListBase

XYZ trajectory file reader and manager.

This class provides functionality to read, manipulate, and write XYZ format trajectory files. It supports both classical XYZ format (element + coordinates) and extended XYZ format (with periodic boundary conditions and additional properties).

Two reading modes are available: - Serial mode: Read frames sequentially (default) - Fast mode: Optimized batch reading assuming all frames have identical columns

Parameters:

filename (Optional[str]) – Path to XYZ trajectory file to load
systems (Optional[List[System]]) – List of System objects to initialize trajectory from memory
fast_mode (bool, default=False) – If True, use optimized reading assuming all frames have the same columns. This mode is significantly faster but requires all frames to have identical column structure.

Raises:

ValueError – If neither filename nor systems is provided

Examples

>>> # Load trajectory in serial mode
>>> traj = XYZTrajectory("trajectory.xyz")

>>> # Load trajectory in fast mode
>>> traj = XYZTrajectory("trajectory.xyz", fast_mode=True)

>>> # Create trajectory from existing System objects
>>> traj = XYZTrajectory(systems=[system1, system2])

>>> # Access frames
>>> print(len(traj))  # Number of frames
>>> frame = traj[0]  # Get first frame
>>> sub_traj = traj[0:10]  # Slice trajectory

>>> # Iterate over frames
>>> for frame in traj:
...     print(frame.N)

>>> # Save trajectory
>>> traj.save("output.xyz")
>>> traj.save("output.xyz", frames=[0, 1, 2])  # Save specific frames

save(filename: str, frames: List[int] | int | None = None, mode: str = 'w', vacuum: float = 0.0) → None[source]¶

Save trajectory to XYZ file.

Parameters:

filename (str) – Output file path
frames (Optional[Union[List[int], int]], default=None) – Frame indices to save. Can be: - None: save all frames - int: save single frame - List[int]: save specified frames
mode (str, default='w') – Writing mode can be: - ‘w’ : write mode - ‘a’ : append mode
vacuum (float, default=0.0) – When > 0, every non-periodic axis of every saved frame is padded by vacuum Å (atoms shifted by vacuum / 2 so they sit centred in the new cell, boundary flipped to periodic). Useful for auto-boxing classical-XYZ frames so downstream MD code sees a well-defined supercell. The in-memory trajectory is not mutated — padding is applied to a per-frame copy at write time.

Examples

>>> traj.save("output.xyz")              # save all frames
>>> traj.save("output.xyz", 0)           # save first frame only
>>> traj.save("output.xyz", [0, 5, 10])  # save specific frames
>>> traj.save("output.xyz", vacuum=200)  # auto-box FFF frames

class mdapy.trajectory.Trajectory(filename: str | None = None, systems: List[System] | None = None, format: str | None = None, fast_mode: bool = False, verbose: bool = True)[source]¶

Bases: _TrajectoryListBase

Multi-frame trajectory container — supports XYZ and LAMMPS dump (read + write).

Parameters:

filename (str, optional) – Path to load (.xyz, .dump, .lammpstrj, optionally .gz).
systems (list of mdapy.System, optional) – Pre-built list of frames.
format ({'xyz', 'dump'}, optional) – Override file-format detection. Defaults to inferring from the filename extension.
fast_mode (bool, default=False) – Use the vectorised XYZ reader (_read_xyz_fast). Amortises per-frame overhead (regex-parsing the comment line, building the column schema, constructing Box) over a single bulk pass — typically 5–7× faster on a long, regular XYZ trajectory. Requires identical column schema across frames AND single-character whitespace separators in the per-atom block; raises a ValueError naming the offending frame otherwise. fast_mode=True is not supported for LAMMPS dump — the dump serial reader already vectorises each frame internally (via pl.read_csv), so a separate “bulk” path adds complexity without measurable speedup. Pass fast_mode=True on a dump file and you get a ValueError saying so.
verbose (bool, default=True) – Print a one-line frame i/N progress update every 200 frames during loading. Pass verbose=False to silence the output (useful for tests / scripts).

Examples

>>> traj = mp.Trajectory("dump.lammpstrj")
>>> for frame in traj: print(frame.N)
>>> traj.save("subset.xyz", frames=[0, 2, 4])
>>> traj.append(other_system)
>>> # boolean / integer-array indexing for filtering frames
>>> hot = traj[traj.get_atoms_count() > 100]      # bool mask
>>> first_few = traj[[0, 5, 7]]                    # integer mask
>>> # silence progress output
>>> traj = mp.Trajectory("big.xyz", fast_mode=True, verbose=False)

save(filename: str, frames: List[int] | int | None = None, mode: str = 'w', format: str | None = None, vacuum: float = 0.0) → None[source]¶

Write the trajectory to disk.

Parameters:

filename (str) – Output path. Format is auto-detected from the suffix.
frames (int or list of int, optional) – Subset of frame indices to save. Default: all frames.
mode ({'w', 'a'}, default 'w') – Open mode (write/truncate vs. append).
format ({'xyz', 'dump'}, optional) – Override format detection.
vacuum (float, default=0.0) – XYZ output only. When > 0, every non-periodic axis of every saved frame is padded by vacuum Å (atoms shifted by vacuum / 2 so they sit centred in the new cell, boundary flipped to periodic). Useful for auto-boxing training-set frames that came in as classical XYZ (PBC = FFF) so downstream MD code sees a well-defined supercell. The in-memory trajectory is not mutated — padding is applied to a per-frame copy at write time. Ignored for the dump format (LAMMPS dumps already require an explicit box).

unwrap() → Trajectory[source]¶: Return a new Trajectory with continuous (unwrapped) particle positions. See mdapy.unwrap_trajectory() for the full algorithm description and edge cases.

mdapy.unwrap_trajectory module¶

Unwrap periodic-boundary trajectories so that particle paths become continuous across frames — equivalent to OVITO’s UnwrapTrajectoriesModifier.

mdapy.unwrap_trajectory.unwrap_trajectory(traj: Trajectory) → Trajectory[source]¶

Make particle trajectories continuous across periodic boundaries.

Equivalent to OVITO’s UnwrapTrajectoriesModifier. The function picks one of three methods, in order of decreasing reliability, by inspecting the columns of frame 0:

Pre-unwrapped coordinates (xu, yu, zu). These come straight from the simulation code and are already continuous, so the function just renames them as the new x, y, z.
Periodic-image flags (ix, iy, iz). LAMMPS writes these whenever dump_modify ... pbc yes is set. We compute \(\vec{r}_\text{unwrapped} = \vec{r} + i_x \vec{a} + i_y \vec{b} + i_z \vec{c}\) using each frame’s own cell vectors, so trajectories that change box (NPT) are handled correctly.
Minimum-image heuristic. Walks the trajectory in time and accumulates per-atom integer PBC shift vectors from fractional- coordinate jumps greater than half a cell edge. This is the same algorithm OVITO uses when no image flags are available.

When an id column is present every frame is sorted by id before unwrapping, so the output rows line up across frames and cross-frame analyses (MSD, drift, …) are a single subtraction.

Parameters:: traj (mdapy.Trajectory) – Multi-frame trajectory. Every frame must contain the same number of atoms and (if an id column is present) the same id set.
Returns:: mdapy.Trajectory – New trajectory with the same number of frames. Each frame keeps only the id, type, element columns (whichever existed in the input) plus the unwrapped x, y, z and the original simulation box. The selected unwrap method is recorded on the returned object as trajectory._unwrap_method.
Raises:: ValueError – If the trajectory is empty, atom counts differ between frames, ids contain duplicates, or the id set differs between frames.

Warning

A RuntimeWarning is issued when boundary flags change across frames (the frame-0 flags are used throughout) and when a likely LAMMPS triclinic-cell flip is detected during the minimum-image scan.

Notes

The minimum-image heuristic cannot distinguish a true periodic- boundary crossing from genuine physical motion that exceeds half a cell edge in a single frame. If your dump frequency is too coarse, the unwrapped trajectory will be silently wrong — re-dump with a smaller stride or have LAMMPS write the ix iy iz image flags so this routine takes the second branch.

mdapy.visualize module¶

class mdapy.visualize.View(system: System)[source]¶

Bases: object

Visualize atomic systems using k3d.

Parameters:: system (System) – MDAPY System object containing atomic positions, box information and optional per-atom properties such as element, type, radius, etc.

system¶

The input atomic system.

Type:: System

plot¶

The k3d canvas used for visualization.

Type:: Plot

atoms¶

Object storing atomic coordinates, colors and radii.

Type:: k3d.points

box¶

Object representing the simulation box.

Type:: k3d.lines

label¶

The text label for colorbar. Only created when colored_by() is used.

Type:: Text2d or None

colored_by_type() → None[source]¶

Color atoms using their type values.

Notes

The type is mapped cyclically into nine predefined colors. Colors are updated in system.data[“color”] and applied to the k3d point object if already initialized.

colored_by_element() → None[source]¶

Color atoms using their element name.

Notes

Element symbols are mapped to integer colors according to ele_dict.

init_plot() → None[source]¶

Initialize the visualization canvas.

Notes

Creates box lines, atomic points.
Initializes default color and radius.
Creates element/type legends using text2d.

display() → None[source]¶: Display the k3d plot in supported environments (e.g., Jupyter).

close() → None[source]¶: Close the k3d plot and release the rendering canvas.

hide_object_by_group_name(name: str, remove: bool = False) → None[source]¶

Hide or remove k3d objects by their group name.

Parameters:

name (str) – The object group to hide/remove.
remove (bool, default False) – If True, remove the object entirely. If False, only hide it visually.

delete_color_bar() → None[source]¶: Remove existing colorbar from the plot.

colored_by(name: str, vmin: float | None = None, vmax: float | None = None, cmap: str = 'rainbow') → None[source]¶

Color atoms based on a given scalar per-atom quantity.

Parameters:

name (str) – Column name in system.data used for coloring.
vmin (float, optional) – Minimum value of the colormap. If None, automatically determined.
vmax (float, optional) – Maximum value of the colormap. If None, automatically determined.
cmap (str, default "rainbow") – Matplotlib colormap name.

Notes

If name is “element” or “type”, discrete coloring is applied.
If name is a structure classifier (“ptm”, “cna”, “aja”, “ids”), structure coloring is used.
Otherwise, continuous colormap coloring is used.
Colorbar is updated accordingly.

mdapy.void_analysis module¶

class mdapy.void_analysis.VoidAnalysis(system: System, rc: float)[source]¶

Bases: object

Detect voids in an atomic system by discretizing the simulation box into cells and locating empty cells.

Parameters:

system (System) – The atomic system containing coordinates and box information.
rc (float) – Cutoff radius that controls the size of the spatial grid cells.

system¶

Input system.

Type:: System

rc¶

Grid cell size parameter.

Type:: float

void_system¶

A System containing the centers of detected void regions. None if no voids are found.

Type:: Optional[System]

void_number¶

Number of distinct void clusters.

Type:: int

void_volume¶

Estimated total void volume, computed as N_void_cells * rc^3.

Type:: float

compute()[source]¶

Perform the void detection.

Compute the number of grid cells ncell along each dimension based on the system thickness and cutoff rc.
Use _neighbor._fill_cell_for_void to assign each cell a value: 0 for empty, 1 for occupied.
If empty cells exist, compute their geometric centers and create a System object containing void positions.
Perform cluster analysis on these void positions.
Filter clusters that contain only one void cell (noise avoidance).
Renumber clusters and compute final void_number and void_volume.

mdapy.voronoi module¶

This module provides efficient Voronoi tessellation analysis for atomic systems, enabling calculation of Voronoi cells, neighbor identification, volumes, and detailed geometric information. The implementation uses OpenMP parallelization through the Voro++ library for high performance.

The module supports both orthogonal and triclinic simulation boxes with periodic boundary conditions.

References

Lu J, Lazar E A, Rycroft C H. An extension to Voro++ for multithreaded computation of Voronoi cells. Computer Physics Communications, 2023, 291: 108832. https://doi.org/10.1016/j.cpc.2023.108832

class mdapy.voronoi.Voronoi(box: Box, data: DataFrame)[source]¶

Bases: object

Voronoi tessellation analysis for atomic systems.

This class provides methods to compute Voronoi cells and their properties for a system of atoms/particles. It supports both orthogonal and triclinic simulation boxes with periodic boundary conditions.

Parameters:

box (Box) – The simulation box object containing box vectors, boundary conditions, and origin information.
data (pl.DataFrame) – A Polars DataFrame containing atomic positions with columns ‘x’, ‘y’, ‘z’.

box¶

The simulation box object.

Type:: Box

data¶

The atomic position data.

Type:: pl.DataFrame

_enlarge_data¶

Internal replicated atomic data when periodic extension is required.

Type:: pl.DataFrame, optional

_enlarge_box¶

The enlarged simulation box corresponding to replicated atoms.

Type:: Box, optional

Notes

The implementation automatically handles small systems by replicating atoms to ensure sufficient neighbors for accurate Voronoi tessellation. For triclinic boxes, automatic rotation to LAMMPS-compatible format is performed when necessary.

get_neighbor(a_face_area_threshold: float = -1.0, r_face_area_threshold: float = -1.0) → Tuple[ndarray, ndarray, ndarray, ndarray][source]¶

Calculate Voronoi neighbors and face properties for all atoms.

This method performs Voronoi tessellation to identify nearest neighbors based on shared Voronoi cell faces. It can filter neighbors based on face area thresholds to exclude insignificant contacts.

Parameters:

a_face_area_threshold (float, optional) – Absolute face area threshold. Faces with area below this value are ignored as neighbors. Default is -1.0 (no filtering).
r_face_area_threshold (float, optional) – Relative face area threshold (fraction of average face area). Faces with relative area below this value are ignored. Default is -1.0 (no filtering).

Returns:

verlet_list (np.ndarray) – 2D array of shape (N, max_neighbors) containing neighbor indices for each atom. Unused entries are filled with -1.
distance_list (np.ndarray) – 2D array of shape (N, max_neighbors) containing distances to each neighbor (in simulation units).
face_area (np.ndarray) – 2D array of shape (N, max_neighbors) containing the area of the Voronoi face shared with each neighbor.
neighbor_number (np.ndarray) – 1D array of shape (N,) containing the coordination number (number of Voronoi neighbors) for each atom.

Notes

If both thresholds are provided, the larger effective threshold is used.
The method automatically handles periodic boundary conditions.
For small systems (<50 atoms), automatic replication is performed to ensure accurate neighbor identification.
OpenMP parallelization is used for performance with large systems.

get_cell_info() → Tuple[List[List[List[int]]], List[List[List[float]]], List[float], List[float], List[List[float]]][source]¶

Retrieve detailed Voronoi polygon information for all atoms.

This method provides comprehensive geometric information about each Voronoi cell, including face vertices, positions, and areas. Due to the detailed nature of the output, this method is computationally intensive and memory-demanding, making it suitable primarily for small systems.

Returns:

face_vertices_indices (List[List[List[int]]]) – For each atom, a list of faces, where each face is defined by a list of vertex indices that form the face polygon.
face_vertices_positions (List[List[List[float]]]) – For each atom, a list of faces, where each face contains the 3D positions (x, y, z) of all vertices forming that face.
volume (List[float]) – The volume of each atom’s Voronoi cell (in cubic simulation units).
radius (List[float]) – The cavity radius for each atom - the distance from the atom to the farthest vertex of its Voronoi cell (in simulation units).
face_areas (List[List[float]]) – For each atom, a list containing the area of each face of its Voronoi cell (in square simulation units).

Raises:

AssertionError – If the box is triclinic (only orthogonal boxes are supported).
AssertionError – If the system contains less than 2 atoms.

Warning

This method is computationally expensive and requires significant memory, especially for large systems. It should primarily be used for small systems.

For routine analysis of large systems, use get_neighbor() or get_volume() instead.

get_volume() → Tuple[ndarray, ndarray, ndarray][source]¶

Calculate Voronoi cell volumes, coordination numbers, and cavity radii.

This method efficiently computes key Voronoi cell properties without storing detailed geometric information, making it suitable for large-scale analysis of atomic systems.

Returns:

volume (np.ndarray) – 1D array of shape (N,) containing the volume of each atom’s Voronoi cell (in cubic simulation units of length).
neighbor_number (np.ndarray) – 1D array of shape (N,) containing the coordination number for each atom. This is the number of faces of the Voronoi cell, which equals the number of nearest neighbors.
cavity_radius (np.ndarray) – 1D array of shape (N,) containing the cavity radius for each atom. This is the distance from the atom to the farthest vertex of its Voronoi cell, representing the radius of the largest empty sphere (containing no other atoms) that touches the atom.

class mdapy.voronoi.Cell(face_vertices: List[List[int]], vertices: ndarray, volume: float, cavity_radius: float, face_areas: ndarray, pos: ndarray)[source]¶

Bases: object

A lightweight container representing a single Voronoi cell and all its geometric properties.

Parameters:

face_vertices (List[List[int]]) – A list of faces, where each face is represented by a list of vertex indices referring to positions in the vertices array.
vertices (np.ndarray) – Array of shape (M, 3) containing the 3D coordinates of all polygon vertices that form the Voronoi cell.
volume (float) – The volume of the Voronoi cell.
cavity_radius (float) – The distance from the particle to the farthest vertex of its Voronoi cell, representing the largest empty-sphere radius.
face_areas (np.ndarray) – Array containing the area of each face of the Voronoi cell.
pos (np.ndarray) – The (x, y, z) coordinates of the particle associated with this Voronoi cell.

Notes

This class stores only geometric results and contains no computational methods. Instances of this class are typically created internally by the Container class.

face_vertices: List[List[int]]¶

vertices: ndarray¶

volume: float¶

cavity_radius: float¶

face_areas: ndarray¶

pos: ndarray¶

class mdapy.voronoi.Container(data: DataFrame | ndarray, box: Box)[source]¶

Bases: object

High-level container that stores Voronoi cell geometry for all atoms in a system.

This class wraps the Voronoi computation and provides Python‑friendly access to each atom’s Voronoi cell, represented by a Cell object.

Parameters:

data (Union[pl.DataFrame, np.ndarray]) –
Atomic coordinates. Accepted formats:
- A Polars DataFrame with columns ‘x’, ‘y’, ‘z’.
- A NumPy array of shape (N, 3).
box (Box) – The simulation box defining boundaries and coordinate transformation.

_data¶

A list containing a Cell instance for each atom in the input system.

Type:: List[Cell]

Examples

>>> container = Container(pos_array, box)
>>> cell_0 = container[0]
>>> print(cell_0.volume)

The object behaves like a Python list of Cell objects.

mdapy.warren_cowley_parameter module¶

class mdapy.warren_cowley_parameter.WarrenCowleyParameter(verlet_list: ndarray, neighbor_number: ndarray, data: DataFrame)[source]¶

Bases: object

Calculate the Warren-Cowley Parameter (WCP) for chemical short-range order analysis.

The Warren-Cowley Parameter quantifies chemical ordering in multi-component systems by measuring the deviation of local atomic arrangements from random mixing. It is particularly useful for analyzing segregation, clustering, and ordering tendencies in alloys.

Parameters:

verlet_list (np.ndarray) – Neighbor list array of shape (N, max_neigh).
neighbor_number (np.ndarray) – Number of neighbors for each atom, shape (N,).
data (pl.DataFrame) – Atomic data containing either ‘element’ or ‘type’ column.

verlet_list¶

Neighbor indices.

Type:: np.ndarray

neighbor_number¶

Neighbor counts.

Type:: np.ndarray

data¶

Input atomic data.

Type:: pl.DataFrame

type_list¶

Atom types (0-indexed) derived from data.

Type:: np.ndarray

Ntype¶

Number of distinct atom types.

Type:: int

ele2type¶

Mapping from element symbols to type indices (if using element column).

Type:: dict, optional

WCP¶

WCP matrix of shape (Ntype, Ntype) after calling compute(). Element (i,j) represents WCP for type-j atoms around type-i atoms.

Type:: np.ndarray

Notes

The Warren-Cowley Parameter α_ij is defined as:

\[\alpha_{ij} = 1 - \frac{P_{ij}}{c_j}\]

where P_ij is the probability of finding a type-j atom around a type-i atom, and c_j is the global concentration of type-j atoms.

Interpretation: - α_ij < 0: Preference for unlike neighbors (ordering) - α_ij = 0: Random mixing - α_ij > 0: Preference for like neighbors (clustering/segregation) - α_ij = 1: Complete segregation (no unlike neighbors)

References

Warren, B. E. (1990). X-ray Diffraction. Dover Publications.

Cowley, J. M. (1950). An approximate theory of order in alloys. Physical Review, 77(5), 669.

compute() → None[source]¶

Compute the Warren-Cowley Parameter matrix.

This method calculates the WCP for all type pairs and stores the result in the WCP attribute as a symmetric matrix.

Notes

After calling this method, the WCP attribute will contain an (Ntype, Ntype) array with WCP values for each type pair.

plot(elements_list: List[str] | None = None, fig: Figure | None = None, ax: Axes | None = None, vmin: float = -2, vmax: float = 1, cmap: str = 'GnBu') → Tuple[Figure, Axes][source]¶

Visualize the Warren-Cowley Parameter matrix as a heatmap.

Parameters:

elements_list (list of str, optional) – Element symbols for axis labels. If None and ele2type exists, uses element symbols from ele2type.
fig (matplotlib.figure.Figure, optional) – Existing figure to plot on. If None, creates new figure.
ax (matplotlib.axes.Axes, optional) – Existing axes to plot on. If None, creates new axes.
vmin (float, default=-2) – Minimum value for colormap scaling.
vmax (float, default=1) – Maximum value for colormap scaling.
cmap (str, default='GnBu') – Matplotlib colormap name.

Returns:

fig (matplotlib.figure.Figure) – Figure object containing the plot.
ax (matplotlib.axes.Axes) – Axes object containing the plot.

Notes

The plot displays WCP values as a colored matrix with numerical annotations. Rows represent central atom types, columns represent neighboring atom types.

mdapy.wigner_seitz_defect module¶

class mdapy.wigner_seitz_defect.WignerSeitzAnalysis(ref: System, affine: bool = False)[source]¶

Bases: object

Wigner-Seitz defect analysis.

This class performs point defect analysis using the Wigner-Seitz cell method, similar to the implementation in OVITO. It identifies vacancies and interstitials by comparing a reference configuration (ideal lattice sites) with a current configuration (displaced atoms). Each atom in the current configuration is assigned to the nearest reference site. Site occupancies are then computed to detect defects: sites with occupancy 0 are vacancies, and sites with occupancy >1 indicate interstitials.

Parameters:

ref (System) – Reference system defining the ideal lattice sites.
affine (bool, optional) – If True, applies an affine transformation to map current positions to the reference cell before assignment (handles cell deformation). Default is False.

compute(current: System) → Dict[str, int | ndarray][source]¶

Perform Wigner-Seitz analysis.

This method assigns each atom in the current configuration to the nearest reference site and computes site occupancies to identify defects.

Parameters:

current (System) – The system to analyze (the “current” or “displaced” configuration).

Returns:

dict –

Dictionary containing:

’site_occupancy’np.ndarray (ref.N,): Occupancy counts per reference site.
’atom_site_index’np.ndarray (current.N,): For each atom, the index of the assigned reference site.
’atom_site_type’np.ndarray (current.N,): For each atom, the type of the assigned reference site.
’atom_occupancy’np.ndarray (current.N,): For each atom, the occupancy of its assigned site.
’vacancy_count’int: Number of vacant sites (occupancy == 0).
’interstitial_count’int: Number of sites with interstitials (occupancy > 1).

Index¶

mdapy.ackland_jones_analysis module¶

mdapy.angular_distribution_function module¶

mdapy.atomic_strain module¶

mdapy.atomic_temperature module¶

mdapy.bond_analysis module¶

mdapy.bond_stiffness module¶

Workflow¶

Citation¶

mdapy.box module¶

mdapy.build_bond module¶

mdapy.build_lattice module¶

mdapy.calculator module¶

Abstract Base Class for Static Calculators¶

Classes¶

mdapy.centro_symmetry_parameter module¶

mdapy.chill_plus module¶

mdapy.cluster_analysis module¶

mdapy.common_neighbor_analysis module¶

mdapy.common_neighbor_parameter module¶

mdapy.create_polycrystal module¶

mdapy.data module¶

mdapy.eam module¶

Embedded Atom Method (EAM) Potential Calculator¶

mdapy.elastic module¶

mdapy.identify_diamond_structure module¶

mdapy.identify_fcc_planar_faults module¶

mdapy.knn module¶

mdapy.lammps_potential module¶

mdapy.lammps_runner module¶

mdapy.lindemann_parameter module¶

mdapy.load_save module¶

mdapy.md_elastic module¶

Finite-temperature elastic constants from MD stress-strain¶

mdapy.mean_squared_displacement module¶

mdapy.minimizer module¶

mdapy.neighbor module¶

mdapy.nep module¶

mdapy.nep4ase module¶

mdapy.orthogonal_cell module¶

Usage¶

mdapy.parallel module¶

Usage¶

Rules¶

mdapy.phonon module¶

mdapy.pigz module¶

mdapy.plotset module¶

Plotting Utilities for Scientific Publications¶

Functions¶

mdapy.polyhedral_template_matching module¶

mdapy.potential_tool module¶

mdapy.qha_elastic module¶

mdapy.radial_distribution_function module¶

mdapy.render module¶

Perspective mode¶

Orthographic mode¶

OVITO coordinate convention¶

All eight available view names¶

mdapy.spatial_binning module¶

mdapy.spline module¶

mdapy.sqs module¶

Quick recipe¶

How to choose cutoffs¶

What the parameters mean¶

Mathematical background¶

Citation¶

mdapy.steinhardt_bond_orientation module¶

mdapy.structure_entropy module¶

mdapy.structure_factor module¶

mdapy.system module¶

Core System Class for MDAPY¶

Classes¶

mdapy.tool_function module¶

mdapy.trajectory module¶

mdapy.unwrap_trajectory module¶

mdapy.visualize module¶

mdapy.void_analysis module¶

mdapy.voronoi module¶

mdapy.warren_cowley_parameter module¶

mdapy.wigner_seitz_defect module¶

How to choose `cutoffs`¶