"""Extract converged magnetic moments from WIEN2k ``.scf`` files.
Layer: **WIEN2k-specific** (backs the ``mag4-moments`` CLI).
"""
from __future__ import annotations
import glob
import math
import os
import re
from typing import Dict, Iterable, List, Optional, Tuple
#: report line written at generation time, e.g.
#: " MAGMOM (S = 0.5, |MAGMOM| = 2·S = 1 μ_B)"
_TARGET_RE = re.compile(r"2[·*]S\s*=\s*([\d.]+)")
_SPIN_RE = re.compile(r"\bS\s*=\s*([\d.]+)")
[docs]
def parse_report_target_moment(report_path: str) -> Optional[float]:
"""Read the intended per-site ``|MAGMOM| = 2·S`` from a mag4 report.
Returns the target moment (μ_B) or ``None`` if the report has no
``MAGMOM (S = …, |MAGMOM| = 2·S = …)`` header (older/other reports).
"""
if not os.path.isfile(report_path):
return None
with open(report_path, encoding="utf-8", errors="replace") as fh:
for line in fh:
if "MAGMOM" not in line:
continue
m = _TARGET_RE.search(line)
if m:
return float(m.group(1))
m = _SPIN_RE.search(line)
if m:
return 2.0 * float(m.group(1))
return None
[docs]
def find_report(base_dir: str = ".") -> Optional[str]:
"""Locate a ``*_report.dat`` in ``base_dir`` (prefers ``four_state_report``)."""
cand = os.path.join(base_dir, "four_state_report.dat")
if os.path.isfile(cand):
return cand
hits = sorted(glob.glob(os.path.join(base_dir, "*_report.dat")))
return hits[0] if hits else None
def _resolve_scf(case_dir: str, scf_basename: Optional[str] = None) -> str:
"""Path of the ``.scf`` file to read inside ``case_dir``.
Default (``scf_basename is None``) follows the WIEN2k convention
``<dirname>/<dirname>.scf``. A ``scf_basename`` (with or without the
``.scf`` suffix) overrides it — ``configN/<basename>.scf`` — matching
``mag4-extract --scffile`` so runs saved side by side as e.g.
``pbe.scf`` / ``r2scan.scf`` can be selected.
"""
if scf_basename:
base = (scf_basename[:-4] if scf_basename.lower().endswith(".scf")
else scf_basename)
else:
base = os.path.basename(case_dir.rstrip("/\\"))
return os.path.join(case_dir, f"{base}.scf")
[docs]
def find_cases(base_dir: str = ".",
scf_basename: Optional[str] = None,
code: str = "wien2k") -> List[str]:
"""Auto-detect config subdirectories with results for ``code`` (sorted).
``code='wien2k'`` → directories with a ``.scf`` file (or exactly
``<scf_basename>.scf`` when given); ``code='vasp'`` → directories with an
``OUTCAR``.
"""
out = []
for d in sorted(glob.glob(os.path.join(base_dir, "*"))):
if not os.path.isdir(d):
continue
if code.lower() == "vasp":
if os.path.isfile(os.path.join(d, "OUTCAR")):
out.append(d)
elif scf_basename:
if os.path.isfile(_resolve_scf(d, scf_basename)):
out.append(d)
elif glob.glob(os.path.join(d, "*.scf")):
out.append(d)
return out
[docs]
def build_table(cases: List[str], atom_indices: List[int],
scf_basename: Optional[str] = None,
target: Optional[float] = None,
tol: float = 0.5,
code: str = "wien2k") -> Tuple[List[str], dict]:
"""Textual table of moments for each case + statistics + a target check.
``code`` selects the parser: ``'wien2k'`` reads per-sphere ``:MMI:`` values
from ``<dir>/<scf_basename>.scf``; ``'vasp'`` reads the per-ion ``tot``
column of the last ``magnetization (x)`` block in ``<dir>/OUTCAR`` (needs
``LORBIT=10/11``). For VASP the indices are 1-based POSCAR ion numbers
(mag4 writes the magnetic atoms first, so they are ``1 … n_mag``).
Returns ``(lines, summary)`` where ``summary`` carries
``n_data`` (configs with a readable SCF), ``n_missing`` (SCF not found),
``failures`` (list of ``(case, reason)`` where a site did not reach the
``target`` moment within ``tol``, or collapsed to ~0), and ``target``.
``target`` is the expected per-site ``|μ|`` (e.g. ``2·S`` = the generated
``|MAGMOM|``). When ``None`` the check only flags collapsed moments
(``|μ| < 0.1 μ_B``), which usually mean the config fell into a
non-magnetic solution rather than the intended spin state.
"""
col_w = 10
lbl_w = 18
total_w = lbl_w + col_w * (len(atom_indices) + 1)
sep, dash = "=" * total_w, "-" * total_w
src_label = ("VASP OUTCAR (per-ion)" if code.lower() == "vasp"
else "WIEN2k SCF (per-sphere)")
lines = ["", sep, f" Magnetic moments (muB) -- {src_label}", sep]
hdr = f"{'Case':<{lbl_w}}"
for i in atom_indices:
hdr += f"{'At.' + str(i):>{col_w}}"
hdr += f"{'M_tot':>{col_w}}"
lines += [hdr, dash]
from mag4.energy import (read_site_moments, read_vasp_magnetization,
read_wien2k_magnetization)
is_vasp = code.lower() == "vasp"
all_data = []
n_missing = 0
failures: List[Tuple[str, str]] = []
for case_dir in cases:
case_name = os.path.basename(case_dir.rstrip("/\\"))
src = (os.path.join(case_dir, "OUTCAR") if is_vasp
else _resolve_scf(case_dir, scf_basename))
if not os.path.isfile(src):
print(f" Warning: {src} not found — skipping.")
n_missing += 1
continue
moments = read_site_moments(case_dir, code, set(atom_indices),
scf_basename)
m_tot = (read_vasp_magnetization(case_dir) if is_vasp
else read_wien2k_magnetization(case_dir, scf_basename))
row = f"{case_name:<{lbl_w}}"
for idx in atom_indices:
val = moments.get(idx)
row += f"{val:>{col_w}.4f}" if val is not None else f"{'N/A':>{col_w}}"
row += f"{m_tot:>{col_w}.4f}" if m_tot is not None else f"{'N/A':>{col_w}}"
lines.append(row)
all_data.append((case_name, moments, m_tot))
lines.append(sep)
# ---- target / consistency checks (run once all configs are collected) ----
# 1. Missing :MMI: line for a requested site -> SCF incomplete/wrong index.
# 2. Consistency: |μ| at a magnetic site must be ~constant across configs
# (only its SIGN flips between configs), so a config whose |μ| deviates
# from that site's cross-config median by > tol is the broken one. This
# is covalency-proof — it uses the ensemble as its own reference, not the
# formal 2·S (which the converged moment always undershoots in oxides).
# 3. Collapse floor (only when a target is known): catches the case where
# EVERY config collapsed (consistency alone would call them "consistent").
# A site is collapsed if |μ| < 0.25·target.
per_case_missing = {}
for name, moments, _ in all_data:
miss = [i for i in atom_indices if moments.get(i) is None]
if miss:
per_case_missing[name] = miss
# site medians over the configs that reported a value
medians = {}
for idx in atom_indices:
vals = sorted(abs(m[idx]) for _n, m, _t in all_data if m.get(idx) is not None)
if vals:
medians[idx] = vals[len(vals) // 2]
floor = 0.25 * target if target is not None else None
for name, moments, _ in all_data:
if name in per_case_missing:
why = ("OUTCAR has no per-ion moment block (set LORBIT=11) or "
"wrong ion index" if is_vasp
else "SCF incomplete or wrong atom index")
failures.append((name, f"no moment parsed for site(s) "
f"{per_case_missing[name]} — {why}"))
continue
dev = [i for i in atom_indices
if i in medians and abs(abs(moments[i]) - medians[i]) > tol]
low = ([i for i in atom_indices if abs(moments[i]) < floor]
if floor is not None else [])
bad = sorted(set(dev) | set(low))
if bad:
worst = max(bad, key=lambda i: abs(abs(moments[i]) - medians.get(i, 0)))
failures.append(
(name, f"|μ| inconsistent at site(s) {bad} "
f"(e.g. site {worst}: |μ|={abs(moments[worst]):.3f} vs "
f"ensemble {medians.get(worst, float('nan')):.3f})"))
if all_data:
lines += ["", " Statistics for |moment| per atom site across all configurations:", dash]
mean_row = f"{'Mean <|mu|>':<{lbl_w}}"
std_row = f"{'Std. Dev.':<{lbl_w}}"
for idx in atom_indices:
vals = [abs(d[1][idx]) for d in all_data if d[1].get(idx) is not None]
if vals:
mean = sum(vals) / len(vals)
variance = sum((x - mean) ** 2 for x in vals) / len(vals)
std = math.sqrt(variance)
mean_row += f"{mean:>{col_w}.4f}"
std_row += f"{std:>{col_w}.4f}"
else:
mean_row += f"{'N/A':>{col_w}}"
std_row += f"{'N/A':>{col_w}}"
lines += [mean_row, std_row, dash]
# ---- verdict -------------------------------------------------------------
lines += ["", " Magnetic-state check", dash]
basis = (f"|μ| consistent across configs within {tol:g} μ_B"
+ (f", none below the collapse floor 0.25·2S = {0.25*target:.3f}"
if target is not None else ""))
if not all_data:
lines.append(" ✗ no SCF file was read — nothing to check "
"(see warnings above).")
elif failures:
lines.append(f" ✗ {len(failures)} of {len(all_data)} config(s) did NOT "
"reach a consistent magnetic state:")
for case, why in failures:
lines.append(f" {case}: {why}")
lines.append(" → their total energies are unreliable for J/Jring "
"extraction; re-run these configs")
lines.append(" (tighter convergence / better MAGMOM seed) before "
"trusting the couplings.")
else:
lines.append(f" ✓ all {len(all_data)} config(s) OK ({basis}).")
lines.append(dash)
summary = {"n_data": len(all_data), "n_missing": n_missing,
"failures": failures, "target": target}
return lines, summary