"""
StructureAnalyser — main class that delegates to submodules.
Usage
-----
from amorphgen.analysis import StructureAnalyser
sa = StructureAnalyser("output_dir/") # default cutoff="auto-rdf"
sa.summary()
sa.plot(output_dir="plots/")
"""
from __future__ import annotations
import os
import glob
import numpy as np
from ase.io import read
from .cutoff import auto_cutoff_minsep, auto_cutoff_rdf
from .structure import (compute_density, compute_coordination,
compute_bond_distances, compute_all_angles,
compute_bond_angle_stats, build_neighbour_dict)
from .rdf import (compute_rdf, compute_structure_factor,
compute_averaged_rdf, DEFAULT_SMEARING, DEFAULT_SQ_SMOOTH)
from .rings import compute_ring_statistics
from .voronoi import compute_voronoi
from .energy import compute_energy_ranking
from .plotting import plot_analysis
from .uncertainty import summarize_structures, summarize_site_groups
[docs]
class StructureAnalyser:
"""
Analyse amorphous structures: density, coordination, distances,
angles, RDF, S(q), rings, Voronoi, voids, oxygen speciation,
crystal-like bond order, elastic moduli, vibrational DOS and energy ranking.
Results retain per-structure observations and ``uncertainty`` summaries
of equal-weight structure means (SEM, Student-t intervals and percentile
bootstrap bounds). Pooled site/bond/angle ``std`` fields describe spread,
not uncertainty. Curve bands are pointwise. Independent structures are
assumed; fewer than two contributing structures leave intervals undefined.
Parameters
----------
source : str or list of str
Path to a structure file, a directory, or a list of file paths.
cutoff : float, dict, or str
- float: single cutoff (A) for all pairs
- dict: pair-specific, e.g. {"Si-O": 2.2}
- "auto-rdf": from RDF first minimum (default — physically correct
first coordination shell, matches the standard convention in
neutron-diffraction analysis of glasses/liquids)
- "auto": from bonding radii (fast; can truncate the first peak
for systems with broad bond distributions such as a-Si, a-HfO2,
chalcogenides — prefer "auto-rdf" unless you have a specific
reason)
"""
_file_list: list # populated by _load; used by log-energy lookup
[docs]
def __init__(self, source, cutoff="auto-rdf"):
"""
Parameters
----------
source : str, list of str, or list of Atoms
Path to a structure file, a directory of structure files,
a list of file paths, or a list of ASE Atoms objects.
cutoff : float, dict, or str
Neighbour cutoff for CN, bond distances, and angles.
- float: single cutoff (A) for all pairs
- str with overrides: ``"In-O=2.6,Zn-O=2.3"`` keeps auto-rdf for
every other pair; ``"auto,In-O=2.6"`` or ``"2.4,In-O=2.6"``
set the base rule explicitly
- dict listing only some pairs: the rest come from auto-rdf
(or from its ``"default"`` entry: a mode name or a number)
- dict: pair-specific, e.g. ``{"Si-O": 2.2, "O-O": 2.8}``
- ``"auto-rdf"`` (default): first minimum of the partial RDF —
the standard physical definition of the first coordination
shell. Robust across material classes; prefer this for
published numbers.
- ``"auto"``: derived from Shannon/Cordero/Goldschmidt radii
(minsep). Fast but can under-count coordination for
systems with broad first-shell distributions (a-Si,
a-HfO2, chalcogenides). Kept for backward compatibility
and for the placement-stage workflows where no RDF is
available yet.
Notes
-----
Changed in v1.0.0: default switched from ``"auto"`` to
``"auto-rdf"`` to match neutron-diffraction analysis convention
and to avoid systematically under-counting coordination for
materials with broad bond-length distributions.
"""
self.atoms_list, self._file_list = self._load(source)
if not self.atoms_list:
raise ValueError("No structures loaded. Check the source path.")
# Warn if structures have different compositions
formulas = set(a.get_chemical_formula(mode="hill")
for a in self.atoms_list)
if len(formulas) > 1:
import warnings
warnings.warn(
f"Structures have different compositions: {formulas}. "
f"Analysis assumes uniform composition; energy/atom and "
f"CN statistics may mix incompatible systems."
)
# Any accepted form: "auto-rdf", "auto", a number, "In-O=2.6,Zn-O=2.3",
# "auto,In-O=2.6", or a dict (optionally with a "default" entry).
# A dict that lists only some pairs is completed from auto-rdf.
from .cutoff import resolve_cutoffs
# Auto cutoffs pool frame RDFs and use the first frame for species
# discovery. Canonicalize that calculation so generation order cannot
# alter extracted descriptors or their convergence curves. Keep the
# caller's ordering for all per-structure results and file alignment.
cutoff_atoms = sorted(self.atoms_list, key=lambda atoms: (
atoms.numbers.tobytes(), np.asarray(atoms.cell).tobytes(),
atoms.positions.tobytes(), atoms.pbc.tobytes()))
self.cutoff, self._cutoff_mode = resolve_cutoffs(cutoff_atoms, cutoff)
self._max_cutoff = (
max(self.cutoff.values()) if isinstance(self.cutoff, dict)
else self.cutoff
)
@staticmethod
def _load(source):
"""Return (atoms_list, file_paths). file_paths may be [] if the
caller passed in-memory Atoms objects rather than disk paths."""
from ase import Atoms
from ..utils.relaxation import read_relaxation_metadata
def load_file(path):
atoms = read(path)
read_relaxation_metadata(path, atoms)
return atoms
if isinstance(source, list):
if source and isinstance(source[0], Atoms):
return source, []
file_list = list(source)
return [load_file(f) for f in file_list], file_list
if os.path.isdir(source):
# One structure per stem: the optimiser writes s_opt.xyz AND
# s_opt.cif (and .traj), which must not count twice. Priority
# xyz > extxyz > vasp > cif (the formats that carry energy first).
by_stem = {}
for ext in ("*.xyz", "*.extxyz", "*.vasp", "*.cif"):
for f in glob.glob(os.path.join(source, ext)):
by_stem.setdefault(os.path.splitext(os.path.basename(f))[0], f)
files = [by_stem[k] for k in sorted(by_stem)]
if not files:
raise FileNotFoundError(f"No structure files in {source}/")
return [load_file(f) for f in files], files
return [load_file(source)], [source]
[docs]
def screen(self, config=True):
"""Label candidates and record exclusion decisions without changing them.
Each screen has an independent ``exclude`` policy (default false).
Coordination uses this analyser's resolved chemical cutoffs; default
allowed sets come from random generation's inferred targets and
tolerance. See :func:`amorphgen.analysis.screen_structures` for settings.
The returned audit starts with zero analysed structures.
"""
from .screening import screen_structures
return screen_structures(self.atoms_list, config, cutoff=self.cutoff,
source_names=self._file_list or None)
[docs]
def screened(self, config=True):
"""Return ``(retained_analyser_or_None, audit)`` without mutating inputs.
Selection keeps the full candidate ensemble's resolved cutoffs fixed
and preserves file alignment. Labels alone do not remove a structure.
Call ``mark_screening_analysed(audit, audit['retained_indices'])`` after
successfully analysing the subset to record its contribution.
"""
from copy import copy
report = self.screen(config)
indices = report["retained_indices"]
if not indices:
return None, report
selected = copy(self)
selected.atoms_list = [self.atoms_list[i] for i in indices]
selected._file_list = ([self._file_list[i] for i in indices]
if self._file_list else [])
return selected, report
def _get_cutoff(self, s1, s2):
if isinstance(self.cutoff, (int, float)):
return self.cutoff
key1 = f"{s1}-{s2}"
key2 = f"{s2}-{s1}"
if key1 in self.cutoff:
return self.cutoff[key1]
if key2 in self.cutoff:
return self.cutoff[key2]
if not hasattr(self, "_warned_pairs"):
self._warned_pairs = set()
if key1 not in self._warned_pairs:
self._warned_pairs.add(key1)
import warnings
warnings.warn(f"no cutoff for pair {key1}; treating it as not bonded "
f"(add it to the cutoff dict to include it)")
return 0.0
def _build_neighbour_dict(self, atoms):
return build_neighbour_dict(atoms, self._max_cutoff, self._get_cutoff)
[docs]
def total_coordination(self, centre=None, partners=None):
"""First-shell coordination of an element counting several partner
types together.
Parameters
----------
centre : str, optional
Element at the centre (``"O"``). Default: every element.
partners : iterable of str, optional
Partner elements to count, such as In and Ga. Default: every
*bonded* partner type (bonded = cation-anion or hetero covalent,
the rule of the "Bonding coordination numbers" table). When given,
the named partners are counted within their pair cutoffs whether
or not the pair is classed as a bond.
Returns
-------
dict
Keyed by element, with ``mean``, ``std``, ``min``, ``max`` and
``distribution`` (CN to percent). For IGZO the default gives the
total O coordination (Ga + In + Zn around O) that the per-pair
O-Ga / O-In / O-Zn entries only give in parts; ``centre='O'`` with
``partners`` set to In and Ga gives the count over the two larger
cations only.
"""
from .structure import is_bonding_pair
from collections import Counter
wanted = set(partners) if partners is not None else None
frames = []
for atoms in self.atoms_list:
frame = {}
nbr, syms = self._build_neighbour_dict(atoms)
elements = Counter(syms)
for i, si in enumerate(syms):
if centre is not None and si != centre:
continue
cn = sum(1 for _, sj, _, _ in nbr[i]
if (sj in wanted if wanted is not None
else is_bonding_pair(si, sj, elements)))
frame.setdefault(si, []).append(cn)
frames.append(frame)
out = {}
for s in sorted({s for frame in frames for s in frame}):
out[s] = summarize_site_groups([frame.get(s, []) for frame in frames])
return out
# ── Core analysis methods ────────────────────────────────────────────
[docs]
def density(self):
"""Compute density (g/cm3) for each structure.
Returns
-------
dict
{"values": list[float], "mean": float, "std": float}
"""
return compute_density(self.atoms_list)
[docs]
def coordination(self, pair=None):
"""Compute coordination numbers across all structures.
Parameters
----------
pair : str, optional
Pair to analyse, e.g. "Si-O". If None, all pairs computed.
Returns
-------
dict
Keyed by pair string. Each value is a dict with
"mean", "std", "min", "max", "distribution", "total_atoms".
"""
return compute_coordination(self.atoms_list, self._max_cutoff,
self._get_cutoff, pair)
[docs]
def cutoff_robustness(self, window=0.1, points=5):
"""Report near-cutoff pair shares and directional coordination curves.
Sweep the resolved cutoffs over +/- ``window`` Angstrom (default
0.10) at ``points`` evenly spaced values (odd and >= 3, default 5).
Auto cutoffs are held fixed, with no RDF refit. Positive cutoffs are
shifted together, lower radii clipped to zero, and zero cutoffs kept
at zero. The analyser and its structures are not modified.
``pairs`` maps unordered element pairs to ``cutoff``, actual
``cutoffs``, ``near_pairs``, ``n_pairs`` and ``near_fraction``.
The fraction divides contacts added across the window by contacts
included at its upper endpoint; it is None for an empty denominator.
Reciprocal contacts count once, distinct periodic images separately.
Each pair's ``coordination`` has one curve per direction: ``mean``
pools central sites, ``ensemble_mean`` weights contributing structures
equally, ``delta`` / ``ensemble_delta`` give upper-minus-lower changes,
and ``per_structure`` retains aligned curves (None for absent centres).
Pair counts also retain aligned ``per_structure`` observations.
Zero-neighbour sites are included. The central point matches
:meth:`coordination`, including its exact distance-boundary rules.
These perturbations measure sensitivity, not sampling uncertainty.
"""
from .robustness import compute_cutoff_robustness
report = compute_cutoff_robustness(
self.atoms_list, self._max_cutoff, self._get_cutoff, window, points)
report["cutoff_mode"] = self._cutoff_mode
return report
[docs]
def dimer_report(self, threshold_frac=0.85):
"""Detect unphysical same-element / anion-anion close contacts.
Flags "dimers" — pairs closer than ``threshold_frac`` times the
radii-derived minimum separation (O-O peroxide below ~1.9 A, Cl2,
metal-metal dimers). These form when cold MLIP relaxation collapses
under-coordinated seeds and mark structures that usually rank higher
in energy; the recommended workflow is to generate an ensemble and
keep the dimer-free, lowest-energy members.
Parameters
----------
threshold_frac : float
Fraction of the pair minsep below which a contact counts as a
dimer (default 0.85).
Returns
-------
dict
``{"pairs": {pair: {"count", "min_distance", "threshold"}},
"per_structure": list[int], "total": int, "n_structures": int,
"threshold_frac": float}``
"""
from .structure import compute_dimers
return compute_dimers(self.atoms_list, threshold_frac=threshold_frac)
[docs]
def bond_order(self, qbar6_threshold=0.3, min_neighbors=4, cutoff=None):
"""Steinhardt q6, Lechner--Dellago qbar6 and ordered clusters.
Uses this analyser's resolved neighbour cutoffs unless ``cutoff``
is provided. Ordered atoms have qbar6 at least ``qbar6_threshold``
and at least ``min_neighbors`` neighbour images. Calibrate the
threshold against the material's crystal and liquid structures.
See :func:`amorphgen.analysis.compute_bond_order` for output fields.
"""
from .bond_order import compute_bond_order
return compute_bond_order(
self.atoms_list, self.cutoff if cutoff is None else cutoff,
qbar6_threshold=qbar6_threshold, min_neighbors=min_neighbors,
)
[docs]
def bond_distances(self, pair=None):
"""Compute bond distance statistics.
Parameters
----------
pair : str, optional
Pair to analyse, e.g. "Si-O". If None, all pairs computed.
Returns
-------
dict
Keyed by pair string. Each value is a dict with
"mean", "std", "min", "max", "count".
"""
return compute_bond_distances(self.atoms_list, self._max_cutoff,
self._get_cutoff, pair)
def _compute_all_angles(self, triplet=None, bonding_only=True):
return compute_all_angles(self.atoms_list, self._max_cutoff,
self._get_cutoff, triplet, bonding_only)
[docs]
def bond_angles(self, triplet=None):
"""Compute bond angle statistics.
Parameters
----------
triplet : str, optional
Triplet to analyse, e.g. "O-Si-O". If None, all triplets computed.
Returns
-------
dict
Keyed by triplet string. Each value is a dict with
"mean", "std", "min", "max", "count" (angles in degrees).
"""
raw = self._compute_all_angles(triplet)
return compute_bond_angle_stats(raw)
[docs]
def angle_distribution(self, triplet=None, bins=90, *, normalise=True,
confidence=0.95, n_bootstrap=1000, seed=0):
"""Per-structure angle histograms and uncertainty of their mean.
Integer ``bins`` spans 0–180 degrees, including linear angles.
Explicit edges must cover that range. Each structure's histogram is
normalized independently before averaging; a structure without the
triplet is missing, not a zero-density sample. Bands are pointwise
percentile intervals from resampling complete structures.
"""
if np.isscalar(bins):
if isinstance(bins, bool) or int(bins) != bins or bins < 1:
raise ValueError("bins must be a positive integer or bin edges")
edges = np.linspace(0.0, 180.0, int(bins) + 1)
else:
edges = np.asarray(bins, dtype=float)
if (edges.ndim != 1 or len(edges) < 2
or not np.isfinite(edges).all() or np.any(np.diff(edges) <= 0)
or edges[0] > 0 or edges[-1] < 180):
raise ValueError("angle bin edges must increase and cover 0–180 degrees")
raw = self._compute_all_angles(triplet)
out = {}
for key in sorted(raw):
values = []
for frame in raw.per_structure:
angles = frame.get(key, [])
if not angles:
values.append([None] * (len(edges) - 1))
continue
hist = np.histogram(angles, bins=edges)[0].astype(float)
if normalise:
hist /= len(angles) * np.diff(edges)
values.append(hist.tolist())
uncertainty = summarize_structures(
values, confidence=confidence, n_bootstrap=n_bootstrap, seed=seed)
out[key] = {"angle": ((edges[:-1] + edges[1:]) / 2).tolist(),
"bin_edges": edges.tolist(), "distribution": uncertainty["mean"],
"normalization": "probability_density" if normalise else "count",
"per_structure": uncertainty["per_structure"],
"uncertainty": uncertainty}
return out
# ── RDF and S(q) ────────────────────────────────────────────────────
[docs]
def rdf(self, pair=None, rmax=None, nbins=200, sigma=DEFAULT_SMEARING, *,
confidence=0.95, n_bootstrap=1000, seed=0):
"""Compute the radial distribution function g(r).
Parameters
----------
pair : str, optional
Pair to analyse, e.g. "Si-O". If None, total RDF.
rmax : float, optional
Maximum radius in A. Auto-detected from cell if None.
nbins : int
Number of histogram bins (default 200).
sigma : float
Gaussian smearing width in A (default ``DEFAULT_SMEARING`` =
0.05). Pass 0.0 for the raw histogram. S(q) methods are
unaffected (they always use the raw g(r)).
Returns
-------
dict
{"r": list[float], "g_r": list[float]}
"""
return compute_rdf(self.atoms_list, pair, rmax, nbins, sigma=sigma,
confidence=confidence, n_bootstrap=n_bootstrap, seed=seed)
[docs]
def structure_factor(self, pair=None, qmax=15.0, nq=300, rmax=None,
weighting="unweighted", *, confidence=0.95,
n_bootstrap=1000, seed=0):
"""Compute the structure factor S(q) from g(r) via Fourier transform.
Parameters
----------
pair : str, optional
Pair to analyse, e.g. ``"Ga-O"``. If ``None``, computes
the total S(q). The ``weighting`` argument only affects the
total case; specific partials are always returned as their
own Faber-Ziman S_AB(q).
qmax : float
Maximum q in inverse-Angstrom (default 15.0).
nq : int
Number of q points (default 300).
rmax : float, optional
Max radius for the underlying g(r). ``None`` = auto (half cell).
weighting : {"unweighted", "xray", "neutron"}, default ``"unweighted"``
How partials are combined into the total. ``"xray"`` uses
q-dependent Waasmaier-Kirfel form factors (Faber-Ziman);
``"neutron"`` uses tabulated coherent scattering lengths.
Match the weighting and normalization of the reference data.
Returns
-------
dict
``{"q": list[float], "s_q": list[float]}``.
"""
return compute_structure_factor(self.atoms_list, pair, qmax, nq,
rmax, weighting=weighting, confidence=confidence,
n_bootstrap=n_bootstrap, seed=seed)
[docs]
def structure_factor_direct(self, qmax=15.0, nq=300,
weighting="xray", sigma_q=0.0,
partials=False, *, confidence=0.95,
n_bootstrap=1000, seed=0):
"""Compute Faber-Ziman S(q) from reciprocal-lattice scattering sums.
Avoids the finite-rmax integral in :meth:`structure_factor`.
The cell size still limits q sampling (q_min ~ 2*pi/L for a cubic
cell); inspect ``n_per_bin`` and check convergence.
Parameters
----------
qmax : float
Maximum q in inverse-Angstrom.
nq : int
Number of q-bins for spherical averaging.
weighting : {"xray", "neutron", "unweighted"}, default ``"xray"``
Per-element scattering factors (q-dependent Waasmaier-Kirfel
form factors for x-rays, Sears scattering lengths for
neutrons). See :func:`compute_structure_factor` for details.
sigma_q : float, default 0.0
Gaussian re-binning width in 1/A, weighted by the number of
q-vectors per shell. Reduces noise but can broaden features;
compare with the raw values. 0 returns raw shell averages.
partials : bool, default False
Also return the Faber-Ziman partial structure factors
``S_ab(q)`` computed from the per-species amplitudes.
Returns
-------
dict
``{"q": list[float], "s_q": list[float], "n_per_bin": list[int]}``,
plus ``"s_q_raw"`` when ``sigma_q > 0`` and ``"partials"``
(``{"A-B": list[float], ...}``) when ``partials=True``.
"""
from .rdf import compute_structure_factor_direct
return compute_structure_factor_direct(self.atoms_list, qmax, nq,
weighting=weighting,
sigma_q=sigma_q,
partials=partials, confidence=confidence,
n_bootstrap=n_bootstrap, seed=seed)
[docs]
def total_correlation(self, weighting="xray", qmin=0.3, qmax=20.0, nq=400,
rmax=10.0, nr=600, window="lorch", *, confidence=0.95,
n_bootstrap=1000, seed=0, sigma_q=DEFAULT_SQ_SMOOTH):
"""Total correlation function T(r) = 4 pi r rho g(r), the curve a
diffraction paper plots beside S(Q).
The g(r) behind it is SCATTERING-WEIGHTED and obtained by Fourier
transforming the weighted S(Q) over the measured Q range, so it is
directly comparable with published data and is NOT the same as
:meth:`rdf` with ``pair=None``, which weights every pair equally.
Set ``qmin``/``qmax``/``window`` to the experiment's own values.
``sigma_q`` is the direct S(q) re-binning width in inverse Angstrom
(default 0.05); use 0 for the unsmoothed shells.
Returns a dict with ``r``, ``g_r``, ``T_r``, ``G_r`` (the reduced PDF),
the ``q``/``s_q`` used, and ``rho``.
"""
from .rdf import compute_total_correlation
return compute_total_correlation(self.atoms_list, weighting=weighting,
qmin=qmin, qmax=qmax, nq=nq, rmax=rmax,
nr=nr, window=window, confidence=confidence,
n_bootstrap=n_bootstrap, seed=seed,
sigma_q=sigma_q)
[docs]
def xrd_pattern(self, wavelength=1.5406, qmax=None, nq=300,
method="direct", rmax=None, sigma_q=0.0, *,
q_batch=4096, confidence=0.95, n_bootstrap=1000, seed=0):
"""Coherent X-ray intensity per atom versus 2theta in degrees.
``wavelength`` is in Angstrom (default Cu K-alpha). Each structure's
X-ray S(q) is converted using its own composition before averaging.
Returns ``q``, ``two_theta``, ``intensity``, ``per_structure`` and
pointwise ``uncertainty``. Missing reciprocal shells stay missing.
This profile includes no instrument, background or sample corrections.
See :func:`amorphgen.analysis.compute_xrd_pattern` for conventions.
"""
from .xrd import compute_xrd_pattern
return compute_xrd_pattern(
self.atoms_list, wavelength=wavelength, qmax=qmax, nq=nq,
method=method, rmax=rmax, sigma_q=sigma_q, q_batch=q_batch,
confidence=confidence, n_bootstrap=n_bootstrap, seed=seed)
[docs]
def compare_experiment(self, experiment, kind=None, *, method="direct",
weighting="xray", calculation_options=None,
load_options=None, x_range=None, confidence=0.95,
n_bootstrap=1000, seed=0):
"""Compare a measured S(q) or T(r) file with this ensemble.
``experiment`` is a path or a mapping returned by ``load_experiment``.
``kind`` is ``"sq"`` (the default for paths) or ``"tr"``. Set
``calculation_options`` to the chosen scattering method's options,
e.g. ``{"qmax": 12, "sigma_q": 0.05}``; for T(r), match the measured
q range and window. ``load_options`` accepts delimiter, columns and
skiprows. Coordinates use inverse Angstrom for S(q), Angstrom for T(r).
Each structure is interpolated onto the measured grid, with no
extrapolation or interpolation across missing bins. The result has
residuals, goodness-of-fit metrics and pointwise mean confidence bands.
No scale, offset or other parameter is fitted. Scattering conventions
and instrument resolution must already match the measured curve.
"""
import inspect
from collections.abc import Mapping
from .experiment import load_experiment, compare_experiment, _kind
if isinstance(experiment, (str, os.PathLike)):
experiment = load_experiment(experiment, kind=kind or "sq",
**(load_options or {}))
elif load_options:
raise ValueError("load_options requires an experimental data path")
if not isinstance(experiment, Mapping):
raise ValueError("experiment must be a path or a mapping from load_experiment")
resolved_kind = _kind(kind or experiment.get("kind", "sq"))
if method not in ("direct", "ft"):
raise ValueError("method must be 'direct' or 'ft'")
if resolved_kind == "tr" and method != "direct":
raise ValueError("T(r) uses the direct S(q) method")
options = dict(calculation_options or {})
for key in ("weighting", "confidence", "n_bootstrap", "seed"):
if key in options:
raise ValueError(f"pass {key} directly, not in calculation_options")
fn = (self.total_correlation if resolved_kind == "tr" else
self.structure_factor_direct if method == "direct" else
self.structure_factor)
calculated = fn(weighting=weighting, confidence=confidence,
n_bootstrap=0, seed=seed, **options)
result = compare_experiment(
calculated, experiment, kind=resolved_kind, x_range=x_range,
confidence=confidence, n_bootstrap=n_bootstrap, seed=seed)
bound = inspect.signature(fn).bind_partial(weighting=weighting, **options)
bound.apply_defaults()
settings = {key: value for key, value in bound.arguments.items()
if key not in ("confidence", "n_bootstrap", "seed")}
if resolved_kind == "sq" and method == "ft":
from .rdf import _shared_rmax
settings["rmax"] = _shared_rmax(self.atoms_list, settings["rmax"])
result["calculation"] = {
"method": method, **settings,
"normalization": ("Faber-Ziman" if resolved_kind == "sq" else
"T(r) = 4*pi*r*rho*g(r)"),
}
result["structure_files"] = list(self._file_list)
return result
[docs]
def averaged_rdf(self, pair=None, rmax=None, nbins=200, *, confidence=0.95,
n_bootstrap=1000, seed=0):
"""Compute RDF per structure with mean and standard deviation.
Parameters
----------
pair : str, optional
Pair to analyse, e.g. "Si-O". If None, total RDF.
rmax : float, optional
Maximum radius in A. Auto-detected from cell if None.
nbins : int
Number of histogram bins (default 200).
Returns
-------
dict
{"r": list, "g_r_mean": list, "g_r_std": list, "n_structures": int}
"""
return compute_averaged_rdf(self.atoms_list, pair, rmax, nbins,
confidence=confidence, n_bootstrap=n_bootstrap, seed=seed)
# ── Advanced analysis ───────────────────────────────────────────────
[docs]
def ring_statistics(self, bond_pair=None, cutoff=None, max_ring=12):
"""Compute ring size statistics for network-forming structures.
Parameters
----------
bond_pair : tuple of str, optional
Bond pair to trace, e.g. ("Si", "O"). Auto-detected if None.
cutoff : float or dict, optional
Positive bond cutoff in Angstrom, or an ASE pair-cutoff mapping.
Uses analyser pair cutoff if None.
max_ring : int
Maximum ring size to search (default 12).
Returns
-------
dict
Size distribution, mean/spread, resolved and unresolved network
edge counts, and per-structure statistics with uncertainty.
Counts (including legacy ``total_rings``) are shortest-cycle
observations per edge, not the number of unique rings. Unresolved
edges have no closure within ``max_ring``. See
:func:`amorphgen.analysis.rings.compute_ring_statistics`.
"""
return compute_ring_statistics(
self.atoms_list, bond_pair, cutoff, max_ring, self._get_cutoff)
[docs]
def polyhedral_connectivity(self, cation=None, anion=None):
"""Corner/edge/face sharing between cation-centred polyhedra.
Two cations sharing one anion neighbour are corner-sharing, two
anions edge-sharing, three or more face-sharing. Returns the
linkage distribution, the percentage of cations in at least one
edge- or face-sharing pair (near 0 in a-SiO2 / a-GeO2; tens of
percent in unrelaxed random placements), per-species means and
per-structure values. See
:func:`amorphgen.analysis.structure.compute_polyhedral_connectivity`.
"""
from .structure import compute_polyhedral_connectivity
return compute_polyhedral_connectivity(
self.atoms_list, self._max_cutoff, self._get_cutoff,
cation=cation, anion=anion)
[docs]
def voronoi(self, element=None):
"""Compute Voronoi tessellation indices (n3, n4, n5, n6).
Parameters
----------
element : str, optional
Element to analyse. If None, all atoms included.
Returns
-------
dict
{"indices": list, "distribution": dict, "top_10": list,
"mean_faces": float, "total_atoms": int}
"""
return compute_voronoi(self.atoms_list, element)
[docs]
def void_distribution(self, n_samples=10000, probe_radius=0.0,
radii=None, nbins=50, seed=0, *, probe_radii=None):
"""Sample periodic free-space clearance and accessible volume.
Distances/radii are in Angstrom. This is a volume-weighted point
clearance distribution; it does not identify connected pores.
``probe_radii`` selects thresholds for an accessible-volume curve
evaluated from the same samples (default: histogram bin edges).
Clearance quantiles describe points accessible to ``probe_radius``.
See :func:`amorphgen.analysis.voids.compute_void_distribution`.
"""
from .voids import compute_void_distribution
return compute_void_distribution(
self.atoms_list, n_samples=n_samples, probe_radius=probe_radius,
radii=radii, nbins=nbins, seed=seed, probe_radii=probe_radii)
[docs]
def oxygen_speciation(self, network_formers=None):
"""Count free/non-bridging/bridging and multiply shared oxygen.
Uses this analyser's pair cutoffs. Defaults to Al, B, Ge, P and Si
present; explicitly select the network formers for other oxides.
See :func:`amorphgen.analysis.oxygen.compute_oxygen_speciation`.
"""
from .oxygen import compute_oxygen_speciation
return compute_oxygen_speciation(
self.atoms_list, network_formers=network_formers, cutoff=self.cutoff)
[docs]
def elastic_moduli(self, calculator=None, strain=0.005, relax=False,
fmax=0.01, steps=200):
"""Finite-strain stress response and Voigt/Reuss/Hill moduli (GPa).
Requires an active ASE stress calculator, supplied explicitly or
attached to each structure. ``relax=True`` relaxes internal atomic
positions at fixed cell. See
:func:`amorphgen.analysis.elasticity.compute_elastic_moduli`.
"""
from .elasticity import compute_elastic_moduli
return compute_elastic_moduli(
self.atoms_list, calculator=calculator, strain=strain,
relax=relax, fmax=fmax, steps=steps)
[docs]
def vibrational_dos(self, calculator=None, displacement=0.01,
npoints=400, sigma=0.1):
"""Harmonic DOS from finite-difference forces, in THz.
Uses Gamma-point normal modes of each supplied cell, with negative
frequencies representing imaginary modes. Requires 6N force calls
per structure. See
:func:`amorphgen.analysis.vibrations.compute_vibrational_dos`.
"""
from .vibrations import compute_vibrational_dos
return compute_vibrational_dos(
self.atoms_list, calculator=calculator, displacement=displacement,
npoints=npoints, sigma=sigma)
[docs]
def energy_ranking(self):
"""Rank structures by potential energy per atom.
Reads energy from atoms.info or attached calculator.
Returns
-------
dict
{"energies_per_atom": dict, "ranking": list, "best": int,
"worst": int, "best_energy": float, "worst_energy": float,
"spread": float}. Returns "warning" key if no energy data.
"""
return compute_energy_ranking(self.atoms_list)
# ── Multi-structure averaging ───────────────────────────────────────
[docs]
def averaged_cn(self, pair=None):
"""Compute per-structure mean CN with overall mean and std.
Parameters
----------
pair : str, optional
Pair to analyse, e.g. "Si-O". If None, all pairs computed.
Returns
-------
dict
Keyed by pair string. Each value is a dict with
"mean_per_structure", "overall_mean", "overall_std", "n_structures".
"""
result = {}
for key, data in self.coordination(pair).items():
uncertainty = data["uncertainty"]
means = [x for x in data["per_structure"] if x is not None]
result[key] = {
"mean_per_structure": data["per_structure"],
"overall_mean": uncertainty["mean"],
"overall_std": float(np.std(means)),
"n_structures": uncertainty["n_structures"],
"uncertainty": uncertainty,
}
return result
# ── Summary and reporting ───────────────────────────────────────────
[docs]
def convergence_report(self, tolerances=None, *, descriptors=None,
confidence=0.95, sizes=None, max_structures=1000000):
"""Report ensemble precision against declared absolute tolerances.
Core descriptor names are ``density`` (g/cm3),
``coordination.Si-O`` and ``total_coordination.Si`` (neighbours),
``bond_distance.O-Si`` (angstrom), and ``bond_angle.O-Si-O`` (degrees),
with the species present in this ensemble replacing these examples.
Each observation is one structure's mean, regardless of its atom count.
``tolerances`` maps descriptor names to positive absolute confidence
half-widths. Undeclared descriptors are reported without a pass/fail
decision. ``descriptors`` can add named per-structure values or existing
uncertainty summaries, for example
``{"rdf.total": self.rdf()["uncertainty"]}``. Additional observations
must align with this ensemble and may not replace core descriptors.
The curve uses full-ensemble variance and Student-t multipliers;
reordering the same observations cannot change it. Projections assume
independent structures with unchanged variance and descriptor
availability. Curve-valued descriptors require every point's half-width
to meet the tolerance; intervals remain pointwise, not simultaneous.
See :func:`amorphgen.analysis.convergence_report` for report fields,
missing-data handling and the distinction from generation-prefix tests.
"""
from collections.abc import Mapping
from .convergence import convergence_report
values = {"density": self.density()["uncertainty"]}
for prefix, results in (
("coordination", self.coordination()),
("total_coordination", self.total_coordination()),
("bond_distance", self.bond_distances()),
("bond_angle", self.bond_angles())):
values.update({f"{prefix}.{key}": result["uncertainty"]
for key, result in results.items()})
core_names = set(values)
if descriptors is not None:
if not isinstance(descriptors, Mapping):
raise ValueError("descriptors must be a mapping of names to observations")
overlap = values.keys() & descriptors.keys()
if overlap:
raise ValueError("Additional descriptors cannot replace core descriptors: "
+ ", ".join(sorted(overlap)))
values.update(descriptors)
report = convergence_report(values, tolerances, confidence=confidence,
sizes=sizes, max_structures=max_structures)
units = {"density": "g/cm^3", "coordination": "neighbours",
"total_coordination": "neighbours", "bond_distance": "angstrom",
"bond_angle": "degrees"}
for name, result in report["descriptors"].items():
prefix = name.split(".", 1)[0]
if name in core_names and prefix in units:
result["units"] = units[prefix]
report["analysis_settings"] = {
"cutoff": self.cutoff, "cutoff_mode": self._cutoff_mode,
}
return report
[docs]
def summary(self, show_angles=True, cutoff_window=0.1):
"""Print and return structural statistics and cutoff sensitivity.
``cutoff_window`` is the positive half-window in Angstrom used by
:meth:`cutoff_robustness` (default 0.10).
"""
from .robustness import format_cutoff_robustness
robustness = self.cutoff_robustness(window=cutoff_window)
lines = []
bar = "=" * 65
n_structs = len(self.atoms_list)
formula = self.atoms_list[0].get_chemical_formula(mode="hill")
n_atoms = len(self.atoms_list[0])
lines.append(f"\n{bar}")
lines.append(f" Structural Analysis: {formula} ({n_atoms} atoms)")
lines.append(f" N structures: {n_structs}")
lines.append(f" Cutoff mode: {self._cutoff_mode}")
if isinstance(self.cutoff, dict):
for pair in sorted(self.cutoff):
lines.append(f" {pair}: {self.cutoff[pair]:.2f} A")
else:
lines.append(f" All pairs: {self.cutoff} A")
lines.append(bar)
d = self.density()
lines.append(f"\n Density: {d['mean']:.2f} g/cm3; structure SD={d['std']:.2f}")
bd = self.bond_distances()
if bd:
lines.append("\n Bond distances: pooled bonds; Std is spread")
lines.append(f" {'Pair':<10} {'Mean (A)':>10} {'Std':>8} "
f"{'Min':>8} {'Max':>8} {'Count':>8}")
lines.append(f" {'-'*54}")
for pair, data in bd.items():
lines.append(
f" {pair:<10} {data['mean']:>10.3f} {data['std']:>7.3f} "
f"{data['min']:>8.3f} {data['max']:>8.3f} "
f"{data['count']:>8}")
cn = self.coordination()
if cn:
from collections import Counter
from .structure import is_bonding_pair
_comp = Counter(self.atoms_list[0].get_chemical_symbols())
bonding_cn = {}
nonbonded_cn = {}
for pair, data in cn.items():
s1, s2 = pair.split("-")
if is_bonding_pair(s1, s2, _comp):
bonding_cn[pair] = data
else:
nonbonded_cn[pair] = data
def _fmt(pair, data):
l1 = (f" {pair}: mean={data['mean']:.1f} (pooled); site SD="
f"{data['std']:.1f} [{data['min']},{data['max']}]")
parts = [f"CN={cn}: {pct:.1f}%"
for cn, pct in sorted(data["distribution"].items())
if pct >= 0.5]
l2 = " Fraction of sites: " + ", ".join(parts)
prevalence = ", ".join(
f"CN={cn}: {100 * fraction:.1f}%"
for cn, fraction in data["fraction_of_structures"].items())
return [l1, l2, " Fraction of structures (any such site): " + prevalence,
" " + _format_uncertainty(data["uncertainty"])]
if bonding_cn:
lines.append(f"\n Bonding coordination numbers:")
for pair, data in bonding_cn.items():
lines.extend(_fmt(pair, data))
# elements bonded to more than one partner type (O in IGZO:
# O-Ga + O-In + O-Zn) also get their total first-shell CN
partners = {}
for pair in bonding_cn:
s1, s2 = pair.split("-")
partners.setdefault(s1, []).append(s2)
multi = {s: ps for s, ps in partners.items() if len(ps) > 1}
if multi:
tot = self.total_coordination()
lines.append(f"\n Total coordination (all bonded partners):")
for s, ps in multi.items():
if s in tot:
lines.extend(_fmt(f"{s}-({'+'.join(ps)})", tot[s]))
if nonbonded_cn:
nb = {k: v for k, v in nonbonded_cn.items()
if v["mean"] > 0.05}
if nb:
lines.append(f"\n Non-bonded contacts:")
for pair, data in nb.items():
lines.extend(_fmt(pair, data))
if show_angles:
ba = self.bond_angles()
if ba:
lines.append("\n Bond angles: pooled angles; Std is spread")
lines.append(f" {'Triplet':<12} {'Mean (deg)':>10} "
f"{'Std':>8} {'Count':>8}")
lines.append(f" {'-'*42}")
for triplet, data in ba.items():
lines.append(
f" {triplet:<12} {data['mean']:>10.1f} "
f"{data['std']:>7.1f} {data['count']:>8}")
lines.append("\n Ensemble means (equal weight per independent structure):")
lines.append(" Density: " + _format_uncertainty(d["uncertainty"]))
for pair, data in bd.items():
lines.append(f" Bond {pair}: " + _format_uncertainty(data["uncertainty"]))
if show_angles:
for triplet, data in ba.items():
lines.append(f" Angle {triplet}: " + _format_uncertainty(data["uncertainty"]))
lines.append(" SD above describes spread; SEM and t intervals describe the mean.")
lines.append(" Intervals assume independent structures; n < 2 is unavailable.")
lines.append(format_cutoff_robustness(robustness))
lines.append(f"\n{bar}\n")
text = "\n".join(lines)
print(text)
return text
def _lookup_log_energies(self) -> dict[int, float]:
"""Locate the random_gen.log next to these structures and return a
``{file_index: e_per_atom_eV}`` map. Empty dict if no log found.
Searches the directory of the first file path, plus its parent
(handles ``random_opt/`` and ``random_initial/`` subdirs in the
v1.0.0rc2 layout). The log's `random_NNNN` indices map to the
sorted file list positionally (i.e. the i-th file → i-th log
entry by index).
"""
if not self._file_list:
return {}
from .energy import rank_from_log
# Search in current dir and one level up
first_dir = os.path.dirname(os.path.abspath(self._file_list[0]))
candidates = [
os.path.join(first_dir, "random_gen.log"),
os.path.join(os.path.dirname(first_dir), "random_gen.log"),
]
log_path = next((p for p in candidates if os.path.isfile(p)), None)
if log_path is None:
return {}
try:
result = rank_from_log(log_path)
except Exception:
return {}
# rank_from_log returns rows of (idx, e_total, e_per_atom, fmax, n, status)
log_by_idx = {row[0]: row[2] for row in result.get("rows", [])}
if not log_by_idx:
return {}
# Map structure-position -> log index by parsing "random_NNNN" out
# of each file name. Fall back to positional if no match.
import re
out = {}
for pos, path in enumerate(self._file_list):
m = re.search(r"random_(\d+)", os.path.basename(path))
log_idx = int(m.group(1)) if m else pos
if log_idx in log_by_idx:
out[pos] = log_by_idx[log_idx]
return out
[docs]
def per_structure_summary(self, cutoff_window=0.1) -> str:
"""
Analyse each structure individually and produce a comparison table.
Appends ensemble cutoff sensitivity over +/- ``cutoff_window`` A.
Returns
-------
str — formatted table (also printed).
"""
from .robustness import format_cutoff_robustness
robustness = self.cutoff_robustness(window=cutoff_window)
lines = []
bar = "=" * 80
formula = self.atoms_list[0].get_chemical_formula(mode="hill")
n_atoms = len(self.atoms_list[0])
lines.append(f"\n{bar}")
lines.append(f" Per-Structure Analysis: {formula} ({n_atoms} atoms)")
lines.append(f" N structures: {len(self.atoms_list)}")
lines.append(bar)
# Determine bonding pairs for CN
from collections import Counter
from .structure import is_bonding_pair
unique = sorted(set(self.atoms_list[0].get_chemical_symbols()))
comp = Counter(self.atoms_list[0].get_chemical_symbols())
bonding_pairs = [f"{s1}-{s2}" for s1 in unique for s2 in unique
if is_bonding_pair(s1, s2, comp)]
# If single element, use same-species
if not bonding_pairs:
bonding_pairs = [f"{unique[0]}-{unique[0]}"]
# Header
cn_headers = [f"CN({p})" for p in bonding_pairs[:3]]
header = f" {'#':<5} {'Density':>8} {'E/atom':>10}"
for h in cn_headers:
header += f" {h:>10}"
lines.append(f"\n{header}")
lines.append(f" {'-'*(len(header)-2)}")
# Per-structure data
all_densities = []
all_energies = []
all_cns = {p: [] for p in bonding_pairs[:3]}
ensemble_cn = self.coordination()
# v1.0.0rc2: if the structure files were produced by --random-gen
# (which writes per-step energies into random_gen.log) and the
# files themselves don't carry energy in their headers (e.g. VASP
# format strips energy on write), scan the sibling log so the
# report's "E/atom" column isn't full of N/A.
log_energies = self._lookup_log_energies()
for i, atoms in enumerate(self.atoms_list):
sa_single = StructureAnalyser([atoms], cutoff=self.cutoff)
# Density
d = sa_single.density()
dens = d["mean"]
all_densities.append(dens)
# Energy (per atom, using THIS structure's atom count)
n_at = len(atoms)
e_str = ""
for key in ['energy', 'Energy', 'potential_energy']:
if key in atoms.info:
e = atoms.info[key] / n_at
all_energies.append(e)
e_str = f"{e:.4f}"
break
if not e_str:
try:
e = atoms.get_potential_energy() / n_at
all_energies.append(e)
e_str = f"{e:.4f}"
except Exception:
# Final fallback: random_gen.log
e = log_energies.get(i)
if e is not None:
all_energies.append(e)
e_str = f"{e:.4f}"
else:
e_str = "N/A"
# CN for bonding pairs
row = f" {i:<5} {dens:>8.2f} {e_str:>10}"
for p in bonding_pairs[:3]:
cn_val = ensemble_cn.get(p, {}).get("per_structure", [None] * len(self.atoms_list))[i]
all_cns[p].append(cn_val)
row += f" {cn_val:>10.1f}" if cn_val is not None else f" {'N/A':>10}"
lines.append(row)
# Summary statistics
lines.append(f" {'-'*(len(header)-2)}")
mean_row = f" {'Mean':<5} {np.mean(all_densities):>8.2f}"
std_row = f" {'Std':<5} {np.std(all_densities):>8.2f}"
if all_energies:
mean_row += f" {np.mean(all_energies):>10.4f}"
std_row += f" {np.std(all_energies):>10.4f}"
else:
mean_row += f" {'N/A':>10}"
std_row += f" {'N/A':>10}"
for p in bonding_pairs[:3]:
values = [value for value in all_cns[p] if value is not None]
mean_row += f" {np.mean(values):>10.1f}" if values else f" {'N/A':>10}"
std_row += f" {np.std(values):>10.1f}" if values else f" {'N/A':>10}"
lines.append(mean_row)
lines.append(std_row)
lines.append(" Uncertainty of ensemble means:")
for name, values in [("Density", all_densities), ("E/atom", all_energies),
*[(f"CN({p})", all_cns[p]) for p in bonding_pairs[:3]]]:
lines.append(f" {name}: " + _format_uncertainty(summarize_structures(values)))
lines.append(format_cutoff_robustness(robustness))
lines.append(f"\n{bar}\n")
text = "\n".join(lines)
print(text)
return text
[docs]
def save_report(self, filepath, text=None, show_angles=True, cutoff_window=0.1):
"""Save the summary report to a text file.
Parameters
----------
filepath : str
Output file path.
text : str, optional
Pre-computed report text. If None, calls summary().
show_angles : bool
Include bond angles in the report (default True).
cutoff_window : float
Half-window in Angstrom for cutoff robustness (default 0.10).
Used only when generating report text here.
"""
if text is None:
text = self.summary(show_angles=show_angles, cutoff_window=cutoff_window)
# Create the parent directory like --save-plot does, so a report path
# in a not-yet-existing folder doesn't abort the run.
parent = os.path.dirname(os.path.abspath(filepath))
os.makedirs(parent, exist_ok=True)
with open(filepath, "w") as f:
f.write(text)
print(f" Report saved: {filepath}")
[docs]
def plot(self, **kwargs):
"""Generate and save analysis plots (RDF, CN distribution, bond angles).
Parameters
----------
output_dir : str
Directory for output files (default ".").
prefix : str
Filename prefix (default "analysis").
rdf_pairs : list of str, optional
Pairs to plot, e.g. ["Si-O", "O-O"]. If None, all pairs.
rmax : float, optional
Maximum radius for RDF. Auto-detected if None.
save_csv : bool
Also save raw data as CSV files (default True).
"""
return plot_analysis(self, **kwargs)
def parse_total_cn_spec(spec):
"""``"O"`` -> ("O", None); ``"O:Ga+In"`` -> ("O", ["Ga", "In"])."""
spec = str(spec).strip()
if ":" in spec:
centre, rest = spec.split(":", 1)
partners = [p.strip() for p in rest.replace(",", "+").split("+") if p.strip()]
return centre.strip(), partners
return spec, None
def format_total_cn(sa, specs):
"""Report block for requested total coordinations (``--total-cn``)."""
lines = ["\n Total coordination (requested):"]
for spec in specs:
centre, partners = parse_total_cn_spec(spec)
tot = sa.total_coordination(centre=centre, partners=partners)
if centre not in tot:
lines.append(f" {centre}: not present in the structure")
continue
d = tot[centre]
label = f"{centre}-({'+'.join(partners)})" if partners else f"{centre}-(all bonded)"
lines.append(f" {label}: mean={d['mean']:.1f} (pooled); site SD={d['std']:.1f} "
f"[{d['min']},{d['max']}]")
parts = [f"CN={cn}: {pct:.1f}%" for cn, pct in sorted(d["distribution"].items())
if pct >= 0.5]
lines.append(" Fraction of sites: " + ", ".join(parts))
lines.append(" Fraction of structures (any such site): " + ", ".join(
f"CN={cn}: {100 * fraction:.1f}%"
for cn, fraction in d["fraction_of_structures"].items()))
lines.append(" " + _format_uncertainty(d["uncertainty"]))
return "\n".join(lines)
def _format_uncertainty(uncertainty):
"""Concise scalar ensemble estimate, keeping unavailable errors explicit."""
u = uncertainty
mean = "n/a" if u["mean"] is None else f"{u['mean']:.4g}"
if u["sem"] is None:
return f"mean={mean}; SEM/t interval unavailable (n={u['n_structures']})"
return (f"mean={mean}; SEM={u['sem']:.3g}; {100 * u['confidence']:g}% t CI "
f"[{u['ci_low']:.4g}, {u['ci_high']:.4g}] (n={u['n_structures']})")