"""16-state method for the four-spin **ring (cyclic) exchange** ``Jring``.
Layer: **common** — the ring analogue of the four-state method.
Where the four-state method flips the two atoms of a *bond* (2² = 4 states) and
combines them to isolate a pairwise ``J``, the 16-state method flips the four
atoms of one square *plaquette* ``(i, j, k, l)`` over a fixed bath (2⁴ = 16
states) and projects with the four-spin parity ``χ = s_i s_j s_k s_l``
(mathematically, the Walsh character of the spin-flip group Z₂⁴)::
Jring·S⁴ = (1 / 16M) · Σ_configs χ(config) · E(config)
By parity (character) orthogonality the projection annihilates **every**
lower-order term — the constant, single-spin (bath couplings), pairwise bonds,
and partial (≤3-spin) plaquettes — leaving only the genuine four-spin
coefficient. ``M`` is the plaquette multiplicity (how many squares close on the
same four atoms under PBC; ``M = 1`` for a fully isolated plaquette).
``χ`` is invariant under both a global spin flip (even number of factors) and
any permutation of the four corners, so symmetry-equivalent configurations carry
the **same** ``χ`` and the **same** energy. The 16 states therefore collapse
(via :func:`mag4.spinconfig.find_unique_configs`) to a handful of inequivalent
ones, and the projection is rewritten over the representatives with integer
weights ``w_u = (group size)·χ_u``::
Jring·S⁴ = (Σ_u w_u · E[#u]) / (16M)
Intended for 2D **square-lattice** ring exchange (cuprates, etc.); the plaquette
is the nearest-neighbour square auto-detected by
:func:`mag4.energy_mapping.select_ring_plaquettes`.
"""
from __future__ import annotations
import itertools
import logging
import os
import re
from dataclasses import dataclass
from typing import Dict, List, Optional, Tuple
import numpy as np
from mag4.constants import CrystalData, FourStateConfig, next_config_idx
from mag4.energy_mapping import (
_cartesian_with_spin,
_image_shell,
config_energy_coefficients,
count_ring,
select_ring_plaquettes,
)
log = logging.getLogger(__name__)
[docs]
@dataclass
class RingPlaquette:
"""One square plaquette selected for the ring term.
``idx`` are the four **cell-atom indices** of the corners (in
``crystal.target_species`` order); ``labels`` their labels; ``edge`` the
square edge length (Å); ``multiplicity`` the number of distinct squares that
close on exactly these four atoms (the ``M`` in the ``16M`` divisor).
"""
idx: Tuple[int, int, int, int]
labels: Tuple[str, str, str, str]
edge: float
multiplicity: int
[docs]
def select_ring_quad(crystal: CrystalData, couplings: list, *,
prefer_label: Optional[str] = None,
prefer_dist: Optional[float] = None,
tol_pos: float = 0.15) -> RingPlaquette:
"""Pick the NN square plaquette for the ring term and return it as four
distinct cell-atom indices.
Raises ``ValueError`` if no coupling shell forms a square, or if every square
folds onto fewer than four distinct atoms (cell too small — enlarge it).
"""
lat = np.asarray(crystal.lattice_matrix, dtype=float)
cell = _cartesian_with_spin(crystal.target_species, lat)
image = _image_shell(cell, lat)
em_coup = [(c[0], c[1]) for c in couplings]
plaq, _n_uni, _folded, dist = select_ring_plaquettes(
em_coup, cell, image, prefer_label=prefer_label,
prefer_dist=prefer_dist, tol_pos=tol_pos)
if not plaq:
raise ValueError(
"no coupling shell forms a square plaquette — the 16-state ring "
"method needs a square lattice (pin one with --jring-coupling / "
"--jring-dist, or check --cutoff).")
idx_of = {a[0]: i for i, a in enumerate(cell)}
quad = None
for corners in plaq:
ids = tuple(idx_of[c] for c in corners)
if len(set(ids)) == 4:
quad = (ids, corners)
break
if quad is None:
raise ValueError(
"every square plaquette folds onto fewer than four distinct atoms "
"in this cell — enlarge the supercell so one plaquette has four "
"independent corners.")
ids, corners = quad
target = set(ids)
mult = sum(1 for c in plaq
if {idx_of[x] for x in c} == target) // 4
return RingPlaquette(idx=tuple(ids), labels=tuple(corners),
edge=float(dist), multiplicity=max(1, mult))
[docs]
def find_ring_supercell(crystal: CrystalData, couplings: list, *,
max_mult: int = 4,
prefer_label: Optional[str] = None,
prefer_dist: Optional[float] = None
) -> Tuple[Tuple[int, int, int], CrystalData, RingPlaquette]:
"""Smallest supercell in which the target square plaquette is **fully
isolated** — multiplicity ``M = 1``: four independent atoms forming exactly
one square, so the four-state-style divisor is exactly ``16`` (no ``16·M``).
This is the ring analogue of the four-state isolation criterion. At ``M > 1``
the four atoms wrap onto their own ``J1`` images through the periodic
boundary, so the plaquette is not independent of its copies and the DFT
energies pick up self-interaction bias — the projection math tolerates it,
but the physical extraction does not.
Searches ``(na, nb, nc)`` from smallest atom count upward (a layered square
plaquette isolates with ``nc = 1``). Returns ``((na, nb, nc), expanded
crystal, plaquette)``. Raises ``RuntimeError`` if no cell within
``max_mult`` per axis reaches ``M = 1``.
"""
from mag4.supercell import expand_supercell # local import: avoid cycle
for cand in sorted(itertools.product(range(1, max_mult + 1), repeat=3),
key=lambda t: (t[0] * t[1] * t[2], max(t), t)):
sc = expand_supercell(crystal, *cand)
try:
quad = select_ring_quad(sc, couplings, prefer_label=prefer_label,
prefer_dist=prefer_dist)
except ValueError:
continue
if quad.multiplicity == 1:
return cand, sc, quad
raise RuntimeError(
f"no supercell up to {max_mult}×{max_mult}×{max_mult} isolates the ring "
f"plaquette to multiplicity M=1. Raise --max-supercell, or check that "
f"the magnetic lattice forms a square (--jring-coupling / --jring-dist "
f"to pin the shell).")
[docs]
def ring_design_records(crystal: CrystalData, couplings: list,
unique_groups: List[List[FourStateConfig]],
quad: RingPlaquette, *,
prefer_label: Optional[str] = None,
prefer_dist: Optional[float] = None,
tol: float = 5e-4) -> List[dict]:
"""Design-matrix rows of the unique ring configs (no extra DFT needed).
For each unique configuration the model energy is exactly
E[#u] = E0 + Σ_k c_k·(J_k·S²) + c_ring·(Jring·S⁴)
with ``c_k = Σ_{⟨ij⟩∈k} s_i s_j`` over the unique bonds of shell ``k`` and
``c_ring`` the four-spin plaquette coefficient. Returns energy-mapping
style records ``{"number", "degeneracy", "sz", "coefficients"}`` where the
ring column is labelled ``"Jring"`` — the raw material for deciding *if and
how* the pairwise ``J_k`` can be extracted from the very configurations
the 16-state method already computes for ``Jring``.
"""
atoms = crystal.target_species
em = [(c[0], c[1]) for c in couplings]
coefs = config_energy_coefficients(
atoms, crystal.lattice_matrix, em,
[g[0].spin_vector for g in unique_groups], tol=tol)
lat = np.asarray(crystal.lattice_matrix, dtype=float)
cell_atoms = _cartesian_with_spin(atoms, lat)
image_atoms = _image_shell(cell_atoms, lat)
plaq, _nu, _nf, _dist = select_ring_plaquettes(
em, cell_atoms, image_atoms,
prefer_label=prefer_label, prefer_dist=prefer_dist)
names = [a[0] for a in cell_atoms]
records: List[dict] = []
for uid, (group, (ck, sz)) in enumerate(zip(unique_groups, coefs), 1):
c_ring = count_ring(plaq, {names[i]: group[0].spin_vector[i]
for i in range(len(names))}) if plaq else 0
coefficients = dict(ck)
coefficients["Jring"] = c_ring
records.append({"number": uid, "degeneracy": len(group), "sz": sz,
"coefficients": coefficients})
return records
[docs]
def ring_chi(config: FourStateConfig, quad: RingPlaquette) -> int:
"""Four-spin parity ``χ = Π_{m∈plaquette} s_m`` (±1) for one configuration."""
prod = 1
for t in quad.idx:
prod *= int(config.spin_vector[t])
return prod
[docs]
def generate_ring_16(ref_bath: list, quad: RingPlaquette) -> List[FourStateConfig]:
"""The 16 spin configurations flipping the four plaquette atoms over the bath.
Only the four plaquette corners differ between configs; the rest stay at
``ref_bath``. ``state_label`` is the ``u``/``d`` pattern of the corners.
"""
configs: List[FourStateConfig] = []
for bits in itertools.product((+1, -1), repeat=4):
spins = list(ref_bath)
for t, s in zip(quad.idx, bits):
spins[t] = s
label = "".join("u" if b > 0 else "d" for b in bits)
configs.append(FourStateConfig(
idx=next_config_idx(),
coupling_name="Jring",
state_label=label,
spin_vector=tuple(spins),
dimer=None,
))
return configs
[docs]
def ring_weights(unique_groups: List[List[FourStateConfig]],
quad: RingPlaquette) -> List[int]:
"""Integer projection weight ``w_u = (group size)·χ_u`` per unique config.
All members of a symmetry group share ``χ`` (permutation/flip invariant), so
the per-group weight is well defined; this is verified defensively.
"""
weights: List[int] = []
for group in unique_groups:
chis = {ring_chi(c, quad) for c in group}
if len(chis) != 1:
raise ValueError(
"a symmetry group mixes χ = ±1 — the plaquette is not invariant "
"under the supplied symmetry operations (internal error).")
weights.append(len(group) * chis.pop())
return weights
# ---------------------------------------------------------------------------
# Report (Method: 16-state — recognised by energy_mapping.detect_report_method)
# ---------------------------------------------------------------------------
[docs]
def write_ring_report(path: Optional[str],
unique_groups: List[List[FourStateConfig]],
atoms: list,
quad: RingPlaquette,
supercell: Tuple[int, int, int],
n_sym_ops: int,
*,
n_supercell_atoms: Optional[int] = None,
ref_bath: Optional[list] = None,
design_records: Optional[List[dict]] = None,
coupling_dists: Optional[Dict[str, float]] = None,
config_sgs: Optional[List[Tuple[int, str]]] = None,
config_msgs: Optional[List] = None,
symmetric_structs: bool = False,
tr_pairs: Optional[List[Tuple[int, int]]] = None
) -> None:
"""Print (and optionally write) the 16-state ring report.
Optional extras (all analysis-only — the written structures stay P1):
* ``design_records`` (from :func:`ring_design_records`) — adds the
PAIRWISE-J ANALYSIS section: the design-matrix rows of the unique
configs and, per coupling, whether/why it is extractable from the very
energies computed for ``Jring``.
* ``coupling_dists`` — ``{label: d(Å)}`` shown next to each verdict.
* ``config_sgs`` — per-config ``(spacegroup_number, symbol)`` of the
spin-decorated structure (spglib), shown as an extra column.
* ``tr_pairs`` (from ``--check-degeneracy``) — config pairs related only
by the antiunitary symmetry (unitary op × time reversal); their DFT
energies must coincide and ``mag4-extract`` reports ``|ΔE|`` per pair.
"""
lines: List[str] = []
def out(msg: str = ""):
print(msg)
lines.append(msg)
weights = ring_weights(unique_groups, quad)
divisor = 16 * quad.multiplicity
labels = [a[0] for a in atoms]
na, nb, nc = supercell
out("=" * 80)
out(" 16-STATE METHOD – FOUR-SPIN RING EXCHANGE (Jring)")
out("=" * 80)
out("Method: 16-state")
out(f" Supercell : {na}×{nb}×{nc}")
out(f" Magnetic atoms : {len(atoms)}")
if n_supercell_atoms is not None:
out(f" Total atoms in supercell : {n_supercell_atoms}")
out(f" Ring plaquette : "
+ " -- ".join(quad.labels) + f" (square, edge = {quad.edge:.4f} Å)")
out(f" Plaquette indices : " + " ".join(f"[{i}]" for i in quad.idx))
out(f" Multiplicity M : {quad.multiplicity}"
+ (" (M=1 ⇒ fully isolated plaquette)"
if quad.multiplicity == 1 else
" (M>1: squares fold onto these atoms; a larger cell gives M=1)"))
out(f" Symmetry operations : {n_sym_ops}")
if ref_bath is not None:
nup = sum(1 for s in ref_bath if s > 0)
kind = ("FM" if nup == len(ref_bath)
else "AFM/mixed" if nup else "FM(↓)")
out(f" Reference bath ({kind}): "
+ " ".join(f"{s:+d}" for s in ref_bath))
out()
if tr_pairs is not None:
out(f" 16 plaquette spin states reduced to {len(unique_groups)} "
f"config(s) by UNITARY crystal symmetry only (--check-degeneracy:")
out(f" {len(tr_pairs)} time-reversal pair(s) kept separate as a "
f"degeneracy test" + ("" if tr_pairs else
" — none needed; the unitary group already closes") + ").")
else:
out(f" 16 plaquette spin states reduced to {len(unique_groups)} "
f"inequivalent config(s) by crystal symmetry + time reversal.")
out(f" χ = product of the four corner spins (the four-spin parity).")
out()
has_sg = config_sgs is not None and len(config_sgs) == len(unique_groups)
has_msg = (config_msgs is not None
and len(config_msgs) == len(unique_groups)
and any(m is not None for m in config_msgs))
sg_hdr = f" {'SG':<16}" if has_sg else ""
msg_hdr = (f" {'MSG (BNS)':<12} {'MPG':<10}") if has_msg else ""
out(f" {'Config #':>8} {'χ':>3} {'size':>4}{sg_hdr}{msg_hdr} Spin vector")
out(f" {'--------':>8} {'---':>3} {'----':>4}"
+ (f" {'--'*8}" if has_sg else "")
+ (f" {'--'*6} {'--'*5}" if has_msg else "") + f" {'--'*20}")
for uid, group in enumerate(unique_groups, 1):
chi = ring_chi(group[0], quad)
sv = " ".join(f"{s:+d}" for s in group[0].spin_vector)
sg_col = ""
if has_sg:
num, sym = config_sgs[uid - 1]
sg_col = f" {f'#{num} ({sym})':<16}"
msg_col = ""
if has_msg:
m = config_msgs[uid - 1]
if m is not None:
msg_col = (f" {f'{m.bns_number}-{m.type_roman}':<12}"
f" {m.mpg_symbol or '—':<10}")
else:
msg_col = f" {'—':<12} {'—':<10}"
out(f" #{uid:<7d} {chi:>+3d} {len(group):>4d}{sg_col}{msg_col} {sv}")
if has_sg or has_msg:
out()
if has_sg:
out(" SG = space group of each spin-decorated configuration (spglib).")
if symmetric_structs:
out(" Each WIEN2k configN.struct carries this symmetry (Wyckoff "
"orbits +")
out(" symmetry ops), always in the ONE shared supercell; VASP "
"POSCARs stay P1.")
else:
out(" Shown for information only — the written structures "
"stay P1.")
if has_msg:
out(" MSG = MAGNETIC (Shubnikov) space group, BNS number-type: primed")
out(" operations (unitary × time reversal) included — the SG column")
out(" colours ↑/↓ as species and sees only the unitary part. Types:")
out(" I colourless · III black-white · IV with anti-translations")
out(" (symbol lookup: Bilbao MGENPOS by BNS number).")
out(" MPG = magnetic POINT group (translations dropped): G colourless,")
out(" G1' GREY (an anti-translation becomes pure time reversal),")
out(" G(H) black-white in Shubnikov notation (H = unitary subgroup).")
if tr_pairs:
out()
out("-" * 80)
out(" TIME-REVERSAL DEGENERACY CHECK (--check-degeneracy)")
out("-" * 80)
out(" Each pair below is related only by an ANTIUNITARY operation")
out(" (unitary crystal op × time reversal): the two DFT energies must")
out(" coincide although no unitary operation relates them (note the")
out(" different total Sz). mag4-extract reports |ΔE| per pair — a")
out(" direct test of the antiunitary symmetry and of the spin model.")
for a, b in tr_pairs:
sz_a = sum(unique_groups[a - 1][0].spin_vector)
sz_b = sum(unique_groups[b - 1][0].spin_vector)
out(f" pair: config{a} == config{b} (Sz {sz_a:+d} / {sz_b:+d})")
out()
out("-" * 80)
out(" EXTRACTION FORMULA (yields Jring·S⁴)")
out(" (each Config # refers to one DFT total energy; bath = reference)")
out("-" * 80)
out(format_ring_formula(weights, divisor, quad.edge))
out(f" where {divisor} = 16 × M(={quad.multiplicity}). The χ-weighted sum "
"cancels all")
out(" ≤3-spin terms (pairwise J, bath, partial plaquettes), leaving Jring.")
out("=" * 80)
if design_records:
rec_labels = [l for l in design_records[0]["coefficients"]
if l != "Jring"] + ["Jring"]
verdicts = analyze_pairwise_extractability(design_records, rec_labels)
width = max(6, max(len(l) for l in rec_labels))
out()
out("-" * 80)
out(" PAIRWISE-J ANALYSIS (same configurations — no extra DFT)")
out("-" * 80)
out(" Each config energy obeys the exact model")
out(" E[#u] = E0 + Σ_k c_k·(J_k·S²) + c_ring·(Jring·S⁴)")
out(" so any coupling whose coefficient VARIES independently across "
"the configs")
out(" below can be least-squares fitted from the SAME energies "
"(mag4-extract does")
out(" this automatically alongside the χ-projection).")
out()
if coupling_dists:
out(" Couplings: " + ", ".join(
f"{l}={coupling_dists[l]:.4f}"
for l in rec_labels if l in coupling_dists))
out(f" Design: {'config':>10} {'deg':>5} {'Sz':>5} "
+ " ".join(f"{l:>{width}}" for l in rec_labels))
for rec in design_records:
out(f" Design: {'config' + str(rec['number']):>10} "
f"{rec['degeneracy']:>5} {rec['sz']:>5} "
+ " ".join(f"{rec['coefficients'].get(l, 0):>{width}d}"
for l in rec_labels))
out()
for lbl in rec_labels:
ok, why = verdicts[lbl]
d = (f" (d={coupling_dists[lbl]:.4f} Å)"
if coupling_dists and lbl in coupling_dists else "")
if ok:
out(f" {lbl}{d}: EXTRACTABLE from these configurations.")
else:
out(f" {lbl}{d}: NOT extractable — {why}.")
out(f" → use --method four-state or energy-mapping "
f"for {lbl}.")
n_ok = sum(1 for l in rec_labels if verdicts[l][0])
out()
out(f" Fit: E0 + {n_ok} coupling(s) from {len(design_records)} "
f"energies "
f"({len(design_records) - 1 - n_ok} degree(s) of freedom left).")
out()
out(" CAVEAT — the pairwise fit assumes NO couplings beyond the shells")
out(" listed above. Any unmodelled shell within the plaquette's reach")
out(" (a longer-range J, a 3-spin term) leaks into the fitted Jₖ; "
"mag4-extract")
out(" flags this via the fit RMS and the fit-vs-projection gap. The")
out(" χ-projection Jring is immune to all of these (parity "
"orthogonality) and")
out(" is the robust deliverable; for reliable pairwise Jₖ prefer "
"--method")
out(" four-state, or widen --cutoff so every relevant shell is in "
"the model.")
out("=" * 80)
if path:
with open(path, "w") as f:
f.write("\n".join(lines) + "\n")
log.info("Ring report written to %s", path)
# ---------------------------------------------------------------------------
# Extraction (consumed by mag4-extract --method 16-state)
# ---------------------------------------------------------------------------
_RING_RE = re.compile(
r"Jring\b.*?=\s*\(\s*(?P<terms>[^)]*?)\s*\)\s*/\s*(?P<div>\d+)\s*$")
_TERM_RE = re.compile(r"(?P<w>[+-]\d+)\s*\*\s*E\[#(?P<id>\d+)\]")
[docs]
def parse_ring_report(report_path: str) -> RingFormula:
"""Parse the ``Jring = ( … ) / D`` line of a 16-state report."""
edge = None
with open(report_path) as f:
for line in f:
m_edge = re.search(r"edge=([\d.]+)", line)
if m_edge and edge is None:
try:
edge = float(m_edge.group(1))
except ValueError:
pass
m = _RING_RE.search(line)
if m:
terms = [(int(t.group("w")), int(t.group("id")))
for t in _TERM_RE.finditer(m.group("terms"))]
if not terms:
raise ValueError(
f"ring formula has no E[#i] terms: {line.strip()!r}")
return RingFormula(terms=terms, divisor=int(m.group("div")),
edge=edge)
raise ValueError(f"no 'Jring = ( … ) / D' formula found in {report_path}")
[docs]
def parse_ring_design(report_path: str) -> Optional[dict]:
"""Parse the PAIRWISE-J ANALYSIS design rows of a 16-state report.
Returns ``{"labels": [...], "records": [{"number", "degeneracy", "sz",
"coefficients"}, ...]}`` (energy-mapping record format, ``Jring`` column
included) or ``None`` when the report has no analysis section (older
reports).
"""
labels: Optional[List[str]] = None
records: List[dict] = []
with open(report_path) as f:
for line in f:
s = line.strip()
if not s.startswith("Design:"):
continue
fields = s[len("Design:"):].split()
if fields and fields[0] == "config": # header row
labels = fields[3:]
continue
if labels is None or not fields[0].startswith("config"):
continue
number = int(fields[0][len("config"):])
deg, sz = int(fields[1]), int(fields[2])
coefficients = {l: int(v) for l, v in zip(labels, fields[3:])}
records.append({"number": number, "degeneracy": deg, "sz": sz,
"coefficients": coefficients})
if labels is None or not records:
return None
return {"labels": labels, "records": records}
[docs]
def compute_jring(formula: RingFormula, energies: Dict[int, object]
) -> Tuple[Optional[float], List[str]]:
"""Apply the ring formula to the per-config energies.
``energies`` maps config id → object with an ``energy_eV`` attribute (an
:class:`mag4.energy.EnergyResult`). Returns ``(Jring·S⁴ in eV, problems)``;
the value is ``None`` when an energy is missing/unparsed.
"""
problems: List[str] = []
total = 0.0
for w, cid in formula.terms:
res = energies.get(cid)
e = getattr(res, "energy_eV", None) if res is not None else None
if e is None:
problems.append(f"config{cid}: no energy")
continue
total += w * float(e)
if problems:
return None, problems
return total / float(formula.divisor), problems