API Reference#

All public functions are available directly from import pyscal (or, equivalently, import pyscal3). Results are stored on the ASE Atoms object with the pyscal_ prefix.

Neighbor Finding#

pyscal3.find_neighbors(atoms: Atoms, method='cutoff', cutoff=0, shell_thickness=0, threshold=2, voroexp=1, padding=1.2, nlimit=6, cells=None, nmax=12, assign_neighbor=True)#

Find neighbors of all atoms.

Parameters:

  • atoms (ase.Atoms) – The atomic structure.

  • method ({'cutoff', 'voronoi', 'number'}) – Neighbor finding algorithm.

  • cutoff (float or str) – Cutoff distance. Use ‘sann’ or ‘adaptive’ for adaptive methods. 0 defaults to adaptive.

  • shell_thickness (float, optional) – If > 0, find neighbors in a shell [cutoff, cutoff+shell_thickness].

  • threshold (float, optional) – Safety multiplier for adaptive methods. Default 2.

  • voroexp (int, optional) – Voronoi face area weight exponent. Default 1.

  • padding (float, optional) – Safety padding for adaptive/number methods. Default 1.2.

  • nlimit (int, optional) – Number of atoms for adaptive cutoff estimation. Default 6.

  • cells (bool or None, optional) – Force cell lists on/off. None = auto (>250 atoms).

  • nmax (int, optional) – Number of neighbors for ‘number’ method. Default 12.

  • assign_neighbor (bool, optional) – Whether to assign neighbors (for ‘number’ method). Default True.

Returns:

Results are stored in-place on the atoms object:

  • atoms.arrays["pyscal_neighbors"] — neighbor indices (natoms, max_neighbors)

  • atoms.arrays["pyscal_neighbordist"] — neighbor distances

  • atoms.arrays["pyscal_theta"] — polar angles to neighbors

  • atoms.arrays["pyscal_phi"] — azimuthal angles to neighbors

  • atoms.info["pyscal_neighbors_found"] — set to True

  • atoms.info["pyscal_neighbor_method"] — the method used

For Voronoi, additional keys include face_vertices, face_perimeters, face_areas, vertex_vectors, and voronoivol.

Return type:

None

pyscal3.get_distance(atoms: Atoms, pos1, pos2, vector=False)#

Get the distance between two positions respecting periodic boundaries.

Parameters:

  • atoms (ase.Atoms) – Structure (used for box/PBC info).

  • pos1 (array-like) – Positions.

  • pos2 (array-like) – Positions.

  • vector (bool, optional) – If True, also return the displacement vector.

Returns:

Distance, and optionally the displacement vector.

Return type:

float or (float, list)

Structural Descriptors#

Steinhardt Parameters#

pyscal3.steinhardt_parameter(atoms: Atoms, l, averaged=False)#

Calculate Steinhardt bond order parameters q_l.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed.

  • l (int or list of int) – Steinhardt parameter order(s), e.g. 6 or [4, 6].

  • averaged (bool, optional) – If True, compute neighbor-averaged q values. Default False.

Returns:

One array per requested l value, each of shape (natoms,).

Return type:

list of numpy arrays

Wigner \(W_l\) Parameters#

pyscal3.wigner_w_parameter(atoms: Atoms, l, averaged=False, normalized=True)#

Calculate the third-order Steinhardt invariant W_l.

W_l is the third-order rotational invariant of the bond-orientational order parameters, constructed by contracting q_lm with Wigner 3j symbols. It distinguishes crystal structures that have similar q_l values (e.g., FCC vs HCP via the sign of W_6).

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed via find_neighbors.

  • l (int or list of int) – Order(s) for the W parameter. Only even values give nonzero results (W_l = 0 for odd l when j1=j2=j3=l).

  • averaged (bool, optional) – If True, compute neighbor-averaged W_l (Lechner-Dellago). Default False.

  • normalized (bool, optional) – If True (default), return the normalized hat{W}_l = W_l / (sum |q_lm|^2)^(3/2). If False, return raw W_l.

Returns:

One array per requested l value, each of shape (natoms,).

Return type:

list of numpy arrays

References

Notes

Known values for hat{W}_6:

  • FCC: −0.01316

  • HCP: −0.01244

  • BCC: +0.01316

  • ICO (Mackay): −0.16975

  • Liquid: ≈ 0

Minkowski Structure Metrics#

pyscal3.minkowski_parameter(atoms: Atoms, l, voroexp=1, averaged=False)#

Calculate Minkowski structure metrics \(q_l^{\mathrm{Mink}}\).

These are Voronoi face-area weighted Steinhardt parameters (Mickel et al., J. Chem. Phys. 138, 044501, 2013). For each atom i and angular-momentum order l the metric is

\[q_l^{\mathrm{Mink}}(i) = \sqrt{\frac{4\pi}{2l+1} \sum_{m=-l}^{l} \left| \sum_j \frac{A_j^\alpha}{\sum_k A_k^\alpha} Y_{lm}(\hat{\mathbf{r}}_{ij}) \right|^2}\]

where A_j is the Voronoi face area between atoms i and j, and α is the exponent voroexp.

This function performs Voronoi neighbor finding internally, so there is no need to call find_neighbors() beforehand (any existing neighbor data is overwritten).

Parameters:

  • atoms (ase.Atoms) – The atomic structure (periodic boundary conditions recommended).

  • l (int or list of int) – Steinhardt parameter order(s), e.g. 6 or [4, 6].

  • voroexp (int or float, optional) – Face-area weight exponent α. Default 1.

  • averaged (bool, optional) – If True, compute neighbor-averaged values. Default False.

Returns:

One array per requested l, each of shape (natoms,). Values are also stored in atoms.arrays["pyscal_q{l}"] (or "pyscal_avg_q{l}" when averaged=True).

Return type:

list of numpy arrays

Notes

Unlike a plain steinhardt_parameter call after Voronoi neighbor finding, this function guarantees that Voronoi neighbors are (re)computed with the requested voroexp so results are reproducible regardless of prior state.

References

W. Mickel, S. C. Kapfer, G. E. Schröder-Turk and K. Mecke, “Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter”, J. Chem. Phys. 138, 044501 (2013). doi:10.1063/1.4774084

Disorder Parameter#

pyscal3.disorder(atoms: Atoms, q=6, averaged=False)#

Calculate the disorder parameter.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed.

  • q (int, optional) – Steinhardt parameter order for disorder calc. Default 6.

  • averaged (bool, optional) – If True, average disorder over neighbors. Default False.

Returns:

Per-atom disorder values.

Return type:

numpy array

Common Neighbor Analysis#

pyscal3.common_neighbor_analysis(atoms: Atoms, lattice_constant=None)#

Calculate Common Neighbor Analysis (CNA) or Adaptive CNA.

Parameters:

  • atoms (ase.Atoms) – Structure.

  • lattice_constant (float, optional) – If given, use conventional CNA. If None, use adaptive CNA.

Returns:

Counts: {“fcc”: n, “hcp”: n, “bcc”: n, “ico”: n, “others”: n}

Return type:

dict

Diamond Structure#

pyscal3.diamond_structure(atoms: Atoms)#

Identify diamond structure using extended CNA.

Parameters:

atoms (ase.Atoms) – Structure.

Returns:

Counts per structure type.

Return type:

dict

Centrosymmetry Parameter#

pyscal3.centrosymmetry(atoms: Atoms, nmax=12)#

Calculate the centrosymmetry parameter.

Parameters:

  • atoms (ase.Atoms) – Structure.

  • nmax (int, optional) – Number of neighbors (must be positive even integer). Default 12.

Returns:

Per-atom centrosymmetry values.

Return type:

numpy array

Voronoi Vector#

pyscal3.voronoi_vector(atoms: Atoms, edge_cutoff=0.05, area_cutoff=0.01)#

Calculate the Voronoi structure identification vector (n3, n4, n5, n6).

Parameters:

  • atoms (ase.Atoms) – Structure (must have Voronoi neighbors computed).

  • edge_cutoff (float, optional) – Minimum edge length fraction. Default 0.05.

  • area_cutoff (float, optional) – Minimum face area fraction. Default 0.01.

Returns:

Voronoi vectors [n3, n4, n5, n6] per atom.

Return type:

numpy array of shape (natoms, 4)

Entropy Parameter#

pyscal3.entropy(atoms: Atoms, rm, sigma=0.2, rstart=0.001, h=0.001, local=False, average=False)#

Calculate the entropy parameter for each atom.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors computed.

  • rm (float) – Cutoff distance for integration.

  • sigma (float, optional) – Broadening parameter. Default 0.2.

  • rstart (float, optional) – Integration start. Default 0.001.

  • h (float, optional) – Integration step (trapezoidal). Default 0.001.

  • local (bool, optional) – Use local density instead of global. Default False.

  • average (bool, optional) – Compute neighbor-averaged entropy. Default False.

Returns:

Per-atom entropy (or averaged entropy) values.

Return type:

numpy array

Short-Range Order#

pyscal3.short_range_order(atoms: Atoms, reference_type=1, compare_type=2, average=True)#

Calculate Warren-Cowley short-range order parameter.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors computed.

  • reference_type (int, optional) – Atomic number of reference type. Default 1.

  • compare_type (int, optional) – Atomic number to compare. Default 2.

  • average (bool, optional) – If True, return system average. Default True.

Returns:

Per-atom SRO values, or system average.

Return type:

numpy array or float

Radial Distribution Function#

pyscal3.radial_distribution_function(atoms: Atoms, rmin=0, rmax=5.0, bins=100)#

Calculate radial distribution function g(r).

Parameters:

  • atoms (ase.Atoms) – Structure.

  • rmin (float) – Distance range.

  • rmax (float) – Distance range.

  • bins (int) – Number of histogram bins.

Returns:

(rdf, r)

Return type:

tuple of numpy arrays

Angular Criteria#

pyscal3.angular_criteria(atoms: Atoms)#

Calculate angular criteria for diamond structure identification.

Parameters:

atoms (ase.Atoms) – Structure with neighbors computed.

Returns:

Per-atom angular parameter A values.

Return type:

numpy array

Chi Parameters#

pyscal3.chi_params(atoms: Atoms, angles=False)#

Calculate chi-parameter vector for structure identification.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors computed.

  • angles (bool, optional) – If True, also return cosine angles. Default False.

Returns:

Chi parameter vectors.

Return type:

numpy array of shape (natoms, 9)

Ackland-Jones Classification#

pyscal3.identify_ackland_jones(atoms: Atoms)#

Classify atomic environments with the Ackland–Jones method.

Uses the chi-parameter histogram (angular signature) to assign each atom a structure type via a decision tree. The chi histogram bins cosines of all pairwise neighbor angles into 9 ranges; the counts in specific bins discriminate FCC, HCP, BCC, and icosahedral (ICO) environments.

The algorithm follows Ackland & Jones, Phys. Rev. B 73, 054104 (2006), adapted for the 9-bin chi scheme used by pyscal3. :param atoms: pyscal3.find_neighbors()). :type atoms: ase.Atoms Structure with neighbors already computed (via

Returns:

  • labels (numpy.ndarray of int, shape (natoms,)) – Per-atom structure label:

    Value

    Structure

    0

    other / unknown

    1

    FCC

    2

    HCP

    3

    BCC

    4

    ICO (icosahedral)

  • names (list of str) – Human-readable name for each atom ("fcc", "hcp", etc.). Also stored in atoms.arrays["pyscal_structure"].

Notes

Chi-parameter bin edges (9 bins, cosine of the angle):

Bin

Cosine range

χ₀

[−1.000, −0.945) ~180°

χ₁

[−0.945, −0.915) ~160°

χ₂

[−0.915, −0.755) ~139°

χ₃

[−0.755, −0.705) ~135°

χ₄

[−0.705, −0.195)

χ₅

[−0.195, +0.195) ~90°

χ₆

[+0.195, +0.245)

χ₇

[+0.245, +0.795) ~60°

χ₈

[+0.795, +1.000] ~0°

Decision tree (ideal counts for perfect structures):

  • BCC (14 neighbours): χ₀ = 7, χ₅ = 12, χ₇ = 36

  • FCC (12 neighbours): χ₀ = 6, χ₅ = 12, χ₇ = 24, χ₁₊₂₊₃ = 0

  • HCP (12 neighbours): χ₀ = 3, χ₂ = 6, χ₇ = 24

  • ICO (12 neighbours): χ₀ = 6, χ₅ = 0, χ₇ = 30

References

G. J. Ackland and A. P. Jones, “Applications of local crystal structure measures in experiment and simulation”, Phys. Rev. B 73, 054104 (2006). doi:10.1103/PhysRevB.73.054104

Coordination Numbers#

pyscal3.coordination_number(atoms: Atoms)#

Return the simple coordination number (integer neighbor count).

Parameters:

atoms (ase.Atoms) – Structure with neighbors already computed.

Returns:

Per-atom coordination number. Also stored in atoms.arrays["pyscal_cn"].

Return type:

numpy.ndarray of int, shape (natoms,)

pyscal3.effective_coordination_number(atoms: Atoms)#

Calculate the effective coordination number (ECoN).

ECoN weights each neighbor by a continuous function of distance relative to the average neighbor distance, converging iteratively:

\[\mathrm{ECoN}_i = \sum_j \exp\!\left[ 1 - \left(\frac{d_{ij}}{d_{\mathrm{av},i}}\right)^6 \right]\]

where \(d_{\mathrm{av},i}\) is the weighted average distance that is itself computed from the weights (self-consistent iteration).

Parameters:

atoms (ase.Atoms) – Structure with neighbors already computed.

Returns:

Per-atom ECoN values. Also stored in atoms.arrays["pyscal_econ"].

Return type:

numpy.ndarray, shape (natoms,)

References

R. Hoppe, “Effective coordination numbers (ECoN) and mean fictive ionic radii (MEFIR)”, Z. Kristallogr. 150, 23 (1979).

pyscal3.generalized_coordination_number(atoms: Atoms, cn_max=None)#

Calculate the generalized coordination number (GCN).

GCN accounts for the coordination of each neighbor, penalizing under-coordinated surface atoms:

\[\mathrm{GCN}_i = \frac{1}{N_{\max}} \sum_{j \in \mathrm{neigh}(i)} \mathrm{CN}_j\]

where \(N_{\max}\) is the maximum (bulk) coordination number.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed.

  • cn_max (int or None) – Bulk coordination number. If None, uses the maximum CN observed in the system.

Returns:

Per-atom GCN values. Also stored in atoms.arrays["pyscal_gcn"].

Return type:

numpy.ndarray, shape (natoms,)

References

F. Calle-Vallejo et al., “Finding optimal surface sites on heterogeneous catalysts by counting nearest neighbors”, Science 350, 185 (2015). doi:10.1126/science.aab3501

pyscal3.local_density(atoms: Atoms)#

Estimate the local atomic number density.

For each atom, the local density is estimated from the mean neighbor distance:

\[\rho_i = \frac{N_i}{\frac{4}{3}\pi \bar{d}_i^3}\]

where \(N_i\) is the coordination number and \(\bar{d}_i\) is the mean neighbor distance.

Parameters:

atoms (ase.Atoms) – Structure with neighbors already computed.

Returns:

Per-atom local density (atoms per unit volume). Also stored in atoms.arrays["pyscal_local_density"].

Return type:

numpy.ndarray, shape (natoms,)

Angular and Bond Length Distributions#

pyscal3.angular_distribution_function(atoms: Atoms, bins=180)#

Calculate the angular distribution function (ADF).

The ADF is the histogram of all bond angles \(\theta_{jik}\) formed by pairs of neighbors (j, k) around each atom i. It characterises the local angular environment independently of a specific order parameter.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed.

  • bins (int) – Number of histogram bins in the angle range [0, 180] degrees. Default 180 (1-degree resolution).

Returns:

(adf, angles) – adf is the normalised probability density of bond angles, and angles is the array of bin left-edges in degrees.

Return type:

tuple of numpy.ndarray

Notes

Angles are computed from the cosine of the angle between displacement vectors, using the same C++ infrastructure as chi_params(). Results are also stored in atoms.info["pyscal_adf"] and atoms.info["pyscal_adf_angles"].

pyscal3.bond_length_distribution(atoms: Atoms, bins=100, rmin=None, rmax=None)#

Calculate the bond-length distribution function (BLDF).

Unlike the full radial distribution function, the BLDF histograms only the bonds defined by the current neighbor list (no shell-volume normalisation).

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed.

  • bins (int) – Number of histogram bins. Default 100.

  • rmin (float or None) – Distance range. If None, determined from the data.

  • rmax (float or None) – Distance range. If None, determined from the data.

Returns:

(bldf, r) – bldf is the normalised probability density, and r is the array of bin left-edges. Results are also stored in atoms.info["pyscal_bldf"] and atoms.info["pyscal_bldf_r"].

Return type:

tuple of numpy.ndarray

Atomic Deformation#

pyscal3.atomic_strain(atoms: Atoms, reference: Atoms)#

Calculate the atomic (Green-Lagrange) strain tensor for each atom.

The local deformation gradient F is computed by minimizing the squared difference between reference and deformed neighbor vectors. The Lagrangian strain is then E = (F^T F - I) / 2. :param atoms: :type atoms: ase.Atoms Deformed configuration with neighbors computed. :param reference: Reference (undeformed) configuration with the same neighbor list. :type reference: ase.Atoms

Returns:

Green-Lagrange strain tensor per atom. Also stored as atoms.arrays["pyscal_strain"].

Return type:

numpy.ndarray of shape (natoms, 3, 3)

Notes

Both configurations must have neighbors computed with identical neighbor lists (same method and cutoff). Reference: Falk & Langer, PRE 57 (1998) 7192 (D^2_min); Shimizu, Ogata, Li, Mat. Trans. 48 (2007) 2923 (atomic strain).

pyscal3.von_mises_strain(atoms: Atoms, reference: Atoms)#

Compute the von Mises shear strain invariant from the atomic strain.

\[\eta^{\text{Mises}} = \sqrt{ \frac{1}{2} \left[ (E_{xx}-E_{yy})^2 + (E_{yy}-E_{zz})^2 + (E_{zz}-E_{xx})^2 \right] + E_{xy}^2 + E_{yz}^2 + E_{xz}^2 }\]
Parameters:

  • atoms (ase.Atoms) – Deformed configuration.

  • reference (ase.Atoms) – Reference configuration.

Returns:

Scalar von Mises strain per atom. Also stored as atoms.arrays["pyscal_von_mises"].

Return type:

numpy.ndarray of shape (natoms,)

pyscal3.d2min(atoms: Atoms, reference: Atoms)#

Compute the D^2_min non-affine displacement (Falk & Langer 1998).

D^2_min is the residual mean-squared displacement after subtracting the best-fit affine deformation.

Parameters:

  • atoms (ase.Atoms) – Deformed configuration with neighbors computed.

  • reference (ase.Atoms) – Reference configuration with neighbors computed.

Returns:

D^2_min per atom (Angstrom^2). Also stored as atoms.arrays["pyscal_d2min"].

Return type:

numpy.ndarray of shape (natoms,)

pyscal3.slip_vector(atoms: Atoms, reference: Atoms)#

Compute the slip vector for each atom.

The slip vector s_i = (1/N_s) sum_j (dX_ij - dx_ij) is the average difference between reference and current displacement vectors, without fitting an affine transformation.

Parameters:

  • atoms (ase.Atoms) – Deformed configuration.

  • reference (ase.Atoms) – Reference configuration.

Returns:

Slip vector per atom. Also stored as atoms.arrays["pyscal_slip_vector"].

Return type:

numpy.ndarray of shape (natoms, 3)

Notes

Ref: Zimmerman, Kelchner, Klein, Hamilton, Foiles, PRL 87 (2001) 165507.

Wigner-Seitz Defect Analysis#

pyscal3.wigner_seitz_analysis(atoms: Atoms, reference: Atoms, affine_mapping: str = 'none', per_type_occupancies: bool = False) → dict#

Wigner-Seitz cell analysis for vacancy/interstitial detection.

This assigns each atom in atoms to the nearest lattice site in reference and counts site occupancies. Sites with occupancy 0 are vacancies; sites with occupancy > 1 contain interstitials.

Parameters:

  • atoms (Atoms) – The configuration to analyze (containing defects).

  • reference (Atoms) – The perfect/reference configuration defining lattice sites.

  • affine_mapping (str, optional) – How to handle cell distortion: - “none” (default): Use positions as-is - “to_reference”: Rescale atoms positions to match reference cell

  • per_type_occupancies (bool, optional) – If True, track occupancy by atom type (for antisite detection).

Returns:

Keys: - “vacancy_count”: Number of sites with occupancy = 0 - “interstitial_count”: Total excess atoms (sum of max(occ-1, 0)) - “occupancy”: (N_ref,) array of site occupancies - “site_index”: (N_atoms,) array mapping each atom to its assigned site - “vacancy_indices”: indices of reference sites with vacancies - “interstitial_sites”: indices of sites with excess atoms

If per_type_occupancies=True, also includes: - “occupancy_by_type”: dict mapping type → (N_ref,) occupancy array

Return type:

dict

Notes

Results are also stored on atoms: - atoms.arrays[“pyscal_ws_site_index”]: site assignment - atoms.arrays[“pyscal_ws_occupancy”]: occupancy of assigned site - atoms.info[“pyscal_ws_vacancy_count”] - atoms.info[“pyscal_ws_interstitial_count”]

The algorithm uses nearest-neighbor search via cKDTree. For periodic systems, reference sites near cell boundaries are replicated to handle atoms that may have wrapped to different periodic images.

Examples

>>> from ase.build import bulk
>>> import pyscal3
>>> # Perfect FCC reference
>>> ref = bulk("Cu", "fcc", cubic=True).repeat(3)
>>> # Create vacancy by deleting atom
>>> defected = ref.copy()
>>> del defected[0]
>>> result = pyscal3.wigner_seitz_analysis(defected, ref)
>>> print(f"Vacancies: {result['vacancy_count']}")  # 1
>>> print(f"Interstitials: {result['interstitial_count']}")  # 0
pyscal3.identify_defect_atoms(atoms: Atoms, reference: Atoms, affine_mapping: str = 'none') → dict#

Identify which atoms are at vacancies, interstitials, or antisites.

This is a convenience wrapper around wigner_seitz_analysis that returns masks for different defect types.

Parameters:

  • atoms (Atoms) – The configuration to analyze.

  • reference (Atoms) – The reference configuration.

  • affine_mapping (str, optional) – Affine mapping mode (see wigner_seitz_analysis).

Returns:

  • “perfect_mask”: bool array, atoms at singly-occupied sites

  • ”interstitial_mask”: bool array, atoms at multiply-occupied sites

  • ”vacancy_positions”: (N_vac, 3) positions of empty reference sites

  • ”defect_summary”: string description of defects found

Return type:

dict

ACE Descriptors#

pyscal3.ace(atoms: Atoms, nmax=4, lmax=4, nu_max=2, cutoff=None, normalize=True)#

Compute Atomic Cluster Expansion (ACE) descriptors.

ACE provides a systematic and complete expansion of atomic environments, with SOAP (nu=2) and bispectrum (nu=3) as special cases. The descriptors are rotationally, translationally, and permutationally invariant.

The implementation follows Drautz (2019) and computes B-basis descriptors by coupling A-functions (single-particle basis) to form rotationally invariant combinations at each correlation order.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors already computed.

  • nmax (int, default 4) – Number of radial basis functions. Higher values capture finer radial resolution but increase computation.

  • lmax (int, default 4) – Maximum angular momentum quantum number. Higher values capture more angular detail. Typically 3-6 for ML potentials.

  • nu_max (int, default 2) – Maximum correlation order: - nu=1: Radial density (neighbor count per shell) - nu=2: Pair correlations (SOAP power spectrum) - nu=3: Triplet correlations (bispectrum) Higher orders rapidly increase descriptor count.

  • cutoff (float, optional) – Neighbor cutoff radius. If None, uses the cutoff from find_neighbors.

  • normalize (bool, default True) – If True, normalize descriptors to unit norm per atom.

Returns:

Dictionary with keys: - ‘nu1’: ndarray (natoms, nmax) - radial density descriptors - ‘nu2’: ndarray (natoms, n2) - power spectrum (if nu_max >= 2) - ‘nu3’: ndarray (natoms, n3) - bispectrum-like (if nu_max >= 3) - ‘full’: ndarray (natoms, n_total) - concatenated descriptors

Return type:

dict

Notes

The descriptor count scales as: - nu=1: O(nmax) - nu=2: O(nmax^2 * lmax) - nu=3: O(nmax^3 * lmax^3) but limited for tractability

References

Examples

>>> from ase.build import bulk
>>> import pyscal3
>>> atoms = bulk("Cu", "fcc", cubic=True).repeat(3)
>>> pyscal3.find_neighbors(atoms, method="cutoff", cutoff=4.0)
>>> desc = pyscal3.ace(atoms, nmax=4, lmax=3, nu_max=2)
>>> print(desc['full'].shape)
>>> print("nu=2 descriptors:", desc['nu2'].shape[1])

Solid/Liquid Classification#

pyscal3.find_solids(atoms: Atoms, bonds=0.5, threshold=0.5, avgthreshold=0.6, cluster=True, q=6, cutoff=0, right=True)#

Distinguish solid and liquid atoms.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors computed.

  • bonds (int or float) – Min solid bonds (int) or min fraction of solid neighbors (float 0-1).

  • threshold (float) – Bond correlation cutoff. Default 0.5.

  • avgthreshold (float) – Average bond cutoff. Default 0.6.

  • cluster (bool) – If True, cluster solid atoms and return largest cluster size.

  • q (int) – Steinhardt parameter order. Default 6.

  • cutoff (float) – Cluster cutoff (0 = use neighbor cutoff).

  • right (bool) – If True, use > comparison. Default True.

Returns:

Largest cluster size if cluster=True.

Return type:

int or None

pyscal3.find_clusters(atoms: Atoms, condition, largest=True, cutoff=0, d=None)#

Cluster atoms based on a boolean condition.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors computed.

  • condition (array-like of bool) – Which atoms to include in clustering.

  • largest (bool) – If True, return largest cluster size.

  • cutoff (float) – Cluster cutoff (0 = use neighbor cutoff).

Returns:

Largest cluster size if largest=True.

Return type:

int

Utilities#

pyscal3.average_over_neighbors(atoms: Atoms, key: str, include_self=True)#

Average a per-atom property over each atom’s neighbors.

Parameters:

  • atoms (ase.Atoms) – Structure with neighbors computed.

  • key (str) – Key in atoms.arrays (with or without ‘pyscal_’ prefix).

  • include_self (bool) – Include the atom itself in the average. Default True.

Return type:

numpy array

Structure Creation#

pyscal3.structures.make_crystal(structure, lattice_constant=1.0, repetitions=None, ca_ratio=1.633, noise=0, element=None, primitive=False)#

Create a crystal structure as an ASE Atoms object.

Parameters:

  • structure (str) – Crystal type: ‘sc’, ‘bcc’, ‘fcc’, ‘hcp’, ‘dhcp’, ‘diamond’, ‘a15’, ‘l12’, ‘b2’.

  • lattice_constant (float) – Lattice constant in Angstroms. Default 1.0.

  • repetitions (int or tuple of 3 ints, optional) – Unit cell repetitions. Default (1,1,1).

  • ca_ratio (float) – c/a ratio for hcp/dhcp. Default 1.633.

  • noise (float) – Standard deviation of Gaussian noise on positions. Default 0.

  • element (str or list of str, optional) – Chemical element(s).

  • primitive (bool) – If True, use primitive cell. Default False.

Return type:

ase.Atoms

pyscal3.structures.make_element(symbol, repetitions=None, primitive=False)#

Create a crystal of a specific element using its known structure.

Parameters:

  • symbol (str) – Element symbol, e.g. ‘Fe’, ‘Cu’, ‘Al’.

  • repetitions (tuple of 3 ints, optional) – Default (1,1,1).

  • primitive (bool) – Use primitive cell. Default False.

Return type:

ase.Atoms

pyscal3.structures.make_general_lattice(positions, types, box, lattice_constant=1.0, repetitions=None, noise=0, element=None)#

Create a custom lattice from fractional positions.

Parameters:

  • positions (list of [x,y,z]) – Fractional coordinates.

  • types (list of int) – Per-atom type indices.

  • box (list of 3 vectors) – Unit cell vectors (before lattice_constant scaling).

  • lattice_constant (float) – Scaling factor. Default 1.0.

  • repetitions (tuple of 3 ints) – Repetitions. Default (1,1,1).

  • noise (float) – Gaussian noise std dev. Default 0.

  • element (str or list of str, optional) – Chemical element(s).

Return type:

ase.Atoms

pyscal3.structures.make_grain_boundary(axis, sigma, gb_plane, structure=None, element=None, lattice_constant=1, repetitions=(1, 1, 1), overlap=0.0)#

Create a symmetric tilt grain boundary as an ASE Atoms object.

Parameters:

  • axis (list of 3 int) – Tilt axis, e.g. [1, 0, 0].

  • sigma (int) – CSL sigma value.

  • gb_plane (list of 3 int) – Grain boundary plane Miller indices.

  • structure (str, optional) – Crystal structure (‘bcc’, ‘fcc’, etc.). Required if element not given.

  • element (str, optional) – Element symbol. Structure and lattice constant auto-determined.

  • lattice_constant (float) – Lattice constant. Default 1.

  • repetitions (tuple of 3 int) – Supercell repetitions. Default (1,1,1).

  • overlap (float) – Atom overlap distance to remove. Default 0.0.

Return type:

ase.Atoms

pyscal3.structures.available_structures()#

Return list of available crystal structure names.

pyscal3.structures.available_elements()#

Return list of available element symbols.

Trajectory#

class pyscal3.Trajectory(filename)#

Lazy reader for LAMMPS dump trajectory files.

Supports block-level random access via indexing and slicing. Blocks are loaded on demand and can be converted to ASE Atoms.

Parameters:

filename (str) – Path to a LAMMPS dump file.

nblocks#

Number of timestep blocks in the file.

Type:

int

natoms#

Number of atoms per block.

Type:

int

Examples

>>> traj = pyscal3.Trajectory("dump.lammpstrj")
>>> ts = traj[0]                    # single block → Timeslice
>>> atoms_list = ts.to_atoms(species=["Cu"])
>>> atoms_list = traj[0:10].to_atoms(species=["Cu"])  # slice
get_block(blockno)#

Get a block from the file as raw data

Parameters:

blockno (int) – number of the block to be read, starts from 0

Returns:

data – list of strings containing data

Return type:

list

load(blockno)#

Load the data of a block into memory as a dictionary of numpy arrays

Parameters:

blockno (int) – number of the block to be read, starts from 0

Return type:

None

Notes

When the data of a block is loaded, it is accessible through Trajectory.data[x]. This data can then be modified. When the block is written out, the modified data is written instead of existing one. But, loaded data is kept in memory until unloaded using unload method.

unload(blockno)#

Unload the data that is loaded to memory using load method

Parameters:

blockno (int) – number of the block to be read, starts from 0

Return type:

None