"""Pretty-printers for distance tables, J-shell summaries and four-state reports."""
from __future__ import annotations
import logging
from typing import Dict, List, Optional, Tuple
from mag4.constants import FourStateConfig
log = logging.getLogger(__name__)
[docs]
def print_mm_distances(results, element, cutoff, structure) -> None:
"""Tabular dump of all M–M pairs found within ``cutoff``."""
n_m = sum(1 for s in structure if s.specie.symbol == element)
W = 66
print()
print("=" * W)
print(f" {element}-{element} distances | cutoff = {cutoff} Å")
print(f" {n_m} {element} sites in unit cell | {len(results)} pairs found")
print("=" * W)
print(f" {'Pair':<18} {'Image (h k l)':<16} {'Distance (Å)':>12}")
print("-" * W)
for d in results:
pair = f"{d['label_i']} - {d['label_j']}"
image = "({:+d} {:+d} {:+d})".format(*d["image"])
print(f" {pair:<18} {image:<16} {d['distance']:>12.4f}")
print("=" * W)
if results:
dmin = results[0]["distance"]
dmax = results[-1]["distance"]
avg = sum(d["distance"] for d in results) / len(results)
print(f" d_min = {dmin:.4f} Å d_max = {dmax:.4f} Å d_avg = {avg:.4f} Å")
print("=" * W)
print()
def _fmt_per_fu(mult: int, z: int) -> str:
"""Format ``mult / z`` with one decimal only when not an integer."""
if z <= 0:
return ""
q = mult / z
return f"{int(q)}" if q == int(q) else f"{q:.2f}"
[docs]
def print_inequivalent(groups, element, tol,
*,
ligand_elements=None,
bond_cutoff: float = 2.5,
angle_tol: float = 1.0,
ll_cutoff: float = 0.0,
z_formula_units: int = 1) -> None:
"""Tabular dump of the inequivalent J shells.
The ``Multip.`` column gives the number of distinct M–M dimers per
primitive cell, taking periodic boundary conditions into account
(computed via pymatgen's :meth:`Structure.get_neighbors` with
canonical ``(min_idx, max_idx, image)`` deduplication). A second
``/f.u.`` column normalises that by ``z_formula_units`` (the number
of formula units in the primitive cell, taken from the reduced
composition).
When ``ligand_elements`` is supplied the shells were produced by
:func:`mag4.geometry.group_inequivalent_geom`; an extra ``M-L-M bridge``
column is shown. Without it, the legacy :func:`mag4.io_cif.group_inequivalent`
output is printed.
"""
geom = ligand_elements is not None
W = 110 if geom else 76
print()
print("=" * W)
print(f" Inequivalent {element}-{element} distances -> magnetic couplings J")
if geom:
ligand_str = "/".join(ligand_elements)
bits = [f"dist tol = {tol} Å",
f"bridge via {ligand_str} within {bond_cutoff:.2f} Å",
f"angle tol = {angle_tol}°"]
if ll_cutoff and ll_cutoff > 0:
bits.append(f"M-L-L-M fallback within {ll_cutoff:.2f} Å")
print(" (" + " | ".join(bits) + ")")
else:
print(f" (grouping tolerance = {tol} Å)")
print(f" Multiplicity counts unique dimers per primitive cell (PBC-aware). "
f"Z = {z_formula_units} f.u./cell.")
print("=" * W)
if geom:
print(f" {'Coupling':<10} {'d_mean (Å)':>10} "
f"{'Multip.':>8} {'/f.u.':>6} {'M-L-M bridge / M-L-L-M dih.':<40} Pairs")
else:
print(f" {'Coupling':<10} {'d_mean (Å)':>10} {'d_min (Å)':>10} "
f"{'d_max (Å)':>10} {'Multip.':>8} {'/f.u.':>6} Pairs")
print("-" * W)
for g in groups:
pair_str = " ".join(
f"{p['label_i']}-{p['label_j']}" for p in g["pairs"])
per_fu = _fmt_per_fu(g["multiplicity"], z_formula_units)
if geom:
fp_str = g.get("fingerprint_str", "")
print(f" {g['label']:<10} {g['distance']:>10.4f} "
f"{g['multiplicity']:>8} {per_fu:>6} "
f"{fp_str:<40} {pair_str}")
else:
print(f" {g['label']:<10} {g['distance']:>10.4f} "
f"{g['dist_min']:>10.4f} {g['dist_max']:>10.4f} "
f"{g['multiplicity']:>8} {per_fu:>6} {pair_str}")
print("=" * W)
print(f"\n => {len(groups)} inequivalent coupling(s): "
+ ", ".join(g["label"] for g in groups))
print()
[docs]
def print_report(all_configs: List[FourStateConfig],
unique_groups: List[List[FourStateConfig]],
atoms: list,
couplings: list,
ref_bath: list,
n_sym_ops: int,
supercell: Tuple[int, int, int],
output_file: Optional[str] = None,
*,
per_config_cells: Optional[List[dict]] = None,
n_supercell_atoms: Optional[int] = None,
multiplicities: Optional[Dict[str, int]] = None,
config_coefficients: Optional[List] = None,
ring_coefficients: Optional[List[int]] = None,
ring_edge: Optional[float] = None,
spectator_shells: Optional[List] = None,
spectator_coefficients: Optional[List] = None,
tr_pairs: Optional[List[Tuple[int, int]]] = None,
config_msgs: Optional[List] = None) -> None:
"""Print (and optionally write) the full four-state analysis report.
``per_config_cells`` (optional) is a list, one dict per unique
config in ``unique_groups`` order, with keys ``cell`` (6-tuple
a,b,c,α,β,γ), ``n_atoms`` (int), ``volume_A3`` (float), ``sg_number``
(int), ``sg_symbol`` (str). When supplied, a one-line cell summary
is emitted after each ``Config #N`` block so ``mag4-extract`` can
rescale per-config energies for the four-state formula.
"""
lines_out: List[str] = []
def out(msg: str = ""):
print(msg)
lines_out.append(msg)
n_atoms = len(atoms)
n_total = len(all_configs)
n_unique = len(unique_groups)
idx_to_uid: Dict[int, int] = {}
for uid, group in enumerate(unique_groups, 1):
for cfg in group:
idx_to_uid[cfg.idx] = uid
na, nb, nc = supercell
out("=" * 80)
out(" FOUR-STATE METHOD – CONFIGURATION ANALYSIS")
out("=" * 80)
out()
out(f" Supercell : {na}×{nb}×{nc}")
out(f" Magnetic atoms : {n_atoms}")
if n_supercell_atoms is not None:
# Total atom count in the supercell (all elements). This is
# the reference cell size the four-state formula is defined
# against — ``mag4-extract`` reads this line and rescales
# each config's WIEN2k :ENE up to this count (DFT extensivity),
# so configs computed in reduced cells contribute correctly
# to the J extraction.
out(f" Total atoms in supercell : {n_supercell_atoms}")
out(f" Couplings : "
+ ", ".join(f"{c[0]}={c[1]:.4f} Å" for c in couplings))
out(f" Symmetry operations : {n_sym_ops}")
out()
labels = [a[0] for a in atoms]
out("Reference magnetic bath:")
out(" " + " ".join(
f"{labels[i]}:{'+' if ref_bath[i] > 0 else '-'}"
for i in range(n_atoms)
))
out()
for entry in couplings:
coup_name = entry[0]
out(f"--- {coup_name} ---")
cfgs = [c for c in all_configs if c.coupling_name == coup_name]
if not cfgs:
out(" No dimer found for this coupling!")
out()
continue
d = cfgs[0].dimer
out(f" Dimer : {d.label_a} -- {d.label_b} (d = {d.distance:.4f} Å)")
out(f" Dimer indices : [{d.idx_a}] [{d.idx_b}]")
out()
out(f" {'State':<6} {'Spin vector':<60} {'Config #':>8}")
out(f" {'------'} {'--'*30} {'--------':>8}")
for c in cfgs:
sv_str = " ".join(f"{s:+d}" for s in c.spin_vector)
out(f" {c.state_label:<6} {sv_str:<60} #{idx_to_uid[c.idx]}")
out()
out("=" * 80)
out(f" UNIQUE CONFIGURATIONS : {n_unique} out of {n_total} total")
out(f" Savings : {n_total - n_unique} redundant DFT calculation(s) avoided")
out("=" * 80)
out()
# --check-degeneracy: pairs the grey group (unitary ops × time reversal)
# proves degenerate but that were kept as separate DFT runs on purpose.
if tr_pairs:
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.")
out(" mag4-extract reports |ΔE| per pair — a direct test of the")
out(" antiunitary symmetry and of the spin model's completeness.")
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()
for uid, group in enumerate(unique_groups, 1):
sv_str = " ".join(f"{s:+d}" for s in group[0].spin_vector)
member_strs = [f"{m.coupling_name}({m.state_label})" for m in group]
out(f" Config #{uid}: {', '.join(member_strs)}")
out(f" Spins: {sv_str}")
# Magnetic (Shubnikov) group of this spin configuration — primed
# (time-reversal-combined) operations included; BNS number is the
# citable label (symbol lookup: Bilbao MGENPOS).
if config_msgs is not None and uid - 1 < len(config_msgs) \
and config_msgs[uid - 1] is not None:
m = config_msgs[uid - 1]
pg = f" · point group {m.point_group()}" if m.mpg_symbol else ""
out(f" Magnetic group: {m.short()}{pg} "
f"[{m.n_unitary} unitary + {m.n_primed} primed ops]")
if per_config_cells is not None and uid - 1 < len(per_config_cells):
meta = per_config_cells[uid - 1]
cell = meta["cell"]
# Stable, parseable format — keys=values separated by ' | '.
# ``n_atoms`` is the count WIEN2k's ``:ENE`` actually refers
# to (= sum of MULTs in the struct = conv cell count /
# centring factor); ``mag4-extract`` reads this value for
# LCM rescaling. ``n_atoms_conv`` is the conv cell count
# for human inspection.
extras = ""
if "n_atoms_conv" in meta and meta["n_atoms_conv"] != meta["n_atoms"]:
extras = (f" | n_atoms_conv={meta['n_atoms_conv']}"
f" | lattice={meta.get('lattice_letter', '?')}")
out(f" Cell: a={cell[0]:.4f} | b={cell[1]:.4f} | "
f"c={cell[2]:.4f} | α={cell[3]:.3f} | β={cell[4]:.3f} | "
f"γ={cell[5]:.3f} | n_atoms={meta['n_atoms']}{extras} | "
f"V={meta['volume_A3']:.3f} ų | "
f"SG=#{meta['sg_number']} ({meta['sg_symbol']})")
out()
out("-" * 80)
out(" EXTRACTION FORMULAS")
out(" (each Config # refers to one DFT total energy)")
out("-" * 80)
seen_formulas: Dict[tuple, str] = {}
formula_warnings: List[str] = []
formula_ids: Dict[str, tuple] = {} # coup_name -> (ids, divisor)
for entry in couplings:
coup_name, coup_dist = entry[0], entry[1]
cfgs = [c for c in all_configs if c.coupling_name == coup_name]
if len(cfgs) != 4:
out(f" {coup_name} (d={coup_dist:.4f} Å): INCOMPLETE -- only {len(cfgs)}"
f" states found. Cell too small -- see supercell warnings above.")
continue
ids = {c.state_label: idx_to_uid[c.idx] for c in cfgs}
formula_key = (ids['uu'], ids['dd'], ids['ud'], ids['du'])
# divisor = 4·M, where M is the dimer multiplicity (how many Jk bonds
# link the two atoms, counting periodic images). M=1 for a geometrically
# isolated dimer; M>1 in a smaller cell where the dimer also couples to
# its own image through the SAME J — the four-state self-coefficient is
# then 4·M, so we divide by 4·M (not the textbook /4).
m_k = (multiplicities or {}).get(coup_name, 1)
div = 4 * m_k
formula_ids[coup_name] = (ids, div)
formula_str = (f"( E[#{ids['uu']}] + E[#{ids['dd']}]"
f" - E[#{ids['ud']}] - E[#{ids['du']}] ) / {div}")
if formula_key in seen_formulas:
prev = seen_formulas[formula_key]
out(f" {coup_name} (d={coup_dist:.4f} Å) = {formula_str}")
warn = (f" *** {coup_name} and {prev} share the same Config numbers: "
f"their dimers produce equivalent spin patterns. "
f"Enlarge the supercell to give them distinct representatives.")
formula_warnings.append(warn)
out(warn)
else:
seen_formulas[formula_key] = coup_name
out(f" {coup_name} (d={coup_dist:.4f} Å) = {formula_str}")
out()
if formula_warnings:
out(" NOTE: couplings sharing Config numbers cannot be extracted")
out(" independently from those DFT runs. A larger supercell is needed.")
out()
has_ring = ring_coefficients is not None
if config_coefficients:
coup_names = [entry[0] for entry in couplings]
out("-" * 80)
out(" MAGNETIC-ORDER ENERGY EQUATIONS (per config, to check the formulas)")
ring_hdr = " + c_ring·(Jring·S⁴)" if has_ring else ""
out(f" E[#i] = E0 + Σ_k c_k·(J_k·S²){ring_hdr}; "
"c_k = Σ_⟨ij⟩ s_i·s_j" + (" c_ring = Σ_□ s_i s_j s_k s_l"
if has_ring else " over the cell"))
out("-" * 80)
for uid, group in enumerate(unique_groups, 1):
if uid - 1 >= len(config_coefficients):
break
coefs, sz = config_coefficients[uid - 1]
terms = " ".join(
f"{coefs.get(nm, 0):+d}·({nm}·S²)" for nm in coup_names)
ring_term = ""
if has_ring and uid - 1 < len(ring_coefficients):
ring_term = f" {ring_coefficients[uid - 1]:+d}·(Jring·S⁴)"
out(f" E[#{uid}] = E0 {terms}{ring_term} [|Sz|={sz}]")
out()
out(" Check: substituting these into an EXTRACTION FORMULA above must")
out(" leave only its own J (all other c_k combine to 0).")
out()
# Accumulate every leak term so the report can close with a single corrected-J
# formula per Jk: coup_name -> [(term_label, leak_norm), …], where the
# reported four-state Jk = Jk_true + Σ leak_norm·term·S² (so J_true subtracts
# them). Populated by the --jring and --jprime blocks below.
corrections: Dict[str, List[Tuple[str, float]]] = {}
# --jring: how much the four-spin ring exchange LEAKS into each bilinear J.
# The four-state combination has no Jring term, so a non-zero leak biases the
# extracted J by (leak/divisor)·Jring·S⁴ — i.e. the reported J (= JS²/S²)
# shifts by (leak/divisor)·Jring·S².
if has_ring and ring_coefficients and formula_ids:
def _r(uid):
return ring_coefficients[uid - 1] if uid - 1 < len(ring_coefficients) else 0
out("-" * 80)
out(" Jring CONTAMINATION CHECK (four-spin ring leaking into each J)")
if ring_edge is not None:
out(f" ring on the {ring_edge:.4f} Å square plaquette; leak = the "
"J-formula's ± pattern applied to c_ring")
out("-" * 80)
any_leak = False
for coup_name, (ids, div) in formula_ids.items():
leak = _r(ids['uu']) + _r(ids['dd']) - _r(ids['ud']) - _r(ids['du'])
if leak == 0:
out(f" {coup_name}: ring leak = 0 → clean (no Jring contamination)")
else:
any_leak = True
norm = leak / div
# Ring enters E as Jring·S⁴; the reported (bilinear) J is the
# four-state quantity ÷ S², so the leak shows up as Jring·S².
corrections.setdefault(coup_name, []).append(("Jring·S²", norm))
out(f" {coup_name}: ring leak = {leak:+d} / {div} → reported "
f"{coup_name} = {coup_name}_true {norm:+g}·Jring·S²")
if any_leak:
out(" ⚠ A non-zero leak means the bilinear four-state J is biased by")
out(" the ring term. Fit Jring explicitly (--method energy-mapping")
out(" --jring, or --method 16-state for Jring alone) for the true J.")
out()
# --jprime: how much each BEYOND-CUTOFF bilinear shell J' contributes to the
# configN energies and leaks into each extracted Jk. The four-state
# combination only cancels couplings that do not fold onto a self-image, so a
# J' with d > --cutoff can survive and bias Jk in a cell that is too small.
if spectator_shells and spectator_coefficients and formula_ids:
def _s(label, uid):
if uid - 1 >= len(spectator_coefficients):
return 0
return spectator_coefficients[uid - 1][0].get(label, 0)
out("-" * 80)
out(" SPECTATOR-COUPLING (J') CONTAMINATION CHECK")
out(" beyond-cutoff bilinear shells: contribution to each configN energy")
out(" (c'[#i]) and leak into each extracted J (the J-formula's ± pattern")
out(" applied to c', ÷ divisor)")
out("-" * 80)
any_leak = False
clean_shells: List[str] = []
for label, dist in spectator_shells:
# Does this shell leak into ANY extracted J?
leaks = {}
for coup_name, (ids, div) in formula_ids.items():
num = (_s(label, ids['uu']) + _s(label, ids['dd'])
- _s(label, ids['ud']) - _s(label, ids['du']))
if num:
leaks[coup_name] = (num, div)
if not leaks:
clean_shells.append(label)
continue
any_leak = True
cper = " ".join(f"#{uid}:{_s(label, uid):+d}"
for uid in range(1, n_unique + 1))
out(f" {label} (d={dist:.4f} Å) c'[#i] = {cper}")
for coup_name, (ids, div) in formula_ids.items():
if coup_name in leaks:
num, div = leaks[coup_name]
# J' is bilinear (enters E as J'·S², same power as Jk), so
# after the four-state ÷ S² the leak is just (num/div)·J' —
# NO S² factor (unlike the four-spin ring term). Tag the
# consolidated term with the shell distance so far (≈ 0)
# spectators are obvious at a glance in the corrected-J line.
corrections.setdefault(coup_name, []).append(
(f"{label}[{dist:.2f}Å]", num / div))
out(f" {coup_name} ← {label}: leak = {num:+d} / {div} → "
f"reported {coup_name} = {coup_name}_true "
f"{num / div:+g}·{label}")
else:
out(f" {coup_name} ← {label}: leak = 0 (clean)")
out()
if clean_shells:
out(" Clean (leak = 0 into every extracted J): "
+ ", ".join(clean_shells))
if any_leak:
out(" ⚠ A non-zero leak means the extracted J is biased by that")
out(" beyond-cutoff shell in THIS cell. Grow the supercell so the")
out(" shell's distance is below the shortest lattice translation (no")
out(" self-image folding), or raise --cutoff to fit it explicitly.")
out()
# Consolidated correction: combine every leak (ring + each spectator J')
# collected above into one corrected-J formula per Jk. J_true SUBTRACTS the
# leaks, so a +norm leak term enters here with the opposite sign. Symbolic:
# the four-state run has no DFT energies yet, so substitute the four-state Jk
# (from mag4-extract) and Jring / J' (from a mapping fit or estimate).
if corrections:
out("-" * 80)
out(" CORRECTED (true) bilinear J — four-state Jk with every leak removed")
out("-" * 80)
for coup_name in [entry[0] for entry in couplings]:
terms = corrections.get(coup_name)
if not terms:
continue
# `lbl` already carries its decoration: "Jring·S²" for the ring
# term; "J5[6.61Å]" for a bilinear spectator (distance tagged so
# negligible far shells are obvious).
rhs = " ".join(f"{-coef:+g}·{lbl}" for lbl, coef in terms)
out(f" {coup_name}_true = {coup_name}_4state {rhs}")
out(" (Jk_4state = the value mag4-extract prints from the four-state")
out(" formula; substitute Jring / J' from a mapping fit or estimate.)")
out()
if output_file:
with open(output_file, "w") as f:
f.write("\n".join(lines_out) + "\n")
log.info("Report written to %s", output_file)