"""CLI: extract J values from a completed four-state DFT campaign.
Usage example::
cd MyCompound/ # (contains configN/ + four_state_report.dat)
mag4-extract --code vasp --spin 3.5
# → prints a table; also writes MyCompound/j_values.dat
"""
from __future__ import annotations
import argparse
import logging
import os
import sys
from typing import List, Optional
import numpy as np
from mag4.constants import KB_MEV
from mag4.energy import RY_TO_EV, read_band_gap, read_magnetization
from mag4.energy_mapping import detect_report_method, parse_energy_mapping_report
from mag4.extract import (
DEFAULT_REPORT_NAME,
JResult,
compute_j_values,
parse_report_config_cells,
parse_report_config_sz,
parse_report_supercell,
parse_report_supercell_atoms,
read_struct_atom_counts,
parse_report_formulas,
parse_report_tr_pairs,
read_all_energies,
write_j_values_dat,
)
from mag4.fit import FitResult, fit_couplings, spin_power
from mag4.ring16 import (
analyze_pairwise_extractability,
compute_jring,
parse_ring_design,
parse_ring_report,
)
logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
log = logging.getLogger(__name__)
[docs]
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="mag4-extract",
description=("Extract Heisenberg J values from completed four-state "
"DFT calculations."),
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog="""\
Workflow
1. Run mag4-magnetic to generate the per-config inputs and the
<compound>/four_state_report.dat file.
2. Run the DFT code (VASP or WIEN2k) inside each <compound>/configN/.
3. Run mag4-extract to walk the configN/ tree, check convergence,
pull the total energies (TOTEN for VASP, :ENE for WIEN2k), and
apply the formulas from four_state_report.dat to obtain JS².
Per-config cell sizes (WIEN2k symmetric path, default)
mag4-magnetic --code wien2k writes each configN.struct in the
primitive cell of its colored magnetic structure, so different
configs can have different cell sizes. mag4-extract reads the
per-config 'Cell: ... n_atoms=N ...' line emitted under each
Config #N block of four_state_report.dat and rescales each
contributing config's total energy to the least-common-multiple
(LCM) atom count of the four configs of each J:
E_X(in N_lcm-atom cell) = E_X(in N_X-atom primitive) * (N_lcm / N_X)
This works because DFT total energy is extensive: scaling a primitive
cell to a k-multiplied supercell scales the total energy by the same
factor k exactly (no approximation). The four-state formula is then
applied to the scaled energies, yielding the same J value you would
get if every config had been computed in a single common cell.
When all configs share a cell (the default --no-reduce layout, or
--code vasp), the LCM is N_X itself and the scaling is a no-op.
Examples
cd MyCompound/
mag4-extract --code vasp # JS² in eV / meV
mag4-extract --code vasp --spin 3.5 # adds J in meV and K
mag4-extract --code wien2k --spin 3.5 # auto-detects per-config cells
mag4-extract --code vasp --compound MyCompound --report custom_report.dat
""",
)
parser.add_argument("--code", choices=["vasp", "wien2k"], required=True,
help="DFT code that produced the output files. "
"'vasp' reads OUTCAR; 'wien2k' reads <case>.scf.")
parser.add_argument("--scffile", default=None, metavar="NAME",
help="WIEN2k only. Basename (with or without the "
"'.scf' suffix) of the SCF file to read inside "
"each configN/ directory. Use this when you "
"save runs side by side as e.g. pbe.scf / "
"r2scan.scf — then '--scffile pbe' picks "
"configN/pbe.scf for every config. "
"Default: configN/configN.scf "
"(WIEN2k 'case == basename(pwd)' convention).")
parser.add_argument("--compound", default=".", metavar="DIR",
help="Compound directory containing configN/ + the "
"four_state_report.dat (default: cwd).")
parser.add_argument("--report", default=None, metavar="PATH",
help=f"Path to the report file "
f"(default: <compound>/{DEFAULT_REPORT_NAME}).")
parser.add_argument("--spin", type=float, default=None,
help="Magnetic spin S (e.g. 3.5 for Gd³⁺). When "
"given, the table also reports J = JS²/S² in "
"meV and J/k_B in Kelvin.")
parser.add_argument("--method",
choices=["auto", "four-state", "energy-mapping", "16-state"],
default="auto",
help="Extraction method. 'auto' (default) reads the "
"'Method:' tag of the report and dispatches; "
"'four-state' applies the difference formulas; "
"'energy-mapping' does the linear least-squares fit "
"over all enumerated orders; '16-state' applies the "
"χ-weighted ring projection to get Jring·S⁴.")
parser.add_argument("--weight-by-deg", action="store_true",
help="energy-mapping: weight each residual by "
"sqrt(degeneracy) in the least-squares fit.")
parser.add_argument("--keep-collinear", action="store_true",
help="energy-mapping: keep every coupling column "
"(skip the automatic drop of constant / aliased "
"columns). Individual J values may be unreliable.")
parser.add_argument("--n-sigma", type=float, default=2.0, metavar="K",
dest="n_sigma",
help="energy-mapping: coverage of the reported error "
"bars in standard deviations (default: 2 ≈ 95%%).")
parser.add_argument("--check-mag", action="store_true",
help="energy-mapping: exclude configurations whose "
"converged total moment is inconsistent with their "
"intended Sz (a sign the SCF settled into the wrong "
"magnetic state). Expected |mag| = μ·|Sz|.")
parser.add_argument("--moment", type=float, default=None, metavar="MU",
help="--check-mag: per-site moment μ_B used for the "
"expected |mag| = MU·|Sz| (overrides --moment-ref).")
parser.add_argument("--moment-ref", choices=["spin", "median"], default="spin",
dest="moment_ref",
help="--check-mag reference per-site moment when --moment "
"is unset: 'spin' (default) = saturated 2·S from "
"--spin (catches a global moment collapse); "
"'median' = median |mag|/|Sz| over the data.")
parser.add_argument("--mag-tol", type=float, default=0.5, metavar="TOL",
dest="mag_tol",
help="--check-mag per-site tolerance (μ_B): a config is "
"excluded when | |mag|/|Sz| − μ | > TOL (Sz≠0), or "
"|mag| > TOL for a singlet (Sz=0). Default 0.5.")
parser.add_argument("--exclude", type=int, nargs="+", default=None,
metavar="N", dest="exclude",
help="energy-mapping: configuration number(s) to drop "
"from the fit, e.g. --exclude 11 12 13 to remove "
"high-energy orders the bilinear model cannot "
"capture. Combined with any --check-mag exclusions; "
"excluded configs still appear in the per-config "
"table marked EXCLUDED.")
parser.add_argument("--couplings", type=str, default=None, metavar="LIST",
help="energy-mapping: fit ONLY these couplings and force "
"the rest to zero, re-solving the least-squares with "
"the reduced design matrix. Comma-separated, e.g. "
"--couplings J1,J2,J3 (drop Jring) or --couplings "
"J1,J2 (drop J3 and Jring). Mutually exclusive with "
"--drop-couplings.")
parser.add_argument("--drop-couplings", type=str, default=None, metavar="LIST",
dest="drop_couplings",
help="energy-mapping: fit all couplings EXCEPT these "
"(force them to zero) and re-solve, e.g. "
"--drop-couplings Jring or --drop-couplings J3,Jring. "
"Mutually exclusive with --couplings.")
parser.add_argument("--no-gap", action="store_true", dest="no_gap",
help="Skip the per-config band-gap check (by default "
"mag4-extract scans each OUTCAR / .scf for the gap "
"and warns about metallic configs).")
parser.add_argument("--gap-tol", type=float, default=0.05, metavar="EV",
dest="gap_tol",
help="Band-gap threshold (eV) below which a config is "
"flagged metallic (default: 0.05). A config with "
"partially occupied bands is metallic regardless.")
parser.add_argument("--output", default=None, metavar="PATH",
help="Where to write the J values text file "
"(default: <compound>/j_values.dat).")
return parser
# ---------------------------------------------------------------------------
# Pretty printing
# ---------------------------------------------------------------------------
def _fmt_or_na(x: Optional[float], width: int, fmt: str = "{:>{w}.4f}") -> str:
return fmt.format(x, w=width) if x is not None else "{:>{w}}".format("N/A", w=width)
[docs]
def print_j_table(results: List[JResult], spin: Optional[float] = None) -> None:
"""Pretty-print the J table to stdout."""
has_spin = spin is not None and spin > 0
print()
print("=" * 96)
print(" mag4-extract — Heisenberg J values from the four-state method")
print("=" * 96)
cols = ["Shell", "d (Å)", "JS² (eV)", "JS² (meV)"]
widths = [6, 10, 14, 12]
if has_spin:
cols += ["J (meV)", "J (K)"]
widths += [12, 10]
cols += ["Status"]
header = " ".join(f"{c:>{w}}" for c, w in zip(cols, widths))
print(" " + header + " Status")
print(" " + "-" * (sum(widths) + 2*len(widths)) + "-" * 8)
for r in results:
row = [
f"{r.name:>6}",
f"{r.distance_A:>10.4f}",
_fmt_or_na(r.JS2_eV, 14, "{:>{w}.6f}"),
_fmt_or_na(r.JS2_meV, 12, "{:>{w}.4f}"),
]
if has_spin:
row += [
_fmt_or_na(r.J_meV, 12, "{:>{w}.4f}"),
_fmt_or_na(r.J_K, 10, "{:>{w}.2f}"),
]
status = "OK" if r.computed else (r.error or "missing")
# Truncate status to keep the table readable
if len(status) > 60:
status = status[:57] + "..."
row.append(status)
print(" " + " ".join(row))
print("=" * 96)
# ---------------------------------------------------------------------------
# Band-gap / metallicity check
# ---------------------------------------------------------------------------
def _print_band_gaps(args, compound: str, config_ids, scf) -> dict:
"""Read each config's band gap and print a summary, warning about metallic
configs (no gap) — for a magnetic insulator these usually mean the SCF
settled into the wrong solution. Returns ``{config_id: BandGap}``."""
gaps = {idx: read_band_gap(os.path.join(compound, f"config{idx}"),
args.code, scf_basename=scf, gap_tol=args.gap_tol)
for idx in sorted(config_ids)}
parsed = {i: g for i, g in gaps.items() if g.gap_eV is not None}
if not parsed:
log.info("Band-gap check: no eigenvalue/`:GAP`/`:BAN` data found "
"(skipped).")
return gaps
W = 64
print()
print("=" * W)
print(" Band gaps (per config)")
print("=" * W)
metals = []
for idx, g in gaps.items():
if g.gap_eV is None:
print(f" config{idx:<4d} gap = n/a ({g.note})")
elif g.metallic:
metals.append(idx)
note = f" — {g.note}" if g.note else ""
print(f" config{idx:<4d} gap = {g.gap_eV:6.3f} eV ⚠ METALLIC{note}")
else:
print(f" config{idx:<4d} gap = {g.gap_eV:6.3f} eV")
print("=" * W)
if metals:
log.warning("%d config(s) are METALLIC (no band gap): %s. For a magnetic "
"insulator this usually signals the wrong SCF solution — "
"inspect those configs (tighter smearing, AFM start, …).",
len(metals), ", ".join(f"config{i}" for i in metals))
return gaps
# ---------------------------------------------------------------------------
# Convergence summary
# ---------------------------------------------------------------------------
def _print_convergence_summary(energies, code: str) -> int:
"""Print a one-line summary of convergence and warn about failures.
Returns the number of non-converged configurations.
"""
n_total = len(energies)
failed = [(idx, er) for idx, er in sorted(energies.items())
if not er.converged]
print()
print(f" Configurations read: {n_total} ({code} outputs)")
if not failed:
print(f" All {n_total} configurations converged. ✓")
else:
print(f" ⚠ {len(failed)} non-converged configuration(s):")
for idx, er in failed:
diag = er.diagnostics[0] if er.diagnostics else "(no diagnostic)"
print(f" config{idx}: {diag}")
if code == "wien2k":
_print_wien2k_verification_table(energies)
return len(failed)
def _print_wien2k_verification_table(energies) -> None:
"""One-line-per-config table showing the last ``:DIS`` and ``:ENE``
values the parser picked up — so the user can verify nothing else
in the SCF file slipped past us.
"""
print()
print(" WIEN2k parser verification (last :DIS and :ENE per config):")
header = (f" {'config':<9} {'iter':>5} {'last :DIS':>12} "
f"{'last :ENE (Ry)':>18} status")
print(header)
print(" " + "-" * (len(header) - 4))
for idx, er in sorted(energies.items()):
dis_s = f"{er.last_dis:.6g}" if er.last_dis is not None else "—"
if er.last_ene_line is not None and er.energy_eV is not None:
ene_ry = er.energy_eV / RY_TO_EV
ene_s = f"{ene_ry:.6f}"
else:
ene_s = "—"
status = "OK" if er.converged else "FAILED"
print(f" config{idx:<3} {er.n_iter:>5} "
f"{dis_s:>12} {ene_s:>18} {status}")
print()
def _print_config_energies(energies, config_sz=None, mags=None,
spin=None) -> None:
"""Per-config total energies (eV) for the four-state path, plus — when
moments are available — the converged ``mag`` vs the configuration's ``Sz``
(expected ``|mag| = 2·S·|Sz|``)."""
mu = 2.0 * spin if (spin and spin > 0) else None
has_mag = mags is not None and any(v is not None for v in mags.values())
print()
print(" Per-configuration energies")
hdr = f" {'config':>8} {'iter':>5} {'energy (eV)':>18}"
if has_mag:
hdr += f" {'Sz':>4} {'mag':>9} {'exp|mag|':>9}"
hdr += " status"
print(hdr)
print(" " + "-" * (len(hdr) - 4))
for idx, er in sorted(energies.items()):
e_s = f"{er.energy_eV:.6f}" if er.energy_eV is not None else "—"
line = f" config{idx:<3} {er.n_iter:>5} {e_s:>18}"
status = "OK" if er.converged else "FAILED"
if has_mag:
sz = config_sz.get(idx) if config_sz else None
mg = mags.get(idx)
exp = abs(mu * sz) if (mu is not None and sz is not None) else None
line += (f" {(sz if sz is not None else '—')!s:>4} "
f"{(f'{mg:.3f}' if mg is not None else '—'):>9} "
f"{(f'{exp:.3f}' if exp is not None else '—'):>9}")
if mg is not None and exp is not None and abs(abs(mg) - exp) > 0.5:
status += " ⚠mag≠Sz"
print(line + f" {status}")
print()
# |ΔE| above this is suspicious for a T-degenerate pair (paper-grade DFT
# agrees to tens of µeV); it signals inconsistent/unconverged runs or a
# symmetry error, not a failure of time reversal itself.
_TR_PAIR_WARN_MEV = 1.0
def _print_tr_pairs(pairs, energies) -> None:
"""|ΔE| per time-reversal-degenerate pair (``--check-degeneracy``).
Each pair is related only by an antiunitary operation (unitary crystal op
× time reversal), so the two DFT energies must coincide — their agreement
directly tests the antiunitary symmetry and the consistency of the runs.
"""
if not pairs:
return
W = 72
print()
print("=" * W)
print(" TIME-REVERSAL DEGENERACY CHECK (antiunitary symmetry)")
print("=" * W)
print(" Each pair must be degenerate: related by unitary op × time")
print(" reversal, which collinear no-SOC DFT obeys exactly.")
print(f" {'pair':>22} {'E_a (eV)':>16} {'E_b (eV)':>16} {'|ΔE| (meV)':>12}")
print(" " + "-" * (W - 2))
worst = None
n_ok = 0
for a, b in pairs:
ea = getattr(energies.get(a), "energy_eV", None)
eb = getattr(energies.get(b), "energy_eV", None)
tag = f"config{a} == config{b}"
if ea is None or eb is None:
missing = [f"config{i}" for i, e in ((a, ea), (b, eb)) if e is None]
print(f" {tag:>22} {'—':>16} {'—':>16} {'—':>12} "
f"(no energy: {', '.join(missing)})")
continue
de_mev = abs(ea - eb) * 1e3
flag = " ⚠" if de_mev > _TR_PAIR_WARN_MEV else ""
print(f" {tag:>22} {ea:>16.6f} {eb:>16.6f} {de_mev:>12.4f}{flag}")
n_ok += 1
if worst is None or de_mev > worst:
worst = de_mev
if worst is not None:
if worst <= _TR_PAIR_WARN_MEV:
print(f" max |ΔE| = {worst:.4f} meV ({worst*1e3:.1f} µeV) — "
"antiunitary degeneracy confirmed.")
else:
print(f" ⚠ max |ΔE| = {worst:.4f} meV exceeds {_TR_PAIR_WARN_MEV:g} "
"meV: the pair members disagree.")
print(" Time reversal is exact in collinear no-SOC DFT, so this "
"points to")
print(" inconsistent runs (k-mesh, smearing, U, convergence) or "
"a different")
print(" SCF solution in one member — re-examine before trusting "
"the J values.")
print("=" * W)
# ---------------------------------------------------------------------------
# Energy-mapping (linear least-squares) path
# ---------------------------------------------------------------------------
def _print_fit_table(res: FitResult, labels: List[str],
spin: Optional[float], n_sigma: float,
title: str = "energy-mapping (linear least-squares)"
) -> None:
has_S = spin is not None and spin > 0
S = spin if has_S else 1.0
ksig = n_sigma
have_se = bool(res.se)
W = 80
print()
print("=" * W)
print(f" mag4-extract — {title}")
print("=" * W)
print(f" Used {res.n_data} configuration(s) | fit = E0 + "
f"{len(res.kept)} coupling(s)"
+ (f" (S = {spin:g})" if has_S else " (no --spin; reporting raw JS^p)"))
print()
print("Fitted parameters" + (f" (value ± {ksig:g}σ)" if have_se else ""))
print("-" * W)
print(f" E0 = {res.e0:16.6f} eV"
+ (f" ± {ksig*res.se_e0*1000:.4f} meV" if res.se_e0 is not None else ""))
for lbl in labels:
if lbl in res.js2:
p = spin_power(lbl)
js2_meV = res.js2[lbl] * 1000.0
se = res.se.get(lbl)
raw_bar = f" ± {ksig*se*1000:.4f}" if se is not None else ""
if has_S:
J_meV = js2_meV / (S ** p)
J_bar = f" ± {ksig*(se / (S ** p))*1000:.4f}" if se is not None else ""
print(f" {lbl:<5} = {J_meV:14.4f}{J_bar} meV "
f"(raw fit {lbl}*S^{p} = {js2_meV:.4f}{raw_bar} meV)")
else:
print(f" {lbl:<5} = {js2_meV:14.4f}{raw_bar} meV "
f"(this is {lbl}*S^{p}; pass --spin to get the bare J)")
else:
print(f" {lbl:<5} = indeterminate ({res.dropped.get(lbl, '')})")
print("-" * W)
print(f" RMS error = {res.rms_eV*1000:11.4f} meV ({res.rms_eV:.4e} eV)")
print(f" RMS (per DOF, n-p) = {res.rms_dof_eV*1000:11.4f} meV")
bad_vif = [(l, v) for l, v in res.vif.items() if v >= 10.0]
if bad_vif:
print()
print(" ⚠ near-degenerate couplings (VIF ≥ 10) — their split is "
"unreliable (a too-small supercell aliases the shells):")
for l, v in bad_vif:
print(f" VIF({l}) = {'inf' if not np.isfinite(v) else f'{v:.1f}'}")
if res.max_corr is not None and res.max_corr[2] > 0.7:
a, b, c = res.max_corr
print(f" max |corr| between couplings = {c:.2f} ({a}↔{b}) — "
f"{'strongly' if c > 0.9 else 'moderately'} correlated.")
def _print_per_config_fit_table(records: list, energies: dict, mags: dict,
res: FitResult, excluded: set) -> None:
"""Per-configuration table: DFT energy vs fitted-model energy + residual."""
has_mag = any(v is not None for v in mags.values())
print()
print("Per-configuration fit")
hdr = (f" {'cfg':>4} {'deg':>5} {'Sz':>4} "
+ (f"{'mag':>9} " if has_mag else "")
+ f"{'DFT (eV)':>16} {'model (eV)':>16} {'resid (meV)':>12} status")
print(hdr)
print(" " + "-" * (len(hdr) - 2))
for r in sorted(records, key=lambda x: x["number"]):
c = r["number"]
if c not in energies:
continue # missing / unconverged — already warned
dft = energies[c]
model = res.e0 + sum(res.js2[l] * r["coefficients"].get(l, 0)
for l in res.kept)
mg = mags.get(c)
mg_s = f"{mg:9.3f} " if mg is not None else f"{'—':>9} "
status = "EXCLUDED" if c in excluded else "fit"
print(f" {c:>4} {r['degeneracy']:>5} {r['sz']:>4} "
+ (mg_s if has_mag else "")
+ f"{dft:>16.6f} {model:>16.6f} {(model - dft)*1000:>12.3f} {status}")
print("=" * 80)
def _write_fit_dat(path: str, res: FitResult, labels: List[str], dist_of: dict,
spin: Optional[float], compound: str, code: str) -> None:
has_S = spin is not None and spin > 0
S = spin if has_S else 1.0
with open(path, "w") as f:
f.write("# mag4-extract — energy-mapping J values\n")
f.write(f"# compound: {compound} code: {code}"
+ (f" spin: {spin:g}" if has_S else "") + "\n")
f.write(f"# E0 = {res.e0:.6f} eV RMS = {res.rms_eV*1000:.4f} meV\n")
# sig(meV) = 1σ standard error of the fitted JS² — downstream tools
# (mag4-magnon) use it to decide whether a coupling is resolved from
# zero (|JS²| < 2σ ⇒ statistically zero).
cols = ("# name d(A) JS2(meV) sig(meV)"
+ (" J(meV) J(K)" if has_S else "") + " status\n")
f.write(cols)
for lbl in labels:
d = dist_of.get(lbl, float("nan"))
if lbl in res.js2:
js2_meV = res.js2[lbl] * 1000.0
se = res.se.get(lbl)
se_s = f"{se*1000.0:>11.6f}" if se is not None else f"{'N/A':>11}"
if has_S:
J_meV = js2_meV / (S ** spin_power(lbl))
f.write(f"{lbl:<6} {d:>10.4f} {js2_meV:>14.6f} {se_s} "
f"{J_meV:>12.6f} {J_meV/KB_MEV:>10.3f} OK\n")
else:
f.write(f"{lbl:<6} {d:>10.4f} {js2_meV:>14.6f} {se_s} OK\n")
else:
f.write(f"{lbl:<6} {d:>10.4f} {'N/A':>14} indeterminate\n")
def _check_mag(args, records: list, energies: dict, mags: dict) -> set:
"""``--check-mag``: return the set of config numbers whose converged moment
disagrees with their Sz (and print the per-site-moment diagnostics). Does
not mutate ``energies``."""
ratios = [abs(mags[r["number"]]) / r["sz"] for r in records
if r["number"] in energies and r["sz"] != 0
and mags.get(r["number"]) is not None]
mu_median = float(np.median(ratios)) if ratios else None
mu_spin = 2.0 * args.spin if (args.spin and args.spin > 0) else None
if args.moment is not None:
mu, mu_src = args.moment, "given (--moment)"
elif args.moment_ref == "spin":
mu = mu_spin
mu_src = (f"2·S (S={args.spin:g})" if mu_spin is not None
else "spin (no --spin given)")
else:
mu, mu_src = mu_median, "median |mag|/|Sz| over the data"
if mu is None:
log.warning("--check-mag: cannot infer a per-site moment (no --spin, "
"and no Sz≠0 config with a moment); only Sz=0 singlets are "
"checked.")
mu = 0.0
print(f" --check-mag: per-site moment μ = {mu:.3f} μ_B ({mu_src}), "
f"tolerance = {args.mag_tol:g} μ_B")
if mu_median is not None and mu_spin and mu_median < 0.5 * mu_spin:
log.warning("median per-site moment %.3f μ_B is far below 2·S = %.3f μ_B "
"— the local moments appear COLLAPSED across the whole "
"dataset; the fitted J's may be unphysical%s.",
mu_median, mu_spin,
" (and --moment-ref median hides it — use 'spin')"
if args.moment_ref == "median" else "")
excluded: set = set()
for r in records:
c, sz = r["number"], r["sz"]
if c not in energies:
continue
mg = mags.get(c)
if mg is None:
continue
if sz == 0:
dev, expected = abs(mg), 0.0
else:
dev, expected = abs(abs(mg) / sz - mu), mu * sz
if dev > args.mag_tol:
log.warning("config%d EXCLUDED: |mag| = %.3f μ_B, expected ≈ %.3f "
"(Sz=%d).", c, abs(mg), expected, sz)
excluded.add(c)
if not excluded:
print(" --check-mag: all configurations consistent with their Sz.")
else:
print(f" --check-mag: excluded {len(excluded)} configuration(s) — "
+ ", ".join(f"config{c}" for c in sorted(excluded)) + ".")
return excluded
def _run_ring16_fit(args, compound: str, report_path: str) -> int:
"""``mag4-extract --method 16-state``: parse the ring report's χ-weighted
formula, read each config energy, and apply the projection to get Jring·S⁴.
All configs share one supercell so no per-config rescaling is needed.
When the report carries the PAIRWISE-J ANALYSIS design rows, the same
energies are additionally least-squares fitted to
``E = E0 + Σ_k c_k·(J_k·S²) + c_ring·(Jring·S⁴)`` — the extractable
pairwise couplings come for free, and the fitted ``Jring·S⁴`` cross-checks
the exact χ-projection."""
try:
formula = parse_ring_report(report_path)
except FileNotFoundError:
log.error("Report not found: %s", report_path)
log.error("Hint: run 'mag4-magnetic --method 16-state' here first.")
return 2
except ValueError as exc:
log.error("%s", exc)
return 2
design = parse_ring_design(report_path) # None on old reports
max_idx = max(formula.config_ids, default=0)
if design:
max_idx = max([max_idx] + [r["number"] for r in design["records"]])
scf = args.scffile if args.code == "wien2k" else None
energies = read_all_energies(compound, args.code, max_index=max_idx,
scf_basename=scf)
if not energies:
log.error("No configN/ subdirectory found under %s", compound)
return 2
n_failed = _print_convergence_summary(energies, args.code)
if not args.no_gap:
_print_band_gaps(args, compound, energies.keys(),
args.scffile if args.code == "wien2k" else None)
_print_tr_pairs(parse_report_tr_pairs(report_path), energies)
jring_eV, problems = compute_jring(formula, energies)
W = 64
print()
print("=" * W)
print(" mag4-extract — four-spin ring exchange (16-state method)")
print("=" * W)
if jring_eV is None:
print(" Jring : N/A — " + "; ".join(problems))
print("=" * W)
return 1 if n_failed == 0 else 3
print(f" Jring·S⁴ = {jring_eV*1e3:>12.4f} meV ({jring_eV:.6f} eV)")
if args.spin: # divide out S⁴ when --spin given
s4 = float(args.spin) ** spin_power("Jring")
if s4 > 0:
jring = jring_eV / s4
print(f" Jring = {jring*1e3:>12.4f} meV "
f"(divided by S⁴ = {s4:g}, S = {args.spin})")
if formula.edge is not None:
print(f" square plaquette edge = {formula.edge:.4f} Å, "
f"divisor = {formula.divisor} (= 16·M)")
# Sensitivity: Jring·S⁴ is a χ-weighted near-cancellation — a small residual
# of the (much larger) spread of config energies. Report the ratio so the
# user sees how tightly the DFT must be converged; a few-meV energy shift at
# ratio ≪ 1 % moves Jring by tens of %.
evs = [getattr(energies.get(cid), "energy_eV", None)
for cid in formula.config_ids]
evs = [e for e in evs if e is not None]
if len(evs) >= 2:
spread_meV = (max(evs) - min(evs)) * 1e3
if spread_meV > 1e-6:
ratio = abs(jring_eV * 1e3) / spread_meV
print(f" sensitivity: |Jring·S⁴| is {ratio*100:.2f}% of the "
f"{spread_meV:.0f} meV config-energy spread")
if ratio < 0.02:
print(" ⚠ Jring is a small residual of near-cancelling "
"energies — it demands tight,")
print(" consistent DFT convergence (same k-mesh, smearing, "
"U, precision across ALL")
print(" configs); a few-meV per-config error moves Jring by "
"tens of percent.")
print("=" * W)
# ---- pairwise-J fit from the same energies (report analysis section) ----
if design:
labels, records = design["labels"], design["records"]
verdicts = analyze_pairwise_extractability(records, labels)
keep = [l for l in labels if verdicts[l][0]]
skipped = [l for l in labels if not verdicts[l][0]]
energies_all = {idx: er.energy_eV for idx, er in energies.items()
if er.converged and er.energy_eV is not None}
for lbl in skipped:
print(f" {lbl}: not fitted — {verdicts[lbl][1]}.")
try:
res = fit_couplings(records, keep, energies_all,
weight_by_deg=args.weight_by_deg)
except ValueError as exc:
log.warning("pairwise-J fit skipped: %s", exc)
return 0 if n_failed == 0 else 3
_print_fit_table(res, keep, args.spin, args.n_sigma,
title="pairwise Js from the 16-state configs "
"(least-squares)")
# The χ-projection Jring is EXACT (parity orthogonality annihilates every
# ≤3-spin and longer-range term); the least-squares pairwise Js are NOT —
# they assume no couplings beyond the shells listed here, and any
# unmodelled shell within the flip's reach leaks into them. The RMS
# residual and the projection-vs-fit gap below are the health flags.
if "Jring" in res.js2 and jring_eV is not None:
diff = (res.js2["Jring"] - jring_eV) * 1e3
print(f" cross-check: fitted Jring·S⁴ = "
f"{res.js2['Jring']*1e3:.4f} meV vs χ-projection "
f"{jring_eV*1e3:.4f} meV (Δ = {diff:+.4f} meV)")
rel = abs(diff) / max(abs(jring_eV) * 1e3, 1e-9)
if rel > 0.05:
print(" ⚠ fit and projection disagree by >5% — the config "
"energies carry couplings NOT in this model")
print(" (shells beyond the ring cutoff, 3-spin terms, …). "
"TRUST the χ-projection Jring above;")
print(" treat the fitted pairwise Js as contaminated — "
"re-extract them with --method four-state")
print(" or a larger --cutoff / supercell.")
if res.rms_eV * 1e3 > 0.05:
print(f" note: fit RMS = {res.rms_eV*1e3:.3f} meV > 0 ⇒ the 3-term "
"model is incomplete for these energies;")
print(" the fitted pairwise Js are a least-squares compromise "
"(the χ-projection Jring is unaffected).")
print(" (pairwise Js fitted from the SAME 16-state configurations — "
"no extra DFT; Jring itself comes from the exact χ-projection "
"above)")
return 0 if n_failed == 0 else 3
def _select_couplings(labels: List[str], args) -> Optional[List[str]]:
"""Resolve which coupling columns to fit from ``--couplings`` /
``--drop-couplings`` (mutually exclusive). Returns the kept labels in the
report's order, or ``None`` on a user error (already logged)."""
inc, exc = args.couplings, args.drop_couplings
if inc and exc:
log.error("--couplings and --drop-couplings are mutually exclusive.")
return None
if not inc and not exc:
return list(labels)
def _parse(s):
return [t.strip() for t in s.split(",") if t.strip()]
known = set(labels)
if inc:
want = set(_parse(inc))
unknown = sorted(want - known)
if unknown:
log.error("--couplings: unknown coupling(s) %s; report has [%s].",
", ".join(unknown), ", ".join(labels))
return None
kept = [l for l in labels if l in want]
else:
drop = set(_parse(exc))
unknown = sorted(drop - known)
if unknown:
log.error("--drop-couplings: unknown coupling(s) %s; report has [%s].",
", ".join(unknown), ", ".join(labels))
return None
kept = [l for l in labels if l not in drop]
if not kept:
log.error("No couplings left to fit after the --couplings/"
"--drop-couplings selection.")
return None
return kept
def _run_energy_mapping_fit(args, compound: str, report_path: str) -> int:
"""``mag4-extract --method energy-mapping``: parse the design-matrix report,
read every config energy (VASP TOTEN / WIEN2k :ENE), and fit JS² by linear
least-squares."""
try:
rep = parse_energy_mapping_report(report_path)
except FileNotFoundError:
log.error("Report not found: %s", report_path)
log.error("Hint: run 'mag4-magnetic --method energy-mapping' here first.")
return 2
except ValueError as exc:
log.error("%s", exc)
return 2
labels, records, dist_of = rep["labels"], rep["records"], rep["dist_of"]
log.info("Energy-mapping report: %d order(s), couplings [%s], "
"supercell %dx%dx%d", len(records), ", ".join(labels), *rep["supercell"])
# Optional model restriction: refit with only a subset of couplings (the
# rest forced to zero). Re-solving with fewer design-matrix columns lets the
# user compare models — e.g. with vs without Jring, or without J3 and Jring.
kept = _select_couplings(labels, args)
if kept is None:
return 2
if kept != labels:
dropped = [l for l in labels if l not in kept]
log.info("Restricting the fit to [%s]; forcing [%s] to zero "
"(their contribution moves into the residual / E0).",
", ".join(kept), ", ".join(dropped))
labels = kept
max_idx = max((r["number"] for r in records), default=0)
energies_res = read_all_energies(
compound, args.code, max_index=max_idx,
scf_basename=args.scffile if args.code == "wien2k" else None)
if not energies_res:
log.error("No configN/ subdirectory found under %s", compound)
return 2
_print_convergence_summary(energies_res, args.code)
if not args.no_gap:
_print_band_gaps(args, compound, energies_res.keys(),
args.scffile if args.code == "wien2k" else None)
_print_tr_pairs(parse_report_tr_pairs(report_path), energies_res)
energies_all = {idx: er.energy_eV for idx, er in energies_res.items()
if er.converged and er.energy_eV is not None}
missing = sorted(r["number"] for r in records if r["number"] not in energies_all)
if missing:
log.warning("%d configuration(s) excluded from the fit (missing or "
"unconverged): %s", len(missing),
", ".join(f"config{i}" for i in missing))
# Total moment per config — used by the per-config table and --check-mag.
scf = args.scffile if args.code == "wien2k" else None
mags = {r["number"]: read_magnetization(
os.path.join(compound, f"config{r['number']}"),
args.code, scf_basename=scf)
for r in records if r["number"] in energies_all}
mag_excluded: set = set()
if args.check_mag:
mag_excluded = _check_mag(args, records, energies_all, mags)
user_excluded: set = set()
if args.exclude:
known = {r["number"] for r in records}
unknown = sorted(set(args.exclude) - known)
if unknown:
log.warning("--exclude: no such configuration(s): %s (report has "
"config1..config%d).",
", ".join(f"config{i}" for i in unknown), max_idx)
user_excluded = {n for n in args.exclude if n in known}
present = sorted(n for n in user_excluded if n in energies_all)
if present:
log.info("--exclude: dropping %d configuration(s) from the fit: %s",
len(present), ", ".join(f"config{i}" for i in present))
excluded = mag_excluded | user_excluded
fit_energies = {k: v for k, v in energies_all.items() if k not in excluded}
try:
res = fit_couplings(records, labels, fit_energies,
weight_by_deg=args.weight_by_deg,
keep_collinear=args.keep_collinear)
except ValueError as exc:
log.error("%s", exc)
return 1
_print_fit_table(res, labels, args.spin, args.n_sigma)
_print_per_config_fit_table(records, energies_all, mags, res, excluded)
out_path = args.output or os.path.join(compound, "j_values.dat")
_write_fit_dat(out_path, res, labels, dist_of, args.spin, compound, args.code)
print(f"\n J values written to: {out_path}")
if res.dropped:
print(f"\n ⚠ {len(res.dropped)} coupling(s) indeterminate in this cell "
f"(see notes above) — enlarge the supercell to resolve them.")
return 1
print(f"\n All {len(labels)} coupling(s) fitted.")
return 0
# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------
class _Tee:
"""Duplicate writes to several streams (console + the persisted log)."""
def __init__(self, *streams):
self._streams = streams
def write(self, data):
for s in self._streams:
s.write(data)
def flush(self):
for s in self._streams:
s.flush()
[docs]
def main(argv=None) -> int:
args = build_parser().parse_args(argv)
compound = os.path.abspath(args.compound)
if not os.path.isdir(compound):
log.error("compound directory not found: %s", compound)
return 2
report_path = args.report or os.path.join(compound, DEFAULT_REPORT_NAME)
# Persist the complete screen output — it carries the convergence,
# gap, moment and degeneracy diagnostics that j_values.dat does not.
tee_path = os.path.join(compound, "mag4_extract.out")
try:
_tee_file = open(tee_path, "w")
except OSError as exc:
log.warning("Cannot write %s (%s); output not persisted.",
tee_path, exc)
_tee_file = None
_old_stdout = sys.stdout
_file_handler = None
if _tee_file is not None:
sys.stdout = _Tee(_old_stdout, _tee_file)
_file_handler = logging.StreamHandler(_tee_file)
_file_handler.setFormatter(logging.Formatter("%(levelname)s: %(message)s"))
logging.getLogger().addHandler(_file_handler)
try:
# Dispatch on the method recorded in the report (or forced via --method).
method = args.method
if method == "auto":
method = detect_report_method(report_path)
if (args.couplings or args.drop_couplings) and method != "energy-mapping":
log.warning("--couplings/--drop-couplings only apply to the "
"energy-mapping fit (method is '%s'); ignoring.", method)
if method == "energy-mapping":
return _run_energy_mapping_fit(args, compound, report_path)
if method == "16-state":
return _run_ring16_fit(args, compound, report_path)
# 1. parse the report
try:
formulas = parse_report_formulas(report_path)
except FileNotFoundError as exc:
log.error("%s", exc)
log.error("Hint: did you run 'mag4-magnetic' in this directory?")
return 2
except ValueError as exc:
log.error("%s", exc)
return 2
log.info("Parsed %d formula(s) from %s", len(formulas), report_path)
# 2. read every config energy
max_idx = max((max(f.config_indices) for f in formulas if not f.incomplete),
default=0)
# --check-degeneracy partner configs sit above the formulas' indices —
# include them so their energies are read for the |ΔE| check.
tr_pairs = parse_report_tr_pairs(report_path)
if tr_pairs:
max_idx = max([max_idx] + [i for pair in tr_pairs for i in pair])
if args.scffile and args.code != "wien2k":
log.warning("--scffile is only used with --code wien2k (ignored).")
energies = read_all_energies(
compound, args.code, max_index=max_idx,
scf_basename=args.scffile if args.code == "wien2k" else None,
)
if not energies:
log.error("No configN/ subdirectory found under %s", compound)
return 2
# Supercell, per-config Sz (from the report) and converged moments — for the
# per-configuration table and the moment-vs-Sz consistency check.
supercell = parse_report_supercell(report_path)
config_sz = parse_report_config_sz(report_path)
scf = args.scffile if args.code == "wien2k" else None
mags = {idx: read_magnetization(os.path.join(compound, f"config{idx}"),
args.code, scf_basename=scf)
for idx in energies}
n_failed = _print_convergence_summary(energies, args.code)
_print_config_energies(energies, config_sz, mags, args.spin)
_print_tr_pairs(tr_pairs, energies)
if not args.no_gap:
_print_band_gaps(args, compound, energies.keys(),
args.scffile if args.code == "wien2k" else None)
# 3. apply formulas
# Per-config atom counts (when mag4-magnetic wrote per-config reduced
# cells) so compute_j_values can rescale energies to a common cell.
report_n_atoms = parse_report_config_cells(report_path)
# Override with the **actual** atom count from each configN.struct
# whenever it disagrees with the report. This catches the case
# where mag4's pre-emptive cell reduction was conservative and
# WIEN2k's sgroup (auto-accepted by job.init) reduced further --
# e.g. CrSb config7 where mag4 wrote 8 atoms but sgroup ran
# SCF on 4. The struct file on disk is the authoritative source.
config_n_atoms: Dict[int, int] = dict(report_n_atoms)
if args.code == "wien2k":
struct_n_atoms = read_struct_atom_counts(compound, len(energies))
for cfg_idx, n in struct_n_atoms.items():
old = config_n_atoms.get(cfg_idx)
if old is not None and old != n:
log.info("config%d: struct file has %d atoms but report "
"says %d -- using struct (post-sgroup reduction).",
cfg_idx, n, old)
config_n_atoms[cfg_idx] = n
# Total supercell atom count — the reference cell the four-state
# formula is defined against. Each config's WIEN2k :ENE gets
# rescaled UP to this count via DFT extensivity.
n_supercell_atoms = parse_report_supercell_atoms(report_path)
if config_n_atoms:
n_uniq = len(set(config_n_atoms.values()))
if n_uniq > 1:
if n_supercell_atoms is not None:
log.info("Per-config cell metadata found: %d configs span %d "
"distinct atom counts %s — energies will be rescaled "
"to the %d-atom supercell.",
len(config_n_atoms), n_uniq,
sorted(set(config_n_atoms.values())),
n_supercell_atoms)
else:
log.info("Per-config cell metadata found: %d configs span %d "
"distinct atom counts %s — energies will be rescaled "
"to per-J LCM cells (no supercell-atoms line in "
"report).",
len(config_n_atoms), n_uniq,
sorted(set(config_n_atoms.values())))
results = compute_j_values(formulas, energies, spin=args.spin,
config_n_atoms=config_n_atoms or None,
n_supercell_atoms=n_supercell_atoms)
# 4. print + persist
print_j_table(results, spin=args.spin)
out_path = args.output or os.path.join(compound, "j_values.dat")
write_j_values_dat(results, out_path, compound_dir=compound,
code=args.code, spin=args.spin, supercell=supercell,
energies=energies, config_sz=config_sz, mags=mags)
print(f"\n J values written to: {out_path}")
# 5. exit code: 1 only if at least one J could not be computed
n_unrecoverable = sum(1 for r in results if not r.computed)
if n_unrecoverable:
print(f"\n ⚠ {n_unrecoverable} J shell(s) could not be computed "
f"(see 'Status' column above).")
return 1
print()
print(f" All {len(results)} J shell(s) computed successfully.")
return 0
finally:
if _file_handler is not None:
logging.getLogger().removeHandler(_file_handler)
sys.stdout = _old_stdout
if _tee_file is not None:
_tee_file.close()
print(f" Full output saved to: {tee_path}")
if __name__ == "__main__": # pragma: no cover
sys.exit(main())