Utilities

User-facing helper functions used across the pipeline.

Common helpers

MD-dynamics builder, trajectory I/O, config merging, density helper (amorphgen.utils.common).

amorphgen.utils.common

Shared helpers used across all pipeline stages: cell manipulation, MD dynamics builder, temperature ramps, logging, trajectory I/O, config merging, and snapshot extraction.

Calculator-related functions are in amorphgen.utils.calculators.

amorphgen.utils.common.compute_density_gcm3(atoms)[source]

Compute density of an Atoms object in g/cm3.

Return type:

float

exception amorphgen.utils.common.DivergenceError[source]

Bases: RuntimeError

Invalid state or an MLIP safety limit exceeded during MD or relaxation.

Almost always a foundation-model MLIP going out-of-distribution in the high-temperature liquid regime, or too large a timestep. Raised eagerly (see assert_finite()) so a NaN/Inf never silently propagates into a saved structure or trajectory — a wrong-but-plausible result is worse than a clear failure.

amorphgen.utils.common.assert_finite(atoms, context='', step=None)[source]

Raise DivergenceError if the current energy/forces are non-finite.

Reads the energy/forces already computed for this step (the MD integrator and the optimiser both evaluate them every step, and ASE caches the result until the atoms change), so it adds no calculator call and is cheap enough to run every step. context and step are woven into the message to pinpoint where the divergence happened.

A calculator that raises (rather than returning NaN) is left alone — that is a different failure and must surface on its own, not be masked here.

Parameters:

context (str)

Return type:

None

amorphgen.utils.common.resolve_device(device)[source]

Resolve device="auto" in priority order: CUDA, MPS, then CPU.

Torch is an optional dependency (pulled in by the MLIP extras), so a torch-free install — random generation, analysis, or classical-potential pipelines — resolves auto to "cpu" instead of crashing on the import. Any explicit device string is passed through unchanged.

Parameters:

device (str)

Return type:

str

amorphgen.utils.common.make_cubic(atoms)[source]

Reshape the cell to a cube of equal volume, rescaling atom positions.

amorphgen.utils.common.cubic_cell_filter(atoms)[source]

Cell filter for cell_filter="cubic": volume relaxes, shape stays.

FrechetCellFilter with hydrostatic strain, in place of the ExpCellFilter that ASE deprecated in 3.23. Under hydrostatic strain the two return the same forces; exp_cell_factor=1 keeps ExpCellFilter’s scale for the cell rows (the virial, |P| V), which the optimisation loops’ max|force| < fmax test includes. Frechet’s default divides them by the number of atoms, which loosens the pressure criterion by that factor, to ~0.1 GPa at fmax = 0.01 eV/A.

amorphgen.utils.common.calculator_supports_stress(calc)[source]

Return True if calc advertises a stress tensor.

Variable-cell operations (NPT barostats, cell-filter optimisation) need the stress. The classical pair potentials (Lennard-Jones, Buckingham) implement only energy + forces, so this returns False for them.

Return type:

bool

amorphgen.utils.common.require_stress(calc, context)[source]

Raise a clear error if calc cannot provide stress for context.

Prevents an opaque PropertyNotImplementedError from surfacing deep inside ASE when a stress-less calculator is used with a barostat or a cell filter.

Parameters:

context (str)

Return type:

None

amorphgen.utils.common.build_md_dynamics(atoms, ensemble='NVT', T=300.0, timestep=1.0, friction=0.01, ttime=25.0, pfactor=None, external_stress=0.0, npt_method='berendsen', taup_factor=10.0, compressibility_GPa=100.0, rng=None, **kwargs)[source]

Create an NVT or NPT ASE dynamics object.

Parameters:

  • atoms (ase.Atoms) – Must already have a calculator attached.

  • ensemble (str) – "NVT" or "NPT".

  • T (float) – Temperature in Kelvin.

  • timestep (float) – Time step in fs.

  • friction (float) – Langevin friction coefficient (for NVT). The Langevin thermostat leaves the centre of mass free (fixcm=False); pass fixcm in kwargs to override.

  • ttime (float) – Thermostat time constant in fs. For "berendsen" it is taut; for "mtk" and "parrinello-rahman" it is the Nose-Hoover-chain time constant (ttime in the ASE API).

  • pfactor (float, optional) – Barostat coupling factor for "mtk" and "parrinello-rahman". If None, defaults to (ttime * taup_factor fs)**2 * compressibility_GPa GPa, giving a barostat taup_factor times slower than the thermostat (same spirit as the Berendsen taup). Ignored by "berendsen".

  • external_stress (float) – External pressure in GPa (for NPT). For "mtk" and "parrinello-rahman" this is converted to an isotropic stress tensor.

  • npt_method ({"berendsen", "mtk", "parrinello-rahman"}) –

    NPT integrator to use when ensemble == "NPT".

    • "berendsen" (default) — weak-coupling Berendsen barostat and thermostat (ase.md.nptberendsen.NPTBerendsen). Robust during 300 K -> 3000 K melt ramps; does not produce true canonical fluctuations, so heat capacities and isothermal compressibilities derived from fluctuations are incorrect. Averages are correct.

    • "mtk" — Martyna-Tobias-Klein Nose-Hoover-chain NPT (ase.md.nose_hoover_chain.IsotropicMTKNPT). Produces true canonical fluctuations. Recommended for the equilibration stages (2, 4, 6); may become unstable during rapid temperature ramps (stages 3, 5).

    • "parrinello-rahman" — Nose-Hoover + Parrinello-Rahman flexible-cell NPT (ase.md.npt.NPT). Allows the cell shape (not just volume) to change; useful for anisotropic glasses but requires upper-triangular cell.

    Ignored when ensemble == "NVT".

  • taup_factor (float, default 10.0) – Ratio of barostat coupling time to thermostat coupling time, i.e. taup = taup_factor * ttime. Larger values give a slower, more stable barostat — useful for damping cell-volume excursions during the 300 K -> 3000 K melt ramp. Applied to the Berendsen taup and to the MTK / Parrinello-Rahman barostat-time defaults.

  • compressibility_GPa (float, default 100.0) – Reference isothermal compressibility used by the Berendsen barostat as 1/(compressibility_GPa * GPa). The default (100 GPa) is intentionally soft and gives liquid-like responsiveness; for stiffer oxides (a-In2O3, a-Ga2O3, a-HfO2, bulk modulus ~150-300 GPa) using 200 GPa gives more realistic and more stable volume control. Ignored by "mtk" and "parrinello-rahman".

  • **kwargs – Extra arguments forwarded to the ASE dynamics class.

Return type:

ASE dynamics object

amorphgen.utils.common.parse_index_spec(spec, n_total=None)[source]

"80-90", "0,5,7-9" or a sequence of ints -> set of indices.

Ranges are inclusive. n_total (if given) bounds the result.

Parameters:

n_total (int | None)

Return type:

set[int]

amorphgen.utils.common.stage_rng(seed, stage, run_index=0)[source]

Per-stage, per-run NumPy Generator derived from the global seed.

None seed -> None (ASE’s default, unseeded, generator). The stream depends only on (seed, stage, run_index), so stage 4 of run 7 draws the same thermostat noise whatever ran before it or on another machine.

Parameters:
amorphgen.utils.common.scoped_run_index(local, scope=None, source='local')[source]

Seed index for one MD run: a local identity inside a per-source band.

local is the run’s own identity (its snapshot_NNNN number, or its position in the loop). source names where the enclosing scope came from, and scope is that scope’s number:

"local" no enclosing scope; the index is just local "batch" an explicit --run-index on a batch / ensemble run "slurm" SLURM_ARRAY_TASK_ID on a batch / ensemble run "pipeline" an explicit --run-index on the single-structure pipeline "pipeline-slurm" SLURM_ARRAY_TASK_ID on the single-structure pipeline

Two runs therefore share a seed stream only when they come from the same source with the same scope and the same local identity. Out-of-range values are refused rather than wrapped: folding snapshot_100003 onto snapshot_0003, or letting --run-index 100000 spill into the next band, would silently put two runs on one stream.

Parameters:
Return type:

int

amorphgen.utils.common.run_index_from_cwd()[source]

Index of a run_NNNN working directory (batch / ensemble modes).

Outside such a directory (a single-snapshot run writes straight into the work dir) the SLURM array task id is used when present, so array jobs sharing one --seed still draw different velocities and noise; else 0.

Return type:

int

amorphgen.utils.common.run_index_for(cfg)[source]

Run index for the MD seed stream: an explicit run_index in the config (set by batch_quench per run, or --run-index) beats the working-directory / SLURM inference.

Parameters:

cfg (dict)

Return type:

int

amorphgen.utils.common.resolve_ramp(T_start, T_end, T_step)[source]

Generate the list of temperatures for a ramp from T_start to T_end.

The ramp direction is taken from the endpoints, so only the magnitude of T_step matters — a mis-signed step (e.g. a positive step for a cooling ramp) can no longer produce an empty list or an infinite loop. Float steps are supported. T_end is always the final entry, even when the span is not an integer multiple of the step, and the ramp never overshoots past T_end. T_start itself is NOT in the list: the system already sits at T_start when the ramp begins, so the segments are the targets T_start+step, T_start+2*step, …, T_end. With steps_per_T MD steps per segment the realised rate then equals the configured one (an extra segment at T_start used to lower it by n/(n+1)).

Raises:

ValueError – If T_step has zero magnitude.

Parameters:
Return type:

list[float]

class amorphgen.utils.common.MDLogger(logfile, mode='w', step_offset=0)[source]

Bases: object

Per-step MD logger that writes to both a file and stdout.

Logs step number, time (ps), temperature (K), potential energy (eV), kinetic energy (eV), total energy (eV), and volume (ų).

Parameters:
log(dyn, atoms)[source]
close()[source]
class amorphgen.utils.common.TrajectoryWriter(filename, fmt='extxyz', append=False)[source]

Bases: object

Unified trajectory output supporting multiple formats.

Wraps ASE’s write() for extxyz/xyz/lammps-dump and ASE’s Trajectory for .traj binary format.

Parameters:
write(atoms=None)[source]
close()[source]
amorphgen.utils.common.attach_outputs(dyn, atoms, logfile, trajfile, fmt='extxyz', interval=100, append=False, step_offset=0, safety=None)[source]

Attach an MDLogger and TrajectoryWriter to dyn.

Parameters:

  • append (bool) – Continue existing log/trajectory files instead of truncating them (frame-level resume). The step-0 observer call of the resumed run is suppressed so the resume point is not written twice — the trajectory stays one-frame-per-interval-steps, which is what read_md_checkpoint() relies on to count elapsed steps.

  • (logger (Returns)

  • later. (traj_writer) so they can be closed)

  • logfile (str)

  • trajfile (str)

  • fmt (str)

  • interval (int)

  • step_offset (int)

amorphgen.utils.common.read_md_checkpoint(trajfile, interval=100)[source]

Last complete frame of an MD trajectory and the MD steps it represents.

The returned frame carries the MD momenta (extxyz stores them), so the stage can continue from it. A trajectory whose LAST frame was torn by a walltime kill is truncated to its complete frames (which are kept) rather than discarded whole.

Returns None when the file is missing, empty or unreadable.

Parameters:

  • trajfile (str)

  • interval (int)

amorphgen.utils.common.resume_md_stage(trajfile, resume, stage_label, legacy_trajfile=None)[source]

Shared frame-level-resume entry point for the MD stages (2-6).

Returns (checkpoint_atoms_or_None, elapsed_steps). Holds the resume invariants in ONE place (see also needs_velocity_init() and ramp_resume_position()) so the three stage modules cannot drift. legacy_trajfile lets ramp stages pick up a trajectory written under the pre-rename default name by an older AmorphGen version.

Parameters:

  • trajfile (str)

  • stage_label (str)

  • legacy_trajfile (str | None)

amorphgen.utils.common.needs_velocity_init(atoms, elapsed)[source]

Should the stage (re)draw Maxwell-Boltzmann velocities?

Fresh runs always do. Resumed runs keep the checkpoint’s momenta — unless the frame carries none (all-zero momenta cannot occur mid-MD, so zeros mean the trajectory format dropped them) and re-initialisation is the only option.

Parameters:

elapsed (int)

Return type:

bool

amorphgen.utils.common.ramp_resume_position(elapsed, steps_per_T, n_temps)[source]

Position in a temperature ramp after elapsed completed steps.

Returns (k0, offset): k0 full segments are done and offset steps of segment k0 — the caller skips segments < k0 and runs steps_per_T - offset for segment k0. When elapsed equals the ramp total, k0 == n_temps and the loop runs nothing: the stage output is then written from the checkpoint frame, which can lag the true final state by up to TRAJ_LOG_INTERVAL - 1 steps (the frames between write intervals are not recoverable) — physically negligible for equilibrium MD, but a resumed run is not byte-identical to an uninterrupted one.

Parameters:

  • elapsed (int)

  • steps_per_T (int)

  • n_temps (int)

amorphgen.utils.common.set_md_temperature(dyn, T)[source]

Change the target temperature of a running ASE dynamics object.

Langevin, NPTBerendsen and the Melchionna/NPT integrators expose set_temperature; IsotropicMTKNPT (npt_method “mtk”) does not, so its thermostat and barostat kT are updated directly. Used by the heating and cooling ramps.

Parameters:

T (float)

Return type:

None

amorphgen.utils.common.merge_config(defaults, overrides)[source]

Deep-merge overrides into a copy of defaults.

Parameters:

  • defaults (dict)

  • overrides (dict | None)

Return type:

dict

amorphgen.utils.common.stage_file(name, work_dir=None)[source]

Path a stage runner writes its file name to.

Relative names stay in the current directory, which is the run’s work dir once MeltQuenchPipeline or batch_quench has changed into it. A stage called on its own is given work_dir instead: name goes inside it (created if missing), unless name is absolute.

Parameters:

name (str)

Return type:

str

amorphgen.utils.common.extract_snapshots(traj_file, n_snapshots=20, select='uniform', output_dir='snapshots', burn_in_frames=None, output_format='extxyz', *, timestep_fs=0.5, frame_stride=100, decorrelation_distance=None, report_path=None, resume=False)[source]

Extract snapshot frames from a trajectory file.

Parameters:

  • traj_file (str) – Path to the trajectory file.

  • n_snapshots (int) – Number of snapshots to extract.

  • select (str) – Selection strategy: "uniform" (evenly spaced), "last" (final n frames), or "decorrelated" (autocorrelation and species diffusion determine the minimum spacing).

  • output_dir (str) – Directory for output files.

  • burn_in_frames (int or None) – Number of leading frames to skip before sampling. Useful for discarding the equilibration period at the start of an MD trajectory. Sampling indices run over the closed interval [burn_in_frames, n_frames - 1]. Raises ValueError if burn_in_frames >= n_frames. None uses automatic burn-in for decorrelated selection and zero for uniform/last selection.

  • output_format (str, default "extxyz") – Output file format. Accepted values: "extxyz", "xyz" (both write extended XYZ with a .xyz extension), "vasp" (POSCAR-style), "cif", "traj".

  • timestep_fs (float, int) – Integration timestep and saved-frame stride for time diagnostics.

  • frame_stride (float, int) – Integration timestep and saved-frame stride for time diagnostics.

  • decorrelation_distance (float or None) – Length scale in Angstrom for the species displacement correlation proxy; inferred from nearest neighbours when omitted.

  • report_path (str or None) – JSON report filename, with a sibling .txt summary. Decorrelated extraction defaults to output_dir/snapshot_sampling.json.

  • resume (bool) – Protect correspondence with existing downstream runs by comparing saved selection settings, indices and source-frame fingerprints before writing. Enable this when resuming existing quenches.

Returns:

Paths to extracted snapshot files.

Return type:

list of str

Radii, classification and auto-derivation

Shannon ionic / Cordero covalent / Goldschmidt metallic radii, bonding-type classification, automatic minsep + density + target-CN derivation, and charge-balance oxidation-state inference (amorphgen.utils.radii). These helpers are exercised in Tutorial 2 (“Zero-config random structure generation”; see Tutorials) and underpin the automatic parameter choices described in Random structure generation.

amorphgen.utils.radii

Atomic radii data, bonding classification, minsep calculation, and density estimation for random structure generation.

This module contains: - Shannon ionic radii (CN=4 and CN=6) - Metallic radii (Goldschmidt CN=12) - Pauling electronegativities - Element classification (nonmetal, metalloid) - Elemental solid densities - Bond type classification - Minsep calculation from radii - Cell volume / density estimation

amorphgen.utils.radii.infer_oxidation_state(sym, composition)[source]

Infer the oxidation state of sym in composition by charge balance.

Sums the charge of the anions (anion_elements(), charges from ANION_CHARGES) and assigns the remaining positive charge equally across cations of the queried element when only one cation type is present, or for any element that appears as the only non-anion type. A nonmetal that acts as a cation here (cation_nonmetals(): P in Li3PO4, S in Li2SO4) is solved like any other; the Li of Li2SO4 is +1, not the +5 that counting S as an anion gave. Returns None for an anion, for Te (a metalloid on the anion table, left to the highest-state default whether it is the anion or the cation), for the covalent hydrogenated networks (a-Si:H, where charge balance against H- made Si50H50 Si+1), and when the inference is ambiguous (multiple cation species, non-integer oxidation state, or no anions at all), letting the caller fall back to a default Shannon entry.

Parameters:
Return type:

int | None

amorphgen.utils.radii.cation_nonmetals(composition)[source]

Nonmetals that act as CATIONS in this compound.

The centre of an oxoanion (C in a carbonate, N in a nitrate, P in a phosphate, S in a sulfate, Se in a selenate, Cl, Br or I in a halate) and the H of a hydroxide or an acid salt bond to the anions like any other cation, although NONMETALS puts them with the anions. Charge balance picks them out: the elements anion_elements() promotes, plus C and P, which that rule never makes anions. They count as cations only where an anion more electronegative than them is present and they balance the charge better as cations than as C4- or P3-: the P of Li3PO4 or Li3PS4 is a cation, the carbide C of an oxycarbide (SiOC) or a carbonitride stays an anion, and so do the C and P of carbides, phosphides and a-C:H. Only elements with a tabulated cation radius qualify.

composition may be a mapping of counts (preferred) or a bare set of symbols, in which case one of each is assumed.

Return type:

frozenset

amorphgen.utils.radii.auto_target_cn(composition)[source]

Auto-detect target coordination numbers and tolerance from composition.

Returns (target_cn, cn_tolerance) tuple.

Rules:

  • Metalloids (Si, Ge, B): CN=4, tolerance=0 (strict tetrahedral)

  • Pure Si/Ge: CN=4, tolerance=0

  • Hydrogenated networks (a-Si:H, a-C:H, a-SiC:H): C, Si, Ge CN=4 and H CN=1 (one X-H bond), tolerance=0

  • Metals in nitrides: CN=4, tolerance=0 (wurtzite tetrahedral)

  • Metals in halides/sulfides: CN=6, tolerance=0 (strict octahedral)

  • Metals in oxides: CN=5, tolerance=1 (flexible 4-6 range)

  • Nonmetal cations (cation_nonmetals()): the ligand count of their oxoanion, 4 in PO4 3-, SO4 2- and ClO4 -, 3 in CO3 2-, NO3 - and IO3 -, and 1 for the H of a hydroxide

Parameters:

composition (dict) – e.g. {“In”: 32, “O”: 48}

Returns:

(target_cn, cn_tolerance) – target_cn: dict or None cn_tolerance: int (0=strict, 1=flexible)

Return type:

tuple

amorphgen.utils.radii.classify_bond(sym_a, sym_b, composition=None)[source]

Classify a pair of elements by bonding type.

Uses an element-type classification rule (nonmetal / metalloid / metal membership) with a Pauling-electronegativity refinement: pairs that the type rules would call “ionic” but whose Pauling electronegativity difference is small (Δχ < 1.0) are reclassified as “covalent”. This catches III–V semiconductors (GaAs Δχ=0.37, AlAs ~0.6) and related compounds where the chemistry is predominantly covalent and Shannon ionic radii give unrealistically small inter-atomic distances.

Returns one of: “ionic”, “covalent”, “metallic”.

With composition (a {symbol: count} mapping), the roles charge balance gives the two elements in that compound come first. A nonmetal cation (cation_nonmetals(): P in a phosphate, S in a sulfate, H in a hydroxide) bonds to the anions like any other cation, so the pair is “ionic”. Two cations of a compound with anions only meet across an anion: “cation-cation” when one is a nonmetal cation, and never “ionic” otherwise, since the Δχ rule is for bonded pairs (Na-B in a borate and K-Si in a silicate are then “covalent”, as Na-Si already is).

Parameters:
Return type:

str

amorphgen.utils.radii.get_ionic_radius(sym, cn=None, oxidation_state=None)[source]

Get Shannon ionic radius for an element.

Parameters:

  • sym (str) – Element symbol.

  • cn (int, optional) – Coordination number. If provided, use CN-specific radius. Defaults to CN=6.

  • oxidation_state (int, optional) – Specific oxidation state to look up (e.g.5 for Sb in Sb_2O_5). If supplied and present in the Shannon table the matching radius is returned. If supplied but absent, falls through to the default selection. When None (default), anions return their most-negative state and cations return their highest positive state - usually the right choice for binary oxides / halides / nitrides where the higher state dominates (TiO_2, V_2O_5, WO_3 etc.).

Return type:

float or None

amorphgen.utils.radii.get_metallic_radius(sym)[source]

Get metallic radius for an element. Returns None if unavailable.

Parameters:

sym (str)

Return type:

float | None

amorphgen.utils.radii.get_effective_radius(sym)[source]

Get the most appropriate radius for volume estimation.

  • Metals -> metallic radius

  • Anion-formers (O, S, Cl, …) -> Shannon ionic radius

  • Metalloids (Si, Ge, …) -> metallic radius if available

  • Fallback -> covalent radius

Parameters:

sym (str)

Return type:

float

amorphgen.utils.radii.default_minsep(symbols, scale=0.85, target_cn=None)[source]

Build a minsep dict for all element pairs.

Uses bonding-type-aware radii. When target_cn is provided, uses CN-specific Shannon radii (e.g. CN=4 for Si in SiO2).

Pairs are classified in the compound (classify_bond with the composition), so a nonmetal that charge balance makes a cation, the P of a phosphate or the C of a carbonate, is kept at bonding distance from its anions (P-O 1.26 A, C-O 1.06 A) and away from the other cations.

A hydrogenated network (a-Si:H, a-C:H) has no anions: see _hydrogenated_network_minsep().

Parameters:

  • symbols (list of str) – One symbol per atom, or a {symbol: count} mapping; the counts decide which nonmetals act as cations.

  • scale (float) – Fallback scale factor (default 0.85).

  • target_cn (dict, optional) – e.g. {“Si”: 4, “O”: 2}. Uses CN-specific radii when available.

Returns:

Minimum separations keyed as “A-B” with A <= B alphabetically.

Return type:

dict

amorphgen.utils.radii.estimate_density(composition, amorphous_factor=0.8)[source]

Estimate density from elemental solid densities.

Returns density in g/cm3, or None if any element lacks density data.

Parameters:

  • composition (dict)

  • amorphous_factor (float)

Return type:

float | None

amorphgen.utils.radii.get_packing_factor(cls, composition=None)[source]

Packing factor for a material class — composition-aware where needed.

The single dispatch point between the static PACKING_FACTORS table and classes whose packing depends on the composition itself: oxyhalide, which interpolates between its halogen-free base class and halide by halogen fraction, and hydrogenated_network, which takes the packing factor of its H-free host. New mixed-anion families (oxynitrides, oxysulfides) should add their interpolation here rather than special-casing estimate_cell_length(). Note the obvious oxynitride interpolation is currently a no-op — metal_oxide and nitride share pf 0.52 — so no oxynitride branch exists yet; add one only with calibration data that motivates it.

Parameters:

  • cls (str)

  • composition (dict | None)

Return type:

float

amorphgen.utils.radii.anion_elements(composition)[source]

Which elements act as ANIONS in this compound, decided by charge balance.

Membership of the anion table is not enough on its own: tellurium is the anion in CdTe and the cation in TeO2, sulfur the anion in ZnS and the cation in a sulfate, hydrogen the anion in LiH and the cation in a hydroxide. So every element on the anion table starts as a candidate anion, and while the cations present cannot supply the positive charge the candidates demand, the least electronegative candidate is promoted to cation. This reproduces the chemistry without a lookup table of exceptions: La2O2S balances with sulfur as an anion and keeps it, H2SO4 does not and promotes first hydrogen then sulfur, leaving the sulfate O as the anion.

A promotion happens only when it brings the compound CLOSER to charge balance, which is what keeps real cells intact without any tolerance setting. In a chalcogen-rich glass or a doped cell, promoting the major anion would overshoot far past neutrality (Ge20S10Se70: a gap of 80 before, 480 after; F-doped silica: 2 before, 254 after), so every anion is kept: the mixed chalcogen glasses keep Ge-Se and Ge-Te, F-doped silica keeps Si-O, an O impurity in NaCl keeps Na-Cl and LiPON keeps P-N. Where a promotion really is the chemistry, it improves the balance and happens: TeO2, a tellurite, a sulfate, a nitrate, a hydroxide.

Known limitation: an element that is a cation in one site of a polyatomic group and an anion in another cannot be both. In a thiosulfate the central sulfur is cation-like and the terminal sulfur is an anion; the rule keeps both as anions, so the S-O bonds are not counted. Rare in amorphous work.

composition may be a mapping of counts (preferred) or a bare set of symbols, in which case one of each is assumed.

Return type:

set

amorphgen.utils.radii.estimate_cell_length(composition, target_density=None, packing_factor=None, density_scale=1.0)[source]

Estimate cubic cell length from composition via sphere packing.

A single class-aware sphere-packing model handles all amorphous compositions:

\[V_{\text{cell}} = \frac{\sum_i N_i \cdot \tfrac{4}{3}\pi r_i^3} {f_{\text{pack}} \cdot s},\]

where \(r_i\) is the radius for element \(i\) chosen according to the compound’s bonding regime (Shannon ionic for ionic compounds, Cordero covalent for covalent compounds, Goldschmidt metallic for alloys; see _radius_for_density()). \(f_{\text{pack}}\) is the class-aware packing fraction stored in PACKING_FACTORS, and \(s\) is the user-supplied density_scale (default 1.0).

If target_density is supplied, the cell is sized directly from it and density_scale / packing_factor are ignored.

Parameters:

  • composition (dict)

  • target_density (float | None)

  • packing_factor (float | None)

  • density_scale (float)

Return type:

float

amorphgen.utils.radii.format_auto_derive_summary(composition, target_cn, minsep, est_density, cell_length)[source]

Single-line summary of the auto-derivation chain for the random_gen log.

Captures the chemistry-informed decisions AmorphGen makes from a bare composition: material class, oxidation states (where inferable), per-cation target coordination, per-pair bond classification + Pauling Δχ + minsep, and the final estimated density / cell length. Designed to fit on one log line for typical 2–4 element systems (e.g. Ga₂O₃ ≈ 200 chars).

Format:
[auto-derive] <formula> → <class>, OS{…}, CN{…},

minsep{pair:val class [Δχ=val] | …}, ρ=<val> g/cm³, L=<val> Å

For ionic pairs (where the Pauling Δχ < 1.0 refinement fires) the Δχ value is shown so the user can see exactly why the pair was classified as ionic vs covalent. For metallic and anion-packing pairs the Δχ is omitted (the type rule alone decides those).

Parameters:
Return type:

str

Format conversion

Format-conversion helper used by amorphgen --convert and the convert-mode YAML config (amorphgen.utils.convert).

amorphgen.utils.convert

Public file-format conversion utility.

Convert one structure file or every ASE-readable file in a directory to a target format (xyz / vasp / cif). VASP outputs are sorted by species so the resulting POSCAR is clean.

Usable from CLI (amorphgen --convert PATH --format vasp), Python API (from amorphgen import convert), or via a convert: block in a YAML config file (amorphgen --config convert.yaml).

amorphgen.utils.convert.convert(input_path, output_format='vasp', output_dir=None, sort=True, verbose=True)[source]

Convert a structure file or directory of files to output_format.

Parameters:

  • input_path (str) – Path to a single ASE-readable structure file, or to a directory containing one or more such files. Directory globs match *.xyz, *.extxyz, *.vasp, *.cif and POSCAR*.

  • output_format (str) – Target format key. Allowed values: "xyz" / "extxyz" (both write ASE extended XYZ to .xyz), "vasp" (POSCAR to .vasp), "cif".

  • output_dir (str, optional) – Directory to write converted files into. If None, defaults to "<input>_<format>" for directory inputs, or the parent directory of input_path for single-file inputs.

  • sort (bool, default True) – For VASP outputs, sort atoms by species so the POSCAR is clean. Ignored for other formats.

  • verbose (bool, default True) – Print one progress line per file plus a summary footer.

Returns:

Paths to the converted output files, in input order.

Return type:

list of str

Raises:

ValueError – If an output would overwrite any input, including through a symbolic or hard link, or if two outputs alias the same file. No files are written when such a collision is found.

Examples

>>> from amorphgen import convert
>>> convert("snapshots/", output_format="vasp",
...         output_dir="snapshots_vasp/")
['snapshots_vasp/snapshot_0000_frame00000.vasp', ...]

MD equilibration convergence analysis

Functions for assessing whether an MD equilibration stage has converged (amorphgen.utils.equilibration). Useful as a post-pipeline quality check: did the high-T equilibration in your melt-quench actually reach steady-state, or did you quench from a still-thermalising liquid?

The high-level entry point is convergence_report, which generates energy-vs-time, block-averaging, MSD, temperature, and time-window-RDF plots in one call. Lower-level functions are exposed for users who want to assemble custom diagnostics.

amorphgen.utils.equilibration

Equilibration convergence analysis for AmorphGen MD trajectories.

Provides functions to assess whether an MD equilibration run has converged:

  1. Energy vs time (running average + drift detection)

  2. Block averaging (quantitative equilibration test)

  3. RDF in time windows (structural convergence)

  4. Mean Square Displacement (diffusion / liquid vs glass)

  5. Coordination number vs time

  6. Temperature vs time

Usage

from amorphgen.utils.equilibration import convergence_report

# Quick all-in-one convergence report
report = convergence_report(
    "stage4_eq.xyz",
    timestep_fs=1.0,
    T_target=3000,
)

# Or individual analyses
fig_energy, drift = plot_energy_convergence(traj, timestep_fs=1.0)
is_eq, block_data = block_average_test(traj, n_blocks=4)
fig_msd, D_dict = plot_msd(traj, timestep_fs=1.0)

Notes

  • Two input modes: trajectory file (.xyz/.traj) or log file (.log). Log files are faster (no atomic positions to load) but only support energy, temperature, and block-average analyses. Trajectory files are needed for MSD, RDF, and CN analyses.

  • Default timestep is taken from AmorphGen’s default_config.py (0.5 fs), so the time axes are right for a default run. Trajectories are written every TRAJ_LOG_INTERVAL MD steps; pass frame_stride if yours differ.

amorphgen.utils.equilibration.parse_md_log(logfile)[source]

Parse an AmorphGen MD stage log file.

Expects whitespace-separated columns:

Step Time(ps) T(K) Epot(eV) Ekin(eV) Etot(eV) Vol(A^3)

Lines starting with ‘Step’, ‘-’, or ‘->’ are skipped as headers/markers.

Returns:

  • dict with keys (‘step’, ‘time_ps’, ‘T_K’, ‘Epot_eV’, ‘Ekin_eV’,)

  • ’Etot_eV’, ‘Vol_A3’ — each a numpy array.

Parameters:

logfile (str)

Return type:

dict

amorphgen.utils.equilibration.running_average(data, window)[source]

Compute running average with given window size.

Parameters:
Return type:

ndarray

amorphgen.utils.equilibration.extract_energies(source, n_atoms=None)[source]

Extract potential energies (eV) and per-atom energies.

Parameters:

  • source (str or list) – Trajectory file, log file (.log), or list of Atoms.

  • n_atoms (int, optional) – Required if source is a .log file.

Returns:

  • energies (np.ndarray — total potential energy (eV))

  • energies_per_atom (np.ndarray — energy per atom (eV/atom))

Return type:

tuple[ndarray, ndarray]

amorphgen.utils.equilibration.plot_energy_convergence(source, timestep_fs=0.5, window_ps=0.5, per_atom=True, n_atoms=None, ax=None, frame_stride=100)[source]

Plot potential energy vs time with running average and linear drift.

Parameters:

  • source (str or list) – Trajectory file, log file (.log), or list of Atoms.

  • timestep_fs (float) – MD timestep in femtoseconds (defaults to AmorphGen’s DEFAULT_CONFIG value, 0.5 fs). Pass your run’s actual timestep if it differs.

  • window_ps (float) – Running average window in picoseconds.

  • n_atoms (int, optional) – Required if source is a .log file.

  • per_atom (bool)

  • frame_stride (int)

Returns:

  • fig (matplotlib Figure)

  • drift_eV_per_ps (float) – Linear drift in energy. Should be ~0 if equilibrated.

amorphgen.utils.equilibration.block_average_test(source, n_blocks=4, discard_fraction=0.1, n_atoms=None)[source]

Split trajectory into blocks and compare mean energies.

A system is considered equilibrated if the block means agree within 2x the standard error of the mean (SEM).

Parameters:

  • source (str or list) – Trajectory file, log file (.log), or list of Atoms.

  • n_blocks (int) – Number of equal blocks to split the production phase into.

  • discard_fraction (float) – Fraction of initial trajectory to discard (thermalisation).

  • n_atoms (int, optional) – Required if source is a .log file.

Returns:

  • is_equilibrated (bool)

  • block_data (dict) – Keys: block_means, overall_mean, overall_std, sem, max_deviation, threshold, is_equilibrated.

Return type:

tuple[bool, dict]

amorphgen.utils.equilibration.plot_block_averages(source, n_blocks=4, discard_fraction=0.1, timestep_fs=0.5, n_atoms=None, ax=None, frame_stride=100)[source]

Visualise block averaging: block means vs overall mean +/- 2*SEM.

Parameters:

  • n_blocks (int)

  • discard_fraction (float)

  • timestep_fs (float)

  • n_atoms (int | None)

  • frame_stride (int)

amorphgen.utils.equilibration.compute_msd(traj, timestep_fs=0.5, by_element=True, frame_stride=100)[source]

Compute MSD from trajectory using unwrapped positions.

Handles both orthorhombic and non-orthorhombic cells via fractional coordinate unwrapping. Displacements are measured relative to the centre of mass, so a drift of the whole system does not count as diffusion.

timestep_fs is the MD integration timestep. AmorphGen writes one trajectory frame every frame_stride MD steps (TRAJ_LOG_INTERVAL, default 100), so the real time between consecutive frames is timestep_fs * frame_stride. Pass frame_stride=1 for a trajectory that stores every step. Getting this wrong rescales the time axis (and hence the fitted diffusion coefficient) by frame_stride.

Parameters:

  • traj (str, list of Atoms, or Trajectory)

  • timestep_fs (float)

  • by_element (bool) – If True, return MSD per element type in addition to total.

  • frame_stride (int)

Returns:

  • time_ps (np.ndarray)

  • msd_dict (dict — keys are element symbols (+ “all”),) – values are MSD arrays (A^2).

Return type:

tuple[ndarray, dict]

amorphgen.utils.equilibration.plot_msd(traj, timestep_fs=0.5, ax=None, frame_stride=100)[source]

Plot MSD vs time per element. Fits diffusion coefficient D.

D is fitted from the linear regime (last 50% of trajectory) using MSD = 6*D*t (3D diffusion).

Returns:

  • fig (matplotlib Figure)

  • D_dict (dict — diffusion coefficients (cm^2/s) per element.) – D > 1e-6 cm^2/s indicates liquid/diffusive behaviour. D < 1e-6 cm^2/s indicates frozen/glass behaviour.

Parameters:

  • timestep_fs (float)

  • frame_stride (int)

amorphgen.utils.equilibration.plot_temperature(source, timestep_fs=0.5, T_target=None, ax=None, frame_stride=100)[source]

Plot instantaneous temperature vs time.

Parameters:

  • source (str or list) – Log file (.log) or trajectory / list of Atoms.

  • T_target (float, optional) – Target temperature (K). If given, shows expected fluctuation band.

  • timestep_fs (float)

  • frame_stride (int)

amorphgen.utils.equilibration.plot_rdf_time_windows(traj, pairs=None, n_windows=4, rmax=4.0, nbins=100, timestep_fs=0.5, frame_stride=100)[source]

Overlay partial RDFs from different time windows.

If the RDFs from all windows overlap, the structure is equilibrated.

Parameters:
amorphgen.utils.equilibration.compute_cn_vs_time(traj, centre, neighbour, cutoff=None, window_size=50, timestep_fs=0.5, frame_stride=100)[source]

Compute average coordination number vs time using a sliding window.

Parameters:
Return type:

tuple[ndarray, ndarray, ndarray]

amorphgen.utils.equilibration.plot_cn_vs_time(traj, pairs, cutoffs=None, window_size=50, timestep_fs=0.5, ax=None, frame_stride=100)[source]

Plot CN vs time for multiple centre-neighbour pairs.

Parameters:

  • pairs (list of (centre, neighbour, expected_cn)) – e.g. [(“Si”, “O”, 4.0), (“O”, “Si”, 2.0)]

  • cutoffs (list of float, optional) – Bond cutoffs per pair. None = auto from covalent radii.

  • window_size (int)

  • timestep_fs (float)

  • frame_stride (int)

amorphgen.utils.equilibration.convergence_report(source, timestep_fs=0.5, T_target=None, n_atoms=None, pairs_rdf=None, pairs_cn=None, cn_cutoffs=None, rmax=4.0, n_blocks=4, output_dir=None, prefix='convergence', frame_stride=100)[source]

Generate a comprehensive convergence report.

Parameters:

  • source (str, list of Atoms, or Trajectory) – Trajectory file (.xyz/.traj), log file (.log), or list of Atoms. Log files are faster but only provide energy/temperature (no MSD, RDF, or CN analysis).

  • timestep_fs (float) – MD timestep in femtoseconds (default 0.5).

  • T_target (float, optional) – Target temperature (K) for temperature check.

  • n_atoms (int, optional) – Required if source is a .log file.

  • pairs_rdf (list of (str, str), optional) – Element pairs for RDF windows. Auto-detected if None.

  • pairs_cn (list of (str, str, float), optional) – Element pairs + expected CN for CN tracking. e.g. [(“Si”, “O”, 4.0), (“O”, “Si”, 2.0)]

  • cn_cutoffs (list of float, optional) – Bond cutoffs for CN pairs. None = auto from covalent radii.

  • rmax (float) – Max distance for RDF (must not exceed half cell length).

  • n_blocks (int) – Number of blocks for block average test.

  • output_dir (str, optional) – If provided, save all plots as PNG and a text report here.

  • prefix (str) – Filename prefix for saved plots.

  • frame_stride (int)

Returns:

report – Keys include: n_frames, n_atoms, total_time_ps, elements, energy_drift_eV_per_atom_per_ps, block_test_passed, block_data, diffusion_coefficients_cm2_s, summary_text.

Return type:

dict

Other utility modules

  • Calculators, multi-backend calculator factory (get_calculator).

  • amorphgen.utils.classical, Lennard-Jones and Buckingham+Coulomb calculators; loaded via get_calculator("lennard-jones" | "buckingham", ...).