Index

mdapy.ackland_jones_analysis module

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

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

[1]

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

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)

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:

1. For each central atom i of type A

2. Find neighbors j of type B within [rc_AB_min, rc_AB_max]

3. Find neighbors k of type C within [rc_AC_min, rc_AC_max]

4. Calculate angle θ between j-i-k bonds

5. 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()

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]

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)

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:

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

2. Calculate deformation gradient tensor:

   \[F = (W \cdot V^{-1})^T\]

3. Compute strain tensor (Green-Lagrange strain):

   \[\varepsilon = \frac{1}{2}(F^T F - I)\]

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

[1]

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)

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)

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

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)

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()

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]

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]

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.box module

class mdapy.box.Box(box: int | float | Iterable[float] | ndarray | Box, boundary: Iterable[int] | ndarray | None = None, origin: Iterable[float] | ndarray | None = None)

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

box

The 3x3 box matrix where each row represents a box vector.

Type:

np.ndarray

origin

The origin position of the box.

Type:

np.ndarray

boundary

Boundary condition flags for each dimension.

Type:

np.ndarray

triclinic

True if the box is triclinic (non-orthogonal), False otherwise.

Type:

bool

volume

The volume of the simulation box.

Type:

float

inverse_box

The inverse of the box matrix for coordinate transformations.

Type:

np.ndarray

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

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

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

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.

align_to_lammps_box() → Tuple[Box, ndarray]

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

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

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

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_lattice module

mdapy.build_lattice.build_crystal(name: str, 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_over_a: float = np.float64(1.632993161855452)) → System

Build crystal structure with optional Miller indices orientation.

This function creates a crystal structure with specified lattice type and orientation. It supports standard crystallographic structures (FCC, BCC, HCP, diamond, graphene) and allows custom orientations via Miller indices for cubic structures.

Parameters:

• name (str) – Element symbol (e.g., ‘Cu’, ‘Al’, ‘Mg’, ‘C’).

• structure (str) –

  Crystal structure type. Supported values:

  • 'fcc': Face-centered cubic

  • 'bcc': Body-centered cubic

  • 'diamond': Diamond cubic

  • 'hcp': Hexagonal close-packed

  • 'graphene': Graphene layer

• a (float) –

  Lattice constant in Angstroms.

  • For cubic structures (FCC/BCC/diamond): edge length of cubic unit cell

  • For HCP: basal plane lattice parameter

  • For graphene: C-C bond length

• miller1 (tuple of int, optional) – First Miller index (h, k, l) defining the x-axis direction. Only applicable for cubic structures. If None, uses standard orientation [1,0,0].

• miller2 (tuple of int, optional) – Second Miller index (h, k, l) defining the y-axis direction. Only applicable for cubic structures. If None, uses standard orientation [0,1,0].

• miller3 (tuple of int, optional) –

  Third Miller index (h, k, l) defining the z-axis direction. Only applicable for cubic structures. If None, uses standard orientation [0,0,1].

  Note

  The three Miller indices must be mutually orthogonal and satisfy the right-hand rule: miller1 × miller2 = miller3.

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

• c_over_a (float, default=sqrt(8/3)) – Ratio of c to a for HCP structure. Default value gives ideal HCP packing. Ignored for other structures.

Returns:

System – A System object containing the built crystal structure with atomic positions and simulation box.

Raises:

ValueError –

• If an unsupported structure type is specified - If Miller indices are provided for non-cubic structures - If Miller indices are not mutually orthogonal - If Miller indices don’t satisfy the right-hand rule

Examples

Create a standard FCC copper crystal:

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

Create an FCC crystal with custom Miller orientation:

>>> cu_rotated = build_crystal(
...     "Cu",
...     "fcc",
...     a=3.615,
...     miller1=(1, -1, 0),
...     miller2=(1, 1, -2),
...     miller3=(1, 1, 1),
...     nx=5,
...     ny=5,
...     nz=5,
... )

Create an HCP magnesium crystal:

>>> mg = build_crystal("Mg", "hcp", a=3.21, c_over_a=1.624, nx=10, ny=10, nz=10)

Create a graphene sheet:

>>> graphene = build_crystal("C", "graphene", a=1.42, nx=20, ny=20, nz=1)

Notes

• For cubic structures with Miller indices, the function automatically finds the minimal periodic cell and aligns it to coordinate axes.

• The resulting structure has an orthogonal simulation box aligned with the Cartesian coordinate system.

• Atomic positions are generated by replicating the unit cell nx×ny×nz times.

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_over_a: float = np.float64(1.632993161855452), random_seed: int | None = None) → System

Build a high-entropy alloy (HEA) crystal structure with an optional orientation defined by Miller indices.

Parameters:

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

• structure (str) –

  Crystal structure type. Supported values are:

  • 'fcc' — Face-centered cubic

  • 'bcc' — Body-centered cubic

  • 'hcp' — Hexagonal close-packed

• a (float) – Lattice constant in Ångströms. - For cubic structures (FCC/BCC): edge length of the cubic unit cell. - For HCP: basal plane lattice parameter.

• miller1 (tuple of int, optional) –

  Miller indices (h, k, l) defining the orientation of x-, y-, and z-axes. Only applicable to cubic structures. If any are None, the standard orthogonal orientation ([1,0,0], [0,1,0], [0,0,1]) is used.

  Note

  The three Miller indices must be orthogonal and follow the right-hand rule: miller1 × miller2 = miller3.

• miller2 (tuple of int, optional) –

  Miller indices (h, k, l) defining the orientation of x-, y-, and z-axes. Only applicable to cubic structures. If any are None, the standard orthogonal orientation ([1,0,0], [0,1,0], [0,0,1]) is used.

  Note

  The three Miller indices must be orthogonal and follow the right-hand rule: miller1 × miller2 = miller3.

• miller3 (tuple of int, optional) –

  Miller indices (h, k, l) defining the orientation of x-, y-, and z-axes. Only applicable to cubic structures. If any are None, the standard orthogonal orientation ([1,0,0], [0,1,0], [0,0,1]) is used.

  Note

  The three Miller indices must be orthogonal and follow the right-hand rule: miller1 × miller2 = miller3.

• nx (int, default=1) – Number of unit cell repetitions along x, y, and z directions.

• ny (int, default=1) – Number of unit cell repetitions along x, y, and z directions.

• nz (int, default=1) – Number of unit cell repetitions along x, y, and z directions.

• c_over_a (float, default=sqrt(8/3)) – The c/a ratio for the HCP structure. The default corresponds to the ideal HCP packing. Ignored for other structures.

• random_seed (int, optional) – Random seed for reproducible element assignment.

Returns:

System – A System object containing the generated crystal structure with atomic positions and simulation box information.

Examples

>>> system = 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

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

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

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.

results

Dictionary storing computed results. Keys typically include:

• ‘energies’: per-atom potential energies

• ‘forces’: per-atom forces

• ‘stress’: system stress tensor

• ‘virials’: per-atom virial tensors

• ‘energy’: system potential energy

Type:

dict

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

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

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

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

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

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)

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

[1]

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

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

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()

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)

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:

1. Adaptive cutoff (default): determines neighbors based on the 14 nearest atoms.

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

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

2. 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()

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)

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

[1]

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

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)

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

Compute and generate the polycrystalline structure.

Main workflow:

1. Generate Voronoi tessellation from seeds

2. Fill each cell with rotated unit cell atoms

3. Optionally add graphene at grain boundaries

4. Remove overlapping atoms

5. 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 Å). Values are based on Alvarez (2013). 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]}

The elemental radius and colors come from OVITO for visualization. The original statements are as below:

https://gitlab.com/stuko/ovito/-/blob/master/src/ovito/particles/objects/ParticleType.cpp#L225-329

// Define default names, colors, and radii for some predefined particle types. // // Van der Waals radii have been adopted from the VMD software, which adopted them from A. Bondi, J. Phys. Chem., 68, 441 - 452, 1964, // except the value for H, which was taken from R.S. Rowland & R. Taylor, J. Phys. Chem., 100, 7384 - 7391, 1996. // For radii that are not available in either of these publications use r = 2.0. // The radii for ions (Na, K, Cl, Ca, Mg, and Cs) are based on the CHARMM27 Rmin/2 parameters for (SOD, POT, CLA, CAL, MG, CES). // // Colors and covalent radii of elements marked with ‘//’ have been adopted from OpenBabel.

Here we time 2 for radius.

Here is the code to extract the radius and color information:

res = [i.split(’, ‘)[:5] for i in context.split(‘ParticleType::PredefinedChemicalType’)[2:]] ele_rgb, ele_radius = {}, {} for i in res:

ele, r, g, b, size = i #print(ele) ele = ele.split(‘(‘)[-1].split(‘”’)[1] if ‘/’ in r:

L, R = r.split(‘(‘)[1].split(‘/’) r = float(L[:-1]) / float(R[:-1]) L, R = g.split(‘/’) g = float(L[:-1]) / float(R[:-1]) L, R = b.split(‘/’) b = float(L[:-1]) / float(R[:-2])

else:

r = float(r.split()[1].replace(‘f’, ‘’)) g = float(g.split()[0].replace(‘f’, ‘’)) b = float(b.split()[0][:-1].replace(‘f’, ‘’))

size = float(size[:-1]) ele_rgb[ele] = [int(r*255), int(g*255), int(b*255)] ele_radius[ele] = size

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:

1. Embedding energy: \(F_i(\rho_i)\) - energy to embed atom i into the electron density \(\rho_i\)

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

[1]

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.

[2]

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)

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.

filename

Path to the EAM potential file.

Type:

str

header

Header lines from the potential file (first 3 lines).

Type:

List[str]

Nelements

Number of elements in the potential.

Type:

int

elements_list

List of element symbols supported by this potential.

Type:

List[str]

nrho

Number of points in the electron density grid.

Type:

int

drho

Spacing of the electron density grid.

Type:

float

nr

Number of points in the distance grid.

Type:

int

dr

Spacing of the distance grid.

Type:

float

rc

Cutoff radius for the potential (in Angstroms).

Type:

float

F_rho

Embedding function F(ρ) for each element. Shape: (Nelements, nrho).

Type:

NDArray[np.float64]

rho_r

Electron density function ρ(r) for each element. Shape: (Nelements, nr).

Type:

NDArray[np.float64]

phi_r

Pair potential φ(r) between element pairs. Shape: (Nelements, Nelements, nr).

Type:

NDArray[np.float64]

r

Distance grid points. Shape: (nr,).

Type:

NDArray[np.float64]

rho

Electron density grid points. Shape: (nrho,).

Type:

NDArray[np.float64]

results

Dictionary storing calculation results (energies, forces, stress, virials).

Type:

Dict[str, NDArray]

write_eam_alloy(output_name: str | None = None)

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

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

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

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

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

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

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

Creates three subplots showing:

1. Embedding function F(ρ) vs electron density ρ

2. Electron density function ρ(r) vs distance r

3. 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)

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)

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

class mdapy.elastic.ElasticConstant(system: System)

Bases: object

A class to compute the elastic constants of a crystal structure using finite strain and stress analysis based on molecular dynamics or energy minimization simulations.

The workflow involves applying small deformations to the simulation box, relaxing the atomic structure, and evaluating the resulting stress to fit the elastic stiffness tensor components.

Parameters:

system (System) – The initial system to be used for deformation and stress calculation.

system

The reference atomic system.

Type:

System

Cij

The calculated elastic constants in Voigt notation.

Type:

np.ndarray or None

static triclinic_deform(u: ndarray) → ndarray

Construct the strain-stress coupling matrix for a triclinic system.

Parameters:

u (np.ndarray) – Strain vector in the order of [uxx, uyy, uzz, uyz, uxz, uxy].

Returns:

np.ndarray – A (6, 18) array representing the linear relationship between strain and stress components for fitting.

static get_pressure(stress: ndarray) → float

Compute the hydrostatic pressure from the stress tensor.

Parameters:

stress (array-like) – Stress components in ASE format: [sxx, syy, szz, syz, sxz, sxy].

Returns:

float – The hydrostatic pressure (negative of mean normal stress).

static relax_struc(system: System, change_cell: bool = False, show_process: bool = False, fmax: float = 0.0001, steps: int = 10000)

Relax the atomic structure using the FIRE minimization algorithm.

Parameters:

• system (System) – The system to relax.

• change_cell (bool, optional) – Whether to optimize the cell shape and volume (default: False).

• show_process (bool, optional) – Whether to print relaxation progress (default: False).

• fmax (float, optional) – Force convergence criterion in eV/Å (default: 1e-4).

• steps (int, optional) – Maximum number of FIRE steps (default: 10000).

get_cart_deformed_cell(base_cryst: System, axis: int = 0, size: float = 1) → System

Generate a deformed structure by applying a small strain to the box.

Parameters:

• base_cryst (System) – The reference (undeformed) system.

• axis (int, optional) – Index of deformation mode: 0–2 for normal strains, 3–5 for shear strains (default: 0).

• size (float, optional) – Magnitude of deformation in percent (default: 1).

Returns:

System – A new system object with deformed box and updated atomic coordinates.

generate_deformations(n: int = 5, d: float = 2) → List[System]

Generate a series of deformed systems for elastic constant fitting.

Parameters:

• n (int, optional) – Number of deformation steps for each mode (default: 5).

• d (float, optional) – Maximum deformation amplitude in percent (default: 2).

Returns:

list of System – List of deformed system objects.

static get_strain(cryst: System, refcell: System) → ndarray

Compute the engineering strain tensor between deformed and reference cells.

Parameters:

• cryst (System) – The deformed system.

• refcell (System) – The undeformed (reference) system.

Returns:

np.ndarray – Strain vector [uxx, uyy, uzz, uyz, uxz, uxy].

compute(n: int = 5, d: float = 2)

Compute the elastic constants by fitting stress–strain data.

Parameters:

• n (int, optional) – Number of deformation steps per mode (default: 5).

• d (float, optional) – Maximum deformation magnitude in percent (default: 2).

Notes

The method performs the following steps:

1. Generate deformed configurations based on a fully relaxed structure.

2. Relax each configuration to minimize energy.

3. Compute the resulting stress tensors.

4. Fit the linear stress–strain relation to extract Cij.

print_Cij(scale: float = 160.2176621)

Print the calculated elastic constants in GPa.

Parameters:

scale (float, optional) – Unit conversion factor (default: 160.2176621, eV/ų → GPa).

Notes

The output follows Voigt notation with components:

C11, C22, C33, C12, C13, C23, C44, C55, C66, C16, C26, C36, C46, C56, C14, C15, C25, C45.

mdapy.identify_diamond_structure module

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

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

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 (N < 500) 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)

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 [1]_, 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

[1]

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

compute()

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)

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()

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

class mdapy.lammps_potential.LammpsPotential(pair_parameter: str, element_list: List[str], units: str = 'metal', centroid_stress: bool = False)

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.

calculate(data: DataFrame, box: Box) → None

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

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

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

Return per-atom forces; compute if necessary.

get_stress(data: DataFrame, box: Box) → Any

Return global stress in Voigt order [xx, yy, zz, yz, xz, xy]; compute if necessary.

get_virials(data: DataFrame, box: Box) → Any

Return per-atom virials (9 components) and compute if necessary.

mdapy.lindemann_parameter module

class mdapy.lindemann_parameter.LindemannParameter(pos_list: ndarray, only_global: bool = False)

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

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]

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

Bases: object

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

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

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]

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]

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]

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]

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]

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

Parameters:

filename (str) – Path to XYZ file (can be .xyz.gz).

Returns:

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

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

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

Parameters:

filename (str) – Path to POSCAR file (can be .gz).

Returns:

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

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

Read system from LAMMPS data file (supports .gz compression).

Parameters:

filename (str) – Path to data file (can be .gz).

Returns:

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

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

Read system from LAMMPS dump file.

Parameters:

filename (str) – Path to dump file.

Returns:

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

class mdapy.load_save.SaveSystem

Bases: object

static to_ase(data: pl.DataFrame, box: Box) → Atoms

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

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

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)

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)

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)

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)

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.

mdapy.mean_squared_displacement module

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

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:

1. 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)\).

2. 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*}

Important:

• 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

[1] 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.

[2] pyFFTW: https://github.com/pyFFTW/pyFFTW

compute() → None

Compute MSD for all particles and time frames.

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

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)

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

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)

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)

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

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)

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.

• 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) – Maximum number of neighbors allowed for each atom. If not provided, the neighbor list will be estimated as a largest safe value, and it will raise some performance overhead, expecially for memory.

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()

Build the neighbor list and compute interatomic distances.

This method determines whether the box is large enough to directly perform neighbor searching or needs to be replicated.

Raises:

AssertionError – If max_neigh is provided but smaller than the actual maximum neighbor count.

mdapy.nep module

class mdapy.nep.NEP(filename: str)

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]

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)

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

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

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

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

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

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

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

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

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

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)

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

implemented_properties

Properties that this calculator can compute: [‘energy’, ‘energies’, ‘forces’, ‘stress’, ‘virials’]

Type:

list

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]

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'])

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

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

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

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

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

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

Compress a file to .gz format using multiprocessing.

This function provides parallel gzip compression for faster processing of large files. It automatically uses all available CPU cores and maintains full gzip format compatibility.

Args:

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

Examples:

>>> # Compress with automatic output name
>>> compress_file("data.txt")
'data.txt.gz'

>>> # Compress with custom output name
>>> compress_file("input.txt", "output.gz")
'output.gz'

Note:

• Small files (<5MB) use single-process compression

• Large files automatically use all CPU cores

• Uses 512KB chunks for optimal parallelism

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")

mdapy.plotset.set_figure(figsize: Tuple[float | int, float | int] = (8.5, 7.0), figdpi: int = 150, nrow: int = 1, ncol: int = 1, color_cycler: List[str] | Tuple[str, ...] | None = None, **kwargs: Any) → tuple[Figure, Axes | List[Axes] | List[List[Axes]]]

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

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)

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

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

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

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

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

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

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

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

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)

Bases: object

Simple PCA implementation (similar to sklearn PCA).

fit_transform(X: ndarray) → ndarray

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

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

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 {‘Al’:0, ‘C’:1}.

• 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.radial_distribution_function module

class mdapy.radial_distribution_function.RadialDistributionFunction(rc: float, nbin: int, box: Box, verlet_list: ndarray[tuple[Any, ...], dtype[int32]], distance_list: ndarray[tuple[Any, ...], dtype[float64]], neighbor_number: ndarray[tuple[Any, ...], dtype[int32]], type_list: ndarray[tuple[Any, ...], dtype[int32]])

Bases: object

Calculate the radial distribution function (RDF) for atomic systems.

The radial distribution function g(r) describes the probability of finding an atom at distance r from a reference atom, normalized by the probability expected for a completely random distribution at the same density.

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

\[g(r) = c_{\alpha}^2 g_{\alpha\alpha}(r) + 2c_{\alpha}c_{\beta}g_{\alpha\beta}(r) + c_{\beta}^2 g_{\beta\beta}(r)\]

where \(c_{\alpha}\) and \(c_{\beta}\) denote the concentration of atom types α and β, respectively, and \(g_{\alpha\beta}(r) = g_{\beta\alpha}(r)\).

The RDF is computed as:

\[g_{\alpha\beta}(r) = \frac{V}{N_{\alpha}N_{\beta}} \sum_{i \in \alpha} \sum_{j \in \beta} \delta(r - r_{ij})\]

where V is the system volume, \(N_{\alpha}\) and \(N_{\beta}\) are the number of atoms of types α and β, and \(r_{ij}\) is the distance between atoms i and j.

Parameters:

• rc (float) – Cutoff distance for RDF calculation (in Angstroms).

• nbin (int) – Number of bins for histogram.

• box (Box) – System box object containing boundary information.

• verlet_list (NDArray[np.int32]) – Shape (N_particles, max_neigh). Verlet neighbor list where verlet_list[i, j] gives the index of the j-th neighbor of particle i. Value of -1 indicates no neighbor at that position.

• distance_list (NDArray[np.float64]) – Shape (N_particles, max_neigh). Distance list where distance_list[i, j] gives the distance between particle i and its j-th neighbor.

• neighbor_number (NDArray[np.int32]) – Shape (N_particles,). Number of neighbors for each particle.

• type_list (NDArray[np.int32]) – Shape (N_particles,). Atom type indices for each particle.

r

Shape (nbin,). Distance values at bin centers.

Type:

NDArray[np.float64]

g_total

Shape (nbin,). Total (averaged) radial distribution function.

Type:

NDArray[np.float64]

Ntype

Number of distinct atom types in the system.

Type:

int

g

Shape (Ntype, Ntype, nbin). Partial RDFs for each pair of atom types. g[i, j, :] gives the RDF between type i and type j atoms.

Type:

NDArray[np.float64]

vol

Volume of the simulation box.

Type:

float

N

Total number of particles in the system.

Type:

int

compute() → None

Compute the radial distribution function.

This method calculates both the total RDF and partial RDFs for all pair combinations of atom types. The results are stored in the g_total, g, and r attributes.

plot(fig: Figure | None = None, ax: Axes | None = None) → Tuple

Plot the total (global) radial distribution function.

Parameters:

• fig (Optional[Figure]) – Existing matplotlib figure.

• ax (Optional[Axes]) – Existing matplotlib axes.

Returns:

• fig (matplotlib.figure.Figure) – The matplotlib figure object.

• ax (matplotlib.axes.Axes) – The matplotlib axes object.

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

Plot partial radial distribution functions for all atom type pairs.

Parameters:

• elements_list (list of str, optional) – List of element symbols corresponding to each atom type, e.g., [‘Al’, ‘Ni’]. If None, numeric labels are used. Length must match the number of atom types (Ntype).

• fig (Optional[Figure]) – Existing matplotlib figure.

• ax (Optional[Axes]) – Existing matplotlib axes.

Returns:

• fig (matplotlib.figure.Figure) – The matplotlib figure object.

• ax (matplotlib.axes.Axes) – The matplotlib axes object containing the partial RDF plot.

Raises:

AssertionError – If length of elements_list does not match Ntype.

mdapy.spatial_binning module

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

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

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

class mdapy.spline.Spline(x: List | Tuple | ndarray, y: List | Tuple | ndarray)

Bases: object

Cubic Spline Interpolation with clamped boundary conditions, where the spline’s endpoint derivatives are set equal to the original derivatives

This class provides a convenient interface for cubic spline interpolation. It supports evaluation of the spline function and its first derivative at arbitrary points within the interpolation range.

Parameters:

• x (array_like) – 1-D array of x-coordinates of data points. Must be strictly increasing and contain at least 2 points. Can be list, tuple, or numpy array.

• y (array_like) – 1-D array of y-coordinates of data points. Must have the same length as x. Can be list, tuple, or numpy array.

x

The x-coordinates of interpolation points (sorted and validated).

Type:

np.ndarray

y

The y-coordinates of interpolation points.

Type:

np.ndarray

Examples

>>> # Basic usage with lists
>>> x = [0, 1, 2, 3, 4]
>>> y = [0, 1, 4, 9, 16]
>>> sp = Spline(x, y)
>>>
>>> # Evaluate at a single point
>>> value = sp.evaluate(2.5)
>>> print(f"f(2.5) = {value}")
>>>
>>> # Evaluate at multiple points
>>> x_new = np.linspace(0, 4, 100)
>>> y_new = sp.evaluate(x_new)
>>>
>>> # Compute derivative
>>> derivative = sp.derivative(2.5)
>>> print(f"f'(2.5) = {derivative}")
>>>
>>> # Using numpy arrays
>>> x = np.array([0.0, 1.0, 2.0, 3.0])
>>> y = np.array([1.0, 2.0, 1.5, 3.0])
>>> sp = Spline(x, y)
>>> sp.evaluate([0.5, 1.5, 2.5])

Notes

• The spline is only valid within the range [x[0], x[-1]]. Attempting to evaluate outside this range will raise an assertion error.

• The spline uses natural boundary conditions (zero second derivatives at endpoints).

evaluate(x: float | int | List | Tuple | ndarray) → float | ndarray

Evaluate the spline at given point(s).

Computes the interpolated value of the spline function at the specified x-coordinate(s). Supports both scalar and array inputs.

Parameters:

x (float, int, list, tuple, or np.ndarray) – Point(s) at which to evaluate the spline. Must be within the range [self.x[0], self.x[-1]]. Can be: - A single number (int or float) - A list or tuple of numbers - A numpy array

Returns:

float or np.ndarray – The interpolated value(s) at x. Returns a float if input is scalar, otherwise returns a numpy array of the same shape as input.

Raises:

AssertionError – If any value in x is outside the interpolation range [self.x[0], self.x[-1]].

Examples

>>> sp = Spline([0, 1, 2, 3], [0, 1, 4, 9])
>>>
>>> # Scalar evaluation
>>> sp.evaluate(1.5)
2.125
>>>
>>> # List evaluation
>>> sp.evaluate([0.5, 1.5, 2.5])
array([0.375, 2.125, 6.125])
>>>
>>> # Numpy array evaluation
>>> x = np.linspace(0, 3, 10)
>>> y = sp.evaluate(x)

derivative(x: float | int | List | Tuple | ndarray) → float | ndarray

Evaluate the first derivative of the spline at given point(s).

Computes the first derivative (slope) of the interpolated spline function at the specified x-coordinate(s). Supports both scalar and array inputs.

Parameters:

x (float, int, list, tuple, or np.ndarray) – Point(s) at which to evaluate the derivative. Must be within the range [self.x[0], self.x[-1]]. Can be: - A single number (int or float) - A list or tuple of numbers - A numpy array

Returns:

float or np.ndarray – The derivative value(s) at x. Returns a float if input is scalar, otherwise returns a numpy array of the same shape as input.

Raises:

AssertionError – If any value in x is outside the interpolation range [self.x[0], self.x[-1]].

Examples

>>> sp = Spline([0, 1, 2, 3], [0, 1, 4, 9])
>>>
>>> # Scalar derivative
>>> sp.derivative(1.5)
3.5
>>>
>>> # Array derivative
>>> x = np.array([0.5, 1.5, 2.5])
>>> derivatives = sp.derivative(x)
>>> print(derivatives)

Notes

The derivative is computed analytically from the cubic spline coefficients, not by numerical differentiation, ensuring high accuracy.

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)

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

[1]

Paul J. Steinhardt. Bond-orientational order in liquids and glasses. Physical Review B, 28(2):784-805, 1983. doi:10.1103/PhysRevB.28.784.

[2]

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

[3]

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

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

[1]

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)

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()

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'] = 'direct')

Bases: object

Compute the static structure factor S(k) of an atomic system.

The static structure factor characterizes the spatial distribution of particles in reciprocal space and is related to the Fourier transform of the pair correlation function.

Two computational methods are available:

1. Direct method: Computes S(k) by directly evaluating the structure factor in reciprocal space using discrete k-space bins:

   \[S(\vec{k}) = \frac{1}{N} \sum_{i=0}^{N} \sum_{j=0}^N e^{i\vec{k} \cdot \vec{r}_{ij}}\]

2. Debye method: Uses the Debye scattering equation, which computes S(k) for specific k values. Better for large systems or when specific k values are needed.

   \[S(k) = \frac{1}{N} \sum_{i,j} \frac{\sin(k r_{ij})}{k r_{ij}}\]

where N is the number of particles, k is the wavenumber magnitude, and r_ij is the distance between particles i and j.

Parameters:

• data (pl.DataFrame) – Polars DataFrame containing atomic positions. Must have ‘x’, ‘y’, ‘z’ columns. For partial structure factors, must also contain ‘element’ or ‘type’ column.

• box (Box) – Simulation box object defining boundary conditions and box dimensions.

• k_min (float) – Minimum wavenumber for S(k) calculation (must be >= 0).

• k_max (float) – Maximum wavenumber for S(k) calculation (must be > k_min).

• nbins (int) – Number of bins for discretizing the k-space (must be > 0). - For ‘direct’ mode: creates nbins+1 bin edges, returns nbins values - For ‘debye’ mode: creates nbins linearly spaced k values

• cal_partial (bool, default=False) – If True, compute partial structure factors S_αβ(k) for each pair of atom types α and β. Requires ‘element’ or ‘type’ column in data.

• atomic_form_factors (bool, default=False) – If True, use atomic form factors f to weigh the atoms’ individual contributions to S(k). Atomic form factors are taken from TU Graz.

• mode ({'direct', 'debye'}, default='direct') – Computational method: - ‘direct’: Direct reciprocal space summation - ‘debye’: Debye scattering equation

k

Array of wavenumber values where S(k) is computed.

Type:

np.ndarray

Sk

Total static structure factor values.

Type:

np.ndarray

Sk_partial

Dictionary of partial structure factors (only if cal_partial=True). Keys are formatted as ‘element1-element2’.

Type:

dict[str, np.ndarray], optional

Sk_partial_xray

Dictionary of partial structure factors weighted by atomic_form (only if cal_partial=True and atomic_form_factors=True). This is helpful to directly compare with experimental results. Keys are formatted as ‘element1-element2’.

Type:

dict[str, np.ndarray], optional

Sk_xray

Total static structure factor values weighted by atomic_form (only if cal_partial=True and atomic_form_factors=True).

Type:

np.ndarray

References

[1]

Frenkel, D. & Smit, B. Understanding Molecular Simulation. Academic Press (2002)

[2]

Hansen, J.-P. & McDonald, I. R. Theory of Simple Liquids. Academic Press (2006)

[3]

freud.diffraction.StaticStructureFactorDirect https://freud.readthedocs.io/en/stable/modules/diffraction.html

compute() → None

Compute the static structure factor S(k).

This method performs the actual calculation of S(k) using the specified computational mode. Results are stored in the Sk attribute, and if cal_partial=True, partial structure factors are stored in Sk_partial.

Raises:

AssertionError – If required columns (‘x’, ‘y’, ‘z’) are missing from data, or if ‘element’ or ‘type’ column is missing when cal_partial=True.

Notes

For partial structure factors in a binary system with species A and B, three partial S(k) are computed: S_AA(k), S_AB(k), and S_BB(k). The total structure factor is:

\[S(k) = S_{AA}(k) + 2S_{AB}(k) + S_{BB}(k)\]

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

Plot the total static structure factor S(k).

Creates a line plot of S(k) versus wavenumber k with markers.

Parameters:

• fig (matplotlib.figure.Figure, optional) – Matplotlib figure object. If None, a new figure is created.

• ax (matplotlib.axes.Axes, optional) – Matplotlib axes object. If None, new axes are created.

Returns:

• fig (matplotlib.figure.Figure) – The matplotlib figure object containing the plot.

• ax (matplotlib.axes.Axes) – The matplotlib axes object containing the plot.

plot_partial(fig: Figure | None = None, ax: Axes | None = None) → Tuple[Figure, Axes]

Plot partial structure factors S_αβ(k) for multi-component systems.

Creates a line plot showing all partial structure factors computed for different atom type pairs. Each partial S(k) is plotted with a different color and labeled accordingly.

Parameters:

• fig (matplotlib.figure.Figure, optional) – Matplotlib figure object. If None, a new figure is created.

• ax (matplotlib.axes.Axes, optional) – Matplotlib axes object. If None, new axes are created.

Returns:

• fig (matplotlib.figure.Figure) – The matplotlib figure object containing the plot.

• ax (matplotlib.axes.Axes) – The matplotlib axes object containing the plot.

Raises:

AttributeError – If compute() was not called with cal_partial=True.

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

class mdapy.system.System(filename: str | None = None, data: DataFrame | None = None, pos: ndarray | None = None, box: int | float | Iterable[float] | ndarray | Box = None, ase_atom: Atoms | None = None, ovito_atom: DataCollection | None = None, format: str | None = None, global_info: Dict[str, Any] | None = {})

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.

data

Polars DataFrame containing per-atom properties. Always includes ‘x’, ‘y’, ‘z’ columns for positions.

Type:

pl.DataFrame

box

Simulation box object containing cell vectors and boundary conditions.

Type:

Box

global_info

Dictionary of global system properties.

Type:

dict

N

Number of atoms in the system.

Type:

int

calc

Attached calculator for energy/force/virial computations.

Type:

CalculatorMP or None

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")

See also

mdapy.box.Box

Simulation box representation

mdapy.calculator.CalculatorMP

Abstract calculator interface

mdapy.load_save.BuildSystem

System construction utilities

mdapy.load_save.SaveSystem

System export utilities

property global_info: Dict[str, Any]

Obtain global system information.

Returns:

dict – Dictionary containing global properties such as total energy, virial tensor, stress tensor, timestep, temperature, etc. Keys and values depend on the data source and attached calculator.

property data: DataFrame

Access the atomic data DataFrame.

Returns:

pl.DataFrame – Polars DataFrame containing per-atom properties. Always includes columns ‘x’, ‘y’, ‘z’ for Cartesian coordinates. May include additional columns such as ‘type’, ‘element’, ‘vx’, ‘vy’, ‘vz’ (velocities), ‘fx’, ‘fy’, ‘fz’ (forces), and various computed properties from analysis methods.

Examples

>>> system = System("config.dump")
>>> print(system.data.columns)
['x', 'y', 'z', 'type', 'vx', 'vy', 'vz']

property N: int

Get the number of atoms in the system.

Returns:

int – Total number of atoms/particles in the system.

Examples

>>> system = System(pos=np.random.rand(100, 3), box=10)
>>> print(system.N)
100

set_element(element: str | List[str] | Tuple[str] | ndarray)

Set element names for atoms in the system.

Parameters:

element (str, list of str, tuple of str, or np.ndarray) –

Element specification:

• str: Assigns the same element to all atoms

• list/tuple/array: Must have length equal to N, specifying element for each atom individually

Examples

Assign same element to all atoms:

>>> system.set_element("Cu")

Assign different elements:

>>> elements = ["Cu"] * 50 + ["Al"] * 50
>>> system.set_element(elements)

Notes

This method updates the ‘element’ column in the data DataFrame. If an ‘element’ column exists, it is replaced.

set_type_by_element(element_list: List[str] | Tuple[str] | ndarray)

Set atom types based on element ordering.

Parameters:

element_list (list of str, tuple of str, or np.ndarray) – Ordered list of unique elements. Atoms are assigned type numbers based on the index of their element in this list (starting from 1).

Raises:

AssertionError – If data does not contain ‘element’ column or if any element in data is not present in the provided element list.

Examples

>>> system.set_element(["Cu", "Cu", "Al", "Al"])
>>> system.set_type_by_element(["Cu", "Al"])
>>> print(system.data.select("element", "type"))
┌─────────┬──────┐
│ element │ type │
├─────────┼──────┤
│ Cu      │ 1    │
│ Cu      │ 1    │
│ Al      │ 2    │
│ Al      │ 2    │
└─────────┴──────┘

Notes

This method creates or updates the ‘type’ column in the data DataFrame. Type numbers start from 1 and correspond to the order in the element list.

get_positions(reduced: bool = False) → DataFrame

Extract atomic positions from the system.

Parameters:

reduced (bool, optional) – If True, return fractional (reduced) coordinates in the box basis. If False (default), return Cartesian coordinates.

Returns:

pl.DataFrame – DataFrame containing position columns:

• reduced=False: Columns ‘x’, ‘y’, ‘z’ (Cartesian coordinates)

• reduced=True: Columns ‘r_x’, ‘r_y’, ‘r_z’ (fractional coordinates)

Examples

>>> system = System("POSCAR")
>>> cart_pos = system.get_positions(reduced=False)
>>> frac_pos = system.get_positions(reduced=True)

Notes

Fractional coordinates are computed using the inverse box matrix: reduced = positions @ box.inverse_box

get_velocities() → DataFrame

Extract atomic velocities from the system.

Returns:

pl.DataFrame – DataFrame with columns ‘vx’, ‘vy’, ‘vz’ containing velocity components for each atom.

Raises:

AssertionError – If the data does not contain velocity columns (‘vx’, ‘vy’, ‘vz’).

Examples

>>> system = System("relax.dump")  # File with velocities
>>> velocities = system.get_velocities()

set_pka(energy: float, direction: ndarray, index: int | None = None, element: str | None = None, factor: float = 1.0)

Set primary knock-on atom (PKA) for radiation damage simulation.

This method assigns initial kinetic energy and direction to a PKA atom, commonly used to initialize cascade simulations in radiation damage studies. We assume that the velocity in the units of A/fs, you need to use factor to convert your velocity to this units.

Parameters:

• energy (float) – Kinetic energy to assign to the PKA (in eV).

• direction (np.ndarray) – 1D array of length 3 specifying the velocity direction vector. Will be normalized automatically.

• index (int, optional) – Index of the atom to set as PKA. If None, selects by nearest center atom.

• element (str, optional) – Element name to select PKA. If provided and index is None, selects by nearest center atom with same element.

• factor (float) – Convert your velocity to the units of A/fs. Defaults to 1.0.

Notes

This method updates the velocity columns (‘vx’, ‘vy’, ‘vz’) in the data DataFrame.

Examples

Set PKA by atom index:

>>> system.set_pka(energy=1000.0, direction=np.array([1, 3, 5]), index=50)

Set PKA by element:

>>> system.set_pka(energy=2000.0, direction=np.array([1, 3, 5]), element="Fe")

get_energy() → float

Calculate total potential energy of the system.

Returns:

float – Total potential energy in the units eV.

Raises:

AssertionError – If no calculator is attached to the system.

Examples

>>> from mdapy.eam import EAM
>>> system = System("config.dump")
>>> system.calc = EAM("potential.eam")
>>> total_energy = system.get_energy()

get_energies() → ndarray

Calculate per-atom potential energies in the units eV.

Returns:

np.ndarray – 1D array of shape (N,) containing potential energy for each atom.

Raises:

AssertionError – If no calculator is attached to the system.

get_force() → ndarray

Calculate forces acting on each atom in the units eV/A.

Returns:

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

Raises:

AssertionError – If no calculator is attached to the system.

Examples

>>> from mdapy.eam import EAM
>>> system = System("config.dump")
>>> system.calc = EAM("potential.eam")
>>> forces = system.get_force()
>>> print(forces.shape)
(1000, 3)

get_stress() → ndarray

Calculate stress tensor of the system in the units eV/A^3.

Returns:

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

Raises:

AssertionError – If no calculator is attached to the system.

Examples

>>> from mdapy.eam import EAM
>>> system = System("config.dump")
>>> system.calc = EAM("potential.eam")
>>> stress = system.get_stress()
>>> print(stress)  # [σ_xx, σ_yy, σ_zz, σ_yz, σ_xz, σ_xy]

get_virials() → ndarray

Calculate per-atom virial tensors in the units eV.

Returns:

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

Raises:

AssertionError – If no calculator is attached to the system.

Examples

>>> from mdapy.eam import EAM
>>> system = System("config.dump")
>>> system.calc = EAM("potential.eam")
>>> virials = system.get_virials()
>>> print(virials.shape)
(1000, 9)

update_data(data: DataFrame, reset_calcolator: bool = False, reset_neighbor: bool = False) → None

Update the atomic data DataFrame.

Parameters:

• data (pl.DataFrame) – New DataFrame containing updated atomic information. Should maintain the same structure (columns) as the original data.

• reset_calcolator (bool, optional) – If True, clears cached results from the attached calculator. Set to True when atomic positions or types/element change. Default is False.

• reset_neighbor (bool, optional) – If True, clear all neighbor information. Set to True when atomic positions change. Default is False.

Examples

>>> # Add a new property column
>>> new_data = system.data.with_columns(
...     temperature=np.random.rand(system.N) * 300
... )
>>> system.update_data(new_data)

>>> # Update positions (reset calculator)
>>> new_data = system.data.with_columns(x=system.data["x"] + 0.1)
>>> system.update_data(new_data, reset_calcolator=True)

Notes

When reset_calcolator=True, any cached energy, force, or stress calculations are invalidated and will be recomputed on next access.

update_box(box: int | float | Iterable[float] | ndarray | Box, scale_pos: bool = False) → None

Update the simulation box.

Parameters:

• box (int, float, Iterable[float], np.ndarray, or Box) – New box specification. See Box for accepted formats.

• scale_pos (bool, optional) – If True, scales atomic positions affinely with the box change. Useful for applying strain or changing lattice constants while maintaining fractional coordinates. Default is False.

Raises:

AssertionError – If scale_pos=True but not all boundaries are periodic.

Examples

Change box without scaling positions:

>>> system.update_box([20, 20, 20])

Apply uniform strain by scaling box and positions:

>>> system.update_box([11, 11, 11], scale_pos=True)

Notes

When scale_pos=True, the transformation is: new_pos = old_pos @ (old_box^-1 @ new_box)

This preserves fractional coordinates in the new box.

wrap_pos() → None

Wrap positions into box for PBC boundary.

replicate(nx: int, ny: int, nz: int) → None

Replicate the system by creating a supercell.

Parameters:

• nx (int) – Number of replications along x-direction.

• ny (int) – Number of replications along y-direction.

• nz (int) – Number of replications along z-direction.

Examples

Create a 2x2x2 supercell:

>>> system = System("unit_cell.poscar")
>>> system.replicate(2, 2, 2)
>>> print(system.N)  # 8 times the original number

Notes

This method modifies the system in place. The box is scaled by [nx, ny, nz] and atoms are replicated accordingly.

to_ovito() → DataCollection

Convert system to OVITO DataCollection object.

Returns:

ovito.data.DataCollection – OVITO data structure containing the atomic configuration.

Raises:

ImportError – If OVITO is not installed.

Examples

>>> data = system.to_ovito()
>>> # Use OVITO's analysis or visualization tools
>>> from ovito.modifiers import CommonNeighborAnalysisModifier
>>> data.apply(CommonNeighborAnalysisModifier())

to_ase() → Atoms

Convert system to ASE Atoms object.

Returns:

ase.Atoms – ASE Atoms object containing the atomic configuration.

Raises:

• ImportError – If ASE is not installed.

• AssertionError – If data does not contain an ‘element’ column.

Examples

>>> atoms = system.to_ase()
>>> # Use ASE's functionality
>>> from ase.io import write
>>> write("output.xyz", atoms)

write_mp(output_name: str) → None

Save system to mp format file.

Parameters:

• nx (int) – Number of replications along x-direction.

• ny (int) – Number of replications along y-direction.

• nz (int) – Number of replications along z-direction.

Raises:

• ImportError – If ASE is not installed.

• AssertionError – If data does not contain an ‘element’ column.

Examples

>>> atoms = system.to_ase()
>>> # Use ASE's functionality
>>> from ase.io import write
>>> write("output.xyz", atoms)

write_xyz(output_name: str, classical: bool = False, compress=False)

write_poscar(output_name: str, reduced_pos: bool = False, selective_dynamics: bool = False, compress=False)

write_data(output_name: str, element_list: List[str] | None = None, num_type: int | None = None, data_format: str = 'atomic', compress=False)

write_dump(output_name: str, timestep: float = 0, compress=False)

build_neighbor(rc: float, max_neigh: int | None = None) → None

Build neighbor list for the system.

This method constructs a Verlet neighbor list within a cutoff radius, which is required for many analysis methods. For small systems, the box may be automatically replicated to ensure sufficient neighbors.

Parameters:

• rc (float) – Cutoff radius for neighbor search (in simulation units).

• max_neigh (int, optional) – Maximum number of neighbors per atom. If None, automatically determined based on the actual neighbor number Increase this if you get warnings about insufficient neighbor slots.

• Set (Attributes)

• --------------

• verlet_list (np.ndarray) – 2D array of shape (N, max_neigh) containing neighbor indices. Unfilled slots contain -1.

• distance_list (np.ndarray) – 2D array of shape (N, max_neigh) containing distances to neighbors.

• neighbor_number (np.ndarray) – 1D array of shape (N,) containing number of neighbors for each atom.

• rc – The cutoff radius used.

• _enlarge_box (Box, optional) – replicated box for small system.

• _enlarge_data (pl.DataFrame, optional) – replicated atom dataframe for small system

Examples

>>> system.build_neighbor(rc=5.0)
>>> print(system.verlet_list.shape)
(1000, 50)  # 1000 atoms, up to 50 neighbors each
>>> print(system.neighbor_number)
[12 12 12 ... 12 12 12]  # FCC coordination

See Neighbor for implementation details.

build_voronoi_neighbor(a_face_area_threshold: float = -1.0, r_face_area_threshold: float = -1.0) → None

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

• Set (Attributes)

• --------------

• voro_verlet_list (np.ndarray) – 2D array of shape (N, max_neigh) containing neighbor indices. Unfilled slots contain -1.

• voro_distance_list (np.ndarray) – 2D array of shape (N, max_neigh) containing distances to neighbors.

• voro_neighbor_number (np.ndarray) – 1D array of shape (N,) containing number of neighbors for each atom.

• voro_face_area (np.ndarray) – 2D array of shape (N, max_neighbors) containing the area of the Voronoi face shared with each neighbor.

• _enlarge_box (Box, optional) – replicated box for small system.

• _enlarge_data (pl.DataFrame, optional) – replicated atom dataframe for small system

Examples

>>> system.build_voronoi_neighbor()
>>> print(system.voro_verlet_list.shape)
(1000, 50)  # 1000 atoms, up to 50 neighbors each

See Voronoi for implementation details.

build_nearest_neighbor(k: int) → None

Calculate k nearest neighbor information.

Parameters:

• k (int) – Maximum number of neighbors to find per atom.

• Set (Attributes)

• --------------

• verlet_list (np.ndarray) – 2D array of shape (N, k) containing neighbor indices.

• distance_list (np.ndarray) – 2D array of shape (N, k) containing distances to neighbors.

• neighbor_number (np.ndarray) – 1D array of shape (N,) containing number of neighbors for each atom.

• _enlarge_box (Box, optional) – replicated box for small system.

• _enlarge_data (pl.DataFrame, optional) – replicated atom dataframe for small system

Examples

>>> system.build_nearest_neighbor(12)
>>> print(system.verlet_list.shape)
(1000, 12)  # 1000 atoms, up to 12 sorted neighbors each

See NearestNeighbor for implementation details.

cal_identify_diamond_structure() → None

Identify diamond structure atoms. This will generate ‘ids’ column in self.data:

• 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)

Notes

See IdentifyDiamondStructure for implementation details.

cal_common_neighbor_parameter(rc: float, max_neigh: int | None = None) → None

Calculate common neighbor parameter for structure analysis. This will save results to self.data[‘cnp’].

Parameters:

• rc (float) – Cutoff radius for neighbor search.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Notes

See CommonNeighborParameter for implementation details.

cal_ackland_jones_analysis() → None

Perform Ackland-Jones structure analysis. This will save results to self.data[‘aja’]:

• 0 = Other (unknown coordination)

• 1 = FCC (face-centered cubic)

• 2 = HCP (hexagonal close-packed)

• 3 = BCC (body-centered cubic)

• 4 = ICO (icosahedral coordination)

Notes

See AcklandJonesAnalysis for implementation details.

cal_warren_cowley_parameter(rc: float, max_neigh: int | None = None) → WarrenCowleyParameter

Calculate Warren-Cowley short-range order parameter.

Parameters:

• rc (float) – Cutoff radius for neighbor search.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Returns:

WarrenCowleyParameter – Object containing Warren-Cowley parameters.

Notes

See WarrenCowleyParameter for implementation details and returned attributes.

Examples

>>> wcp = system.cal_warren_cowley_parameter(rc=3.0)
>>> print(wcp.wcp)  # Access Warren-Cowley parameters

cal_atomic_temperature(rc: float, factor: float = 1.0, max_neigh: int | None = None) → None

Calculate the local atomic temperature for each atom by analyzing the velocity fluctuations of the atom and its neighbors. We assume velocity in units of Å/fs, you can use factor to convert your velocity. This will save results to self.data[‘atomic_temp’] in the unit of K.

Parameters:

• rc (float) – Cutoff radius for local averaging.

• factor (float) – Scaling factor for velocities (e.g., for unit conversion). Default is 1.0.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Notes

See AtomicTemperature for implementation details.

cal_steinhardt_bond_orientation(llist: ndarray | List[int], use_voronoi: bool = False, nnn: int = 0, rc: float = -1.0, average: bool = False, use_weight: bool = False, weight: ndarray | None = None, wl: bool = False, wlhat: bool = False, a_face_area_threshold: float = -1, r_face_area_threshold: float = -1, identify_liquid: bool = False, threshold: float = 0.7, n_bond: int = 7, max_neigh: int | None = None) → None

Calculate Steinhardt bond orientation order parameters. One can also use it to identify the solid or liquid atom.

Parameters:

• llist (np.ndarry or List[int]) – Degrees list, such as [4, 6, 8], should be positive, even integer.

• use_voronoi (bool) – If True, compute using Voronoi neighbor. Default is False.

• nnn (int) – If use_voronoi is False and this value is postive, compute using nnn nearest neighbors. Default is 0.

• rc (float) – If use_voronoi is False and nnn is zero, make sure this parameter is postive, and use it to build neighbor. Default is -1.0.

• average (bool) – If True, qlm will be averaged. Default is False.

• use_weight (bool) – If True, the neighbor will be weighed by weight array, otherwise the weight of each neighbor is 1.0. Default is False.

• weight (np.ndarray, optional) – If it is not None, it should be has the same shape with the self.verlet_list. If it is None and the use_voronoi is True, then the voronoi_face_area will be used as weight. Default is None.

• wl (bool) – If True, compute third-order invariant \(w_l\) parameters. Default is False.

• wlhat (bool) – If True, compute normalized third-order invariant \(\hat{w}_l\) parameters. Default is False.

• a_face_area_threshold (float) – 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) – Relative face area threshold (fraction of average face area). Faces with relative area below this value are ignored. Default is -1.0 (no filtering).

• 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. Default to 7.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Notes

See SteinhardtBondOrientation for implementation details. The results will add to self.data[‘ql{l}’] for l in llist. If set wl is True, also will add self.data[‘wl{l}’]. If set wlhat is True, also will add self.data[‘wlh{l}’]. If set identify_liquid, it will add self.data[‘solidliquid’] and self.data[‘nbond’] results.

cal_polyhedral_template_matching(structure='fcc-hcp-bcc', rmsd_threshold=0.1, return_ordering=False, return_rmsd=False, return_atomic_distance=False, return_orientation=False, identify_fcc_planar_faults=False, identify_esf=True)

Identify crystal structures using polyhedral template matching, one can also further distiguish the planar defects in FCC structures.

For PTM, structure results will add as self.data[‘ptm’] and self.data[‘ordering’]

• Structure type (integer, 0=Other, 1=FCC, 2=HCP, 3=BCC, 4=ICO, 5=SC, 6=DCUB, 7=DHEX, 8=Graphene)

• Ordering type (interger, 0=Other, 1=L10, 2=L12 (A-site), 3=L12 (B-site), 4=B2, 5=zincblende / wurtzite)

For planar defects, structure results will add as self.data[‘pft’]

• 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)

Parameters:

• structure (str) – String specifying the structure types to consider, separated by hyphens (Defaults is “fcc-hcp-bcc”). Supported values: “fcc”, “hcp”, “bcc”, “ico”, “sc”, “dcub”, “dhex”, “graphene”. Special values: “all” (all types), “default” (fcc, hcp, bcc).

• rmsd_threshold (float) – Maximum RMSD for a valid structure match. Particles exceeding this are classified as “Other”. Default is 0.1.

• return_ordering (bool) – If True, return the ordering type. Default is False.

• return_rmsd (bool) – If True, return the RMSD value. Default is False.

• return_atomic_distance (bool) – If True, return the interatomic distance. Default is False.

• return_orientation (bool) – If True, return the orientation value. Default is False.

• identify_fcc_planar_faults (bool) – If True, further detect the plannar defects. Default is False.

• identify_esf (bool) – If True, further defect the ESF. Default is True.

Notes

See PolyhedralTemplateMatching and class:~mdapy.identify_fcc_planar_faults.IdentifyFccPlanarFaults for implementation details.

cal_centro_symmetry_parameter(N: int)

Calculate centro-symmetry parameter for defect identification. The results will be added as self.data[‘csp’].

Parameters:

N (int) – Number of nearest neighbors to consider, should be a positive, even integer. 12 for FCC and 8 for BCC.

Notes

See CentroSymmetryParameter for implementation details.

cal_common_neighbor_analysis(rc: float | None = None, max_neigh: int | None = None)

Perform common neighbor analysis for structure identification. This will generate ‘cna’ column in self.data:

• 0 = Other (unknown coordination)

• 1 = FCC (face-centered cubic)

• 2 = HCP (hexagonal close-packed)

• 3 = BCC (body-centered cubic)

• 4 = ICO (icosahedral coordination)

Parameters:

• rc (float, optional) – Cutoff radius. If not given, it will use adaptive mode to do CNA.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Notes

See CommonNeighborAnalysis for implementation details.

cal_structure_factor(k_min: float, k_max: float, nbins: int, cal_partial: bool = False, atomic_form_factors: bool = False, mode: str = 'direct') → StructureFactor

Calculate static structure factor S(k).

Parameters:

• k_min (float) – Minimum wave vector magnitude.

• k_max (float) – Maximum wave vector magnitude.

• nbins (int) – Number of bins for k-space averaging.

• cal_partial (bool, optional) – If True, calculate partial structure factors. Default is False.

• atomic_form_factors (bool, default=False) –

  If True, use atomic form factors f to weigh the atoms’ individual contributions to S(k). Atomic form factors are taken from TU Graz.

• mode (str, optional) – Calculation mode: ‘direct’ or ‘debye’. Default is ‘direct’.

Returns:

StructureFactor – Object containing structure factor data.

Notes

See StructureFactor for implementation details and returned attributes.

cal_bond_analysis(rc: float, nbin: int, max_neigh: int | None = None) → BondAnalysis

Analyze bond length and angle distributions.

Parameters:

• rc (float) – Cutoff radius for bond search.

• nbin (int) – Number of bins for histograms.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Returns:

BondAnalysis – Object containing bond distribution data.

Notes

See BondAnalysis for implementation details and returned attributes.

Examples

>>> ba = system.cal_bond_analysis(rc=5.0, nbin=100)
>>> print(ba.bond_length_distribution)

cal_angular_distribution_function(rc_dict: Dict[str, List[float]], nbin: int, max_neigh: int | None = None) → AngularDistributionFunction

Calculate angular distribution function (ADF).

Parameters:

• 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 angular bins.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Returns:

AngularDistributionFunction – Object containing ADF data.

Notes

See AngularDistributionFunction for implementation details and returned attributes.

cal_radial_distribution_function(rc: float, nbin: int = 100, max_neigh: int | None = None) → RadialDistributionFunction

Calculate radial distribution function g(r).

Parameters:

• rc (float) – Maximum distance for RDF calculation.

• nbin (int, optional) – Number of distance bins. Default is 100.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Returns:

RadialDistributionFunction – Object containing RDF data and partial RDFs if applicable.

Notes

See RadialDistributionFunction for implementation details and returned attributes.

average_by_neighbor(average_rc: float, property_name: str, include_self: bool = True, output_name: str | None = None, max_neigh: int | None = None) → None

Average a property over local neighborhoods.

Parameters:

• average_rc (float) – Cutoff radius for local averaging.

• property_name (str) – Name of the column in data to average.

• include_self (bool) – Whether to include the central atom in averaging. Default is True.

• output_name (str, optional) – Name for the output column. If None, uses ‘{property_name}_ave’.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Raises:

AssertionError – If property_name is not a column in data.

cal_cluster_analysis(rc: float | int | Dict[str, float] = 5.0, max_neigh: int | None = None) → None

Perform cluster analysis to identify connected atomic groups. The results will be added as self.data[‘cluster_id’].

Parameters:

• rc (float, int, or dict, optional) –

  Cutoff radius for cluster connectivity. Can be:

  • float/int: Single cutoff for all atom pairs

  • dict: Type-specific cutoffs, e.g., {‘1-1’: 1.5, ‘1-2’: 1.3}

  Default is 5.0.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Notes

See ClusterAnalysis for implementation details.

Examples

Single cutoff:

>>> system.cal_cluster_analysis(rc=3.0)

Type-specific cutoffs:

>>> rc_dict = {"1-1": 2.8, "1-2": 3.0, "2-2": 3.2}
>>> system.cal_cluster_analysis(rc=rc_dict)
>>> # cluster_id column added to system.data

cal_structure_entropy(rc: float, sigma: float, use_local_density: bool = False, average_rc: float = 0.0, max_neigh: int | None = None) → None

Calculate structural entropy for disorder quantification.

Parameters:

• rc (float) – Cutoff radius for neighbor search.

• sigma (float) – Width parameter for Gaussian kernel.

• use_local_density (bool, optional) – If True, use local density normalization. Default is False.

• average_rc (float, optional) – Radius for local averaging of entropy. If 0, no averaging. Default is 0.0.

• max_neigh (int, optional) – Maximum number of neighbors per atom.

Notes

See StructureEntropy for implementation details.

Examples

>>> system.cal_structure_entropy(rc=5.0, sigma=0.5)
>>> # entropy column added to system.data

With local averaging:

>>> system.cal_structure_entropy(rc=5.0, sigma=0.5, average_rc=3.0)
>>> # entropy and entropy_ave columns added to system.data

cal_voronoi_volume() → None

Calculate Voronoi cell volumes and related properties.

This method computes the Voronoi tessellation and adds three columns to the data DataFrame:

• volume: Atomic volume from Voronoi cell

• neighbor_number: Coordination number from Voronoi neighbors

• cavity_radius: 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.

Notes

See Voronoi for implementation details.

Examples

>>> system.cal_voronoi_volume()
>>> # volume, neighbor_number, cavity_radius columns added to data
>>> avg_volume = system.data["volume"].mean()

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

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)

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

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]

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)

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)

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)

Bases: object

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') → None

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

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

append(system: System) → None

Append a frame to the trajectory.

Parameters:

system (System) – System object to append

Raises:

TypeError – If system is not a System object

extend(systems: List[System]) → None

Extend trajectory with multiple frames.

Parameters:

systems (List[System]) – List of System objects to append

Raises:

TypeError – If any element is not a System object

insert(index: int, system: System) → None

Insert a frame at specified position.

Parameters:

• index (int) – Position to insert at

• system (System) – System object to insert

Raises:

TypeError – If system is not a System object

pop(index: int = -1) → System

Remove and return frame at specified position.

Parameters:

index (int, default=-1) – Position to pop from (default is last frame)

Returns:

System – The removed System object

remove(indices: int) → None

Remove frame at specified index.

Parameters:

indices (int) – Frame index to remove

get_atoms_count() → List[int]

Get atom count for each frame.

Returns:

List[int] – List of atom counts

concatenate(other: XYZTrajectory) → XYZTrajectory

Concatenate two trajectories.

Parameters:

other (XYZTrajectory) – Another trajectory to concatenate

Returns:

XYZTrajectory – New trajectory containing frames from both trajectories

Examples

>>> traj1 = XYZTrajectory("file1.xyz")
>>> traj2 = XYZTrajectory("file2.xyz")
>>> combined = traj1.concatenate(traj2)

mdapy.visualize module

class mdapy.visualize.View(system: System)

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

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

Color atoms using their element name.

Notes

Element symbols are mapped to integer colors according to ele_dict.

init_plot() → None

Initialize the visualization canvas.

Notes

• Creates box lines, atomic points.

• Initializes default color and radius.

• Creates element/type legends using text2d.

display() → None

Display the k3d plot in supported environments (e.g., Jupyter).

close() → None

Close the k3d plot and release the rendering canvas.

hide_object_by_group_name(name: str, remove: bool = False) → None

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

Remove existing colorbar from the plot.

colored_by(name: str, vmin: float | None = None, vmax: float | None = None, cmap: str = 'rainbow') → None

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)

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()

Perform the void detection.

1. Compute the number of grid cells ncell along each dimension based on the system thickness and cutoff rc.

2. Use _neighbor._fill_cell_for_void to assign each cell a value: 0 for empty, 1 for occupied.

3. If empty cells exist, compute their geometric centers and create a System object containing void positions.

4. Perform cluster analysis on these void positions.

5. Filter clusters that contain only one void cell (noise avoidance).

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

[1]

Lu J, Lazar E A, Rycroft C H. An extension to Voro++ for multithreaded computation of Voronoi cells[J]. Computer Physics Communications, 2023, 291: 108832. https://doi.org/10.1016/j.cpc.2023.108832

class mdapy.voronoi.Voronoi(box: Box, data: DataFrame)

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]

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

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]

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)

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)

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)

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

[1]

Warren, B. E. (1990). X-ray Diffraction. Dover Publications.

[2]

Cowley, J. M. (1950). An approximate theory of order in alloys. Physical Review, 77(5), 669.

compute() → None

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]

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)

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]

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