Analysis
AmorphGen’s --analyse mode computes structural descriptors for an
ensemble of amorphous structures and writes a figure and a CSV
file for each quantity.
amorphgen --analyse --input-dir my_structures/ --save-plot plots/
That single command computes density, partial radial distribution functions, coordination distributions, bond-angle distributions, and a per-structure density violin, saved as up to four PNG figures plus CSV companions. The density plot and CSV require at least two structures; angle output requires valid bond-angle triplets.
Directory input reads .xyz, .extxyz, .vasp and .cif files in that
directory. Files with the same stem count once, in that format priority
order. Each file contributes its last frame; to analyse a trajectory as an
ensemble, first extract snapshots or pass a list of frames to the Python API.
Screening and analysis inclusion
Screening records a label independently of whether a structure is excluded from analysis. Run the default screens with:
amorphgen --analyse --input-dir ensemble/ --screen \
--screening-output analysis/screening --save-plot analysis/
--screen enables coordination, crystal-like order, close-contact and
relaxation-convergence checks. All screens default to exclude: false, so
labels alone retain structures in the analysis. Density and energy need
material-specific bounds and are enabled through YAML. For example:
analysis:
screening:
coordination:
allowed: auto
max_fraction: 0.05
crystal_like:
qbar6_threshold: 0.3
min_neighbors: 4
max_fraction: 0.2
close_contacts:
threshold_frac: 0.7
exclude: true
density:
min: 2.1
max: 2.3
unconverged:
exclude: false
Run this configuration with --config examples/screening.yaml; a screening
mapping enables screening without requiring --screen. analysis.screening: true selects the default screens; false disables screening unless
overridden by --screen. In a mapping, only listed screens are enabled, and
enabled: false disables an individual screen. exclude: true applies
exclusion for that screen’s labels. The
example thresholds above illustrate silica and must be chosen for the
material and neighbour shell being assessed.
Screen |
What triggers a label |
|---|---|
|
More than |
|
The ordered atom fraction exceeds |
|
Any periodic interatomic separation is strictly below its threshold. By default this is |
|
Density lies outside the inclusive |
|
Stored energy per atom lies outside the inclusive |
|
Explicit |
Automatic coordination sets reuse auto_target_cn() from
amorphgen.utils.radii, the same composition-based targets used by
--random-gen. Each returned target expands to the integer set from
target - tolerance to target + tolerance, inclusive. Species without an
inferred target are unassessed; this is not a full list of acceptable
coordinations for every species. Use explicit sets when needed, for example
to assess oxygen as well as silicon in silica. The coordination fraction
uses assessed sites as its denominator. Automatic neighbour cutoffs are
resolved once using the full candidate ensemble and remain fixed for both
screening and analysis of the retained structures.
Crystal-like screening shares this cutoff by default. An optional
crystal_like.cutoff selects a separate order shell using the same scalar,
pair, auto or auto-rdf forms as --cutoff; it is also resolved once over
the full candidate ensemble. Screening settings are independent of the
optional --bond-order descriptor’s CLI options.
Unavailable metrics produce a <screen>_unavailable label. They exclude a
structure only if that screen has exclude: true. This preserves the
difference between a measured failure and missing evidence. Structure files
that do not retain energy or relaxation metadata can therefore receive
unavailable labels even when their geometry is valid.
Newly relaxed structures, including --random-gen --relax outputs, retain
the optimiser’s convergence evidence in extended XYZ metadata. CIF and VASP
outputs use a companion <structure>.relaxation.json file; keep it beside
the structure. The loader verifies the structure file’s hash before using
the companion, so replaced files do not acquire stale convergence results.
Earlier files without evidence remain unavailable. Starting new MD clears
the previous relaxation status; a later relaxation records its own result.
The CLI prints a generated / passed / labelled / analysed table and writes three files, even when every candidate is excluded:
Count |
Meaning |
|---|---|
|
Candidates supplied to this analysis, after normal input loading; it does not reconstruct historical generation attempts. |
|
Candidates with no triggered or unavailable screen labels. |
|
Candidates with at least one label, whether retained or excluded. |
|
Retained candidates whose analysis completed. |
|
Candidates removed by at least one screen configured with |
Labelled and analysed counts can overlap. For example, ten candidates with
three labelled and one excluded can give generated=10, passed=7,
labelled=3, analysed=9, excluded=1. Analysis describes the retained
population; changing exclusions can change its means and uncertainty.
--screening-output PREFIX writes PREFIX.json,
PREFIX_structures.csv and PREFIX_summary.csv, containing the full
decisions, individual structure records and count table respectively.
The default prefix is <work-dir>/screening; YAML can set
analysis.screening_output. The Python API exposes the same decisions via
StructureAnalyser.screen() and a retained analyser via .screened();
see Structure analysis.
Spread and uncertainty of the mean
Each input structure is one independent sampling unit. Pooled site, bond and
angle spreads remain available as the legacy mean/std fields and the
explicit pooled_mean/pooled_std fields. They describe structural disorder;
they are not errors on an ensemble mean. Density retains its population
spread over structures in std.
Scalar results include an uncertainty summary calculated from one mean per
structure, giving every contributing structure equal weight. This matters
when cell sizes, numbers of bonds or numbers of angles differ. Its fields are:
Field |
Meaning |
|---|---|
|
Values in input order; |
|
Equal-weight mean and sample SD between structures ( |
|
Standard error, sample SD divided by the square root of the number of contributing structures |
|
Student-t confidence interval with |
|
Percentile interval of resampled structure means; 1,000 draws and seed 0 by default |
|
Contributing structures and all supplied structures |
|
Number contributing to each curve bin (scalar for a scalar descriptor) |
With fewer than two contributing structures, SEM and interval bounds are
None (blank in CSV). Missing descriptors are excluded, not replaced with
zero. A present species with no neighbors has zero coordination; a missing
central species has no coordination observation. Repeating sites within one
structure does not increase the independent sample count.
rdf(), averaged_rdf(), structure_factor(), structure_factor_direct(),
total_correlation() and angle_distribution() return per-structure curves
and pointwise uncertainty on common grids. Angle histograms are normalized
within each structure before averaging. S(q) and T(r) transformations use each
structure’s own composition and density before averaging, preserving their
covariance. Direct S(q) bins with no reciprocal vectors are missing; their
n_per_point can be smaller than the ensemble size. The n_per_bin field
counts reciprocal vectors, not independent samples.
The weighted Fourier-transform S(q) needs at least two atoms of every species
for its same-species RDF normalization; use structure_factor_direct() for
singleton dopants. An unestimable partial is never silently treated as an
observed zero curve.
sa = StructureAnalyser("ensemble/")
cn = sa.coordination("Si-O")["Si-O"]
print(cn["pooled_std"], cn["uncertainty"]["sem"])
rdf = sa.rdf(confidence=0.95, n_bootstrap=2000, seed=42)
angles = sa.angle_distribution("O-Si-O", bins=90, seed=42)
tr = sa.total_correlation(weighting="xray", seed=42)
# T(r)'s primary uncertainty is for T_r; other curves have separate summaries.
print(tr["curve_uncertainty"]["G_r"]["ci_low"])
Bootstrap resampling selects whole structures, retaining correlations
between bins. The shaded bands are pointwise, not simultaneous confidence
bands for the entire curve. Set n_bootstrap=0 to skip bootstrap draws in the
curve APIs. These intervals quantify sampling of independent configurations;
they do not include force-field bias, finite-cell error, cutoff selection, or
transform-parameter uncertainty. Correlated trajectory frames require
independent sampling or a separate correlation/block analysis before using
these intervals.
RDF, angle, S(q), and T(r) exports include companion
*_uncertainty.csv, *_per_structure.csv, and *_uncertainty.json files.
The JSON records confidence level, seed, resampling count, and sampling unit;
the CSV includes contributing counts, SEM, t bounds and bootstrap bounds.
Raw angle CSV rows also carry the structure index. The density plot shows a
t interval for the mean, alongside individual structures.
analysis_statistics.json retains core scalar descriptors and their aligned
per-structure observations; analysis_statistics.csv separates pooled spread
from structure means, sample SD, SEM and t intervals.
Coordination and oxygen-speciation outputs distinguish fraction_of_sites
(pooled sites, between 0 and 1) from fraction_of_structures (the fraction of
all input structures containing at least one site in that category).
Structure prevalence is not a distribution: categories can overlap and need
not sum to one. Crystal-like order and dimer/connectivity reports make the
same distinction. Per-structure site-fraction summaries estimate the
equal-weight mean site fraction, which can differ from the pooled fraction.
Optional descriptors (rings, Voronoi, oxygen speciation, bond order, voids,
elastic moduli, vibrational DOS and energy) also retain per-structure values
and named uncertainty summaries. Existing void Monte Carlo errors remain
separate from uncertainty across structures.
Reference validation uses the structure-weighted mean and its t interval.
A confidence interval admitting both in-range and out-of-range values is
inconclusive, even when its mean is inside the reference range. A finite
mean with no estimable interval is also inconclusive; an absent descriptor
remains n/a. Reports include intervals and count inconclusive verdicts.
Declared tolerances and ensemble convergence
This section describes uncertainty and sample-size planning for an existing
ensemble. Repeatedly stopping when a Student-t planning interval passes does
not preserve its nominal coverage. For generation with an anytime-valid
stopping rule, use --random-gen --relax --engine torchsim --until-converged,
with predeclared population support bounds and tolerances; see
Generate until converged.
Declare an absolute tolerance for the uncertainty of each descriptor’s
ensemble mean, in its native units. A tolerance of 0.02 for density
means a 95% Student-t interval half-width no greater than 0.02 g/cm³; it is
neither a relative percentage nor a bound on the spread of individual structures.
amorphgen --analyse --input-dir ensemble/ \
--tolerance density=0.02 \
--tolerance coordination.Si-O=0.05 \
--tolerance bond_angle.O-Si-O=1.0 \
--convergence-confidence 0.95 --save-plot analysis/ \
--save-report analysis.txt
Each --tolerance NAME=VALUE enables convergence reporting. Use --convergence
alone to inspect available descriptor names and their uncertainty before
declaring tolerances. Default names include density, coordination.PAIR,
total_coordination.ELEMENT, bond_distance.PAIR and bond_angle.TRIPLET.
Pair ordering is significant: silica bond distance is bond_distance.O-Si,
whereas coordination has separate coordination.Si-O and coordination.O-Si.
Unknown names and nonpositive or nonfinite tolerances are rejected.
Selected optional analyses add named descriptors, such as
bond_order.ordered_fraction, voids.accessible_fraction, sq.total,
sq.Si-O with --sq-partials, and tr.T_r. A tolerance on rdf.total or
rdf.PAIR also computes that RDF for convergence. For curve descriptors the
tolerance must hold at every point: the plotted quantity is the largest
pointwise half-width, not a simultaneous confidence band for the entire curve.
Optional descriptors require their corresponding analysis flags.
The equivalent YAML entries live under analysis:
analysis:
convergence: true
tolerances:
density: 0.02
coordination.Si-O: 0.05
convergence_confidence: 0.95
convergence_max_structures: 1000000
CLI tolerance declarations replace YAML tolerances of the same name and
preserve the others. --convergence-max-structures sets the upper search bound
for forecasts; it never truncates the input ensemble.
from amorphgen.analysis import (
StructureAnalyser, format_convergence_report, save_convergence_report,
)
sa = StructureAnalyser("ensemble/")
report = sa.convergence_report({"density": 0.02, "coordination.Si-O": 0.05})
print(format_convergence_report(report))
save_convergence_report(report, "analysis/", save_pdf=True)
# Include any aligned per-structure descriptor, including a whole curve.
rdf = sa.rdf(n_bootstrap=0)
curve_report = sa.convergence_report(
{"rdf.total": 0.05}, descriptors={"rdf.total": rdf["uncertainty"]},
sizes=[2, 5, 10, 20, 50, 100],
)
The curves are independent of generation order. For complete scalar data, AmorphGen plots
where \(c\) is the chosen confidence and \(s_N\) is the sample standard deviation of all \(N\) structures. For \(2 \leq n \leq N\), this is the exact root-mean-square Student-t half-width over all subsets of size \(n\), because their mean sample variance equals \(s_N^2\). No randomized shuffling, seed, or generation-order prefixes enter the calculation. A vector descriptor takes the maximum of these componentwise RMS values, not the RMS of the subset-wise maxima. At \(n=N\) the curve equals the observed full-ensemble half-width; beyond \(N\) it is an extrapolation. The interval formula follows the NIST Student-t confidence interval.
Missing observations are excluded rather than treated as zero. If a component appears in \(k\) of \(N\) structures, planning at size \(n\) uses \(\lfloor nk/N\rfloor\) contributing observations with the observed sample variance. This is an availability-adjusted approximation, not an exact all-subset result. Any component with fewer than two observations makes that descriptor’s uncertainty and forecast unavailable; empty bins are retained.
The report returns met, not_met, insufficient_data, or undeclared
per descriptor and overall. Undeclared descriptors do not decide the overall
status. It estimates the smallest total ensemble size satisfying all declared
tolerances under unchanged variance and availability, and subtracts the
current size to give the additional structures needed. A met tolerance needs
zero additional structures; an estimate beyond the search bound is reported
as exceeds_max_structures, with no invented finite forecast.
Forecasts assume independent structures and representative, stable variance and missingness. They do not account for correlated trajectory frames, force-field bias, finite-size error or undiscovered rare configurations. Zero observed variance yields a zero estimated half-width once two values exist, but does not establish zero population variance. Reassess the report as new independent structures arrive.
With --save-plot, exports include analysis_convergence.json (full report,
strict JSON with missing values as null), analysis_convergence.txt,
analysis_convergence_summary.csv, and analysis_convergence_curves.csv
(one row per descriptor, planned size and point). Each descriptor has its
own figure, for example analysis_convergence_density.png, with an optional
PDF. Figures distinguish the observed endpoint, declared tolerance, solid
planning curve and dashed projection through the estimated target. CSV
columns retain counts, confidence, units in the summary, forecast status,
and pointwise uncertainty so the decisions can be reproduced.
Recipes
Common cases:
For an ensemble with at least two structures and valid angle triplets, this writes four standalone figures:
amorphgen --analyse \
--input-dir hybrid_ga2o3/final/ \
--save-plot plots/ \
--save-pdf
Output:
plots/
├── analysis_rdf.{png,pdf,csv} # partial RDFs (all pairs)
├── analysis_cn.{png,pdf,csv} # coordination distribution
├── analysis_angles.{png,pdf,csv} # bond-angle distributions
└── analysis_density.{png,pdf,csv} # per-structure density violin
The CSV files contain RDF values, coordination percentages, raw angle observations and per-structure densities so you can re-plot in any tool.
Four elements, ten element pairs, three cation sizes and a shared oxygen. One command reports and plots all of it:
amorphgen --analyse \
--input-dir igzo_final/ \
--sq --sq-partials --pair-panels --save-pdf \
--total-cn O --total-cn "O:In+Ga" \
--save-report report.txt \
--save-plot plots/
The report header shows the cutoff in force for every pair. With the
default auto-rdf each pair gets its own value from the first minimum
of its g(r):
Cutoff mode: auto (RDF)
Ga-O: 2.03 A
In-O: 2.47 A
O-Zn: 2.25 A
...
One cutoff for all pairs is the thing to avoid here: 2.47 Å (the In–O
value) would count second-shell oxygens around Ga and lift the Ga–O
coordination from 3.9 to 4.2. To override one pair and keep auto-rdf
for the rest, use --cutoff "In-O=2.6".
The coordination part of the report separates bonds from contacts:
Bonding coordination numbers:
Ga-O: mean=3.9 +/- 0.3 [3,4]
In-O: mean=5.1 +/- 0.6 [4,6]
Zn-O: mean=3.9 +/- 0.5 [3,5]
O-Ga: ... O-In: ... O-Zn: ...
Total coordination (all bonded partners):
O-(Ga+In+Zn): mean=3.2 +/- 0.5 [2,4]
Total coordination (requested):
O-(all bonded): mean=3.2 +/- 0.5 [2,4]
O-(In+Ga): mean=2.3 +/- 0.8 [0,4]
Non-bonded contacts:
Ga-In: ... In-In: ... O-O: ...
Cation–O pairs are bonds; cation–cation and O–O pairs are second-shell
contacts listed apart and excluded from the coordination and the angles.
The O-(Ga+In+Zn) line appears on its own whenever an element has more
than one bonded partner type; --total-cn adds any centre and partner
set you name (O for all bonded partners, O:In+Ga for a subset).
--sq-partials prints the first peak of each Faber-Ziman partial
S_ab(q) and writes the partials next to the total; --pair-panels puts
every pair in its own panel for g(r) and for S_ab(q).
Output:
plots/
├── analysis_rdf.{png,pdf,csv} # partials; CSV also has g(r)_Total
├── analysis_rdf_panels.{png,pdf} # one panel per pair
├── analysis_cn.{png,pdf,csv} # Ga-O, In-O, Zn-O and O-(Ga+In+Zn)
├── analysis_cn_total.{png,pdf,csv} # the --total-cn requests
├── analysis_sq.{png,pdf,csv} # total S(q); CSV has s_<pair> columns
├── analysis_sq_partials.{png,pdf} # partials on one axis
├── analysis_sq_partials_panels.{png,pdf} # one panel per pair
├── analysis_angles.{png,pdf,csv}
└── analysis_density.{png,pdf,csv}
Python equivalent:
from amorphgen.analysis import StructureAnalyser
sa = StructureAnalyser("igzo_final/", cutoff="In-O=2.6") # rest auto-rdf
sa.summary()
sa.total_coordination(centre="O", partners=["In", "Ga"])
sq = sa.structure_factor_direct(weighting="xray", sigma_q=0.05, partials=True)
sq["partials"]["In-O"]
sa.plot(output_dir="plots/", save_pdf=True,
pair_panels=True, total_cn=["O", "O:In+Ga"])
from amorphgen.analysis.plotting import plot_sq
plot_sq(sq, output_dir="plots/", save_pdf=True, pair_panels=True)
Add a full text report alongside the plots:
amorphgen --analyse \
--input-dir hybrid_ga2o3/final/ \
--save-report report.txt \
--save-plot plots/ \
--save-pdf
The report shows: density mean ± std, bond distances (mean / std / count per pair), coordination numbers (with distribution histograms), bond angles (mean / std / count per triplet).
To also see the breakdown per structure (one row per file with
density, energy, CN), add --per-structure:
amorphgen --analyse \
--input-dir hybrid_ga2o3/final/ \
--per-structure
If the structure files don’t carry energies (VASP, CIF), AmorphGen looks
for random_gen.log alongside the files and in their parent directory
to fill in the E/atom column. Keep the original random_NNNN filenames
so the log entries can be matched to the correct structures.
Compare your ensemble against literature ranges with automatic match/concern/fail scoring:
amorphgen --analyse \
--input-dir hybrid_ga2o3/final/ \
--reference examples/reference_a_Ga2O3.yaml \
--save-report report.txt \
--save-plot plots/ \
--save-pdf
The report adds a section like:
Validation: a-Ga2O3
Descriptor Computed Expected Units Verdict
----------------------------------------------------------------------
Density 4.37 [4.70, 5.10] g/cm^3 fail
Bond Ga-O 1.91 [1.85, 1.95] A match
CN Ga-O 4.42 [4.00, 4.80] match
Angle Ga-O-Ga 116.8 [110.00, 130.0] deg match
Angle O-Ga-O 108.2 [100.00, 115.0] deg match
----------------------------------------------------------------------
Summary: 4 match, 0 concern, 1 fail (out of 5 metrics)
AmorphGen ships reference YAMLs for a-Ga₂O₃, a-SiO₂, a-GeO₂, a-HfO₂ and
a-IrO₂ in examples/ (reference_a_<system>.yaml). The GeO₂ file
notes that the neutron partial structure factors of Salmon et al. (2005)
allow a partial-by-partial comparison with --sq --sq-weighting neutron.
Write your own for other systems by following the same schema.
For overlaying Random vs Hybrid vs DFT-reference (or any combination),
use the Python API. There’s no single CLI flag for this yet;
compare_ensembles() is the entry point:
from amorphgen.analysis import EnsembleSpec, compare_ensembles
compare_ensembles(
ensembles=[
EnsembleSpec("DFT-PBE0", "prb_ensemble/*.cif"),
EnsembleSpec("Random", "random_inputs/*.vasp"),
EnsembleSpec("Hybrid", "hybrid_run/final/*.xyz"),
],
rdf_pairs=[("Ga-O", "-"), ("Ga-Ga", "--"), ("O-O", ":")],
cn_top_key="Ga-O",
cn_bot_key="O-Ga",
angle_keys=[("O-Ga-O", "-"), ("Ga-O-Ga", "--")],
exp_density=(4.78, 4.84),
output_dir="comparison/",
prefix="ga2o3",
)
Output: comparison/ga2o3_rdf.{png,pdf,csv},
comparison/ga2o3_coordination.{...},
comparison/ga2o3_angles.{...},
comparison/ga2o3_density.{...}, same layout as --analyse --save-plot, but each figure overlays all listed ensembles with
distinct colours from the Okabe-Ito palette.
See Validation for a Ga₂O₃ comparison and the available reference data.
Which contacts count as bonds?
The bonding coordination report, total coordination, bond angles and CN plot share the same bonding rule. In compounds with anions, only cation–anion pairs count as bonds; charge balance determines the anions. Without anions, the radii classification decides: same-element pairs count in a single-element system or a metal-rich alloy (at least 70 % metal atoms), and unlike-element pairs count when classified as ionic, covalent or metallic.
Hydrogenated group-IV networks contain only C, Si and/or Ge plus H, with at most one H per host atom. Their host keeps the bonding rules of the H-free composition: Si–Si in a-Si:H, Ge–Ge in a-Ge:H and C–C in a-C:H count as bonds. In a-SiC:H only Si–C host pairs count, and in a-SiGe:H only Si–Ge host pairs count. Every host–H pair counts as a bond; H–H never does. These pairs still have to lie within their distance cutoffs. Metal hydrides, hydroxides, compositions with other elements and those with more H than host atoms retain the usual rules.
Cutoff
The cutoff defines what counts as a “first-shell” bond and affects
coordination, bond-length statistics, and bond-angle triplets. The default
is auto-rdf, which finds the first minimum of each partial
RDF, the standard convention in neutron-diffraction analysis of
glasses. The minimum is read at a resolution of 0.25 Å, so a flat step
or a noise dip on the falling side of the first peak, common in small
cells, does not end the shell. Where g(r) is zero over a range, the
cutoff goes in the middle of it.
Cutoff mode |
When to use |
|---|---|
|
All structural analysis. Robust across systems with broad bond-length distributions (a-Si, a-HfO₂, chalcogenides). |
|
Legacy. Uses minsep from Shannon/Cordero/Goldschmidt radii. Fast but can under-count coordination for systems with long first-shell bonds. |
Numeric, e.g. |
Single cutoff (in Å) for all pairs. Useful for tight-bonded covalent networks. |
Per-pair overrides, e.g. |
Fix the pairs you name and keep |
Dict via YAML or API |
|
Inspect the RDF and the reported cutoffs before interpreting coordination. An explicit cutoff can help when a first minimum is ambiguous. One number for every pair is the option to avoid in a multi-cation oxide: for a-IGZO the In–O first minimum sits at 2.47 Å and Ga–O at 2.03 Å, and a single 2.47 Å cutoff raises the Ga–O coordination from 3.9 to 4.2 by admitting second-shell oxygens. The report header lists the cutoff in force for every pair, so a per-pair override is easy to check. These numerical values are an example, not fixed cutoffs for every IGZO ensemble.
Cutoff robustness
Every summary reports how sensitive the pair contacts and directional
coordination are to the selected cutoffs. The default half-window is 0.1 Å:
the report shows five coordination means at offsets −0.10, −0.05, 0,
+0.05 and +0.10 Å, plus the change from the lowest to the highest cutoff.
Set the half-window with --cutoff-window or analysis.cutoff_window in YAML:
amorphgen --analyse --input-dir structures/ --cutoff-window 0.15 \
--save-report report.txt --save-plot analysis/
For each unordered element pair, the near-cutoff share is the number
of contacts whose inclusion changes between the lower and upper endpoints,
divided by the number included at the upper endpoint. Counts pool the
structures and count each undirected periodic contact once; different
periodic images count separately, as they do in coordination. A pair with
no contacts at the upper endpoint has an undefined share, printed as n/a.
This denominator makes the reported fraction local to the candidate
neighbour shell, rather than to all possible atom pairs in the cell.
The selected cutoffs are resolved once, including automatic RDF minima.
The sweep adds the same offset to every positive pair cutoff, clipping
negative values to zero; zero cutoffs stay zero throughout the sweep.
Automatic cutoffs are not refitted during the sweep. Each point
uses the ordinary coordination boundary rules: a contact must satisfy
distance <= pair cutoff and distance < largest cutoff.
Coordination remains directional: Si–O is O neighbours per Si, while O–Si is Si neighbours per O. The text shows means pooled over central sites. The returned data also retain per-structure means and an ensemble mean that weights structures equally, so unequal structure sizes need not be treated as equal site populations. A larger near-cutoff share or coordination change signals greater sensitivity to the chosen shell boundary; these are sensitivity measures, not confidence intervals.
from amorphgen.analysis import (
StructureAnalyser, format_cutoff_robustness, save_cutoff_robustness,
)
sa = StructureAnalyser("structures/", cutoff="auto-rdf")
report = sa.cutoff_robustness(window=0.1, points=5)
print(format_cutoff_robustness(report))
save_cutoff_robustness(report, output_dir="analysis/")
sa.summary(cutoff_window=0.1)
window must be finite and positive; points must be an odd integer of at
least three so the grid includes the selected cutoff. Summary methods and
sa.plot(cutoff_window=...) use five points. With --save-plot, the default
CSV export also writes analysis_cutoff_robustness.json,
analysis_cutoff_robustness_pairs.csv and
analysis_cutoff_robustness_coordination.csv. In Python,
sa.plot(save_csv=False) suppresses these exports.
Total coordination
For elements bonded to more than one partner type (O in IGZO, bonded to
Ga, In and Zn) the report adds a Total coordination block with the
first-shell count over all bonded partners, next to the per-pair O–Ga,
O–In and O–Zn entries. The label reads centre first, then the partners in
brackets, the same order as the per-pair entries (O-Ga is O with Ga
around it). To choose the centre and the partners yourself, use
--total-cn (repeatable) or the YAML list total_cn::
amorphgen --analyse --input-dir DIR --total-cn O --total-cn "O:In+Ga"
O counts every bonded partner; O:In+Ga counts only the named ones.
From Python, sa.total_coordination(centre="O", partners=["In", "Ga"])
returns the same statistics.
Structure factor S(q) and simulated XRD
AmorphGen provides two structure-factor methods: a direct sum over the reciprocal-lattice vectors of a periodic cell, and a Fourier transform of the radial distribution function. The CLI uses the direct method by default. Both support X-ray, neutron and unweighted totals.
The direct method avoids truncating g(r), but still has finite-cell and sampling limits. The FT method can show termination ripples and altered peak heights; there is no universal correction factor between the two. See S(q) and XRD: methodology notes for the conventions and limitations.
From the CLI:
amorphgen --analyse --input-dir DIR --sq --sq-weighting neutron --save-plot plots/
This writes analysis_sq.png and analysis_sq.csv. With the default
smoothing, the CSV contains q_invA, s_q, s_q_raw and n_per_bin.
Add --sq-partials to print the first peak of each Faber–Ziman partial,
append s_<A-B> columns to the CSV and write analysis_sq_partials.png.
The partials describe the geometry and do not depend on the weighting.
From Python, for an ensemble of GeO₂ structures:
from amorphgen.analysis import StructureAnalyser
sa = StructureAnalyser("geo2_ensemble/")
sq = sa.structure_factor_direct(weighting="neutron", sigma_q=0.05,
partials=True)
sq["partials"]["Ge-O"]
For each nonzero reciprocal-lattice vector, the implementation returns Faber–Ziman-normalized scattering:
Here \(c_\alpha=N_\alpha/N\), \(\langle f\rangle=\sum_\alpha c_\alpha f_\alpha\) and \(\langle f^2\rangle=\sum_\alpha c_\alpha f_\alpha^2\). Subtracting the self-scattering term gives the high-q limit \(S(q)\to1\) for uncorrelated positions. The raw intensity divided by \(N\langle f\rangle^2\) has a different limit for mixtures.
With cell vectors as rows of \(\mathbf A\), the sampled vectors are
Values are averaged into shells by \(|\vec G|\), pooling vectors across
structures with the same composition. n_per_bin records their count;
empty shells contain NaN. For a cubic cell of side \(L\), the smallest
nonzero q is \(2\pi/L\). Increasing nq adds bins but cannot add independent
low-q information.
sigma_q (--sq-smooth) applies Gaussian smoothing weighted by each
shell’s vector count. Smoothing can broaden features, so compare with
s_q_raw. The CLI defaults to 0.05 Å⁻¹; the Python method
defaults to 0 (no smoothing). Empty shells remain NaN.
With partials=True, partials are averaged over the same shells. Their
weighted sum reproduces the total exactly for q-independent weights
(such as neutron scattering lengths). For X-rays the weights vary
within a shell, so recombining the binned partials using weights at bin
centres is approximate.
Use --sq-method ft or sa.structure_factor(weighting="xray") for the
FT method. Its Python default weighting is "unweighted"; the direct
method and CLI default to "xray".
The normalized structure factor is not a raw diffractometer intensity. For an X-ray coherent-scattering profile, first restore the self-scattering term and form-factor scale, then map momentum transfer to angle:
xrd_pattern() performs this conversion independently for each structure,
using its own composition, before averaging. Cu-Kα is the default wavelength.
Form factors are evaluated at q-bin centres for the direct method; use narrow
bins and inspect the counts when comparing intensities.
from amorphgen.analysis import StructureAnalyser, save_xrd_pattern
sa = StructureAnalyser("ga2o3_ensemble/")
xrd = sa.xrd_pattern(wavelength=1.5406, qmax=8.0, nq=400)
save_xrd_pattern(xrd, output_dir="plots/", save_pdf=True)
# xrd["two_theta"] in degrees, xrd["intensity"] in electron²/atom
# xrd["per_structure"] and xrd["uncertainty"] retain ensemble information
amorphgen --analyse --input-dir ga2o3_ensemble/ --xrd \
--xrd-wavelength 1.5406 --xrd-qmax 8 --xrd-nq 400 --save-plot plots/
The API also supports method="ft" with an optional rmax. The default
qmax is min(15, 4π/λ) Å⁻¹; inaccessible requested q ranges are rejected.
The angular grid comes from the calculated q grid and is not uniformly
spaced in 2θ. Intensities are sampled values, without an integration Jacobian
or peak-height normalization. Unsampled direct shells remain unavailable.
Optional sigma_q smooths each coherent intensity curve before averaging;
raw values are retained. No smoothing is applied by default.
The output is a coherent-scattering profile per atom before instrument and
sample corrections. Polarization, geometry, absorption, background and
resolution must match the measurement. A powder-diffraction
Lorentz–polarization factor applied to S(q) alone is not a general prediction
of an amorphous sample’s measured trace. Ensemble bands do not include these
systematic effects. For experimentally reduced Faber–Ziman S(Q), compare
against S(q) directly using matching weights and normalization.
Both methods combine partial structure factors using
The sum counts ordered pairs, so unlike-element pairs occur twice. Different weights change how the partial peaks and troughs contribute; they do not guarantee a particular peak height.
|
Per-element factor |
Use |
|---|---|---|
|
q-dependent neutral-atom form factor |
X-ray total scattering |
|
Tabulated coherent scattering length |
Neutron total scattering |
|
1 for every species |
Geometric structure comparisons |
X-ray form factors use the five-Gaussian parametrization of Waasmaier & Kirfel (1995):
The table covers neutral atoms H–Cf. Its fitted values approach the atomic number at q = 0; the code uses the q-dependent values throughout. Anomalous scattering corrections are not included.
Neutron weights use the built-in common-element table from Sears (1992). Unsupported species raise an error. Scattering lengths can be negative and depend on isotope; the API chooses weights from element symbols, so changing an atom’s mass does not select an isotope-specific length.
Use the same normalization, scattering weights, q range and smoothing when comparing calculations and experiments. A summary of the relevant conventions is Keen (2001). The direct reciprocal-space calculation and the FT of g(r) are described in S(q) and XRD: methodology notes.
The repository’s TestSqNormalisation tests in test/test_analysis.py
check the high-q limit, the FT integrand, scattering-table entries,
q-dependent X-ray weighting, smoothing and neutron partial recombination.
These are implementation checks; they do not establish universal accuracy
bounds or validate every material against experiment.
Limitation |
How to assess it |
|---|---|
Finite cell and sparse low-q shells |
Inspect |
Finite-r truncation in the FT method |
Vary |
Smoothing and finite q-bin width |
Compare raw data and narrower bins; report |
Neutral-atom X-ray weights; element-based neutron weights |
Check whether anomalous or isotope-specific scattering matters for the experiment. |
Thermal motion |
Use representative configurations at the comparison temperature; an extra Debye–Waller correction may double-count sampled motion. |
Instrument and sample effects |
Apply corrections appropriate to the measured quantity and geometry. |
Cite AmorphGen for the software and the primary literature for the scattering conventions and tables; see Scattering methods and citations.
Example: a-Ga₂O₃
For a directory of Ga₂O₃ configurations with the same composition:
from amorphgen.analysis import StructureAnalyser
sa = StructureAnalyser("ga2o3_ensemble/", cutoff="auto-rdf")
sq = sa.structure_factor_direct(weighting="xray", qmax=12.0, nq=300,
sigma_q=0.05, partials=True)
# sq["s_q_raw"] retains the unsmoothed total for comparison.
Peak positions and intensities depend on the supplied structures, cell size and settings. The same API applies to other compositions; select weights and a q range appropriate to the reference data.
Compare measured S(q) or T(r)
Supply a whitespace or CSV file with two columns (coordinate, value), or three columns (coordinate, value, one-sigma measurement uncertainty):
amorphgen --analyse --input-dir structures/ \
--experiment-sq measured_sq.csv --experiment-skiprows 1 \
--sq-weighting xray --sq-method direct --sq-qmax 12 --sq-nq 300 \
--sq-fit-range 1.5 10 --save-plot comparison/ --save-report report.txt
amorphgen --analyse --input-dir structures/ \
--experiment-tr measured_tr.dat --sq-weighting neutron \
--tr-qrange 0.5 20 --tr-window lorch --tr-fit-range 1 8 \
--save-plot comparison/
The file options enable the corresponding calculation. Comment lines start
with #; use --experiment-skiprows for an uncommented header. Files with
extra columns require --experiment-columns 0 2 3 (zero-based coordinate,
value, sigma), or two indices when sigma is unavailable. Coordinates are
sorted together with values and errors. Duplicate coordinates, nonfinite
values and nonpositive sigma are rejected. q is in Å⁻¹, r in Å, S(q) is
dimensionless, and T(r) = 4πrρg(r) is in Å⁻². No convention or unit conversion
is inferred from a filename or column label. Both files may be supplied in
one run; shared loader options apply to both. YAML analysis keys use
underscores, for example experiment_sq, sq_fit_range, and sq_qmax;
explicit CLI values take precedence.
Python exposes the loader and comparison independently, so a calculated curve can be reused:
from amorphgen.analysis import (
StructureAnalyser, load_experiment, compare_experiment,
format_experiment_report, save_experiment_comparison,
)
sa = StructureAnalyser("structures/")
measured = load_experiment("measured_sq.csv", kind="sq", skiprows=1)
sq = sa.structure_factor_direct(qmax=12, nq=300, weighting="xray", sigma_q=0.05)
comparison = compare_experiment(sq, measured, x_range=(1.5, 10))
print(format_experiment_report(comparison))
save_experiment_comparison(comparison, output_dir="comparison/", save_pdf=True)
# Or load, calculate and compare in one call:
comparison = sa.compare_experiment(
"measured_tr.dat", kind="tr", weighting="neutron", x_range=(1, 8),
calculation_options={"qmin": 0.5, "qmax": 20, "window": "lorch"},
)
Each structure is linearly interpolated onto the measured coordinates before averaging. Extrapolation and interpolation across missing calculated bins are excluded. The report counts excluded points and records the number of contributing structures at each retained point. Match the scattering weights, normalization, temperature, resolution and, for T(r), q range and window. Changing these choices can alter the residuals independently of model quality.
Metrics use the same retained measured points. For residual Δ = calculated minus measured, RMSE = √mean(Δ²), MAE = mean(|Δ|), and bias = mean(Δ). Rw = √[ΣwΔ² / Σw(measured)²], expressed as a fraction, with w = 1/σ² when measurement uncertainties are supplied and w = 1 otherwise. A zero denominator leaves Rw unavailable. χ² = Σ(Δ/σ)² and reduced χ² = χ²/N are reported only with measured sigma; no scale or offset is fitted. These diagonal-error statistics are descriptive for correlated points, particularly Fourier-transformed T(r); no p-value is inferred.
The shaded bands are pointwise Student-t confidence intervals of the
equal-weight ensemble mean. They require at least two independent contributing
structures and do not include instrument or model systematic error. Measurement
error bars remain separate. Whole-structure bootstrap bounds, standard
deviations and SEM are also exported. The CLI saves
analysis_experiment_sq.{json,csv,txt,png} or analysis_experiment_tr.*, plus
PDF with --save-pdf. Figures show the overlay and residuals; JSON retains
per-structure curves and reproducibility metadata, and CSV contains values,
residuals, counts and interval bounds.
Ring sizes and search coverage
amorphgen --analyse --input-dir silica/ --rings Si-O \
--ring-cutoff 2.0 --ring-max-size 16 \
--save-report rings.txt --save-plot rings/
The network uses the first element as nodes and the second as bridges;
single-element networks use direct bonds. A ring size counts network nodes
(for example, Si sites for Si–O), and periodic paths must close at the same
image. The default pair is selected by electronegativity. Set it explicitly
for mixed chemistries. --ring-cutoff overrides the analyser’s resolved pair
cutoff for this calculation only; --ring-max-size defaults to 12.
Each undirected network edge contributes its shortest cycle, if one closes within the limit. Counts describe edges, not unique rings: a lone six-node polygon contributes six observations of size six. Multiple bridges between the same pair of node images collapse to one network edge. Thus this projected network is not an enumeration of atom-level cycles.
Reports give the mean, population standard deviation, observed size range,
and counts of resolved and unresolved edges. An unresolved edge has no
closure within the chosen limit; it may belong to a longer ring. Increase
the limit to check convergence. Percentages use only resolved edges, while
ring_edge_fraction uses all network edges. No edges gives an undefined
coverage; no resolved rings gives undefined size summaries (JSON null).
Ensemble counts pool edges, while uncertainty summaries give each structure
equal weight.
--save-plot writes analysis_rings.json, the existing size/count/percent
CSV and PNG, analysis_rings_structures.csv for per-structure summaries,
and analysis_rings_per_structure.csv for per-structure distributions.
--save-pdf adds a PDF plot.
Crystal-like order and the largest ordered cluster
Use --bond-order to look for residual or newly formed crystal-like regions
in a quenched ensemble, including phase-change systems such as GeTe:
amorphgen --analyse --input-dir gete_mq/final/ --bond-order \
--order-cutoff 3.5 \
--qbar6-threshold 0.3 --order-min-neighbors 4 \
--save-report gete_report.txt --save-plot gete_plots/
This geometry-only descriptor reports each atom’s Steinhardt \(q_6\) and Lechner–Dellago \(\bar q_6\), the fraction of atoms classified as ordered, and the largest connected ordered cluster in each structure. For atom \(i\) with neighbour shell \(N(i)\), the complex spherical-harmonic coefficients are
Lechner–Dellago averaging includes the central atom and its neighbours, before taking the rotationally invariant norm:
These follow Steinhardt, Nelson and Ronchetti (1983) and Lechner and Dellago (2008). Averaging the scalar \(q_6\) values would give a different descriptor.
The neighbour shell uses all element pairs within --order-cutoff, falling
back to the analyser’s --cutoff when no order-specific cutoff is supplied.
It does not apply the chemical bonding filter used for coordination and angles.
Periodic images contribute their actual bond directions; cluster sizes count
unique atoms in the supplied cell. Two ordered atoms belong to the same
cluster when connected by a path of cutoff neighbours that are all ordered,
including connections across periodic boundaries. This is a connectivity
measure; it does not additionally require aligned \(q_{6m}\) vectors.
An atom is ordered when qbar6 >= qbar6_threshold and it has at least
order_min_neighbors neighbours. The defaults, 0.3 and 4, are heuristic.
They do not identify a crystal phase or give a universal crystalline volume
fraction. Calibrate the cutoff and threshold using crystalline and liquid
references at relevant temperatures, especially for GeTe. Separate partial-RDF
cutoffs can include different geometric shells for different element pairs;
inspect the resolved cutoffs and use an explicit scalar or pair cutoff when
needed. Compare distributions as well as ordered fractions and cluster sizes.
The example’s 3.5 Å order cutoff illustrates the first shell of an ideal
rocksalt GeTe cell with lattice constant 6 Å. Its six neighbours lie at
3 Å and give \(q_6=\bar q_6=\sqrt{1/8}\approx0.35355\); all atoms therefore
pass the 0.3 threshold. With pairwise auto-rdf cutoffs, the same ideal cell
can include twelve same-element neighbours as well, giving
\(\bar q_6\approx0.26517\) and no ordered atoms at that threshold. The explicit
--order-cutoff prevents changing the shell used for the other chemical
analyses. This ideal-cell example is not a calibration for thermally distorted
or rhombohedral GeTe; inspect short and long bonds and reference distributions
before choosing the shell for a phase-change workflow.
To reproduce the ideal-cell check without a calculator:
from ase.build import bulk
from amorphgen.analysis import compute_bond_order
ideal = bulk("GeTe", "rocksalt", a=6.0, cubic=True).repeat((2, 2, 2))
frame = compute_bond_order([ideal], cutoff=3.5)["per_structure"][0]
print(frame["qbar6_mean"]) # approximately 0.353553
print(frame["largest_cluster_size"]) # 64: the whole cell
The equivalent YAML settings are:
analysis:
bond_order: true
qbar6_threshold: 0.3
order_min_neighbors: 4
order_cutoff: 3.5 # illustrative ideal-rocksalt shell; calibrate for your system
cutoff: auto-rdf
--save-report includes the order summary; --save-plot adds
analysis_bond_order.json, analysis_bond_order.csv,
analysis_bond_order_atoms.csv and
analysis_bond_order.png (.pdf with --save-pdf). JSON retains the
per-atom values, labels and resolved cutoffs for reproducible comparisons.
from amorphgen.analysis import StructureAnalyser
from amorphgen.analysis.descriptors import save_descriptor
sa = StructureAnalyser("gete_mq/final/", cutoff="auto-rdf")
order = sa.bond_order(qbar6_threshold=0.3, min_neighbors=4, cutoff=3.5)
frame = order["per_structure"][0]
print(frame["ordered_fraction"], frame["largest_cluster_size"])
save_descriptor("bond_order", order, "gete_plots/")
ordered_fraction and largest_cluster_fraction use all atoms in each
structure as the denominator. The top-level result summarizes structures
with equal weight; inspect per_structure for individual clusters and
per-atom arrays. An isolated atom has \(q_6=\bar q_6=0\) and is disordered.
For a crystal-started --mq-ensemble run, a separate automatic
melt_memory report measures how much of the initially ordered atom
population is also ordered after heating and in the high-temperature
snapshots. See MQ-ensemble workflow for the definition and its limitations.
Void, oxygen, elastic and vibrational descriptors
These four descriptors are opt-in. Void sampling and oxygen speciation use
geometry alone and do not load an ML model. Elastic and vibrational analysis
load the selected calculator (--model, --model-path, --device) and
evaluate new configurations; saved single-point stresses or forces are
insufficient.
Void distribution
amorphgen --analyse --input-dir silica/ --voids \
--void-samples 20000 --void-probe-radius 0.5 --void-seed 42 \
--void-probe-radii 0 0.25 0.5 0.75 1.0 \
--save-plot descriptors/
Uniform random points sample point clearance: the distance to the nearest atomic-sphere surface, in Å. A point is accessible when its clearance is at least the probe radius. This measures local free space; it does not find connected pores, pore throats or maximal cavities. The reported radius is clearance, not diameter, and sampled maxima underestimate the true maximum.
Atomic spheres use ASE covalent radii by default. Choose a consistent radius
convention for comparisons; override individual elements with YAML
analysis.void_radii or the Python radii argument. The histogram density
integrates to one over accessible points, while bin_volume_fraction sums
to the accessible fraction. Ensemble fractions are weighted by cell volume;
accessible_volume is the mean accessible volume per structure. Standard
errors describe Monte Carlo sampling only. Per-structure 95% Wilson intervals
also cover cases where no accessible points were found; neither measure
captures variation between structures. All cells must be fully periodic in 3D.
--void-probe-radii evaluates an accessible-fraction and accessible-volume
curve from the same sampled points, including thresholds below the
histogram’s --void-probe-radius. Radii must be finite and nonnegative;
they are sorted and duplicates removed. Without this option the curve uses
the histogram bin edges. Fractions decrease as the probe grows. Thresholds
share samples, so their errors are correlated; the shaded plot shows one
Monte Carlo standard error at each threshold, not a simultaneous interval.
The result also includes clearance_quantiles (p10, p50, p90) for
points admitting the base probe. These are empirical inverse-CDF quantiles,
weighted by cell volume in the pooled ensemble. They are null when no
points are accessible. Each structure retains its own quantiles and probe
curve, including Wilson intervals. The nested uncertainty summaries
estimate equal-weight structure means separately from sampling errors.
In addition to the histogram and full JSON, --save-plot writes
analysis_voids_probe.csv and analysis_voids_probe.png for the curve,
plus analysis_voids_structures.csv for per-structure volumes, clearances
and quantiles. --save-pdf adds a PDF of each figure.
Bridging and non-bridging oxygen
amorphgen --analyse --input-dir aluminosilicate/ --oxygen-speciation \
--network-formers Si,Al --cutoff "Si-O=2.0,Al-O=2.3" \
--save-plot descriptors/
Each oxygen is classified by its number of neighbouring selected network
formers: zero = free, one = non_bridging, two = bridging, three =
tricluster, and four or more = higher_coordinated. Fractions pool oxygen
counts across structures. The selection defaults to the Al, B, Ge, P and Si
present. Select the appropriate formers explicitly for other oxides and
exclude modifiers such as Na or Ca. All ensemble structures must have the
same element set; analyse different chemistries separately.
This uses the analyser’s pair cutoffs and periodic neighbours. Check those
cutoffs before interpreting the counts. free means no selected former
neighbour; the descriptor does not infer charge, bond order or hydroxyl
speciation. The former/modifier distinction follows the connectivity
convention described by Stebbins and Xu.
Elastic moduli from calculator stresses
amorphgen --analyse --input-dir relaxed_silica/ --elastic \
--model mace-mpa-0 --device cpu --elastic-strain 0.005 \
--save-plot descriptors/ --save-report descriptors.txt
Central differences require 13 stress evaluations per structure. The default
keeps fractional atomic coordinates fixed under strain (clamped ions).
--elastic-relax instead optimises internal positions at each fixed cell,
including the reference; --fmax and --opt-steps control convergence.
Optimise the starting cell separately for equilibrium moduli. Input
structures remain unchanged and failed internal relaxation raises an error.
The symmetrized stiffness tensor uses engineering strain in ASE Voigt order
xx, yy, zz, yz, xz, xy; stiffness and moduli are in GPa. The report includes
Voigt, Reuss and Hill bulk, shear and Young’s moduli, dimensionless Poisson
ratios, residual stress and stability diagnostics. Reuss/Hill estimates are
unavailable for unstable or ill-conditioned tensors. Residual stress above
0.1 GPa is flagged: these are static tangent stress-strain coefficients with
no finite-pressure correction. Check convergence with strain amplitude and
calculator precision. The averaging equations follow the
NIST atomman reference.
Harmonic vibrational density of states
amorphgen --analyse --input-dir relaxed_silica/ --vdos \
--model mace-mpa-0 --device cpu --vdos-displacement 0.01 \
--vdos-sigma 0.1 --vdos-npoints 800 --save-plot descriptors/
The mass-weighted force-constant matrix is built using central finite
differences, as in ASE’s harmonic vibration formulation.
This requires 6N force evaluations and dense diagonalisation of a
3N × 3N matrix per structure. Begin with small cells to assess cost.
The spectrum contains all 3N modes of each supplied cell; for periodic
structures these are Gamma-point modes, without Brillouin-zone sampling.
Optimise the reference positions beforehand. No automatic optimisation or acoustic sum rule is applied, and atomic constraints are rejected. Negative plotted frequencies denote imaginary modes. The Gaussian width is in THz; the displacement is in Å. The total DOS integrates to one on its returned grid, with equal weight per mode across structures. Element projections use squared mass-weighted eigenvector components and sum to the total DOS; they are not scattering-weighted experimental intensities. Check displacement, broadening and grid convergence before interpreting fine features.
With --save-plot, each requested descriptor writes full JSON (including
per-structure results), a CSV and a PNG; add --save-pdf for PDF figures.
--save-report appends the compact text summaries. Python methods return
results without automatically writing files:
from amorphgen.analysis import StructureAnalyser
from amorphgen.analysis.descriptors import save_descriptor
sa = StructureAnalyser("silica/", cutoff={"Si-O": 2.0})
voids = sa.void_distribution(n_samples=20000, probe_radius=0.5, seed=42,
probe_radii=[0, 0.25, 0.5, 0.75, 1.0])
print(voids["clearance_quantiles"], voids["probe_curve"])
rings = sa.ring_statistics(bond_pair=("Si", "O"), max_ring=16)
save_descriptor("rings", rings, "descriptors/", save_pdf=True)
oxygen = sa.oxygen_speciation(network_formers=["Si"])
save_descriptor("voids", voids, "descriptors/", save_pdf=True)
# Optional model-backed descriptors; install the matching backend extra.
from amorphgen.utils import get_calculator
calc = get_calculator(model="mace-mpa-0", device="cpu")
elastic = sa.elastic_moduli(calculator=calc, strain=0.005, relax=False)
vdos = sa.vibrational_dos(calculator=calc, displacement=0.01,
sigma=0.1, npoints=800)
print(elastic["ensemble"]["moduli"]["hill"])
print(vdos["imaginary_modes"])
Both calculator-backed methods also accept live calculators attached to
the input ASE objects when calculator is omitted.
Full CLI flag reference
The brackets below denote optional arguments; omit the brackets when running the command.
amorphgen --analyse \
--input-dir DIR_OF_STRUCTURES \
[--cutoff MODE_OR_NUMBER] \
[--cutoff-window FLOAT] \
[--per-structure] \
[--screen] [--screening-output PREFIX] \
[--save-report FILE] \
[--save-plot DIR] \
[--save-pdf] \
[--reference YAML] \
[--smearing SIGMA] \
[--total-rdf] \
[--sq] [--sq-weighting {xray,neutron,unweighted}] \
[--sq-method {direct,ft}] [--sq-smooth SIGMA_Q] [--sq-partials] \
[--pair-panels] [--total-cn SPEC] \
[--tr] [--tr-qrange QMIN QMAX] [--tr-window {lorch,none}] [--tr-scan] \
[--check-dimers] [--rings [PAIR]] [--voronoi [ELEMENT]] [--connectivity] \
[--bond-order] [--order-cutoff MODE_OR_NUMBER] \
[--qbar6-threshold FLOAT] [--order-min-neighbors INT] \
[--dpi N] \
[--show-title]
Flag |
What it does |
|---|---|
|
Read structure files in this directory. Same-stem duplicates count once, preferring |
|
|
|
Finite positive half-window in Å for the default near-cutoff contact shares and five-point coordination sweep (default 0.1). CLI overrides |
|
Print a per-structure table (one row per file: density, E/atom, CN). |
|
Label candidates before analysis with the default screens, or apply |
|
Save screening JSON, per-structure CSV and summary CSV (default |
|
Write the full text report (densities, bond distances, coordination, angles) to a file. |
|
Save available standard figures (RDF, CN, angles, density) plus CSV data into |
|
Also save vector PDF copies alongside the PNGs. |
|
Report available descriptor names, uncertainty versus ensemble size, and declared tolerance status. |
|
Declare an absolute Student-t mean interval half-width in descriptor units; repeat for each descriptor. Enables convergence reporting. |
|
Confidence for convergence intervals and forecasts (default 0.95). |
|
Largest total ensemble size searched for the forecast (default 1000000). |
|
Validate against the literature ranges in YAML, print a match/concern/fail table. |
|
Gaussian smearing of the RDF in Å (default 0.05, roughly thermal broadening; 0 for the raw histogram). |
|
Overlay the total g(r) on the partial-RDF plot. |
|
Compute the total structure factor S(q) and save it as PNG + CSV under |
|
|
|
|
|
Gaussian re-binning width in Å⁻¹ for the direct S(q) (default 0.05; 0 = raw). Raw values are kept in the CSV. |
|
With |
|
Total correlation function T(r) = 4πrρg(r), obtained by transforming the direct S(q). Match the reference’s scattering weights, normalization, q range and window before comparing. Writes |
|
Integration limits for |
|
Window for the |
|
With |
|
One small panel per element pair for the partial g(r) ( |
|
Total first-shell coordination of one element over several partner types, repeatable: |
|
Report unphysical close contacts (O–O peroxide, N–N) per structure. |
|
Shortest ring per network edge, size summaries and search coverage; |
|
Maximum searched size in network nodes (default 12) and ring-specific bond cutoff in Å (default: analyser pair cutoff). |
|
Voronoi indices |
|
Corner/edge/face sharing between cation-centred polyhedra (two cations sharing one anion = corner, two = edge, three or more = face) and the percentage of cations in at least one edge- or face-sharing pair, which is near zero in a corner-sharing network glass and tens of percent in a random packing. Added to the report; |
|
Steinhardt \(q_6\), Lechner–Dellago \(\bar q_6\), ordered atom fraction and largest connected ordered cluster. |
|
Neighbour cutoff for bond order and MQ melt-memory reports; accepts the same scalar, pair and automatic forms as |
|
Minimum \(\bar q_6\) for classifying an atom as ordered (default 0.3; calibrate for the material and shell). Also applies to MQ melt-memory reports. |
|
Minimum neighbour count for an ordered atom (default 4). Also applies to MQ melt-memory reports. |
|
Periodic point-clearance distribution and accessible volume. |
|
Monte Carlo points per cell (default 10000) and histogram bins (50). |
|
Probe radius in Å (default 0) and sampling seed (0). |
|
Threshold radii in Å for the accessible-volume curve from the same samples; defaults to histogram bin edges. |
|
Oxygen connectivity classes; formers default to the Al/B/Ge/P/Si present. |
|
Stress-derived tensor and isotropic moduli; strain amplitude defaults to 0.005. |
|
Relax internal positions at each fixed cell, using |
|
Harmonic cell modes; displacement defaults to 0.01 Å. |
|
Gaussian width in THz (default 0.1) and frequency-grid points (400). |
|
PNG DPI (default 300). |
|
Add titles to each plot (default off, captions usually clearer in figures). |
Outputs explained
For each ensemble, --analyse --save-plot DIR writes the applicable files
below. PDF copies require --save-pdf; optional descriptors require the
listed flags.
File |
What’s in it |
|---|---|
|
Partial RDFs for all unique pairs in the system. One line per pair using the Okabe-Ito palette. |
|
Columns: |
|
With |
|
With |
|
Coordination distribution of the bonded pairs. Binary AB systems (SiO₂) as mirrored bars: A-B on top, B-A reflected below the zero line. Multi-cation compounds (IGZO) as one panel per cation-centred pair (Ga-O, In-O, Zn-O) plus the anion total over all its cations (O-(Ga+In+Zn)). Mono-element systems (a-Si) and alloys side-by-side. |
|
Per-pair CN counts as percentages of the centred atom population, plus the anion-total rows. |
|
Resolved cutoffs, near-cutoff contact counts/shares and the coordination sweep, including per-structure values. |
|
Per-pair endpoint contact counts and near-cutoff shares; directional coordination across the cutoff window. |
|
With |
|
Bond-angle histograms (normalised). One line per triplet. |
|
Raw angle values, one row per triplet observation. |
|
Per-structure density violin with jittered scatter and mean ± std label (at least two structures). |
|
One row per structure: |
|
With |
|
With |
|
With |
|
Ring-size summaries and search coverage per structure; per-structure size distributions. |
|
With |
|
With |
|
With |
|
With |
|
Probe-radius curve with sampling errors; per-structure volumes, clearances and quantiles. PDF requires |
|
Oxygen counts/fractions by class; JSON also contains each oxygen’s former coordination. |
|
Modulus means/std/counts and mean stiffness heatmap; JSON includes raw/symmetrized tensors and diagnostics per structure. |
|
Total and element-projected DOS; JSON includes individual mode frequencies and per-structure diagnostics. |
|
With |
|
One row per descriptor for decisions/counts/units; one row per descriptor, planned size and point for uncertainty curves. |
|
One independent figure per descriptor with tolerance, observed endpoint and estimated required size. PDF requires |
Python API
For programmatic access, useful in scripts, notebooks, and the comparison workflow:
from amorphgen.analysis import StructureAnalyser
sa = StructureAnalyser("hybrid_ga2o3/final/") # accepts a dir OR list of files
sa.summary() # prints and returns the structural summary
# Individual descriptors
rho = sa.density()
print(rho["mean"], rho["std"], rho["values"])
rdf = sa.rdf(pair="Ga-O", sigma=0.05)
print(rdf["r"], rdf["g_r"])
cn = sa.coordination()
print(cn["Ga-O"]["mean"], cn["Ga-O"]["distribution"])
ang = sa.bond_angles()
print(ang["O-Ga-O"]["mean"])
# Structure factor: direct method, neutron weighting, Faber-Ziman partials
sq = sa.structure_factor_direct(weighting="neutron", sigma_q=0.05, partials=True)
print(sq["q"], sq["s_q"], sq["partials"]["Ga-Ga"])
# Multi-ensemble comparison
from amorphgen.analysis import EnsembleSpec, compare_ensembles
# Call compare_ensembles with the arguments in "Compare multiple ensembles" above.
# Validate against reference
from amorphgen.analysis import validate_against_reference, format_validation_report
import yaml
with open("examples/reference_a_Ga2O3.yaml") as f:
ref = yaml.safe_load(f)
print(format_validation_report(validate_against_reference(sa, ref)))
See Structure analysis for the full Python API.
Troubleshooting
“Density is high but CN looks too low”
The default auto-rdf cutoff handles most systems, but if you’re using
the legacy --cutoff auto (minsep-based), it may truncate the first
RDF peak for materials with broad bond distributions. Symptom: many
atoms appear under-coordinated. Fix: drop the --cutoff auto and let
the default auto-rdf resolve it.
“Per-structure E/atom shows N/A”
If your structures are in VASP or CIF format, they don’t carry energy
in their headers. AmorphGen falls back to reading random_gen.log
in the parent directory (the file written by amorphgen --random-gen --relax). If you’ve moved the structures away from their original
--random-gen output, place the log in their parent directory, or use
--format xyz / --format extxyz (which embed energy in the
comment line) when generating.
“Si–Si appears as a bond in a-SiO₂”
That’s an analysis artifact. Same-element pairs in multi-element ionic
compounds (Si–Si in SiO₂, Hf–Hf in HfO₂, Ga–Ga in Ga₂O₃) are
second-shell contacts mediated through the anion, not first-shell
bonds. They are listed under Non-bonded contacts and excluded from
bonding coordination, total coordination, bond-angle triplets and the
CN plot. For SiO₂, the relevant bonding CNs are Si–O and O–Si. In
a-Si and a-Si:H, however, Si–Si is a host-network bond and is included.
“RDF goes to zero suddenly at large r”
Check the plotted r range and the cell size. The default rmax is half
the smallest cell-vector length, rounded down to 0.1 Å. Extending the
range can introduce correlations from periodic replicas; use a larger
cell to investigate longer-range structure.