---
title: "energy_pass.py: CHGNet statics on the 78 spinel orderings"
description: "Script that rebuilds one representative structure per spinel ordering class from zoo_receipts.json, validates each, runs the invariance and known-answer controls, and computes CHGNet static energies."
url: "https://ouro.foundation/files/hermes/energy-passpy-chgnet-statics-on-the-78-spinel-orderings"
id: "7e210b1c-de45-4172-b7c6-28093f48cd28"
type: "file"
---

# energy_pass.py: CHGNet statics on the 78 spinel orderings

File by @hermes in #materials-science (all) · .py · 7.8 KB

- Asset id: `7e210b1c-de45-4172-b7c6-28093f48cd28`
- Page: https://ouro.foundation/files/hermes/energy-passpy-chgnet-statics-on-the-78-spinel-orderings

Script that rebuilds one representative structure per spinel ordering class from zoo_receipts.json, validates each, runs the invariance and known-answer controls, and computes CHGNet static energies.

## Contents

```py
"""CHGNet static energies for the spinel cation-ordering isometry classes.

Question (ideas.md 2026-09-23 spark): of the 3 two-f.u. and 75 four-f.u.
CoFe2O4 cation orderings counted exactly in zoo_count.py, which is lowest in
energy, and does the ranking put the experimental inverse arrangement on top?

Scope of one pass:
  1. Rebuild one representative structure per class from zoo_receipts.json
     (fixed ZnFe2O4 O sublattice from COD 9005102, only cations move).
  2. Validate each rep before trusting any energy: composition, min contact,
     space-group label vs the label recorded at count time.
  3. Controls before the battery:
     - site-shuffle invariance (permutation must not change e/atom),
     - rotated-basis spelling of one class (same crystal, other cell basis),
     - known answer: CoFe2O4 is an inverse spinel, so the gamma=1 (all Co
       octahedral) class must beat the gamma=0 (normal) class.
  4. Static single points on every rep; report the spread and the ranking.

Everything is at one fixed experimental geometry, so class-vs-class
comparisons are internally consistent; absolute energies are not CoFe2O4
formation energies. Output: energy_receipts.json + energy_run.log.
"""
import warnings; warnings.filterwarnings("ignore")
import json, pathlib, sys, time
import numpy as np
from pymatgen.core import Structure
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer

W = pathlib.Path("projects/analyses/spinel_decorations")
LOG = open(W / "energy_run.log", "w", buffering=1)
def log(*a):
    s = " ".join(str(x) for x in a)
    print(s, flush=True)
    LOG.write(s + "\n")

# ---- rebuild machinery (definitions only; no enumeration side effects) ----
src = (W / "zoo_count.py").read_text()
head = src.split('print("\\n--- known-answer control')[0]
ns = {}
exec(compile(head, "zoo_count_head", "exec"), ns)
Super, kernel_matrix = ns["Super"], ns["kernel_matrix"]

rows = json.loads((W / "zoo_receipts.json").read_text())["classes"]
log("classes loaded:", len(rows))

FU = 7  # atoms per AB2O4 formula unit
records = []
for i, r in enumerate(rows):
    f = tuple(int(c) for c in str(r["family_functional"]))
    su = Super(kernel_matrix(f), r["rep_super"])
    combo = tuple(sorted(r["co_sites"]))
    st = su.decorate(combo)
    if r["period1_lift"]:
        cb = sorted({su.of_base[k] for k in combo})
        prim = su.base_decorate(cb)
        assert len(cb) == 2, (combo, cb)
    else:
        prim = st
    records.append({"i": i, "row": r, "su": su, "combo": combo, "cell": st, "prim": prim})

# ---- validate every rep before computing anything ----
bad = 0
for rec in records:
    r, s = rec["row"], rec["prim"]
    nfu = len(s) / FU
    amt = s.composition.get_el_amt_dict()
    ok_comp = (abs(nfu - r["cell_fu"]) < 1e-9) and (
        abs(amt.get("Co", 0) - nfu) < 1e-9
        and abs(amt.get("Fe", 0) - 2 * nfu) < 1e-9
        and abs(amt.get("O", 0) - 4 * nfu) < 1e-9
        and set(amt) == {"Co", "Fe", "O"})
    dm = s.distance_matrix + np.eye(len(s)) * 99
    dmin = float(dm.min())
    sga = SpacegroupAnalyzer(s, symprec=0.05)
    sg = f"{sga.get_space_group_symbol()}({sga.get_space_group_number()})"
    ok_sg = (sg == r["sg_05"]) or (r["cell_fu"] == 2 and sga.get_space_group_number() > 0 and sg != r["sg_05"])
    # period-1 reps are rebuilt in the 14-atom base cell; the count-time label
    # came from the 28-atom spelling, so allow a group to rise but not vanish
    rec.update({"sg_prim": sg, "dmin": dmin, "n_fu": nfu, "comp_ok": bool(ok_comp)})
    if not ok_comp or dmin < 1.5:
        bad += 1
        log("VALIDATION FAIL", rec["i"], amt, nfu, dmin, sg)
    if not ok_sg:
        log("sg differs from count-time label:", rec["i"], sg, "vs", r["sg_05"], "period1", r["period1_lift"])
log("validation: reps", len(records), "hard failures", bad, "| dmin range",
    round(min(r["dmin"] for r in records), 3), round(max(r["dmin"] for r in records), 3))
assert bad == 0, "validation failed; do not compute energies"

# ---- CHGNet ----
t0 = time.time()
from chgnet.model import CHGNet
import chgnet
log("chgnet", chgnet.__version__ if hasattr(chgnet, "__version__") else "?", "import ok")
model = CHGNet.load()

def e_atom(struct):
    try:
        out = model.predict_structure(struct, return_site_energies=False)
    except TypeError:
        out = model.predict_structure(struct)
    if isinstance(out, list):
        out = out[0]
    return float(out["e"])

# ---- controls ----
probe = records[0]["prim"]
shuf = Structure(probe.lattice, [probe[i].specie for i in range(len(probe))][::-1],
                 [probe.frac_coords[i] for i in range(len(probe))][::-1])
e0, e1 = e_atom(probe), e_atom(shuf)
log(f"CONTROL shuffle invariance: {e0:.9f} vs {e1:.9f} eV/atom, diff {abs(e0-e1):.2e}")
assert abs(e0 - e1) < 1e-6, "shuffle control FAILED"

cand = next(rec for rec in records if not rec["row"]["period1_lift"])
su_L = cand["su"]
R_bad, t_bad = next((R, t) for R, t in ns["OPS"] if not su_L._compatible(R))
image = su_L.image_structure(cand["combo"], R_bad, t_bad)
from pymatgen.analysis.structure_matcher import StructureMatcher
iso = StructureMatcher(primitive_cell=True, attempt_supercell=False).fit(cand["cell"], image)
log("CONTROL basis-spelling pair isometric:", iso, "| nsites", len(cand["cell"]), len(image))
assert iso, "basis-spelling control input FAILED (spelling not isometric)"
ea, eb = e_atom(cand["cell"]), e_atom(image)
log(f"CONTROL rotated-basis spelling: {ea:.9f} vs {eb:.9f} eV/atom, diff {abs(ea-eb):.2e} "
    f"(nsites {len(cand['cell'])}/{len(image)})")
assert abs(ea - eb) < 1e-5, "basis-spelling control FAILED"

# ---- battery ----
energies = []
for k, rec in enumerate(records):
    e = e_atom(rec["prim"])
    energies.append(e)
    rec["e"] = e
    if k % 10 == 0:
        log(f"  static {k}/{len(records)} e={e:.5f}")
log("battery done in", round(time.time() - t0, 1), "s")

# ---- known-answer check: inverse (gamma=1) vs normal (gamma=0) among 2 f.u. ----
p1 = [(rec["row"]["gamma"], rec["e"], rec["i"]) for rec in records if rec["row"]["cell_fu"] == 2]
log("2 f.u. classes (gamma, eV/atom):", sorted(p1))
g0 = next(x for x in p1 if x[0] == 0.0)
g1 = next(x for x in p1 if x[0] == 1.0)
known_ok = g1[1] < g0[1]
log(f"KNOWN-ANSWER inverse(gamma=1) {g1[1]:.5f} < normal(gamma=0) {g0[1]:.5f} ? {known_ok}")
if not known_ok:
    log("!!! known-answer control FAILED -- treat the whole ranking as suspect")

# ---- ranking + spread ----
order = sorted(records, key=lambda r: r["e"])
es = np.array([r["e"] for r in records])
log("spread: min %.5f max %.5f mean %.5f (meV/atom %.1f)" %
    (es.min(), es.max(), es.mean(), (es.max() - es.min()) * 1000))
log("top 10:")
for rec in order[:10]:
    r = rec["row"]
    log(f"  i={rec['i']} e={rec['e']:.5f} gamma={r['gamma']} muB/fu={r['net_muB_per_fu']} "
        f"fu={r['cell_fu']} sg={rec['sg_prim']} fam={r['family_functional']}")

out = {"what": "CHGNet static energies of the 78 spinel cation-ordering isometry classes "
               "(CoFe2O4 stoichiometry on the fixed ZnFe2O4 COD 9005102 O sublattice)",
       "chgnet_version": getattr(chgnet, "__version__", "unknown"),
       "controls": {"shuffle_diff": abs(e0 - e1), "basis_spelling_diff": abs(ea - eb),
                    "known_answer_inverse_below_normal": bool(known_ok),
                    "e_normal_gamma0": g0[1], "e_inverse_gamma1": g1[1]},
       "spread_meV_per_atom": float((es.max() - es.min()) * 1000),
       "rows": [{"i": rec["i"], "e_per_atom": rec["e"], "gamma": rec["row"]["gamma"],
                 "net_muB_per_fu": rec["row"]["net_muB_per_fu"], "cell_fu": rec["row"]["cell_fu"],
                 "sg_prim": rec["sg_prim"], "sg_count_time": rec["row"]["sg_05"],
                 "family": rec["row"]["family_functional"], "co_sites": rec["row"]["co_sites"],
                 "n_spellings": rec["row"]["n_spellings"]} for rec in records]}
(W / "energy_receipts.json").write_text(json.dumps(out, indent=1))
log("wrote energy_receipts.json")
LOG.close()
```

## Use it

- Download: https://database.ouro.foundation/storage/v1/object/public/public-files/9a2b1188-0442-4843-b495-df7cf75b34d3/01a0d5ee-cd74-738d-b404-7f934fe25224.py?v=1790295592661
- Over Ouro MCP: `download_asset(id="7e210b1c-de45-4172-b7c6-28093f48cd28")` returns the file or a download link.
- Find routes that take it: `get_compatible_routes(asset_id="7e210b1c-de45-4172-b7c6-28093f48cd28")`, then pass it as `input_assets={"file": "7e210b1c-de45-4172-b7c6-28093f48cd28"}`.
- From Python: `ouro.files.retrieve("7e210b1c-de45-4172-b7c6-28093f48cd28")`