"""CLI: multi-sublattice LSWT magnon dispersion + Tyablikov RPA T_c."""
from __future__ import annotations
import argparse
import logging
import sys
from typing import Dict, NamedTuple, Optional
from pymatgen.core import Structure
from mag4.io_cif import get_distance_shells
from mag4.magnon.dispersion import compute_magnon_dispersion
from mag4.magnon.exchange import (
build_jq_matrices,
compute_eigenvalues_all_q,
)
from mag4.magnon.kpath import get_kpath
from mag4.magnon.ordering import analyse_ordering
from mag4.magnon.plot import plot_magnon_bands
from mag4.magnon.qmesh import auto_qmesh, build_qmesh
from mag4.magnon.report import print_report
from mag4.magnon.sublattice import get_sublattice_bonds
from mag4.magnon.tc import compute_tc, compute_tc_lswt
logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
log = logging.getLogger(__name__)
[docs]
def parse_j_values(s: str) -> Dict[str, float]:
"""Parse ``"J1:0.26,J2:0.26"`` → ``{"J1": 0.26, "J2": 0.26}``."""
result = {}
for item in s.split(","):
name, val = item.strip().split(":")
result[name.strip()] = float(val.strip())
return result
[docs]
def parse_j_values_file(path: str) -> ExtractedJs:
"""Read a ``j_values.dat`` written by mag4-extract.
Handles both the four-state format (``name d JS2(eV) JS2(meV) J(meV)
J(K) status``) and the energy-mapping format (``name d JS2(meV)
[sig(meV)] …``). ``j_values`` are the **JS² values in meV**
(mag4-magnon's native input); a ``Jring`` row's JS2 column holds the
fitted Jring·S⁴ and is returned separately in ``ring``. The per-coupling
standard errors (``sig(meV)`` column, newer files) and the fit RMS from
the header feed the statistical-significance test used by the
dimensionality guard.
"""
import os
import re
if os.path.isdir(path):
path = os.path.join(path, "j_values.dat")
spin = None
fourstate: bool = True
has_sig = False
rms_meV = None
jv: Dict[str, float] = {}
dists: Dict[str, float] = {}
sigmas: Dict[str, float] = {}
skipped = []
ring_info = None
with open(path) as fh:
for line in fh:
s = line.strip()
if not s:
continue
if s.startswith("#"):
if "energy-mapping J values" in s:
fourstate = False
if "sig(meV)" in s:
has_sig = True
m = re.search(r"spin:\s*([\d.]+)", s)
if m:
spin = float(m.group(1))
m = re.search(r"RMS\s*=\s*([\d.eE+-]+)\s*meV", s)
if m:
try:
rms_meV = float(m.group(1))
except ValueError:
pass
continue
cols = s.split()
if len(cols) < 3 or not cols[0][0].isalpha():
continue
name = cols[0]
try:
d = float(cols[1])
except ValueError:
continue
if "ring" in name.lower():
# JS2 column of a Jring row holds the fitted Jring·S⁴.
try:
ring = {"label": name, "dist": d,
"js4_meV": (float(cols[3]) if fourstate
else float(cols[2]))}
except (ValueError, IndexError):
ring = None
if ring is not None and (s.endswith(" OK")
or s.endswith("\tOK")):
ring_info = ring
skipped.append((name, "four-spin coupling — a bilinear LSWT "
"model cannot include it directly; "
"pass --jring to fold it into "
"effective bilinear couplings "
"(LSWT ring renormalisation)"))
continue
if not s.endswith(" OK") and not s.endswith("\tOK"):
skipped.append((name, f"status not OK ({cols[-1]})"))
continue
try:
js2_meV = float(cols[3]) if fourstate else float(cols[2])
except (ValueError, IndexError):
skipped.append((name, "unparseable row"))
continue
jv[name] = js2_meV
dists[name] = d
if has_sig and not fourstate and len(cols) >= 4:
try:
sigmas[name] = float(cols[3])
except ValueError:
pass
return ExtractedJs(jv, spin, dists, skipped, ring_info, sigmas, rms_meV)
def _network_dimensionality(bond_data, j_tol: float = 1e-9,
exclude_values: frozenset = frozenset()) -> int:
"""Dimensionality of the exchange network: rank of the span of the bond
displacement vectors carrying a resolved, non-zero J (3 = 3D network).
``exclude_values`` holds the J values of couplings classified as
*statistically zero* by the caller (|JS²| below its fit error bar / the
fit RMS / an explicit ``--j-significance`` floor) — physical smallness
alone is never a reason to exclude a coupling, since any genuinely
finite interlayer J orders a quasi-2D magnet."""
import numpy as np
vecs = [v for lst in bond_data["bonds"].values()
for (v, J) in lst
if abs(J) > j_tol and J not in exclude_values]
if not vecs:
return 0
s = np.linalg.svd(np.asarray(vecs, dtype=float), compute_uv=False)
return int((s > 1e-6 * s[0]).sum()) if s[0] > 0 else 0
def _flat_axes_at_q0(evals, q_frac, tol_rel: float = 1e-3):
"""Axes along which the acoustic branch is flat at the LT minimum.
Takes the mesh line through the global λ_min point varying one
fractional coordinate at a time; an axis whose λ_min spread is below
``tol_rel`` of the full bandwidth cannot pin that component of q₀
(e.g. an interlayer coupling ≈ 0 leaves qz undetermined). Returns
``[(axis_letter, spread_meV), …]``.
"""
import numpy as np
lam = evals.min(axis=1)
total = float(lam.max() - lam.min())
if total <= 0:
return []
i0 = int(np.argmin(lam))
q0 = q_frac[i0]
flat = []
for ax in range(3):
others = [i for i in range(3) if i != ax]
mask = np.ones(len(q_frac), dtype=bool)
for i in others:
d = np.abs(((q_frac[:, i] - q0[i]) + 0.5) % 1.0 - 0.5)
mask &= d < 1e-6
if mask.sum() > 1:
spread = float(lam[mask].max() - lam[mask].min())
if spread < tol_rel * total:
flat.append(("abc"[ax], spread))
return flat
[docs]
def apply_ring_renormalisation(j_values: Dict[str, float],
shells: list,
edge_label: str,
js4_meV: float,
*, tol: float = 0.02):
"""Fold a four-spin ring coupling into effective bilinear JS² (LSWT).
Harmonic (LSWT) reduction of the cyclic four-spin operator about the
Néel state (cf. Coldea et al., PRL 86, 5377 (2001)): each pair operator
inside the quartic term is multiplied by the classical value of its
partner bond, so the magnon Hamiltonian is EXACTLY that of a bilinear
model with, per bond and in native JS² units (Jring·S⁴ enters directly —
no explicit S factors),
edge shell (2 plaquettes/bond): JS² → JS² − 2·JringS⁴
diagonal shell (1 plaquette/bond): JS² → JS² − 1·JringS⁴
``edge_label`` names the shell the square plaquette is built on; the
diagonal shell is located geometrically (distance = √2 × edge, within
``tol`` relative). If the diagonal shell carries no fitted J it is
*induced* with JS² = −JringS⁴. Returns ``(new_j_values, lines)`` where
``lines`` document the transformation for the report.
"""
import math
shell_d = {sh["label"]: sh["distance"] for sh in shells}
if edge_label not in shell_d:
raise ValueError(f"--jring: edge shell '{edge_label}' not among the "
f"CIF distance shells within --cutoff "
f"({', '.join(shell_d)}).")
d_edge = shell_d[edge_label]
d_diag = math.sqrt(2.0) * d_edge
diag_label = None
for sh in shells:
if abs(sh["distance"] - d_diag) <= tol * d_diag:
diag_label = sh["label"]
break
out = dict(j_values)
lines = [f" --jring: LSWT ring renormalisation applied "
f"(JringS⁴ = {js4_meV:.4f} meV on the {d_edge:.4f} Å plaquette;",
" harmonic reduction of the cyclic four-spin operator about the "
"Néel state,",
" cf. Coldea et al., PRL 86, 5377 (2001)):"]
old_e = out.get(edge_label, 0.0)
out[edge_label] = old_e - 2.0 * js4_meV
lines.append(f" edge shell {edge_label} (d={d_edge:.4f} Å, "
f"2 plaquettes/bond): JS² {old_e:.4f} → "
f"{out[edge_label]:.4f} meV")
if diag_label is not None:
old_d = out.get(diag_label, 0.0)
out[diag_label] = old_d - js4_meV
induced = "" if diag_label in j_values else " [induced — no fitted J]"
lines.append(f" diagonal shell {diag_label} "
f"(d={shell_d[diag_label]:.4f} Å, 1 plaquette/bond): "
f"JS² {old_d:.4f} → {out[diag_label]:.4f} meV{induced}")
else:
lines.append(f" diagonal shell (d≈{d_diag:.4f} Å) is beyond "
"--cutoff: its induced −1·JringS⁴ contribution cannot "
"be represented — raise --cutoff to include it.")
return out, lines
_LOW_DIM_WARNING = """\
⚠ DIMENSIONALITY — the exchange network spans only {dim}D.
An isotropic Heisenberg model cannot order at T > 0 in two or fewer
dimensions (Mermin–Wagner theorem), so the temperatures above are NOT a
Néel / Curie temperature: they are the mean-field / RPA energy scale of
the short-range correlations in this {dim}D network. (On a finite q-mesh
the RPA value is mesh-dependent and → 0 as the mesh is refined — the RPA
manifestation of Mermin–Wagner.) Real ordering in layered / 2D magnets
is enabled by what this model omits: magnetic anisotropy (single-ion,
XXZ/Ising exchange anisotropy) and/or the interlayer coupling — the true
T_N/T_C is set by those small terms (see e.g.
https://en.wikipedia.org/wiki/Magnetic_2D_materials). Either add the
interlayer J shell (raise --cutoff so it enters the model) or quote the
printed values as an in-plane energy scale, not an ordering temperature."""
[docs]
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="mag4-magnon",
description=("Multi-sublattice magnon LSWT — magnon dispersion via "
"eigenvalues of J(q) in the primitive cell."),
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog="""\
Convention: H = +sum_ij J_ij S_i.S_j (AFM: J>0, FM: J<0)
Input (one of):
--JS2 fitted JS2 = +(E_uu+E_dd-E_ud-E_du)/4 (native; the JS2
column of j_values.dat)
--J bare J = JS2/S^2 (the J(meV) column); x S^2 applied
--from-extract read j_values.dat directly (couplings + spin)
Magnon dispersion: omega_k(q) = 2S * (lambda_k(J(q)) - lambda_min(q0))
Always evaluated in the primitive cell -> n_sub branches.
Goldstone mode (E=0) at k=q0 by construction of the shift.
Works for any q0 (Gamma or not) without supercell.
Example (trigonal, fitted JS2 in meV):
mag4-magnon --cif Gd2O2S.cif --element Gd \\
--JS2 "J1:0.26,J2:0.26,J3:0.12,J4:0.01,J5:0.39" --spin 3.5 \\
--plot --kpath "Gamma M K Gamma A L H"
Example (fitted JS2 in eV):
mag4-magnon --cif structure.cif --element Gd \\
--JS2 "J1:0.00147,J2:-0.00061" --spin 3.5 --unit eV
Example (straight from mag4-extract — no hand-copying, spin from the file):
mag4-magnon --cif structure.cif --element Gd --from-extract MyCompound/
""",
)
parser.add_argument("--cif", required=True,
help="Primitive/conventional CIF file")
parser.add_argument("--element", required=True,
help="Magnetic element symbol, e.g. Gd")
parser.add_argument("--JS2", dest="JS2", default=None,
help="The fitted JS² values, exactly as mag4-extract "
"prints them (the JS2 column of j_values.dat / "
"'raw fit Jk*S^2'): 'J1:1.47,J2:-0.61,...' "
"(AFM > 0). Used as-is — this is mag4-magnon's "
"native input.")
parser.add_argument("--J", dest="J", default=None,
help="The bare J values, i.e. J = JS²/S² — the J(meV) "
"column of j_values.dat: 'J1:5.88,...' (AFM > 0)."
" Multiplied by S² internally (needs --spin); "
"the conversion is logged. Give ONE of --J, "
"--JS2 or --from-extract.")
parser.add_argument("--from-extract", dest="from_extract", default=None,
metavar="DIR|FILE",
help="Read the couplings straight from mag4-extract's "
"output: a compound directory (containing "
"j_values.dat) or the file itself. Takes the "
"JS² column (mag4-magnon's native input) and the "
"spin recorded in the header (--spin overrides); "
"raises --cutoff to cover the longest coupling "
"if needed. Non-OK rows and the four-spin Jring "
"are skipped with a note. Mutually exclusive "
"with --J.")
parser.add_argument("--spin", type=float, default=None,
help="Spin quantum number S (e.g. 3.5 for Gd³⁺). "
"Required unless --from-extract provides it.")
parser.add_argument("--jring", nargs="?", const="FROM_EXTRACT",
default=None, metavar="LABEL:JRINGS4",
help="Fold a four-spin ring coupling into effective "
"bilinear couplings via the LSWT (harmonic) "
"mapping about the Néel state: edge shell "
"JS² −= 2·JringS⁴, plaquette-diagonal shell "
"JS² −= 1·JringS⁴ (native JS² units — JringS⁴ "
"is exactly the Jring row's JS2 column in "
"j_values.dat). Bare --jring with "
"--from-extract takes the file's Jring row; "
"explicitly, give the edge shell and the fitted "
"JringS⁴, e.g. --jring J1:2.0433 (unit follows "
"--unit). The transformation is logged and "
"recorded in the report.")
parser.add_argument("--j-significance", dest="j_significance", type=float,
default=None, metavar="FLOOR",
help="Statistical floor (meV·S²) below which a "
"coupling counts as ZERO for the Mermin-Wagner "
"dimensionality check. Default (unset): decide "
"per coupling from the extraction itself — "
"|JS²| < 2σ of its fit error bar (sig column of "
"j_values.dat), else |JS²| < the fit RMS; with "
"no statistics available nothing is excluded. "
"Pass 0 to disable all exclusion. Physical "
"smallness alone never excludes a coupling — "
"only being unresolved from zero does.")
parser.add_argument("--cutoff", type=float, default=7.0,
help="M-M distance cutoff in Angstrom (default: 7.0)")
parser.add_argument("--qmesh", type=str, nargs="+",
metavar="N", default=["auto"],
help="q-mesh: 3 integers NA NB NC, or 'auto' (default)")
parser.add_argument("--target", type=int, default=125000,
help="Target q-points for auto mesh (default: 125000)")
parser.add_argument("--unit", choices=["meV", "eV"], default="meV",
help="Unit of --J values (default: meV)")
parser.add_argument("--name", default="magnon_sw",
help="Base name for output files (default: magnon_sw)")
parser.add_argument("--plot", action="store_true",
help="Plot magnon band structure (requires matplotlib)")
parser.add_argument("--kpoints", type=int, default=100,
help="k-points per segment (default: 100)")
parser.add_argument("--kpath", type=str, default=None,
help="Custom path, e.g. \"Gamma M K Gamma | A L H\". "
"Use '|' for segment breaks. "
"If omitted, the automatic pymatgen path is used.")
return parser
[docs]
def main(argv=None) -> int:
parser = build_parser()
args = parser.parse_args(argv)
n_sources = sum(x is not None for x in (args.J, args.JS2, args.from_extract))
if n_sources != 1:
log.error("Give exactly one of --JS2 (fitted JS² values), "
"--J (bare J = JS²/S²) or --from-extract.")
return 2
ring_file = None
dists: Dict[str, float] = {}
j_sigmas: Dict[str, float] = {}
fit_rms_meV: Optional[float] = None
if args.from_extract is not None:
try:
xj = parse_j_values_file(args.from_extract)
except (OSError, ValueError) as exc:
log.error("Cannot read mag4-extract output (%s): %s",
args.from_extract, exc)
return 2
j_values, file_spin, dists = xj.j_values, xj.spin, xj.dists
skipped, ring_file = xj.skipped, xj.ring
j_sigmas, fit_rms_meV = xj.sigmas, xj.rms_meV
if not j_values:
log.error("No usable coupling rows in %s.", args.from_extract)
return 2
for name, reason in skipped:
log.warning(" %s skipped: %s", name, reason)
log.info("Couplings read from mag4-extract output (%s): %s",
args.from_extract,
", ".join(f"{k}={v:.4f} meV.S2" for k, v in j_values.items()))
if args.spin is None:
if file_spin is None:
log.error("No spin recorded in the j_values.dat header — "
"pass --spin.")
return 2
args.spin = file_spin
log.info("Spin S = %g taken from the j_values.dat header.",
file_spin)
if args.unit == "eV":
log.warning("--unit eV ignored: j_values.dat carries meV.")
d_max = max(dists.values(), default=0.0)
if d_max > args.cutoff:
log.info("Raising --cutoff %.2f -> %.2f Å to cover the longest "
"extracted coupling.", args.cutoff, d_max + 0.05)
args.cutoff = d_max + 0.05
else:
if args.spin is None:
log.error("--spin is required with --J/--JS2.")
return 2
j_values = parse_j_values(args.JS2 if args.JS2 is not None else args.J)
if args.unit == "eV":
j_values = {k: v * 1000.0 for k, v in j_values.items()}
log.info("Values converted eV -> meV.")
if args.J is not None:
# --J takes the bare J = JS²/S²; the magnon pipeline is JS²-native.
s2 = args.spin ** 2
j_values = {k: v * s2 for k, v in j_values.items()}
log.info("--J takes the bare J: multiplied by S² = %g to build "
"JS² (%s).", s2,
", ".join(f"{k}={v:.4f}" for k, v in j_values.items()))
log.info("JS2 (meV.S2): %s",
", ".join(f"{k}={v:.4f}" for k, v in j_values.items()))
log.info("Loading CIF: %s", args.cif)
structure = Structure.from_file(args.cif)
log.info("Formula: %s Space group: %s",
structure.composition.reduced_formula,
structure.get_space_group_info())
log.info("Distance shells up to %.2f A ...", args.cutoff)
shells = get_distance_shells(args.cif, args.element, args.cutoff)
for sh in shells:
lbl = sh["label"]
J_str = f"{j_values[lbl]:.4f} meV.S2" if lbl in j_values else "not provided"
log.info(" %s d=%.4f A n_bonds=%d JS2=%s",
lbl, sh["distance"], sh["n_bonds"], J_str)
# --from-extract: the CIF/cutoff must reproduce the extraction's shell
# labels — a same-label/different-distance mismatch would silently
# misassign couplings.
if args.from_extract is not None:
shell_d = {sh["label"]: sh["distance"] for sh in shells}
for lbl, d in dists.items():
if lbl in shell_d and abs(shell_d[lbl] - d) > 0.05:
log.error("Shell label mismatch: %s is %.4f Å in this CIF but "
"%.4f Å in j_values.dat — a different cell or "
"cutoff than the extraction used. Aborting.",
lbl, shell_d[lbl], d)
return 2
# --jring: fold the four-spin ring coupling into effective bilinear JS².
ring_lines: list = []
if args.jring is not None:
if args.jring == "FROM_EXTRACT":
if args.from_extract is None:
log.error("Bare --jring needs --from-extract; otherwise give "
"--jring LABEL:JRINGS4 (e.g. --jring J1:2.0433).")
return 2
if ring_file is None:
log.error("--jring: no usable Jring row in the j_values.dat — "
"fit it with 'mag4-magnetic --method energy-mapping "
"--jring' or '--method 16-state'.")
return 2
js4 = ring_file["js4_meV"]
edge_label = None
for sh in shells:
if abs(sh["distance"] - ring_file["dist"]) \
<= 0.02 * sh["distance"]:
edge_label = sh["label"]
break
if edge_label is None:
log.error("--jring: no CIF shell at the ring distance "
"%.4f Å (within --cutoff).", ring_file["dist"])
return 2
else:
if ":" not in args.jring:
log.error("--jring takes LABEL:JRINGS4 (e.g. J1:2.0433), or "
"bare with --from-extract.")
return 2
edge_label, _, val = args.jring.partition(":")
edge_label = edge_label.strip()
try:
js4 = float(val)
except ValueError:
log.error("--jring: cannot parse '%s' as a number.", val)
return 2
if args.unit == "eV":
js4 *= 1000.0
try:
j_values, ring_lines = apply_ring_renormalisation(
j_values, shells, edge_label, js4)
except ValueError as exc:
log.error("%s", exc)
return 2
for ln in ring_lines:
log.info("%s", ln.strip())
log.info("Effective JS2 (meV.S2): %s",
", ".join(f"{k}={v:.4f}" for k, v in j_values.items()))
log.info("Building sublattice bond table ...")
bond_data = get_sublattice_bonds(args.cif, args.element, args.cutoff,
shells, j_values)
if len(args.qmesh) == 1 and args.qmesh[0].lower() == "auto":
na, nb, nc = auto_qmesh(structure, target_points=args.target)
elif len(args.qmesh) == 3:
na, nb, nc = (int(args.qmesh[0]), int(args.qmesh[1]), int(args.qmesh[2]))
else:
sys.exit("--qmesh requires 'auto' or 3 integers (NA NB NC)")
log.info("Building %d x %d x %d q-mesh ...", na, nb, nc)
q_cart, q_frac = build_qmesh(structure, na, nb, nc)
log.info("Computing J_ab(q) and diagonalising (%d points, %d x %d) ...",
na*nb*nc, bond_data["n_sub"], bond_data["n_sub"])
Jq = build_jq_matrices(bond_data, q_cart)
evals = compute_eigenvalues_all_q(Jq)
log.info("J(q) eigenvalue range: [%.4f, %.4f] meV.S2",
float(evals.min()), float(evals.max()))
log.info("Analysing ground-state ordering (Luttinger-Tisza) ...")
ordering_lt = analyse_ordering(Jq, evals, q_frac, q_cart, structure)
log.info(" LT q0 = (%+.4f, %+.4f, %+.4f) frac", *ordering_lt["q0_frac"])
log.info(" LT lambda_min = %.4f meV.S2 | %s",
ordering_lt["lam0"], ordering_lt["character"])
# Flat-axis diagnostic: a component of q0 along an axis with (near-)zero
# dispersion is numerical noise, not a physical spiral — e.g. an
# interlayer coupling fitted as 0 within its error bar leaves qz free.
flat_axes = _flat_axes_at_q0(evals, q_frac)
flat_note = ""
if flat_axes:
axes = ", ".join(a for a, _ in flat_axes)
spreads = ", ".join(f"{s:.4g} meV" for _, s in flat_axes)
flat_note = (
f" note: λ_min(q) is flat along {axes} within {spreads} "
f"(< 10⁻³ of the bandwidth) —\n"
f" the q₀ component(s) along {axes} are UNDETERMINED (the "
f"corresponding\n"
f" coupling is zero within precision); do not read an "
f"incommensurate\n"
f" modulation along {axes} into q₀.")
for ln in flat_note.splitlines():
log.warning("%s", ln.strip())
log.info("Computing T_c (Luttinger-Tisza RPA) ...")
tc_lt = compute_tc(evals, args.spin, ordering_lt["lam0"])
if not tc_lt["note"]:
log.info(" T_c^LT-MF = %.1f K", tc_lt["Tc_MF"])
log.info(" T_c^LT-RPA = %.1f K", tc_lt["Tc_RPA"])
log.info("Computing T_c (LSWT ordered-state RPA) ...")
tc_sw = compute_tc_lswt(bond_data, args.spin, ordering_lt, q_cart, evals)
if not tc_sw["note"]:
log.info(" T_c^%s = %.1f K", tc_sw.get("method", "RPA"), tc_sw["Tc_RPA"])
# The --jring provenance is printed INSIDE the report, directly under the
# exchange-parameter table, so effective couplings are never shown
# without their origin (console and report file alike).
print_report(shells, j_values, bond_data, evals,
ordering_lt, tc_lt, tc_sw,
args.spin, (na, nb, nc), structure,
f"{args.name}_report.dat",
j_values_note="\n".join(ring_lines) if ring_lines else None)
if flat_note:
print()
print(flat_note)
try:
with open(f"{args.name}_report.dat", "a") as fh:
fh.write("\n" + flat_note + "\n")
except OSError: # pragma: no cover
pass
# Mermin-Wagner guard. A coupling counts as ZERO for the dimensionality
# check only when it is STATISTICALLY unresolved from zero — never for
# physical smallness alone (any genuinely finite interlayer J orders a
# quasi-2D magnet). Criterion per coupling, in order of preference:
# explicit --j-significance floor > |JS²| < 2σ of its own fit error bar
# > |JS²| < the fit RMS > no exclusion.
insignificant: Dict[str, str] = {}
for k, v in j_values.items():
if args.j_significance is not None:
if args.j_significance > 0 and abs(v) < args.j_significance:
insignificant[k] = (f"|JS²| < --j-significance "
f"{args.j_significance:g} meV")
elif k in j_sigmas:
if abs(v) < 2.0 * j_sigmas[k]:
insignificant[k] = (f"|JS²| = {abs(v):.4g} < 2σ = "
f"{2.0 * j_sigmas[k]:.4g} meV — "
"consistent with zero")
elif fit_rms_meV is not None:
if abs(v) < fit_rms_meV:
insignificant[k] = (f"|JS²| = {abs(v):.4g} < fit RMS = "
f"{fit_rms_meV:.4g} meV — below the "
"fit's resolution")
exclude = frozenset(j_values[k] for k in insignificant)
dim = _network_dimensionality(bond_data, exclude_values=exclude)
if dim < 3:
warning = _LOW_DIM_WARNING.format(dim=dim)
if insignificant:
head = (" Couplings treated as ZERO for the dimensionality "
"check (STATISTICALLY\n unresolved from zero — override "
"with --j-significance 0):\n"
+ "".join(f" {k}: {why}\n"
for k, why in insignificant.items()))
warning = head + warning
print()
print(warning)
try:
with open(f"{args.name}_report.dat", "a") as fh:
fh.write("\n" + warning + "\n")
except OSError: # pragma: no cover - report file missing/read-only
pass
if args.plot:
log.info("Computing LSWT magnon band structure ...")
try:
kpath = get_kpath(structure, n_points=args.kpoints,
custom_labels=args.kpath)
magnon_bands, kpath = compute_magnon_dispersion(
bond_data, kpath, args.spin, ordering_lt)
out_png = plot_magnon_bands(
magnon_bands, kpath, args.spin,
bond_data["n_sub"], args.name)
if out_png:
print(f"\n Magnon band structure -> {out_png}")
except Exception as exc:
log.error("Magnon band plot failed: %s", exc)
import traceback
traceback.print_exc()
return 1
return 0
if __name__ == "__main__": # pragma: no cover
sys.exit(main())