Structure Analysis Tutorial¶

This notebook walks through every analysis tool available in amorphouspy using a 25 mol% Na₂O-SiO₂ sodium silicate glass as the working example.

The glass network is built from SiO₄ tetrahedra linked by bridging oxygens (BO). Sodium acts as a network modifier: it breaks Si-O-Si linkages, creating non-bridging oxygens (NBO) and lowering the network connectivity. All tools below probe different facets of this local and medium-range order.

Tools covered:

Tool	What it measures
`compute_rdf`	Pair radial distribution functions g(r) + cumulative CN(r)
`compute_coordination`	Coordination number histogram for a chosen element
`compute_qn`	Q⁰-Q⁴ distribution (bridging oxygens per former)
`compute_network_connectivity`	Weighted average Q^n
`compute_angles`	Bond angle distribution for a chosen triplet
`compute_structure_factor`	Total and partial S(q) — neutron or X-ray
`compute_guttmann_rings`	Guttman primitive ring-size distribution
`compute_cavities`	Void size, surface area, and shape statistics
`visualize_cavities`	Interactive 3-D cavity rendering
`analyze_structure`	One-call workflow that runs all of the above
`plot_analysis_results_plotly`	Interactive dashboard of all results

Input file: data/SiONa_25.extxyz — 25 mol% Na₂O, 75 mol% SiO₂, ~600 atoms, periodic cubic box.

0. Imports and structure loading¶

In [1]:

Copied!





import matplotlib.pyplot as plt
import numpy as np
from ase.io import read

from amorphouspy.analysis.bond_angle_distribution import compute_angles
from amorphouspy.analysis.cavities import compute_cavities, visualize_cavities
from amorphouspy.analysis.qn_network_connectivity import compute_network_connectivity, compute_qn
from amorphouspy.analysis.radial_distribution_functions import compute_coordination, compute_rdf
from amorphouspy.analysis.rings import compute_guttmann_rings, generate_bond_length_dict
from amorphouspy.analysis.structure_factor import compute_structure_factor
from amorphouspy.workflows.structural_analysis import analyze_structure, plot_analysis_results_plotly

import plotly.io as pio

pio.renderers.default = "notebook"
import matplotlib.pyplot as plt
import numpy as np
from ase.io import read

from amorphouspy.analysis.bond_angle_distribution import compute_angles
from amorphouspy.analysis.cavities import compute_cavities, visualize_cavities
from amorphouspy.analysis.qn_network_connectivity import compute_network_connectivity, compute_qn
from amorphouspy.analysis.radial_distribution_functions import compute_coordination, compute_rdf
from amorphouspy.analysis.rings import compute_guttmann_rings, generate_bond_length_dict
from amorphouspy.analysis.structure_factor import compute_structure_factor
from amorphouspy.workflows.structural_analysis import analyze_structure, plot_analysis_results_plotly

import plotly.io as pio

pio.renderers.default = "notebook"

In [2]:

Copied!





# Load the structure — ASE reads extxyz natively
atoms = read("data/SiONa_25.extxyz")

# Atomic numbers present in this system
# O=8, Si=14, Na=11
O_TYPE = 8
SI_TYPE = 14
NA_TYPE = 11
TYPE_MAP = {O_TYPE: "O", SI_TYPE: "Si", NA_TYPE: "Na"}

# Quick sanity check
cell = atoms.get_cell().array
volume_A3 = atoms.get_volume()
mass_amu = atoms.get_masses().sum()
density = (mass_amu / 6.022140857e23) / (volume_A3 * 1e-24)

print(f"Atoms     : {len(atoms)}")
print(f"Cell (Å)  : {np.diag(cell).round(3)}")
print(f"Volume    : {volume_A3:.1f} Å³")
print(f"Density   : {density:.3f} g/cm³")
print(f"Elements  : { {sym: int((atoms.numbers == z).sum()) for z, sym in TYPE_MAP.items()} }")
# Load the structure — ASE reads extxyz natively
atoms = read("data/SiONa_25.extxyz")

# Atomic numbers present in this system
# O=8, Si=14, Na=11
O_TYPE = 8
SI_TYPE = 14
NA_TYPE = 11
TYPE_MAP = {O_TYPE: "O", SI_TYPE: "Si", NA_TYPE: "Na"}

# Quick sanity check
cell = atoms.get_cell().array
volume_A3 = atoms.get_volume()
mass_amu = atoms.get_masses().sum()
density = (mass_amu / 6.022140857e23) / (volume_A3 * 1e-24)

print(f"Atoms     : {len(atoms)}")
print(f"Cell (Å)  : {np.diag(cell).round(3)}")
print(f"Volume    : {volume_A3:.1f} Å³")
print(f"Density   : {density:.3f} g/cm³")
print(f"Elements  : { {sym: int((atoms.numbers == z).sum()) for z, sym in TYPE_MAP.items()} }")

Atoms     : 10020
Cell (Å)  : [51.156 51.156 51.156]
Volume    : 133875.2 Å³
Density   : 2.509 g/cm³
Elements  : {'O': 5845, 'Si': 2505, 'Na': 1670}

1. Radial Distribution Functions¶

compute_rdf returns three objects:

r — bin centres in Å, shape (n_bins,)
rdfs — dict[(Z1, Z2) → ndarray] — normalised g(r) for every unique unordered pair found in the structure. By convention the key is (min_Z, max_Z).
cn_cumulative — dict[(Z1, Z2) → ndarray] — mean number of Z2 neighbours within radius r around a Z1 atom. The key is ordered: (Z1, Z2) gives CN of Z2 around Z1, and (Z2, Z1) gives the reverse.

Tip: pass type_pairs=[(8, 14), (8, 11)] to restrict to specific pairs and save memory.

In [3]:

Copied!

r, rdfs, cn = compute_rdf(atoms, r_max=10.0, n_bins=500)

# Keys in rdfs use (min_Z, max_Z) ordering
print("RDF pairs computed:", list(rdfs.keys()))

# Keys in cn_cumulative are ordered (Z_ref, Z_neighbour)
print("CN  pairs computed:", list(cn.keys()))
r, rdfs, cn = compute_rdf(atoms, r_max=10.0, n_bins=500)

# Keys in rdfs use (min_Z, max_Z) ordering
print("RDF pairs computed:", list(rdfs.keys()))

# Keys in cn_cumulative are ordered (Z_ref, Z_neighbour)
print("CN  pairs computed:", list(cn.keys()))

RDF pairs computed: [(8, 8), (8, 11), (8, 14), (11, 11), (11, 14), (14, 14)]
CN  pairs computed: [(8, 8), (8, 11), (11, 8), (8, 14), (14, 8), (11, 11), (11, 14), (14, 11), (14, 14)]

In [4]:

Copied!





fig, axes = plt.subplots(1, 3, figsize=(10.5, 3.5 / 1.333), dpi=150, layout="constrained")

pairs = [
    ((8, 8), (8, 8), "O-O", "C0"),
    ((8, 14), (14, 8), "Si-O", "C1"),
    ((8, 11), (11, 8), "Na-O", "C2"),
]

for ax, (rdf_key, cn_key, label, color) in zip(axes, pairs):
    ax.plot(r, rdfs[rdf_key], color=color, lw=1.5, label="g(r)")
    ax.plot(r, cn[cn_key], color=color, lw=1.5, ls="--", label="CN(r)")
    ax.axhline(1, color="grey", lw=0.6, ls=":")
    ax.set_xlim(0, 4)
    ax.set_ylim(0, 25)
    ax.set_xlabel(r"$r\,(\mathrm{\AA})$")
    ax.set_ylabel(r"$g(r),\;CN(r)$")
    ax.set_title(label)
    ax.legend(fontsize=8)

plt.show()

# Read off the Si-O first-shell coordination number at ~2.5 Å
r_cutoff = 2.5
cn_SiO_at_cutoff = float(cn[(SI_TYPE, O_TYPE)][np.searchsorted(r, r_cutoff)])
print(f"Si-O cumulative CN at {r_cutoff} Å: {cn_SiO_at_cutoff:.2f}")
fig, axes = plt.subplots(1, 3, figsize=(10.5, 3.5 / 1.333), dpi=150, layout="constrained")

pairs = [
    ((8, 8), (8, 8), "O-O", "C0"),
    ((8, 14), (14, 8), "Si-O", "C1"),
    ((8, 11), (11, 8), "Na-O", "C2"),
]

for ax, (rdf_key, cn_key, label, color) in zip(axes, pairs):
    ax.plot(r, rdfs[rdf_key], color=color, lw=1.5, label="g(r)")
    ax.plot(r, cn[cn_key], color=color, lw=1.5, ls="--", label="CN(r)")
    ax.axhline(1, color="grey", lw=0.6, ls=":")
    ax.set_xlim(0, 4)
    ax.set_ylim(0, 25)
    ax.set_xlabel(r"$r\,(\mathrm{\AA})$")
    ax.set_ylabel(r"$g(r),\;CN(r)$")
    ax.set_title(label)
    ax.legend(fontsize=8)

plt.show()

# Read off the Si-O first-shell coordination number at ~2.5 Å
r_cutoff = 2.5
cn_SiO_at_cutoff = float(cn[(SI_TYPE, O_TYPE)][np.searchsorted(r, r_cutoff)])
print(f"Si-O cumulative CN at {r_cutoff} Å: {cn_SiO_at_cutoff:.2f}")

No description has been provided for this image

Si-O cumulative CN at 2.5 Å: 4.00

2. Coordination Numbers¶

compute_coordination(structure, target_type, cutoff, neighbor_types) counts the number of neighbours of neighbor_types within cutoff of each atom of target_type.

It returns:

coord_dist — dict[int → int]: histogram mapping coordination number → atom count.
per_atom_cn — dict[atom_id → int]: per-atom coordination (useful for further analysis).

Good cutoff choices are read from the first minimum of the RDF (valley between first and second shells).

In [5]:

Copied!





# Cutoffs from the first RDF minima (valley after the first peak)
SI_O_CUTOFF = 2.0  # Å — Si-O first shell
NA_O_CUTOFF = 3.1  # Å — Na-O first shell
O_CUTOFF = 2.0  # Å — O neighbours among formers only

# Si coordination by O
si_coord_dist, si_per_atom = compute_coordination(atoms, SI_TYPE, SI_O_CUTOFF, [O_TYPE])

# Na coordination by O
na_coord_dist, na_per_atom = compute_coordination(atoms, NA_TYPE, NA_O_CUTOFF, [O_TYPE])

# O coordination by formers (Si only here)
o_coord_dist, o_per_atom = compute_coordination(atoms, O_TYPE, O_CUTOFF, [SI_TYPE])


def _print_dist(name, dist):
    total = sum(dist.values())
    for cn_val, count in sorted(dist.items()):
        print(f"  {name}_{cn_val}: {count:4d}  ({100 * count / total:.1f} %)")


print("Si coordination by O:")
_print_dist("Si", si_coord_dist)

print("\nNa coordination by O:")
_print_dist("Na", na_coord_dist)

print("\nO coordination by Si (former only):")
_print_dist("O", o_coord_dist)
# Cutoffs from the first RDF minima (valley after the first peak)
SI_O_CUTOFF = 2.0  # Å — Si-O first shell
NA_O_CUTOFF = 3.1  # Å — Na-O first shell
O_CUTOFF = 2.0  # Å — O neighbours among formers only

# Si coordination by O
si_coord_dist, si_per_atom = compute_coordination(atoms, SI_TYPE, SI_O_CUTOFF, [O_TYPE])

# Na coordination by O
na_coord_dist, na_per_atom = compute_coordination(atoms, NA_TYPE, NA_O_CUTOFF, [O_TYPE])

# O coordination by formers (Si only here)
o_coord_dist, o_per_atom = compute_coordination(atoms, O_TYPE, O_CUTOFF, [SI_TYPE])


def _print_dist(name, dist):
    total = sum(dist.values())
    for cn_val, count in sorted(dist.items()):
        print(f"  {name}_{cn_val}: {count:4d}  ({100 * count / total:.1f} %)")


print("Si coordination by O:")
_print_dist("Si", si_coord_dist)

print("\nNa coordination by O:")
_print_dist("Na", na_coord_dist)

print("\nO coordination by Si (former only):")
_print_dist("O", o_coord_dist)

Si coordination by O:
  Si_4: 2505  (100.0 %)

Na coordination by O:
  Na_2:    2  (0.1 %)
  Na_3:   38  (2.3 %)
  Na_4:  336  (20.1 %)
  Na_5:  620  (37.1 %)
  Na_6:  459  (27.5 %)
  Na_7:  174  (10.4 %)
  Na_8:   38  (2.3 %)
  Na_9:    2  (0.1 %)
  Na_10:    1  (0.1 %)

O coordination by Si (former only):
  O_1: 1670  (28.6 %)
  O_2: 4175  (71.4 %)

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

In [6]:

Copied!





fig, axes = plt.subplots(1, 3, figsize=(10.5, 3.5 / 1.333), dpi=150, layout="constrained")

datasets = [
    (si_coord_dist, "Si-O coordination", "C1"),
    (na_coord_dist, "Na-O coordination", "C2"),
    (o_coord_dist, "O-Si coordination", "C0"),
]

for ax, (dist, title, color) in zip(axes, datasets):
    total = sum(dist.values())
    cn_vals = sorted(dist.keys())
    pct = [100 * dist[n] / total for n in cn_vals]
    ax.bar(cn_vals, pct, color=color, edgecolor="black", linewidth=0.8)
    ax.set_xlabel("Coordination number")
    ax.set_ylabel("Fraction (%)")
    ax.set_title(title)
    ax.set_xticks(cn_vals)

plt.show()
fig, axes = plt.subplots(1, 3, figsize=(10.5, 3.5 / 1.333), dpi=150, layout="constrained")

datasets = [
    (si_coord_dist, "Si-O coordination", "C1"),
    (na_coord_dist, "Na-O coordination", "C2"),
    (o_coord_dist, "O-Si coordination", "C0"),
]

for ax, (dist, title, color) in zip(axes, datasets):
    total = sum(dist.values())
    cn_vals = sorted(dist.keys())
    pct = [100 * dist[n] / total for n in cn_vals]
    ax.bar(cn_vals, pct, color=color, edgecolor="black", linewidth=0.8)
    ax.set_xlabel("Coordination number")
    ax.set_ylabel("Fraction (%)")
    ax.set_title(title)
    ax.set_xticks(cn_vals)

plt.show()

3. $Q^n$ Distribution and Network Connectivity¶

The Q^n notation describes how many of the oxygen neighbours of a network-forming tetrahedron are bridging (i.e., shared with another former). In pure SiO₂ all Si are Q⁴; adding Na₂O introduces Q³ and Q² species.

compute_qn(structure, cutoff, former_types, o_type) takes:

cutoff — the Si-O bond cutoff used to identify both former-O bonds and bridging oxygens.
former_types — list of atomic numbers treated as network formers.
o_type — atomic number of oxygen.

Returns (total_qn, partial_qn):

total_qn — dict[int → int] summed over all formers.
partial_qn — dict[Z_former → dict[int → int]] broken down by former type.

compute_network_connectivity(qn_dist) returns the weighted average ⟨n⟩.

In [7]:

Copied!





# cutoff is the Si-O bond length cutoff (same as used for coordination)
qn_total, qn_partial = compute_qn(
    atoms,
    cutoff=SI_O_CUTOFF,
    former_types=[SI_TYPE],
    o_type=O_TYPE,
)

nc = compute_network_connectivity(qn_total)

print("Q^n distribution (total):")
total_si = sum(qn_total.values())
for n, count in sorted(qn_total.items()):
    if count > 0:
        print(f"  Q{n}: {count:4d}  ({100 * count / total_si:.1f} %)")

print(f"\nNetwork connectivity ⟨n⟩ = {nc:.4f}")
print("(Pure SiO₂ → 4.0, maximum depolymerisation → 0.0)")

# Partial Qn per former type (Si=14)
print("\nPartial Q^n for Si (Z=14):")
si_qn = qn_partial[SI_TYPE]
for n, count in sorted(si_qn.items()):
    if count > 0:
        print(f"  Q{n}: {count}")
# cutoff is the Si-O bond length cutoff (same as used for coordination)
qn_total, qn_partial = compute_qn(
    atoms,
    cutoff=SI_O_CUTOFF,
    former_types=[SI_TYPE],
    o_type=O_TYPE,
)

nc = compute_network_connectivity(qn_total)

print("Q^n distribution (total):")
total_si = sum(qn_total.values())
for n, count in sorted(qn_total.items()):
    if count > 0:
        print(f"  Q{n}: {count:4d}  ({100 * count / total_si:.1f} %)")

print(f"\nNetwork connectivity ⟨n⟩ = {nc:.4f}")
print("(Pure SiO₂ → 4.0, maximum depolymerisation → 0.0)")

# Partial Qn per former type (Si=14)
print("\nPartial Q^n for Si (Z=14):")
si_qn = qn_partial[SI_TYPE]
for n, count in sorted(si_qn.items()):
    if count > 0:
        print(f"  Q{n}: {count}")

Q^n distribution (total):
  Q1:   29  (1.2 %)
  Q2:  267  (10.7 %)
  Q3: 1049  (41.9 %)
  Q4: 1160  (46.3 %)

Network connectivity ⟨n⟩ = 3.3333
(Pure SiO₂ → 4.0, maximum depolymerisation → 0.0)

Partial Q^n for Si (Z=14):
  Q1: 29
  Q2: 267
  Q3: 1049
  Q4: 1160

In [8]:

Copied!





fig, ax = plt.subplots(figsize=(4, 3), dpi=150, layout="constrained")

ns = sorted(qn_total.keys())
counts = [qn_total[n] for n in ns]
total = sum(counts)
pct = [100 * c / total for c in counts]

bars = ax.bar(ns, pct, color="C1", edgecolor="black", linewidth=0.8)
ax.set_xlabel(r"$n$ in $Q^n$")
ax.set_ylabel("Fraction (%)")
ax.set_title(rf"$Q^n$ distribution — $\langle n \rangle = {nc:.2f}$")
ax.set_xticks(ns)

# Annotate non-zero bars
for bar, p in zip(bars, pct):
    if p > 0.5:
        ax.text(bar.get_x() + bar.get_width() / 2, bar.get_height(), f"{p:.1f}%", ha="center", va="bottom", fontsize=8)
plt.show()
fig, ax = plt.subplots(figsize=(4, 3), dpi=150, layout="constrained")

ns = sorted(qn_total.keys())
counts = [qn_total[n] for n in ns]
total = sum(counts)
pct = [100 * c / total for c in counts]

bars = ax.bar(ns, pct, color="C1", edgecolor="black", linewidth=0.8)
ax.set_xlabel(r"$n$ in $Q^n$")
ax.set_ylabel("Fraction (%)")
ax.set_title(rf"$Q^n$ distribution — $\langle n \rangle = {nc:.2f}$")
ax.set_xticks(ns)

# Annotate non-zero bars
for bar, p in zip(bars, pct):
    if p > 0.5:
        ax.text(bar.get_x() + bar.get_width() / 2, bar.get_height(), f"{p:.1f}%", ha="center", va="bottom", fontsize=8)
plt.show()

4. Bond Angle Distribution¶

compute_angles(structure, center_type, neighbor_type, cutoff, bins=180) computes the distribution of neighbor-center-neighbor angles.

center_type — scalar atomic number of the central atom.
neighbor_type — scalar atomic number of the two flanking atoms.
cutoff — bond cutoff in Å.

Returns (bin_centers, angle_hist) where angle_hist is normalised (area ≈ 1 over 0-180°).

The two signature angles for sodium silicate are:

O-Si-O ≈ 109.5° (intra-tetrahedron, close to ideal tetrahedral)
Si-O-Si ≈ 145-150° (inter-tetrahedron bridging angle)

In [9]:

Copied!





# O-Si-O: center = Si, neighbors = O
angles_OSiO, hist_OSiO = compute_angles(
    atoms,
    center_type=SI_TYPE,
    neighbor_type=O_TYPE,
    cutoff=SI_O_CUTOFF,
    bins=180,
)

# Si-O-Si: center = O, neighbors = Si
angles_SiOSi, hist_SiOSi = compute_angles(
    atoms,
    center_type=O_TYPE,
    neighbor_type=SI_TYPE,
    cutoff=SI_O_CUTOFF,
    bins=180,
)

fig, axes = plt.subplots(1, 2, figsize=(7, 3.5 / 1.333), dpi=150, layout="constrained")

axes[0].plot(angles_OSiO, hist_OSiO, color="C0", lw=1.5)
axes[0].axvline(109.47, color="grey", lw=0.8, ls="--", label="ideal tetrahedral")
axes[0].set_xlim(60, 160)
axes[0].set_ylim(0, None)
axes[0].set_xlabel(r"Angle (°)")
axes[0].set_ylabel("Frequency (a.u.)")
axes[0].set_title("O-Si-O")
axes[0].legend(fontsize=8)

axes[1].plot(angles_SiOSi, hist_SiOSi, color="C1", lw=1.5)
axes[1].set_xlim(80, 180)
axes[1].set_ylim(0, None)
axes[1].set_xlabel(r"Angle (°)")
axes[1].set_ylabel("Frequency (a.u.)")
axes[1].set_title("Si-O-Si")

plt.show()

# Peak positions
peak_OSiO = float(angles_OSiO[np.argmax(hist_OSiO)])
peak_SiOSi = float(angles_SiOSi[np.argmax(hist_SiOSi)])
print(f"O-Si-O peak: {peak_OSiO:.1f}°  (ideal tetrahedral: 109.5°)")
print(f"Si-O-Si peak: {peak_SiOSi:.1f}°")
# O-Si-O: center = Si, neighbors = O
angles_OSiO, hist_OSiO = compute_angles(
    atoms,
    center_type=SI_TYPE,
    neighbor_type=O_TYPE,
    cutoff=SI_O_CUTOFF,
    bins=180,
)

# Si-O-Si: center = O, neighbors = Si
angles_SiOSi, hist_SiOSi = compute_angles(
    atoms,
    center_type=O_TYPE,
    neighbor_type=SI_TYPE,
    cutoff=SI_O_CUTOFF,
    bins=180,
)

fig, axes = plt.subplots(1, 2, figsize=(7, 3.5 / 1.333), dpi=150, layout="constrained")

axes[0].plot(angles_OSiO, hist_OSiO, color="C0", lw=1.5)
axes[0].axvline(109.47, color="grey", lw=0.8, ls="--", label="ideal tetrahedral")
axes[0].set_xlim(60, 160)
axes[0].set_ylim(0, None)
axes[0].set_xlabel(r"Angle (°)")
axes[0].set_ylabel("Frequency (a.u.)")
axes[0].set_title("O-Si-O")
axes[0].legend(fontsize=8)

axes[1].plot(angles_SiOSi, hist_SiOSi, color="C1", lw=1.5)
axes[1].set_xlim(80, 180)
axes[1].set_ylim(0, None)
axes[1].set_xlabel(r"Angle (°)")
axes[1].set_ylabel("Frequency (a.u.)")
axes[1].set_title("Si-O-Si")

plt.show()

# Peak positions
peak_OSiO = float(angles_OSiO[np.argmax(hist_OSiO)])
peak_SiOSi = float(angles_SiOSi[np.argmax(hist_SiOSi)])
print(f"O-Si-O peak: {peak_OSiO:.1f}°  (ideal tetrahedral: 109.5°)")
print(f"Si-O-Si peak: {peak_SiOSi:.1f}°")

O-Si-O peak: 109.5°  (ideal tetrahedral: 109.5°)
Si-O-Si peak: 150.5°

5. Structure Factor S(q)¶

compute_structure_factor computes the 1-D isotropic structure factor via the Faber-Ziman sine transform of the partial RDFs:

$$S_{\alpha\beta}(q) = 1 + \frac{4\pi\rho}{q} \int_0^{r_\text{max}} r\bigl[g_{\alpha\beta}(r)-1\bigr] M(r)\sin(qr)\,\mathrm{d}r$$

where $M(r) = \operatorname{sinc}(r/r_\text{max})$ is the Lorch modification function.

Key parameters:

radiation="neutron" or "xray" — selects weighting factors (coherent scattering lengths vs. form factors).
lorch_damping=True — suppresses termination ripples from finite r_max.
q_min, q_max, n_q — momentum-transfer grid.

Returns (q, sq_total, sq_partials) where sq_partials is a dict[(Z1, Z2) → ndarray].

In [10]:

Copied!





# --- Neutron S(q) ---
q_n, sq_n_total, sq_n_partials = compute_structure_factor(
    atoms,
    q_min=0.1,
    q_max=20.0,
    n_q=600,
    r_max=10.0,
    n_bins=2000,
    radiation="neutron",
    lorch_damping=True,
)

# --- X-ray S(q) ---
q_x, sq_x_total, sq_x_partials = compute_structure_factor(
    atoms,
    q_min=0.1,
    q_max=20.0,
    n_q=600,
    r_max=10.0,
    n_bins=2000,
    radiation="xray",
    lorch_damping=True,
)

print("Partial pairs:", list(sq_n_partials.keys()))
# --- Neutron S(q) ---
q_n, sq_n_total, sq_n_partials = compute_structure_factor(
    atoms,
    q_min=0.1,
    q_max=20.0,
    n_q=600,
    r_max=10.0,
    n_bins=2000,
    radiation="neutron",
    lorch_damping=True,
)

# --- X-ray S(q) ---
q_x, sq_x_total, sq_x_partials = compute_structure_factor(
    atoms,
    q_min=0.1,
    q_max=20.0,
    n_q=600,
    r_max=10.0,
    n_bins=2000,
    radiation="xray",
    lorch_damping=True,
)

print("Partial pairs:", list(sq_n_partials.keys()))

Partial pairs: [(8, 8), (8, 11), (8, 14), (11, 11), (11, 14), (14, 14)]

In [11]:

Copied!





fig, axes = plt.subplots(2, 2, figsize=(9, 6), dpi=150, layout="constrained")

# --- Total S(q): neutron vs X-ray ---
for ax, (q, sq, label) in zip(
    axes[0],
    [
        (q_n, sq_n_total, "Neutron"),
        (q_x, sq_x_total, "X-ray"),
    ],
):
    ax.plot(q, sq, color="black", lw=1.0)
    ax.axhline(1, color="grey", lw=0.6, ls="--")
    ax.set_xlim(0.5, 20)
    ax.set_xlabel(r"$q\,(\mathrm{\AA}^{-1})$")
    ax.set_ylabel(r"$S(q)$")
    ax.set_title(f"Total {label} S(q)")


# --- Partial S(q): neutron ---
partial_styles = {(8, 8): ("C0", "O-O"), (8, 14): ("C1", "Si-O"), (8, 11): ("C2", "Na-O")}
for pair, (color, label) in partial_styles.items():
    if pair in sq_n_partials:
        axes[1][0].plot(q_n, sq_n_partials[pair], color=color, lw=1.2, label=label)
axes[1][0].axhline(1, color="grey", lw=0.6, ls="--")
axes[1][0].set_xlim(0.5, 20)
axes[1][0].set_xlabel(r"$q\,(\mathrm{\AA}^{-1})$")
axes[1][0].set_ylabel(r"$S_{\alpha\beta}(q)$")
axes[1][0].set_title("Partial neutron S(q)")
axes[1][0].legend(fontsize=8, frameon=False)

# --- Partial S(q): X-ray ---
for pair, (color, label) in partial_styles.items():
    if pair in sq_x_partials:
        axes[1][1].plot(q_x, sq_x_partials[pair], color=color, lw=1.2, label=label)
axes[1][1].axhline(1, color="grey", lw=0.6, ls="--")
axes[1][1].set_xlim(0.5, 20)
axes[1][1].set_xlabel(r"$q\,(\mathrm{\AA}^{-1})$")
axes[1][1].set_ylabel(r"$S_{\alpha\beta}(q)$")
axes[1][1].set_title("Partial X-ray S(q)")
axes[1][1].legend(fontsize=8, frameon=False)

plt.show()
fig, axes = plt.subplots(2, 2, figsize=(9, 6), dpi=150, layout="constrained")

# --- Total S(q): neutron vs X-ray ---
for ax, (q, sq, label) in zip(
    axes[0],
    [
        (q_n, sq_n_total, "Neutron"),
        (q_x, sq_x_total, "X-ray"),
    ],
):
    ax.plot(q, sq, color="black", lw=1.0)
    ax.axhline(1, color="grey", lw=0.6, ls="--")
    ax.set_xlim(0.5, 20)
    ax.set_xlabel(r"$q\,(\mathrm{\AA}^{-1})$")
    ax.set_ylabel(r"$S(q)$")
    ax.set_title(f"Total {label} S(q)")


# --- Partial S(q): neutron ---
partial_styles = {(8, 8): ("C0", "O-O"), (8, 14): ("C1", "Si-O"), (8, 11): ("C2", "Na-O")}
for pair, (color, label) in partial_styles.items():
    if pair in sq_n_partials:
        axes[1][0].plot(q_n, sq_n_partials[pair], color=color, lw=1.2, label=label)
axes[1][0].axhline(1, color="grey", lw=0.6, ls="--")
axes[1][0].set_xlim(0.5, 20)
axes[1][0].set_xlabel(r"$q\,(\mathrm{\AA}^{-1})$")
axes[1][0].set_ylabel(r"$S_{\alpha\beta}(q)$")
axes[1][0].set_title("Partial neutron S(q)")
axes[1][0].legend(fontsize=8, frameon=False)

# --- Partial S(q): X-ray ---
for pair, (color, label) in partial_styles.items():
    if pair in sq_x_partials:
        axes[1][1].plot(q_x, sq_x_partials[pair], color=color, lw=1.2, label=label)
axes[1][1].axhline(1, color="grey", lw=0.6, ls="--")
axes[1][1].set_xlim(0.5, 20)
axes[1][1].set_xlabel(r"$q\,(\mathrm{\AA}^{-1})$")
axes[1][1].set_ylabel(r"$S_{\alpha\beta}(q)$")
axes[1][1].set_title("Partial X-ray S(q)")
axes[1][1].legend(fontsize=8, frameon=False)

plt.show()

6. Ring Statistics (Guttman Algorithm)¶

The Guttman primitive ring criterion finds the shortest closed path through the network graph that passes through each node at most once. In silica networks, rings are counted in terms of the number of Si atoms they contain.

generate_bond_length_dict(atoms, specific_cutoffs) builds the per-pair cutoff dict expected by the ring finder. Use specific_cutoffs to override individual pairs; all other pairs get a conservative default.

compute_guttmann_rings(structure, bond_lengths, max_size, n_cpus) returns:

rings_dist — dict[int → int]: number of rings of each size.
mean_ring_size — weighted average ring size.

The ring search is the most expensive step. Use n_cpus > 1 for large structures.

In [12]:

Copied!





# Build the bond-length dict — only the Si-O cutoff is critical for silicate rings
bond_lengths = generate_bond_length_dict(
    atoms,
    specific_cutoffs={("Si", "O"): SI_O_CUTOFF},
)
print("Bond length cutoffs used:")
for pair, cutoff_val in sorted(bond_lengths.items()):
    print(f"  {pair}: {cutoff_val:.2f} Å")
# Build the bond-length dict — only the Si-O cutoff is critical for silicate rings
bond_lengths = generate_bond_length_dict(
    atoms,
    specific_cutoffs={("Si", "O"): SI_O_CUTOFF},
)
print("Bond length cutoffs used:")
for pair, cutoff_val in sorted(bond_lengths.items()):
    print(f"  {pair}: {cutoff_val:.2f} Å")

Bond length cutoffs used:
  ('Na', 'Na'): -1.00 Å
  ('Na', 'Si'): -1.00 Å
  ('O', 'Na'): -1.00 Å
  ('O', 'O'): -1.00 Å
  ('O', 'Si'): 2.00 Å
  ('Si', 'Si'): -1.00 Å

In [13]:

Copied!





rings_dist, mean_ring_size = compute_guttmann_rings(
    atoms,
    bond_lengths=bond_lengths,
    max_size=40,
    n_cpus=1,  # increase for large cells
)

sizes = np.array(sorted(rings_dist.keys()))
counts = np.array([rings_dist[s] for s in sizes])
total_rings = counts.sum()

print(f"Total rings found: {total_rings}")
print(f"Mean ring size   : {mean_ring_size:.2f} Si atoms")
print("\nRing size distribution:")
for s, c in zip(sizes, counts):
    print(f"  n={s:2d}: {c:4d}  ({100 * c / total_rings:.1f} %)")
rings_dist, mean_ring_size = compute_guttmann_rings(
    atoms,
    bond_lengths=bond_lengths,
    max_size=40,
    n_cpus=1,  # increase for large cells
)

sizes = np.array(sorted(rings_dist.keys()))
counts = np.array([rings_dist[s] for s in sizes])
total_rings = counts.sum()

print(f"Total rings found: {total_rings}")
print(f"Mean ring size   : {mean_ring_size:.2f} Si atoms")
print("\nRing size distribution:")
for s, c in zip(sizes, counts):
    print(f"  n={s:2d}: {c:4d}  ({100 * c / total_rings:.1f} %)")

Finding rings: 100%|██████████| 4175/4175 [00:05<00:00, 767.04edge/s]

Total rings found: 1869
Mean ring size   : 6.14 Si atoms

Ring size distribution:
  n= 3:    3  (0.2 %)
  n= 4:  215  (11.5 %)
  n= 5:  620  (33.2 %)
  n= 6:  482  (25.8 %)
  n= 7:  244  (13.1 %)
  n= 8:  114  (6.1 %)
  n= 9:   79  (4.2 %)
  n=10:   37  (2.0 %)
  n=11:   35  (1.9 %)
  n=12:   15  (0.8 %)
  n=13:    7  (0.4 %)
  n=14:    9  (0.5 %)
  n=15:    9  (0.5 %)

In [14]:

Copied!





fig, ax = plt.subplots(figsize=(5, 3.5), dpi=150, layout="constrained")

ax.bar(
    sizes,
    counts / total_rings,
    width=0.6,
    color="C1",
    edgecolor="black",
    linewidth=0.8,
    align="center",
)
ax.axvline(mean_ring_size, color="C3", lw=1.2, ls="--", label=f"mean = {mean_ring_size:.1f}")
ax.set_xlim(sizes.min() - 0.5, sizes.max() + 0.5)
ax.set_xticks(sizes)
ax.set_xlabel("Si atoms per ring")
ax.set_ylabel("Normalised count")
ax.set_title("Guttman ring-size distribution")
ax.legend(fontsize=9)

plt.show()
fig, ax = plt.subplots(figsize=(5, 3.5), dpi=150, layout="constrained")

ax.bar(
    sizes,
    counts / total_rings,
    width=0.6,
    color="C1",
    edgecolor="black",
    linewidth=0.8,
    align="center",
)
ax.axvline(mean_ring_size, color="C3", lw=1.2, ls="--", label=f"mean = {mean_ring_size:.1f}")
ax.set_xlim(sizes.min() - 0.5, sizes.max() + 0.5)
ax.set_xticks(sizes)
ax.set_xlabel("Si atoms per ring")
ax.set_ylabel("Normalised count")
ax.set_title("Guttman ring-size distribution")
ax.legend(fontsize=9)

plt.show()

7. Cavity Analysis¶

compute_cavities(structure, resolution, cutoff_radii) uses a voxel-based algorithm to detect and characterise empty voids in the glass network.

Each cavity is described by five scalar properties:

volumes — cavity volume (Å³)
surface_areas — cavity surface area (Å²)
asphericities — η ∈ [0, 1]: deviation from a sphere (0 = perfect sphere)
acylindricities — c ∈ [0, 1]: deviation from a cylinder
anisotropies — κ ∈ [0, 1]: relative anisotropy

resolution sets the voxel grid size (64 = 64³ grid); higher values resolve smaller voids but cost more memory.

visualize_cavities(structure, ...) returns an interactive Plotly figure with isosurface-rendered cavities.

In [15]:

Copied!





cavity_results = compute_cavities(structure=atoms, resolution=64)

print("Cavities detected:", len(cavity_results["volumes"]))
print(f"Total void volume : {sum(cavity_results['volumes']):.1f} Å³")
print(f"Largest cavity    : {max(cavity_results['volumes']):.2f} Å³")
print(f"Mean asphericity  : {np.mean(cavity_results['asphericities']):.3f}")
cavity_results = compute_cavities(structure=atoms, resolution=64)

print("Cavities detected:", len(cavity_results["volumes"]))
print(f"Total void volume : {sum(cavity_results['volumes']):.1f} Å³")
print(f"Largest cavity    : {max(cavity_results['volumes']):.2f} Å³")
print(f"Mean asphericity  : {np.mean(cavity_results['asphericities']):.3f}")

/var/folders/6p/qdp8kbgd4m98vyhvp2nptys00000gn/T/ipykernel_42373/4216481361.py:1: UserWarning:

6 percolating cavity/cavities span the simulation cell and will be excluded from the output. Consider using larger cutoff radii.

Cavities detected: 286
Total void volume : 1141.9 Å³
Largest cavity    : 39.83 Å³
Mean asphericity  : 0.428

In [16]:

Copied!





attrs = ["volumes", "surface_areas", "asphericities", "acylindricities", "anisotropies"]
labels = ["Volume (Å³)", "Surface area (Å²)", "Asphericity η", "Acylindricity c", "Anisotropy κ"]

fig, axes = plt.subplots(2, 3, figsize=(10.5, 5), dpi=150, layout="constrained")

for ax, attr, label in zip(axes.flat, attrs, labels):
    data = cavity_results[attr]
    ax.hist(data, bins=30, density=False, edgecolor="black", linewidth=0.6, color="steelblue")
    ax.set_xlabel(label)
    ax.set_ylabel("Count")
    ax.set_title(f"Cavity {label.split('(')[0].strip()}")

# Hide the unused 6th subplot
axes.flat[-1].set_visible(False)

plt.show()
attrs = ["volumes", "surface_areas", "asphericities", "acylindricities", "anisotropies"]
labels = ["Volume (Å³)", "Surface area (Å²)", "Asphericity η", "Acylindricity c", "Anisotropy κ"]

fig, axes = plt.subplots(2, 3, figsize=(10.5, 5), dpi=150, layout="constrained")

for ax, attr, label in zip(axes.flat, attrs, labels):
    data = cavity_results[attr]
    ax.hist(data, bins=30, density=False, edgecolor="black", linewidth=0.6, color="steelblue")
    ax.set_xlabel(label)
    ax.set_ylabel("Count")
    ax.set_title(f"Cavity {label.split('(')[0].strip()}")

# Hide the unused 6th subplot
axes.flat[-1].set_visible(False)

plt.show()

7.1 Interactive 3-D cavity visualisation¶

visualize_cavities renders the cavities as isosurfaces on top of the atom positions. Use excluded_cavities to filter out voids smaller than a given volume threshold (Å³) to avoid visual clutter from tiny interstitial spaces.

In [17]:

Copied!





fig_cavities = visualize_cavities(
    atoms,
    resolution=64,
    show_atoms=False,  # set True to overlay atom positions
    opacity=0.5,
    excluded_cavities=300,  # ignore voids smaller than 300 Å³
)
fig_cavities.show()
fig_cavities = visualize_cavities(
    atoms,
    resolution=64,
    show_atoms=False,  # set True to overlay atom positions
    opacity=0.5,
    excluded_cavities=300,  # ignore voids smaller than 300 Å³
)
fig_cavities.show()

8. Full Workflow: `analyze_structure` + `plot_analysis_results_plotly`¶

analyze_structure(atoms) runs all of the above in a single call. It auto-detects element roles (formers, modifiers, oxygen) and selects cutoffs from the RDF minima. The result is a typed StructureData Pydantic model with logical sub-objects:

StructureData
├── density          float (g/cm³)
├── coordination     CoordinationData
│   ├── oxygen       dict[str, int]
│   ├── formers      dict[str, dict[str, int]]
│   └── modifiers    dict[str, dict[str, int]]
├── network          NetworkData
│   ├── Qn_distribution         dict[str, float]
│   ├── Qn_distribution_partial dict[str, dict[str, float]]
│   └── connectivity            float
├── distributions    StructuralDistributions
│   ├── bond_angles  dict[str, (list[float], list[float])]
│   └── rings        dict  (keys: "distribution", "mean_size")
├── rdfs             RadialDistributionData
│   ├── r                     list[float]
│   ├── rdfs                  dict[str, list[float]]
│   └── cumulative_coordination dict[str, list[float]]
└── elements         ElementInfo
    ├── formers   list[str]
    ├── modifiers list[str]
    └── cutoffs   dict[str, float]

In [18]:

Copied!





results = analyze_structure(atoms)

# --- Summary printout ---
print(f"Density           : {results.density:.3f} g/cm³")
print(f"Network formers   : {results.elements.formers}")
print(f"Modifiers         : {results.elements.modifiers}")
print(f"Auto cutoffs (Å)  : {results.elements.cutoffs}")
print(f"Network connectivity ⟨n⟩: {results.network.connectivity:.4f}")
print("\nQ^n distribution:")
for n, frac in sorted(results.network.Qn_distribution.items(), key=lambda x: int(x[0])):
    total_q = sum(results.network.Qn_distribution.values())
    print(f"  Q{n}: {100 * frac / total_q:.1f} %")

print(f"\nRing mean size: {results.distributions.rings.get('mean_size', 'n/a')}")
print(f"O coordination : {results.coordination.oxygen}")
print(f"Si coordination: {results.coordination.formers}")
results = analyze_structure(atoms)

# --- Summary printout ---
print(f"Density           : {results.density:.3f} g/cm³")
print(f"Network formers   : {results.elements.formers}")
print(f"Modifiers         : {results.elements.modifiers}")
print(f"Auto cutoffs (Å)  : {results.elements.cutoffs}")
print(f"Network connectivity ⟨n⟩: {results.network.connectivity:.4f}")
print("\nQ^n distribution:")
for n, frac in sorted(results.network.Qn_distribution.items(), key=lambda x: int(x[0])):
    total_q = sum(results.network.Qn_distribution.values())
    print(f"  Q{n}: {100 * frac / total_q:.1f} %")

print(f"\nRing mean size: {results.distributions.rings.get('mean_size', 'n/a')}")
print(f"O coordination : {results.coordination.oxygen}")
print(f"Si coordination: {results.coordination.formers}")

Finding rings: 100%|██████████| 4175/4175 [00:05<00:00, 782.41edge/s]

Density           : 2.509 g/cm³
Network formers   : ['Si']
Modifiers         : ['Na']
Auto cutoffs (Å)  : {'O': 1.85, 'Na': 3.3100000000000005, 'Si': 1.85}
Network connectivity ⟨n⟩: 3.3333

Q^n distribution:
  Q0: 0.0 %
  Q1: 1.2 %
  Q2: 10.7 %
  Q3: 41.9 %
  Q4: 46.3 %
  Q5: 0.0 %
  Q6: 0.0 %

Ring mean size: 6.141787051899412
O coordination : {'1': 1670, '2': 4175}
Si coordination: {'Si': {'4': 2505}}

In [19]:

Copied!

# Interactive 3×3 dashboard — coordination, Q^n, bond angles, ring stats, RDFs
fig_dashboard = plot_analysis_results_plotly(results)
fig_dashboard.show()
# Interactive 3×3 dashboard — coordination, Q^n, bond angles, ring stats, RDFs
fig_dashboard = plot_analysis_results_plotly(results)
fig_dashboard.show()

9. Quick Reference¶

Function	Module	Signature summary
`compute_rdf`	`analysis.radial_distribution_functions`	`(atoms, r_max, n_bins, type_pairs) → (r, rdfs, cn)`
`compute_coordination`	`analysis.radial_distribution_functions`	`(atoms, target_type, cutoff, neighbor_types) → (dist, per_atom)`
`compute_qn`	`analysis.qn_network_connectivity`	`(atoms, cutoff, former_types, o_type) → (total, partial)`
`compute_network_connectivity`	`analysis.qn_network_connectivity`	`(qn_dist) → float`
`compute_angles`	`analysis.bond_angle_distribution`	`(atoms, center_type, neighbor_type, cutoff, bins) → (angles, hist)`
`compute_structure_factor`	`analysis.structure_factor`	`(atoms, q_min, q_max, n_q, r_max, n_bins, radiation, lorch_damping) → (q, sq, partials)`
`generate_bond_length_dict`	`analysis.rings`	`(atoms, specific_cutoffs) → dict`
`compute_guttmann_rings`	`analysis.rings`	`(atoms, bond_lengths, max_size, n_cpus) → (dist, mean)`
`compute_cavities`	`analysis.cavities`	`(atoms, resolution, cutoff_radii) → dict`
`visualize_cavities`	`analysis.cavities`	`(atoms, resolution, show_atoms, opacity, excluded_cavities) → Figure`
`analyze_structure`	`workflows.structural_analysis`	`(atoms) → StructureData`
`plot_analysis_results_plotly`	`workflows.structural_analysis`	`(StructureData) → Figure`

Key conventions:

target_type, center_type, neighbor_type are always scalar atomic numbers (int), not lists.
former_types, neighbor_types arguments that expect lists are explicitly typed list[int].
RDF pair keys are canonical: (min_Z, max_Z). CN keys are ordered: (Z_ref, Z_neighbour).
Cutoff for compute_qn is a single float — the Si-O (former-O) bond cutoff.