"""CLI: full four-state pipeline (CIF or POSCAR → spin configurations + VASP/WIEN2k inputs)."""
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 DIST_EQUIV_TOL, FourStateConfig
from mag4.dimers import check_coupling_spacing, find_dimers
from mag4.energy_mapping import (
_cartesian_with_spin,
_image_shell,
build_unique_groups,
config_energy_coefficients,
count_ring,
enumerate_orders,
select_ring_plaquettes,
write_energy_mapping_report,
)
from mag4.io_cif import get_mm_distances, group_inequivalent
from mag4.io_poscar import read_poscar
from mag4.lattice import get_cartesian_coords
from mag4.pipeline import run_cif_analysis
from mag4.reports import print_report
from mag4.ring16 import (
find_ring_supercell,
generate_ring_16,
ring_design_records,
select_ring_quad,
write_ring_report,
)
from mag4.spinconfig import (
build_reference_bath,
build_reference_from_couplings,
build_reference_from_kvector,
build_reference_from_magmom,
find_unique_configs,
time_reversal_pairs,
generate_four_states,
parse_kvector,
)
from mag4.supercell import (
check_four_state_isolation,
expand_supercell,
find_optimal_supercell,
find_supercell_for_enumeration,
shortest_periodic_image,
)
from mag4.symmetry import find_symmetry_ops
from mag4.vasp import build_magmom_lines, create_config_dirs
from mag4.wien2k import (
write_case_inst_files,
write_case_struct_files,
write_job_scripts,
write_symmetric_case_files,
)
logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")
log = logging.getLogger(__name__)
[docs]
def parse_couplings(s: str) -> list:
"""Parse ``"J1:3.94,J2:6.10"`` → ``[["J1", 3.94], ["J2", 6.10]]``."""
result = []
for item in s.split(","):
name, dist = item.strip().split(":")
result.append([name.strip(), float(dist.strip())])
return result
[docs]
def parse_ref_state(s: str) -> List[int]:
"""Parse ``"1,-1,1"`` → ``[1, -1, 1]`` with ±1 validation."""
values = []
for tok in s.split(","):
v = int(tok.strip())
if v not in (+1, -1):
raise argparse.ArgumentTypeError(
f"ref-state values must be +1 or -1, got '{tok}'"
)
values.append(v)
return values
[docs]
def build_parser() -> argparse.ArgumentParser:
parser = argparse.ArgumentParser(
prog="mag4-magnetic",
description="CIF/POSCAR → J extraction + four-state configuration analysis.",
formatter_class=argparse.RawDescriptionHelpFormatter,
epilog="""\
The optimal supercell is defined FOR A GIVEN --cutoff distance between
the magnetic centres: only J shells with d <= --cutoff are extracted and
the auto search picks the smallest cell that gives each of them a
distinct intra-cell pair. Changing --cutoff changes the set of J shells
and therefore the optimal supercell.
Recommended workflow:
1. Pick a trial --cutoff and run with --dry-run to inspect the J shells
and the supercell mag4 would build (no files written).
2. Adjust --cutoff until the J set matches what you want to extract.
3. Re-run the same command WITHOUT --dry-run to write the configN/ tree.
Mode A – CIF (recommended, no POSCAR needed):
# Step 1: dry-run at a tentative cutoff → see J set + optimal supercell:
mag4-magnetic --cif CsVI3-lowT.cif --element V --cutoff 7.0 --dry-run
# Step 2: same cutoff, write inputs:
mag4-magnetic --cif CsVI3-lowT.cif --element V --cutoff 7.0
# Override the auto search with an explicit supercell:
mag4-magnetic --cif CsVI3-lowT.cif --element V --supercell 1 1 2
Mode B – POSCAR + manual couplings:
mag4-magnetic --poscar POSCAR --species 1 \\
--couplings "J1:3.94,J2:5.88"
Reference bath options (mutually exclusive):
(default) # fully FM (all +1)
--ref-couplings "J1:AFM,J2:AFM" # AFM bath respecting chosen J signs
--ref-couplings "J1:AFM,J2:FM" # mixed: J1 antiparallel, J2 parallel
--ref-kvector "1/2 1/2 0" # bath from a propagation vector
--ref-magmom "32*0 1 1 1 -1 -1 -1 -1 1 -1 -1 -1 1 1 1 1 -1 64*0"
# bath as the SUPERCELL MAGMOM string
""",
)
# Options are organised into argument groups purely for --help readability;
# every flag behaves identically wherever it is registered.
g_in = parser.add_argument_group(
"structure input", "What to read and which atoms are magnetic.")
g_shell = parser.add_argument_group(
"J-shell detection",
"Which coupling shells are extracted (distances, tolerances).")
g_cell = parser.add_argument_group(
"supercell & isolation",
"The working cell and the four-state isolation criterion.")
g_ref = parser.add_argument_group(
"reference bath (mutually exclusive)",
"The fixed spin background the probed atoms are flipped over "
"(default: FM, all +1).")
g_method = parser.add_argument_group(
"method & diagnostics",
"four-state (default) / energy-mapping / 16-state, plus the "
"contamination and degeneracy checks.")
g_run = parser.add_argument_group(
"run setup & output",
"DFT code, spin quantum number, dry-run, output naming.")
g_vasp = parser.add_argument_group(
"VASP settings", "Only used when --code includes vasp.")
g_wien = parser.add_argument_group(
"WIEN2k settings", "Only used when --code includes wien2k.")
g_lig = parser.add_argument_group(
"geometry-aware J shells (--ligand)",
"Split same-distance shells by their M-L-M bridging geometry.")
src = g_in.add_mutually_exclusive_group(required=True)
src.add_argument("--cif", default=None,
help="CIF file — auto-extracts J couplings, no POSCAR needed")
src.add_argument("--poscar", default=None,
help="VASP POSCAR file (requires --couplings and --species)")
g_in.add_argument("--element", default=None,
help="Magnetic element symbol, e.g. Gd (required with --cif)")
g_shell.add_argument("--cutoff", type=float, default=7.0,
help="Magnetic-centre distance cutoff in Å (default: 7.0). "
"Only J shells with d <= --cutoff are extracted, "
"AND the auto-supercell search (see --supercell / "
"--no-auto-supercell) is defined relative to this set "
"of J shells. Raising --cutoff may grow the optimal "
"supercell; lowering it may shrink it.")
g_shell.add_argument("--equiv-tol", type=float, default=DIST_EQUIV_TOL,
help=f"Tolerance to merge equivalent distances in Å "
f"(default: {DIST_EQUIV_TOL})")
g_in.add_argument("--species", type=int, default=None,
help="1-based species index in POSCAR (required with --poscar)")
g_in.add_argument("--couplings", type=str, default=None,
help="Manual couplings 'J1:3.94,J2:6.10,...' (required "
"with --poscar). NOTE: unrelated to mag4-extract's "
"--couplings, which restricts the fit.")
g_cell.add_argument("--supercell", type=int, nargs=3,
metavar=("NA", "NB", "NC"), default=None,
help="Supercell expansion along a, b, c. When omitted, "
"mag4 auto-finds the smallest cell in which every "
"coupling with d <= --cutoff has a distinct "
"intra-cell pair (disable with --no-auto-supercell).")
g_cell.add_argument("--no-auto-supercell", dest="auto_supercell",
action="store_false", default=True,
help="Disable the auto-supercell search and keep the "
"input cell as-is when --supercell is not given "
"(the search is on by default).")
# Deprecated positive form: the search is already the default, so the flag
# is a no-op. Kept working (hidden) for old scripts; main() warns.
parser.add_argument("--auto-supercell", dest="auto_supercell",
action="store_true", help=argparse.SUPPRESS)
g_cell.add_argument("--max-supercell", type=int, default=4, metavar="N",
help="Maximum per-axis multiplier explored by the "
"auto-supercell search (default: 4).")
g_cell.add_argument("--full-isolation", action="store_true",
help="Demand a FULLY isolated dimer (multiplicity 1, "
"textbook divisor 4) for every shell, instead of the "
"default algebraic criterion that also accepts "
"compact cells where the dimer wraps onto its own "
"image (divisor 4·M). Yields a larger but "
"noise-safer cell — e.g. La₂CuO₄ 2×2×1 → 3×3×1, "
"CsVI₃ 1×1×2 → 1×1×3 — equivalent to the old "
"geometric '3 copies of d' rule. Shorthand for "
"--max-multiplicity 1.")
g_cell.add_argument("--max-multiplicity", type=int, default=None, metavar="M",
help="Cap the per-shell dimer multiplicity accepted by "
"the auto-supercell search (and enforced against an "
"explicit --supercell). The default algebraic "
"criterion allows any M (divisor 4·M); M=1 is full "
"isolation; intermediate values (e.g. 2) trade cell "
"size against the largest divisor. Must be ≥ 1.")
g_run.add_argument("--dry-run", action="store_true",
help="Analyse only: print couplings, chosen supercell, "
"report, and MAGMOM lines to stdout — do NOT create "
"the compound directory, four_state_report.dat, or "
"any per-config directories. Use this to tune "
"--cutoff before committing to a particular J set "
"and supercell.")
g_shell.add_argument("--tol", type=float, default=0.05,
help="Distance tolerance for dimer search in Å (default: 0.05)")
# Deprecated: superseded by --ref-magmom (same spins, MAGMOM-style, n*v
# repeats, plus the validated full-supercell form). Hidden; main() warns.
parser.add_argument("--ref-state", type=str, default=None,
metavar="S1,S2,...", help=argparse.SUPPRESS)
g_ref.add_argument("--ref-couplings", type=str, default=None,
metavar="J1:AFM,J2:FM",
help="Build the reference bath from requested coupling "
"signs (AFM=antiparallel, FM=parallel), e.g. "
"'J1:AFM,J2:AFM'. Uses exact signed 2-colouring, "
"falling back to energy minimisation if the request "
"is frustrated. Mutually exclusive with the other "
"--ref-* options.")
g_ref.add_argument("--ref-kvector", type=str, default=None,
metavar="KX KY KZ",
help="Build the reference bath from a propagation vector "
"in parent reciprocal-lattice units, e.g. "
"'1/2 1/2 0' (s_i = sign(cos 2π k·R_i)). Warns and "
"suggests a supercell if the cell is incommensurate "
"with k. Mutually exclusive with the other --ref-* "
"options.")
g_ref.add_argument("--ref-magmom", type=str, default=None,
dest="ref_magmom", metavar="MAGMOM",
help="Reference bath given directly as the VASP "
"MAGMOM-style string of the WORKING SUPERCELL "
"(POSCAR atom order, n*v repetitions allowed), "
"e.g. '32*0 1 1 1 -1 -1 -1 -1 1 -1 -1 -1 1 1 1 1 "
"-1 64*0' for a 2×2×1 La₂CuO₄ √2 cell. Only the "
"SIGNS of the magnetic entries are used; "
"non-magnetic entries must be 0 (validated "
"against the species layout). Shorthand: give "
"just the N_mag magnetic values. Mutually "
"exclusive with the other --ref-* options.")
g_method.add_argument("--no-symmetry", action="store_true",
help="Disable crystal symmetry (use time-reversal only)")
g_method.add_argument("--check-degeneracy", action="store_true",
dest="check_degeneracy",
help="Keep the configurations that only TIME REVERSAL "
"(the antiunitary part of the magnetic grey "
"group) proves degenerate, instead of merging "
"them: each such pair is computed twice and "
"mag4-extract reports |ΔE| per pair — a direct "
"test of the antiunitary symmetry (the pairs "
"typically differ in total Sz, so no unitary "
"operation relates them) and of the spin model's "
"completeness. Applies to all three methods; "
"costs the extra DFT runs it emits.")
g_run.add_argument("--name", default="four_state",
help="Base name for output files (default: four_state)")
g_run.add_argument("--spin", type=float, default=3.5,
help="Spin quantum number S of the magnetic ion "
"(default: 3.5 for Gd³⁺ S=7/2). Examples: "
"0.5 for Cu²⁺, 1 for Ni²⁺ (S=1) or V³⁺, "
"1.5 for Cr³⁺, 2 for Mn³⁺, 2.5 for Fe³⁺ or "
"Mn²⁺, 3.5 for Gd³⁺. Matches the meaning "
"of --spin in mag4-extract / mag4-magnon. "
"The on-site moment magnitude |MAGMOM| (or "
"the equivalent WIEN2k case.inst initial "
"moment) is set to 2·S (g≈2 spin-only Landé "
"factor). E.g. --spin 0.5 → |MAGMOM|=1, "
"--spin 3.5 → |MAGMOM|=7.")
g_vasp.add_argument("--kspacing", type=float, default=None, metavar="S",
help="Reciprocal-space k-point spacing in Å⁻¹, à la "
"VASP's KSPACING (https://www.vasp.at/wiki/KSPACING): "
"n_i = max(1, ceil(|b_i| / S)) on the 2π-scaled "
"reciprocal lattice. Smaller = denser. Default "
"when neither flag is given: 0.3 Å⁻¹ (denser than "
"VASP's own 0.5; good for magnetic insulators / "
"semiconductors). Metals need a finer mesh — use "
"~0.15–0.2 Å⁻¹ (with appropriate smearing).")
# Deprecated: legacy real-space k-density; --kspacing (VASP convention)
# is the supported control. Hidden; still honoured; main() warns.
parser.add_argument("--kdens", type=float, default=None,
help=argparse.SUPPRESS)
# ── Functional / DFT+U tuning (VASP INCAR; also feeds WIEN2k job.init) ──
g_vasp.add_argument(
"--functional", "--xc",
dest="functional",
type=lambda s: s.upper(),
choices=["PBE", "PBE+U", "PBE0", "R2SCAN", "METAGGA"],
default="PBE",
help="Exchange-correlation functional used in the per-config "
"INCAR (default: PBE). ``PBE+U`` requires --ldaul and "
"--ldauu. ``PBE0`` enables LHFCALC=.TRUE. AEXX=0.25 "
"(no range separation). ``R2SCAN`` forces METAGGA=R2SCAN; "
"``METAGGA`` accepts any flavour via --metagga (SCAN, "
"RSCAN, R2SCAN, R2SCANL, ...). Case-insensitive. The "
"historic alias ``--xc`` is still accepted.",
)
g_vasp.add_argument(
"--ldaul", "--ldau-l",
dest="ldaul", type=int, default=None, metavar="L",
help="Orbital quantum number ℓ of the magnetic element for "
"DFT+U (e.g. 2 for d-electrons, 3 for f-electrons). "
"Required when --functional=PBE+U.",
)
g_vasp.add_argument(
"--ldauu", "--ldau-u",
dest="ldauu", type=float, default=None, metavar="U",
help="Hubbard U (eV) for the magnetic element. Required when "
"--functional=PBE+U. Written to the INCAR as LDAUU in eV; for "
"--code wien2k it is converted to Ry (÷ 13.6057) in job.init's "
"init_orb block, which is the unit WIEN2k expects.",
)
g_vasp.add_argument(
"--ldauj", "--ldau-j",
dest="ldauj", type=float, default=0.0, metavar="J",
help="Hund's exchange J (eV) for the magnetic element "
"(default: 0.0; combined with U this gives the effective "
"Liechtenstein U_eff = U − J). Written to the INCAR as LDAUJ "
"in eV; for --code wien2k it is converted to Ry (÷ 13.6057) "
"alongside U.",
)
g_vasp.add_argument(
"--aexx", type=float, default=0.25, metavar="X",
help="Fraction of exact exchange used by the PBE0 hybrid "
"(default: 0.25). Ignored unless --functional=PBE0.",
)
g_vasp.add_argument(
"--encut", type=float, default=500.0, metavar="EV",
help="Plane-wave cut-off in eV written to ENCUT in every "
"INCAR (default: 500).",
)
g_vasp.add_argument(
"--metagga", default="R2SCAN", metavar="FLAVOUR",
help="Meta-GGA flavour written for the METAGGA tag when "
"--functional=METAGGA (default: R2SCAN; typical other "
"values: SCAN, RSCAN, R2SCAN, R2SCANL). When "
"--functional=R2SCAN the flavour is always forced to "
"R2SCAN.",
)
def _lmaxtau_arg(v: str):
if str(v).strip().lower() == "auto":
return "auto"
try:
iv = int(v)
except (TypeError, ValueError):
raise argparse.ArgumentTypeError(
"--lmaxtau must be 'auto' or an even integer in {0,2,4,6,8}")
if iv not in (0, 2, 4, 6, 8):
raise argparse.ArgumentTypeError(
"--lmaxtau integer must be one of 0,2,4,6,8 "
"(s→2, p→4, d→6, f→8)")
return iv
g_vasp.add_argument(
"--lmaxtau", type=_lmaxtau_arg, default="auto",
metavar="{auto,0,2,4,6,8}",
help="Meta-GGA LMAXTAU. 'auto' (default): 8 if a "
"lanthanide/actinide is present, else VASP's default 6 is "
"used. Force 6 explicitly for f-in-core POTCARs (RE_3, "
"La_s) where the 4f is not in the valence "
"(s→2, p→4, d→6, f→8).",
)
# ── Method selection ────────────────────────────────────────────────────
g_method.add_argument("--method",
choices=["four-state", "energy-mapping", "16-state"],
default="four-state",
help="J-extraction method these inputs are prepared for: "
"'four-state' (default) writes the uu/dd/ud/du "
"configurations and the per-bond difference "
"formulas; 'energy-mapping' enumerates every "
"inequivalent collinear order of the (super)cell "
"for a linear least-squares fit (consumed by "
"'mag4-extract --method energy-mapping'); "
"'16-state' flips the four atoms of one square "
"plaquette (2⁴=16 states, symmetry-reduced) to "
"isolate the four-spin RING exchange Jring on a 2D "
"square lattice (consumed by 'mag4-extract --method "
"16-state'; pin the square with --jring-coupling / "
"--jring-dist).")
g_method.add_argument("--max-atoms", type=int, default=20, metavar="N",
dest="max_atoms",
help="Safety cap on the magnetic-atom count for "
"--method energy-mapping (the enumeration is 2^N; "
"default: 20).")
g_method.add_argument("--jring", action="store_true", default=None,
help="Account for the four-spin cyclic (ring) exchange "
"Jring (Σ s_i s_j s_k s_l over nearest-neighbour "
"square plaquettes), for square-lattice systems "
"(e.g. cuprates). --method energy-mapping: adds a "
"fitted 'Jring' column (opt-in — it changes the "
"fitted model). --method four-state: adds "
"the Jring·S⁴ term to the per-config energy "
"equations and a CONTAMINATION CHECK showing how "
"much Jring leaks into each bilinear J (the "
"four-state formulas have no ring term, so a "
"non-zero leak biases J). This four-state check "
"is a pure diagnostic — it changes nothing in the "
"generated inputs — and therefore runs BY DEFAULT "
"(disable with --no-jring). The square shell is "
"auto-detected (need not be the shortest M-M "
"distance; pin with --jring-coupling/--jring-dist).")
g_method.add_argument("--no-jring", dest="jring", action="store_false",
help="Disable the default-on four-state Jring "
"contamination diagnostic.")
g_method.add_argument("--jring-coupling", default=None, metavar="LABEL",
dest="jring_coupling",
help="Force the ring term onto a named coupling shell "
"(e.g. J2) instead of auto-detecting the square "
"lattice.")
g_method.add_argument("--jring-dist", type=float, default=None, metavar="ANG",
dest="jring_dist",
help="Force the ring-term square edge length in Å "
"(overrides --jring-coupling and auto-detection).")
g_method.add_argument("--jprime", nargs="?", type=float, const=0.0,
default=None, metavar="RADIUS",
help="--method four-state: check BEYOND-CUTOFF bilinear "
"shells J' for contamination. The four-state "
"formula only cancels couplings that do not fold "
"onto a self-image, so a J' shell with d > --cutoff "
"can leak into the extracted Jk in a small cell "
"(e.g. NiO 1×1×2: J4 folds into J2). Detects every "
"M–M shell out to RADIUS (default: the larger of "
"2×--cutoff and the shortest lattice translation), "
"and for each J' beyond "
"--cutoff prints its contribution to every configN "
"energy and its leak into each Jk (the J-formula's ± "
"pattern applied to c', ÷ divisor). A pure "
"diagnostic — it changes nothing in the generated "
"inputs — and therefore runs BY DEFAULT when --cif "
"is given (disable with --no-jprime).")
g_method.add_argument("--no-jprime", dest="no_jprime", action="store_true",
default=False,
help="Disable the default-on four-state beyond-cutoff "
"J' contamination diagnostic.")
# ── DFT-code selection ──────────────────────────────────────────────────
g_run.add_argument("--code", choices=["vasp", "wien2k", "both"],
default="vasp",
help="Which DFT-code inputs to generate per config dir: "
"'vasp' = POSCAR + INCAR + KPOINTS (default), "
"'wien2k' = POSCAR + configN.inst + configN.struct "
"(struct file written directly in Python — no manual "
"xyz2struct step needed; pass --no-struct-gen to "
"fall back to the legacy POSCAR -> xyz2struct path), "
"'both' = all of the above.")
g_wien.add_argument("--instgen", default=None, metavar="PATH",
help="Path to the WIEN2k 'instgen' script. If omitted, "
"auto-detected via $WIEN2k_INSTGEN, $PATH, $WIENROOT, "
"then ~/bin/instgen. Required when --code includes wien2k.")
# ── WIEN2k struct-file generation (replaces xyz2struct) ────────────────
# Deprecated: job.init runs `setrmt`, which overrides whatever RMT mag4
# writes, and the standard 781-point radial mesh needs no knob. Hidden;
# ignored (struct files use the built-in defaults); main() warns.
parser.add_argument("--rmt", type=float, default=None,
help=argparse.SUPPRESS)
parser.add_argument("--npt", type=int, default=None,
help=argparse.SUPPRESS)
g_wien.add_argument("--no-struct-gen", dest="gen_struct",
action="store_false", default=True,
help="Skip the in-Python configN.struct generation "
"(useful if you prefer to run WIEN2k's xyz2struct "
"yourself from the POSCAR). Implies "
"--no-symmetrize-struct (the symmetric path "
"needs to write the struct).")
g_wien.add_argument("--no-symmetrize-struct", dest="symmetrize_struct",
action="store_false", default=True,
help="Disable spglib-based symmetry detection for "
"configN.struct. The legacy path is used "
"instead: every atom emitted as MULT=1 with the "
"identity-only symmetry block. configN.inst is "
"generated separately (POSCAR atom order).")
g_wien.add_argument("--symprec", type=float, default=1e-3, metavar="E",
help="Distance tolerance (Å) passed to spglib when "
"detecting the symmetry of each spin-decorated "
"configuration (default: 1e-3). Loosen if "
"DFT-relaxed coordinates drift slightly.")
g_wien.add_argument("--no-job-scripts", dest="gen_job_scripts",
action="store_false", default=True,
help="Skip writing the WIEN2k driver scripts "
"(job.init and job.run) at the compound directory.")
g_wien.add_argument("--reduce", dest="reduce_cell",
action="store_true", default=False,
help="Reduce each configN.struct to the primitive "
"cell of its own spin pattern (colored "
"structure), so the Wyckoff multiplicities "
"match the parent SG's standard table. "
"Different configs may then end up in "
"different cell sizes — mag4-extract handles "
"this transparently by rescaling each "
"contributing config's total energy to the "
"least-common-multiple (LCM) atom count of "
"the four configs of each J before applying "
"the four-state formula (DFT energy is "
"extensive, so the scaling is exact). "
"Default: OFF — every configN.struct keeps "
"mag4's one working supercell, so all configs "
"share strictly identical cell and numerics.")
g_wien.add_argument("--no-reduce", dest="reduce_cell",
action="store_false",
help="Keep every configN.struct in the working "
"supercell (this is the default; the flag "
"exists for explicitness and for comparing "
"against --reduce).")
# Deprecated spelling of --no-reduce (which is now the default anyway).
# Hidden; still honoured; main() warns.
parser.add_argument("--no-reduce-cell", dest="reduce_cell",
action="store_false", help=argparse.SUPPRESS)
# ── Geometry-aware J shells (bridging-ligand fingerprint) ──────────────
g_lig.add_argument("--ligand", default=None, metavar="X[,Y,...]",
help="Bridging-ligand element symbol(s) for "
"geometry-aware J-shell discovery (e.g. 'O' for "
"CuO). When set, pairs at the same distance "
"but different M-L-M bridging angle are split "
"into sub-shells (J1a, J1b, ...). Default: "
"disabled (distance-only grouping).")
g_lig.add_argument("--bond-cutoff", type=float, default=2.5, metavar="R",
help="Maximum M-L bond distance (Å) for a ligand to "
"count as a bridge (default: 2.5). Only used "
"when --ligand is set.")
g_lig.add_argument("--angle-tol", type=float, default=1.0, metavar="DEG",
help="Tolerance (°) for grouping M-L-M angles and "
"M-L-L-M dihedrals into the same fingerprint "
"(default: 1.0). Only used when --ligand is set.")
g_lig.add_argument("--ll-cutoff", type=float, default=3.0, metavar="R",
help="Maximum ligand-ligand distance (Å) for an "
"M-L-L-M super-super-exchange path to count as "
"a bridge (default: 3.0). Used **only** for "
"pairs with no direct M-L-M bridge, so shells "
"with a direct bridge keep their existing "
"fingerprint. Set 0 to disable the dihedral "
"fallback. Only used when --ligand is set.")
# Deprecated: exact equivalent of --ll-cutoff 0. Hidden; main() warns.
parser.add_argument("--no-dihedral", dest="use_dihedral",
action="store_false", default=True,
help=argparse.SUPPRESS)
g_lig.add_argument("--signed-dihedral", dest="signed_dihedral",
action="store_true", default=False,
help="Keep the IUPAC sign of M-L-L-M dihedrals so "
"chiral mirror-image super-super-exchange paths "
"are treated as distinct J sub-shells. Use "
"this when Dzyaloshinskii-Moriya / antisymmetric "
"exchange matters. Default: dihedrals folded "
"to |θ| in [0°, 180°] so chiral "
"counterparts merge — the right choice for "
"plain Heisenberg J.")
return parser
def _run_energy_mapping(args, crystal, couplings, supercell,
moment_magnitude) -> int:
"""Energy-mapping path: enumerate every inequivalent collinear order of the
chosen (super)cell, write one config directory each plus the design-matrix
report consumed by ``mag4-extract --method energy-mapping``.
"""
na, nb, nc = supercell
try:
records, labels, ring_dist = enumerate_orders(
crystal, couplings, tol=args.tol, max_atoms=args.max_atoms,
ring=bool(args.jring), ring_coupling=args.jring_coupling,
ring_dist=args.jring_dist)
except ValueError as exc:
log.error("%s", exc)
return 1
# The report lists every coupling column with a distance; the ring term's
# "distance" is the square edge length it was built on.
report_couplings = list(couplings)
if "Jring" in labels and ring_dist is not None:
report_couplings = report_couplings + [["Jring", ring_dist]]
# Self-image folding diagnostic: a coupling at/above the shortest lattice
# translation has a constant coefficient and is indeterminate in the fit.
d_short = shortest_periodic_image(crystal.lattice_matrix)
folded = [name for name, d in [(c[0], c[1]) for c in couplings] if d >= d_short]
if folded:
log.warning("Coupling(s) %s reach/exceed the shortest lattice "
"translation (%.3f Å) — they fold onto self-images, so "
"their coefficient is constant and will be dropped as "
"indeterminate in the fit. Enlarge the supercell.",
", ".join(folded), d_short)
# --check-degeneracy: emit the time-reversal partner (global spin flip) of
# every order whose flip no UNITARY crystal operation reaches. The two
# orders share the model coefficients, so their DFT energies must coincide
# (collinear, no SOC) — mag4-extract reports |ΔE| per pair as a direct
# test of the antiunitary symmetry.
tr_pairs: List[tuple] = []
if args.check_degeneracy:
mag_cart = get_cartesian_coords(crystal.target_species,
crystal.lattice_matrix)
ops = find_symmetry_ops(mag_cart, crystal.lattice_vectors, tol=0.05,
crystal=crystal)
partners = []
next_num = max((r["number"] for r in records), default=0) + 1
for rec in records:
sv = tuple(rec["spin_vector"])
flip = tuple(-s for s in sv)
if any(tuple(sv[op[i]] for i in range(len(sv))) == flip
for op in ops):
continue # a unitary op already reaches the flip
partners.append({
"number": next_num,
"degeneracy": rec["degeneracy"],
"sz": rec["sz"],
"coefficients": dict(rec["coefficients"]),
"spin_vector": [-s for s in rec["spin_vector"]],
})
tr_pairs.append((rec["number"], next_num))
next_num += 1
if partners:
records = records + partners
log.info("--check-degeneracy: %d time-reversal partner order(s) "
"added — each pair must be degenerate; mag4-extract "
"checks |ΔE|.", len(partners))
else:
log.info("--check-degeneracy: every order's global flip is "
"already unitary-equivalent; no partner orders needed.")
unique_groups = build_unique_groups(records)
magmom_lines = build_magmom_lines(unique_groups, crystal, moment_magnitude)
# Magnetic (Shubnikov) group per order — BNS number + type, primed
# operations included (informational).
from mag4.wien2k import analyze_config_magnetic_group
config_msgs = [analyze_config_magnetic_group(
crystal, g, spin_magnitude=moment_magnitude) for g in unique_groups]
compound = (os.path.splitext(os.path.basename(args.cif))[0]
if args.cif else args.name)
n_mag = len(crystal.target_species)
W = 80
print()
print("=" * W)
n_partner = len(tr_pairs)
print(f" ENERGY-MAPPING — {len(records) - n_partner} inequivalent "
f"order(s) of the {na}×{nb}×{nc} supercell"
+ (f" (+{n_partner} time-reversal partner(s))" if n_partner else ""))
print(f" ({n_mag} magnetic atoms → 2^{n_mag} = {2**n_mag:,} collinear orders)")
print(f" Couplings: "
+ ", ".join(f"{c[0]}={c[1]:.4f} Å" for c in report_couplings)
+ (" [Jring = 4-spin ring on the square lattice]"
if "Jring" in labels else ""))
print("=" * W)
print(f" {'config':>9} {'deg':>5} {'Sz':>4} "
+ " ".join(f"{l:>6}" for l in labels)
+ f" {'MSG (BNS)':>12} {'MPG':>10}")
for i, rec in enumerate(records):
coefs = " ".join(f"{rec['coefficients'][l]:>6d}" for l in labels)
m = config_msgs[i] if i < len(config_msgs) else None
msg = f"{m.bns_number}-{m.type_roman}" if m is not None else "—"
mpg = (m.mpg_symbol or "—") if m is not None else "—"
print(f" config{rec['number']:<3d} {rec['degeneracy']:>5} "
f"{rec['sz']:>4} {coefs} {msg:>12} {mpg:>10}")
if tr_pairs:
print("-" * W)
print(" TIME-REVERSAL DEGENERACY CHECK — pairs related only by the")
print(" ANTIUNITARY symmetry (global spin flip); their DFT energies")
print(" must coincide. mag4-extract reports |ΔE| per pair.")
for a, b in tr_pairs:
print(f" pair: config{a} == config{b}")
print("=" * W)
n_param = 1 + len(labels)
if len(records) < n_param:
log.warning("Only %d inequivalent configuration(s) for %d parameters "
"(E0 + %d coupling(s)); the linear fit will be "
"rank-deficient. Use a larger supercell to generate more "
"orders.", len(records), n_param, len(labels))
if args.dry_run:
print(" --dry-run : skipping the report and configN/ directories.")
print("=" * W)
return 0
os.makedirs(compound, exist_ok=True)
report_path = os.path.join(compound, f"{args.name}_report.dat")
write_energy_mapping_report(report_path, records, labels, report_couplings,
(na, nb, nc), args.name, tr_pairs=tr_pairs,
config_msgs=config_msgs)
log.info("Design-matrix report written: %s", report_path)
write_vasp = args.code in ("vasp", "both")
write_wien = args.code in ("wien2k", "both")
# POSCAR is always written; INCAR + KPOINTS only when a VASP run is wanted.
create_config_dirs(unique_groups, crystal, magmom_lines,
args.name, args.kdens, compound,
write_vasp_extras=write_vasp,
kspacing=args.kspacing,
functional=args.functional,
magnetic_element=args.element,
ldaul=args.ldaul,
ldauu=args.ldauu,
ldauj=args.ldauj,
aexx=args.aexx,
encut=args.encut,
metagga=args.metagga,
lmaxtau=args.lmaxtau)
if write_wien:
# Energy-mapping fits every order in the SAME supercell, so each
# configN.struct must be that one cell — only then are the total
# energies directly comparable in the linear fit. The four-state
# symmetric / reduce-cell path writes per-config reduced primitives
# (correct there only because mag4-extract rescales them by atom count),
# so it is intentionally bypassed here.
if args.gen_struct and args.symmetrize_struct:
log.info("energy-mapping: writing every configN.struct in the full "
"%d×%d×%d supercell (the symmetric / reduce-cell path is "
"bypassed so all energies share one cell for the fit).",
na, nb, nc)
try:
write_case_inst_files(unique_groups, crystal,
parent_dir=compound,
instgen_path=args.instgen,
spin_magnitude=moment_magnitude)
except RuntimeError as exc:
log.error("WIEN2k case.inst generation failed: %s", exc)
return 1
if args.gen_struct:
try:
write_case_struct_files(crystal, len(unique_groups), compound)
except (ValueError, OSError) as exc:
log.error("WIEN2k configN.struct generation failed: %s", exc)
return 1
if args.gen_job_scripts:
try:
write_job_scripts(compound, functional=args.functional,
kspacing=args.kspacing,
element=args.element, ldaul=args.ldaul,
ldauu=args.ldauu, ldauj=args.ldauj)
except OSError as exc:
log.error("WIEN2k job script generation failed: %s", exc)
return 1
log.info("Energy-mapping inputs written under %s/", compound)
return 0
def _run_ring16(args, crystal, couplings, supercell, moment_magnitude,
*, base_crystal=None, safe_tol=0.05) -> int:
"""16-state path: flip the four atoms of one square plaquette (2⁴ = 16
states), reduce by crystal symmetry + time reversal, and write the ring
report + configN/ directories consumed by ``mag4-extract --method 16-state``.
All configs share the one supercell so their energies are directly
comparable in the χ-weighted projection.
The plaquette must be fully ISOLATED (multiplicity M=1). The auto-supercell
search guarantees this; an explicit --supercell that leaves M>1 is a hard
error (the four atoms would wrap onto their own J1 images).
"""
na, nb, nc = supercell
atoms = crystal.target_species
n = len(atoms)
try:
quad = select_ring_quad(crystal, couplings,
prefer_label=args.jring_coupling,
prefer_dist=args.jring_dist)
except ValueError as exc:
log.error("%s", exc)
return 1
log.info("Ring plaquette: %s (square edge %.4f Å, multiplicity %d)",
" -- ".join(quad.labels), quad.edge, quad.multiplicity)
# Mandatory isolation: M must be 1 (the ring analogue of four-state dimer
# isolation). Only reachable via an explicit --supercell; recommend the
# smallest isolating cell of the input structure.
if quad.multiplicity > 1:
rec = None
if base_crystal is not None:
try:
rec, _sc, _q = find_ring_supercell(
base_crystal, couplings,
max_mult=max(args.max_supercell, max(na, nb, nc) + 1),
prefer_label=args.jring_coupling, prefer_dist=args.jring_dist)
except (RuntimeError, ValueError):
rec = None
msg = [
"",
"!" * 78,
f" 16-STATE RING ISOLATION ERROR — supercell {na}×{nb}×{nc} leaves "
f"the square plaquette at multiplicity M={quad.multiplicity}; no "
f"inputs were written.",
" The four plaquette atoms fold onto their own J1 images through the "
"periodic boundary, so they are not independent.",
"!" * 78,
]
if rec is not None:
msg += ["", " Smallest fully isolated (M=1) cell:",
f" --supercell {rec[0]} {rec[1]} {rec[2]}",
" Or drop --supercell to let mag4 find it automatically."]
else:
msg += ["", " Raise --max-supercell, or drop --supercell to let "
"mag4 search for an isolated cell automatically."]
log.error("\n".join(msg))
return 2
# Reference bath. Only the four plaquette atoms are flipped over it, and the
# χ-projection is exact for ANY fixed bath — so the bath does not change the
# extracted Jring, only how close the 16 configs sit to the ground state.
# Default is FM, but --ref-couplings / --ref-state / --ref-kvector let the
# 180° superexchange be set AFM (J1 for La₂CuO₄, J2 for SrFeO₂), which DFT
# converges to far more readily.
lat = crystal.lattice_matrix
ref_opts = [args.ref_state, args.ref_couplings, args.ref_kvector,
args.ref_magmom]
if sum(o is not None for o in ref_opts) > 1:
log.error("--ref-state, --ref-couplings, --ref-kvector and "
"--ref-magmom are mutually exclusive.")
return 1
if args.ref_magmom is not None:
# Bath given directly as the supercell MAGMOM string (signs used).
try:
ref_bath = build_reference_from_magmom(args.ref_magmom, crystal)
except ValueError as exc:
log.error("%s", exc)
return 1
elif args.ref_kvector is not None:
try:
kfrac = parse_kvector(args.ref_kvector)
except ValueError as exc:
log.error("Invalid --ref-kvector: %s", exc)
return 1
ref_bath, _rep = build_reference_from_kvector(atoms, kfrac, (na, nb, nc))
elif args.ref_couplings is not None:
try:
ref_bath, rep = build_reference_from_couplings(
atoms, lat, couplings, args.ref_couplings, safe_tol)
except ValueError as exc:
log.error("Invalid --ref-couplings: %s", exc)
return 1
for name, (sat, total) in rep["couplings"].items():
sign = "AFM" if rep["wanted"][name] == -1 else "FM"
if total and sat == total:
log.info(" bath %s (%s): %d/%d bonds satisfied.",
name, sign, sat, total)
elif total:
log.warning(" bath %s (%s): only %d/%d bonds satisfied "
"(frustrated in this cell — still a valid bath, but "
"a commensurate cell converges better).",
name, sign, sat, total)
else:
user_ref = None
if args.ref_state is not None:
try:
user_ref = parse_ref_state(args.ref_state)
except (ValueError, argparse.ArgumentTypeError) as exc:
log.error("Invalid --ref-state: %s", exc)
return 1
try:
ref_bath = build_reference_bath(n, user_ref)
except ValueError as exc:
log.error("%s", exc)
return 1
log.info("Reference bath: %s", " ".join(f"{s:+d}" for s in ref_bath))
if args.no_symmetry:
log.info("Crystal symmetry disabled — using time-reversal only.")
sym_ops = [tuple(range(n))]
else:
log.info("Computing crystal symmetry operations (full crystal) …")
mag_cart = get_cartesian_coords(atoms, crystal.lattice_matrix)
sym_ops = find_symmetry_ops(mag_cart, crystal.lattice_vectors, tol=0.05,
crystal=crystal)
configs = generate_ring_16(ref_bath, quad)
# --check-degeneracy: group by the UNITARY ops only, so configurations the
# antiunitary part (unitary × time reversal) proves degenerate stay
# separate DFT runs; the recorded pairs let mag4-extract verify |ΔE| ≈ 0.
unique_groups = find_unique_configs(
configs, sym_ops, time_reversal=not args.check_degeneracy)
tr_pairs = (time_reversal_pairs(unique_groups, sym_ops)
if args.check_degeneracy else None)
log.info("16 plaquette states reduced to %d inequivalent config(s)%s.",
len(unique_groups),
(f" (unitary only; {len(tr_pairs)} time-reversal pair(s) kept "
"for the degeneracy check)") if args.check_degeneracy else "")
n_supercell_atoms = sum(len(g) for g in crystal.all_species)
compound = (os.path.splitext(os.path.basename(args.cif))[0]
if args.cif else args.name)
# Pairwise-J analysis over the SAME configs: design-matrix rows + per-J
# extractability verdicts (no extra DFT; the section is analysis-only).
design_records = ring_design_records(
crystal, couplings, unique_groups, quad,
prefer_label=args.jring_coupling, prefer_dist=args.jring_dist,
tol=safe_tol)
coupling_dists = {c[0]: float(c[1]) for c in couplings}
coupling_dists["Jring"] = quad.edge
# Space group of each spin-decorated config (spglib) — informational only,
# the written structures stay P1. Alongside the ordinary (unitary-only)
# SG, classify the full MAGNETIC (Shubnikov) group — primed operations
# included — as the BNS number + type.
from mag4.wien2k import analyze_config_magnetic_group, analyze_config_symmetry
config_sgs = []
config_msgs = []
for group in unique_groups:
sym = analyze_config_symmetry(crystal, group,
spin_magnitude=moment_magnitude)
config_sgs.append((sym.spacegroup_number, sym.spacegroup_symbol))
config_msgs.append(analyze_config_magnetic_group(
crystal, group, spin_magnitude=moment_magnitude))
report_path = None
if not args.dry_run:
os.makedirs(compound, exist_ok=True)
report_path = os.path.join(compound, f"{args.name}_report.dat")
symmetric_structs = (args.code in ("wien2k", "both")
and args.gen_struct and args.symmetrize_struct)
write_ring_report(report_path, unique_groups, atoms, quad, (na, nb, nc),
len(sym_ops), n_supercell_atoms=n_supercell_atoms,
ref_bath=ref_bath,
design_records=design_records,
coupling_dists=coupling_dists,
config_sgs=config_sgs,
config_msgs=config_msgs,
symmetric_structs=symmetric_structs,
tr_pairs=tr_pairs)
if args.dry_run:
print("\n --dry-run : no configN/ directories or report file written.")
return 0
magmom_lines = build_magmom_lines(unique_groups, crystal, moment_magnitude)
# Record the target moment (|MAGMOM| = 2·S) in the report so mag4-moments
# can auto-check the configs reached it — same header as the four-state path.
if report_path is not None:
with open(report_path, "a") as f:
f.write("\n" + "=" * 80 + "\n")
f.write(f" MAGMOM (S = {args.spin:g}, |MAGMOM| = 2·S = "
f"{moment_magnitude:g} μ_B)\n")
f.write("=" * 80 + "\n")
for uid, (group, mm) in enumerate(zip(unique_groups, magmom_lines), 1):
f.write(f" Config #{uid}\n MAGMOM = {mm}\n\n")
write_vasp = args.code in ("vasp", "both")
write_wien = args.code in ("wien2k", "both")
create_config_dirs(unique_groups, crystal, magmom_lines, args.name,
args.kdens, compound, write_vasp_extras=write_vasp,
kspacing=args.kspacing, functional=args.functional,
magnetic_element=args.element, ldaul=args.ldaul,
ldauu=args.ldauu, ldauj=args.ldauj, aexx=args.aexx,
encut=args.encut, metagga=args.metagga,
lmaxtau=args.lmaxtau)
if write_wien:
if args.gen_struct and args.symmetrize_struct:
# Symmetry-detected structs (non-P1): each configN.struct carries
# its spin space group (Wyckoff orbits + symmetry-op tail). The
# cell itself is NEVER reduced here — all configs must share the
# one supercell so their energies are directly comparable in the
# χ-projection and the pairwise fit (no per-config rescaling).
if args.reduce_cell:
log.info("16-state: symmetric configN.struct in the full "
"%d×%d×%d supercell (per-config cell reduction is "
"bypassed so all energies share one cell for the "
"projection/fit).", na, nb, nc)
try:
write_symmetric_case_files(
unique_groups, crystal,
parent_dir=compound,
instgen_path=args.instgen,
spin_magnitude=moment_magnitude,
symprec=args.symprec,
reduce_cell=False,
)
except RuntimeError as exc:
log.error("WIEN2k symmetric generation failed: %s", exc)
return 1
else:
try:
write_case_inst_files(unique_groups, crystal,
parent_dir=compound,
instgen_path=args.instgen,
spin_magnitude=moment_magnitude)
except RuntimeError as exc:
log.error("WIEN2k case.inst generation failed: %s", exc)
return 1
if args.gen_struct:
try:
write_case_struct_files(crystal, len(unique_groups),
compound)
except (ValueError, OSError) as exc:
log.error("WIEN2k configN.struct generation failed: %s",
exc)
return 1
if args.gen_job_scripts:
try:
write_job_scripts(compound, functional=args.functional,
kspacing=args.kspacing,
element=args.element, ldaul=args.ldaul,
ldauu=args.ldauu, ldauj=args.ldauj)
except OSError as exc:
log.error("WIEN2k job script generation failed: %s", exc)
return 1
log.info("16-state ring inputs written under %s/", compound)
return 0
def _suggest_ref_supercell(base_crystal, conv_lat_mat, conv_targets, couplings,
ref_spec, safe_tol, current, max_mult, current_frac,
*, max_multiplicity=None):
"""Smallest **isolating** supercell larger than ``current`` (by magnetic-atom
count) that reduces the ``--ref-couplings`` frustration. Returns
``(cell, frac)`` — the chosen multipliers and the broken-bond fraction it
achieves (``frac == 0`` means the order is hosted *fully*) — or
``(None, None)`` if none within ``max_mult`` improves on ``current_frac``.
Isolation alone fixes the supercell for the four-state J extraction, but the
requested reference order can still be *frustrated* in that cell (e.g. a Néel
AFM needs even multipliers — 3×3×1 frustrates it, 4×4×1 hosts it). This
proposes the smallest cell that both isolates and carries the reference with
the fewest broken bonds.
"""
import itertools
cur_atoms = current[0] * current[1] * current[2]
best, best_key, best_frac = None, None, None
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)):
if cand[0] * cand[1] * cand[2] <= cur_atoms:
continue # only larger cells
iso = check_four_state_isolation(conv_lat_mat, cand, couplings,
target_species=conv_targets,
max_multiplicity=max_multiplicity)
if not iso["ok"]:
continue
sc = expand_supercell(base_crystal, *cand)
try:
_, rep = build_reference_from_couplings(
sc.target_species, sc.lattice_matrix, couplings, ref_spec, safe_tol)
except ValueError:
continue
broken = sum(total - sat for sat, total in rep["couplings"].values())
total = sum(total for _, total in rep["couplings"].values())
if total == 0:
continue
frac = broken / total # fraction, not raw count
if frac < current_frac - 1e-9: # strictly less frustrated
# tie-break: lowest frustration, then most symmetric (smallest spread
# of the expanded axes — prefers 4×4×1 over an asymmetric 3×4×1), then
# fewest atoms.
active = [m for m in cand if m > 1]
spread = (max(active) - min(active)) if active else 0
key = (round(frac, 6), spread, cand[0] * cand[1] * cand[2], max(cand))
if best_key is None or key < best_key:
best, best_key, best_frac = cand, key, frac
return best, best_frac
_DEPRECATED_FLAGS = {
# flag → what to use instead (soft deprecation: still works, hidden from
# --help, one-line warning here so old scripts keep running unchanged).
"--ref-state": "use --ref-magmom (same spins, MAGMOM-style, "
"n*v repeats, validated full-supercell form)",
"--kdens": "use --kspacing (VASP-style, Å⁻¹; default 0.3)",
"--no-dihedral": "use --ll-cutoff 0 (exact equivalent)",
"--auto-supercell": "redundant — the auto search is already the default "
"when --supercell is not given",
"--rmt": "removed (ignored) — job.init runs `setrmt`, which "
"overrides any RMT mag4 writes; tune radii there "
"(e.g. `setrmt -a Cu:2.1,O:1.7`)",
"--npt": "removed (ignored) — the struct files use the "
"standard 781-point radial mesh",
"--no-reduce-cell": "keeping the working supercell is now the default; "
"use --reduce to opt back into per-config "
"primitive-cell reduction",
}
[docs]
def main(argv=None) -> int:
parser = build_parser()
args = parser.parse_args(argv)
raw_argv = list(argv) if argv is not None else sys.argv[1:]
for flag, hint in _DEPRECATED_FLAGS.items():
if any(a == flag or a.startswith(flag + "=") for a in raw_argv):
log.warning("%s is deprecated: %s.", flag, hint)
# Effective multiplicity cap for the four-state isolation check.
# --full-isolation is shorthand for --max-multiplicity 1; if both are
# given the most restrictive (smallest) cap wins. None ⇒ unbounded
# (the default algebraic criterion, divisor 4·M).
if args.max_multiplicity is not None and args.max_multiplicity < 1:
log.error("--max-multiplicity must be >= 1 (got %d).",
args.max_multiplicity)
return 2
_caps = [c for c in (1 if args.full_isolation else None,
args.max_multiplicity) if c is not None]
max_mult_cap = min(_caps) if _caps else None
# Human-friendly name of the flag that set the cap (for error messages).
cap_flag = ("--full-isolation" if (max_mult_cap == 1 and args.full_isolation)
else f"--max-multiplicity {max_mult_cap}")
# ``--spin`` is the spin quantum number S of the magnetic ion
# (matching mag4-extract / mag4-magnon semantics). The downstream
# MAGMOM / case.inst code wants a magnitude in Bohr magnetons; we
# use the spin-only Landé factor g≈2, so |MAGMOM| = 2·S.
moment_magnitude = 2.0 * args.spin
# Early-fail on incomplete PBE+U options so the user sees the error
# even in --dry-run mode (no point computing the supercell if the run
# can't write a valid INCAR / job.init init_orb block).
if args.functional == "PBE+U":
missing = []
if args.ldaul is None:
missing.append("--ldaul")
if args.ldauu is None:
missing.append("--ldauu")
if missing:
log.error(
"--functional PBE+U requires %s. Example for a Cu-d "
"system: `--functional PBE+U --ldaul 2 --ldauu 4.0 "
"--ldauj 0.0`.",
" and ".join(missing),
)
return 2
# Parse --ligand once (comma-separated → list). None when not set.
ligand_elements = (
[s.strip() for s in args.ligand.split(",") if s.strip()]
if args.ligand else None
)
# --no-dihedral wins over --ll-cutoff: effective LL cutoff is 0 when
# the user opts out, otherwise the user-supplied (or default 3.0 Å).
eff_ll_cutoff = args.ll_cutoff if args.use_dihedral else 0.0
# ---- Step 1: structure + couplings ----------------------------------------
if args.cif:
element = args.element
if element is None:
element = input("Magnetic element symbol (e.g. Gd): ").strip()
crystal, couplings = run_cif_analysis(
args.cif, element, args.cutoff, args.equiv_tol,
ligand_elements=ligand_elements,
bond_cutoff=args.bond_cutoff,
angle_tol=args.angle_tol,
ll_cutoff=eff_ll_cutoff,
signed_dihedral=args.signed_dihedral,
)
else:
if args.couplings is None:
args.couplings = input("Couplings (e.g. J1:3.94,J2:6.10): ").strip()
couplings = parse_couplings(args.couplings)
species_idx = args.species
if species_idx is None:
species_idx = int(input("Species index in POSCAR (1-based): "))
crystal = read_poscar(args.poscar, species_idx)
log.info("Cell: a=%.3f b=%.3f c=%.3f α=%.1f° β=%.1f° γ=%.1f°", *crystal.cell)
# Stash the conv-cell lattice matrix AND the magnetic-atom list
# BEFORE supercell expansion; ``check_four_state_isolation`` needs
# them both to derive per-J dimer directions and supercell
# recommendations below.
conv_lat_mat = np.asarray(crystal.lattice_matrix, dtype=float)
conv_target_species = list(crystal.target_species)
base_crystal = crystal # conv cell, kept for the ref-state search
# Compute the safe dimer-search tolerance once — both the auto-supercell
# search and the final dimer search below need it, and we want the
# "tol shrunk" warning to fire at most once.
safe_tol = check_coupling_spacing(couplings, args.tol)
# ---- Step 2: supercell ----------------------------------------------------
# Three modes:
# (a) --supercell explicitly given → use as-is.
# (b) auto-supercell (default) → search for the smallest cell
# that gives every J a distinct
# intra-cell pair.
# (c) --no-auto-supercell → fall back to 1×1×1 and let
# find_dimers print its legacy
# per-axis warning.
if args.supercell is not None:
na, nb, nc = args.supercell
log.info("Using user-specified supercell %d×%d×%d "
"(auto search skipped).", na, nb, nc)
if (na, nb, nc) != (1, 1, 1):
print(f"\nExpanding to {na}×{nb}×{nc} supercell …")
crystal = expand_supercell(crystal, na, nb, nc)
elif args.auto_supercell:
cutoff_str = (f"cutoff = {args.cutoff:g} Å" if args.cif
else "user-supplied couplings")
if args.method == "16-state":
# The 16-state ring method needs the target square plaquette fully
# ISOLATED (multiplicity M=1: four independent corners, divisor 16)
# — the ring analogue of four-state dimer isolation.
print(f"\nSearching for the smallest supercell with a fully isolated "
f"(M=1) ring plaquette for {cutoff_str} …")
try:
(na, nb, nc), crystal, _q = find_ring_supercell(
crystal, couplings, max_mult=args.max_supercell,
prefer_label=args.jring_coupling, prefer_dist=args.jring_dist)
except (RuntimeError, ValueError) as exc:
log.error("%s", exc)
return 1
elif args.method == "energy-mapping":
# Energy-mapping needs a cell with no self-image folding (every
# coupling below the shortest lattice translation) — a different
# criterion from four-state dimer isolation.
print(f"\nSearching for the smallest supercell with no self-image "
f"folding for {cutoff_str} …")
(na, nb, nc), crystal = find_supercell_for_enumeration(
crystal, couplings, max_mult=args.max_supercell)
else:
print(f"\nSearching for the smallest supercell that resolves and "
f"isolates every J within {cutoff_str} …")
try:
(na, nb, nc), crystal = find_optimal_supercell(
crystal, couplings, safe_tol,
max_mult=args.max_supercell,
ligand_elements=ligand_elements,
bond_cutoff=args.bond_cutoff,
angle_tol=args.angle_tol,
ll_cutoff=eff_ll_cutoff,
signed_dihedral=args.signed_dihedral,
max_multiplicity=max_mult_cap,
)
except RuntimeError as exc:
log.error("%s", exc)
return 1
if (na, nb, nc) == (1, 1, 1):
print(f" Optimal supercell : 1×1×1 (input cell already sufficient)")
else:
print(f" Optimal supercell : {na}×{nb}×{nc}")
print(f" (defined for {cutoff_str} — change --cutoff to alter the J set)")
else:
na, nb, nc = 1, 1, 1
if (na, nb, nc) != (1, 1, 1):
print(f" New cell : a={crystal.cell[0]:.4f} b={crystal.cell[1]:.4f} "
f"c={crystal.cell[2]:.4f} Å")
print(f" Magnetic atoms in supercell: {len(crystal.target_species)}")
# ---- Energy-mapping branch (alternative to the four-state method) --------
# Everything below this point is four-state-specific (reference bath,
# dimers, isolation, uu/dd/ud/du configs). Energy-mapping instead
# enumerates every inequivalent collinear order of the chosen (super)cell.
if args.method == "energy-mapping":
return _run_energy_mapping(args, crystal, couplings, (na, nb, nc),
moment_magnitude)
if args.method == "16-state":
return _run_ring16(args, crystal, couplings, (na, nb, nc),
moment_magnitude, base_crystal=base_crystal,
safe_tol=safe_tol)
# ---- Four-state isolation check ------------------------------------------
# Algebraic criterion: for each chosen Jk dimer (i, j), no periodic image of
# j may couple to i through a *different* Jm — then the four-state
# combination returns only Jk (with multiplicity M = #Jk images, divided out
# as 4·M). See mag4.supercell.check_four_state_isolation.
# NB: the four-state isolation check is moved BELOW (after the
# greedy dimer selection in Step 6), so it can use the actual
# dimer vectors mag4 picks rather than a generic best-effort
# search in the conv cell.
atoms = crystal.target_species
n = len(atoms)
log.info("Magnetic atoms in working cell: %d", n)
# ---- Step 3: reference bath ----------------------------------------------
ref_opts = [args.ref_state, args.ref_couplings, args.ref_kvector,
args.ref_magmom]
if sum(o is not None for o in ref_opts) > 1:
log.error("--ref-state, --ref-couplings, --ref-kvector and "
"--ref-magmom are mutually exclusive.")
return 1
ref_suggestion = None # (na,nb,nc) of a larger cell hosting the ref order
ref_sugg_frac = None # broken-bond fraction it achieves (0 = full)
if args.ref_magmom is not None:
# Bath given directly as the supercell MAGMOM string (signs used) —
# POSCAR atom order, so the user can paste the intended INCAR MAGMOM.
try:
ref_bath = build_reference_from_magmom(args.ref_magmom, crystal)
except ValueError as e:
log.error("%s", e)
return 1
log.info("Reference bath (--ref-magmom): %s",
" ".join(f"{s:+d}" for s in ref_bath))
elif args.ref_kvector is not None:
# Build a collinear bath from a propagation vector.
try:
kfrac = parse_kvector(args.ref_kvector)
except ValueError as e:
log.error("Invalid --ref-kvector: %s", e)
return 1
ref_bath, ref_report = build_reference_from_kvector(
atoms, kfrac, (na, nb, nc))
log.info("Reference bath (--ref-kvector, %s): %s",
ref_report["method"],
" ".join(f"{s:+d}" for s in ref_bath))
if ref_report["incommensurate"]:
sug = ref_report["suggest_supercell"]
log.warning("Cell is INCOMMENSURATE with k along axis/axes %s — the "
"bath will have a seam. Re-run with a commensurate "
"supercell, e.g. --supercell %d %d %d.",
", ".join(ref_report["incommensurate"]), *sug)
if ref_report["node_atoms"]:
log.warning("%d atom(s) sit on a node of cos(2π k·R) (sign "
"arbitrary, set +1): %s",
len(ref_report["node_atoms"]),
", ".join(ref_report["node_atoms"][:8])
+ (" …" if len(ref_report["node_atoms"]) > 8 else ""))
elif args.ref_couplings is not None:
# Build an AFM/FM bath that respects the requested coupling signs.
try:
ref_bath, ref_report = build_reference_from_couplings(
atoms, crystal.lattice_matrix, couplings,
args.ref_couplings, safe_tol)
except ValueError as e:
log.error("Invalid --ref-couplings: %s", e)
return 1
log.info("Reference bath (--ref-couplings, %s): %s",
ref_report["method"],
" ".join(f"{s:+d}" for s in ref_bath))
for name, (sat, total) in ref_report["couplings"].items():
sign = "AFM" if ref_report["wanted"][name] == -1 else "FM"
if total == 0:
log.warning(" %s (%s): no bonds found in this cell.", name, sign)
elif sat == total:
log.info(" %s (%s): %d/%d bonds satisfied.", name, sign, sat, total)
else:
log.warning(" %s (%s): only %d/%d bonds satisfied (frustrated).",
name, sign, sat, total)
conflict = ref_report["conflict"]
if conflict is not None:
kind = conflict[0]
if kind == "self-image":
log.warning("Frustration: atom %d couples AFM to its own periodic "
"image — the cell is too small/incommensurate along "
"that direction. Enlarge the supercell.", conflict[1])
else:
log.warning("Frustration: the requested coupling signs cannot all "
"be satisfied simultaneously (odd cycle). Used the "
"lowest-energy collinear compromise instead.")
# If the reference is frustrated, look for a larger ISOLATING cell that
# hosts the requested order with fewer broken bonds (e.g. an even cell
# for a Néel AFM) and suggest it at the end of the report.
cur_broken = sum(total - sat
for sat, total in ref_report["couplings"].values())
cur_total = sum(total for _, total in ref_report["couplings"].values())
if cur_broken > 0 and cur_total > 0:
cap = max(args.max_supercell, max(na, nb, nc) + 1)
ref_suggestion, ref_sugg_frac = _suggest_ref_supercell(
base_crystal, conv_lat_mat, conv_target_species, couplings,
args.ref_couplings, safe_tol, (na, nb, nc), cap,
cur_broken / cur_total, max_multiplicity=max_mult_cap)
else:
user_ref = None
if args.ref_state is not None:
try:
user_ref = parse_ref_state(args.ref_state)
except (ValueError, argparse.ArgumentTypeError) as e:
log.error("Invalid --ref-state: %s", e)
return 1
try:
ref_bath = build_reference_bath(n, user_ref)
except ValueError as e:
log.error("%s", e)
return 1
log.info("Reference bath (%s): %s",
"user-defined" if user_ref else "FM (all +1)",
" ".join(f"{s:+d}" for s in ref_bath))
# ---- Step 4: dimers -------------------------------------------------------
log.info("Searching for dimers (tol=±%.5f Å) …", safe_tol)
sc_struct = None
if ligand_elements:
from mag4.geometry import crystal_to_pymatgen
sc_struct = crystal_to_pymatgen(crystal)
dimers = find_dimers(
atoms, crystal.lattice_matrix, couplings, safe_tol,
structure=sc_struct,
ligand_elements=ligand_elements,
bond_cutoff=args.bond_cutoff,
angle_tol=args.angle_tol,
ll_cutoff=eff_ll_cutoff,
signed_dihedral=args.signed_dihedral,
)
# ---- Step 5: symmetry -----------------------------------------------------
if args.no_symmetry:
log.info("Crystal symmetry disabled — using time-reversal only.")
sym_ops = [tuple(range(n))]
else:
log.info("Computing crystal symmetry operations (full crystal) …")
mag_cart = get_cartesian_coords(atoms, crystal.lattice_matrix)
sym_ops = find_symmetry_ops(mag_cart, crystal.lattice_vectors, tol=0.05,
crystal=crystal)
if ligand_elements and sc_struct is not None:
from mag4.geometry import filter_symops_by_fingerprint
kept = filter_symops_by_fingerprint(
sym_ops, sc_struct, n_mag=n,
ligand_elements=ligand_elements,
bond_cutoff=args.bond_cutoff,
angle_tol=args.angle_tol,
ll_cutoff=eff_ll_cutoff,
signed_dihedral=args.signed_dihedral,
)
log.info("Bridging-ligand filter kept %d / %d symmetry op(s) "
"(dropped ops that mix sub-shells with different "
"M-L-M fingerprints).",
len(kept), len(sym_ops))
sym_ops = kept
# ---- Step 6: greedy unique-dimer assignment ------------------------------
W = 72
used_pairs: set = set()
all_configs: List[FourStateConfig] = []
chosen_dimers: dict = {} # coup_name → Dimer (for isolation check below)
for entry in couplings:
coup_name = entry[0]
candidates = dimers.get(coup_name, [])
if not candidates:
log.warning("No dimer found for %s — skipping!", coup_name)
continue
chosen = None
for d in candidates:
key = (d.idx_a, d.idx_b)
if key not in used_pairs:
chosen = d
break
if chosen is None:
print()
print("!" * W)
print(f" ERROR: {coup_name} has no unclaimed dimer pair.")
print(f" All {len(candidates)} candidate pair(s) are already used")
print(f" by an earlier coupling. The cell is too small to give")
print(f" {coup_name} a distinct representative dimer.")
print(f" Expand the supercell so more distinct pairs appear.")
print("!" * W)
print()
continue
used_pairs.add((chosen.idx_a, chosen.idx_b))
chosen_dimers[coup_name] = chosen
log.info("Selected dimer for %s: %s -- %s (d=%.4f Å)",
coup_name, chosen.label_a, chosen.label_b, chosen.distance)
all_configs.extend(generate_four_states(ref_bath, chosen, coup_name))
if not all_configs:
log.error("No configurations generated. Check coupling distances and --tol.")
return 1
# ---- Step 6b: four-state isolation check ---------------------------------
# For each chosen dimer, verify that the supercell stretches at
# least 3·d along the dimer's direction. Uses the *actual* dimer
# vector (Ni1_000 -- Ni1_001 for J2 axial in NiO, etc.), not a
# generic conv-cell search, so the per-J direction matches what
# mag4 will write into ``four_state_report.dat``.
def _dimer_vector(d) -> Optional[np.ndarray]:
"""Cartesian r_b − r_a in the supercell for a chosen ``Dimer``."""
try:
atom_a = crystal.target_species[d.idx_a]
atom_b = crystal.target_species[d.idx_b]
except (IndexError, AttributeError):
return None
frac_a = np.array([float(atom_a[1]), float(atom_a[2]),
float(atom_a[3])])
frac_b = np.array([float(atom_b[1]), float(atom_b[2]),
float(atom_b[3])])
# Take the shortest periodic image of the (b − a) vector under
# the supercell lattice — that's the *physical* dimer.
sc_mat = np.asarray(crystal.lattice_matrix, dtype=float)
dfrac = frac_b - frac_a
# Wrap each component to (−0.5, +0.5]:
dfrac -= np.round(dfrac)
v = sc_mat @ dfrac
return v
iso_couplings = [(name, chosen_dimers[name].distance)
for name, _ in [(c[0], c[1]) for c in couplings]
if name in chosen_dimers]
iso_dimer_vecs = {name: _dimer_vector(chosen_dimers[name])
for name, _ in iso_couplings}
# Multiplicity of each CHOSEN dimer: how many Jk bonds link the selected pair
# across supercell images. The four-state self-coefficient is 4·M (the
# report's divisor), and it is exactly what the isolation cap
# (--max-multiplicity / --full-isolation) must be measured against — so
# compute it here, hand it to the isolation check, and reuse it for the
# report below. Without it the check would judge the best-case symmetry-
# equivalent direction, which need not be the one greedily selected.
_sc_lat = np.asarray(crystal.lattice_matrix, dtype=float)
_shift3 = [_sc_lat @ np.array([float(i), float(j), float(k)])
for i in (-1, 0, 1) for j in (-1, 0, 1) for k in (-1, 0, 1)]
multiplicities: dict = {}
for name, dim in chosen_dimers.items():
ra = _sc_lat @ np.array([float(atoms[dim.idx_a][c]) for c in (1, 2, 3)])
rb = _sc_lat @ np.array([float(atoms[dim.idx_b][c]) for c in (1, 2, 3)])
m = sum(1 for T in _shift3
if abs(float(np.linalg.norm(rb + T - ra)) - dim.distance) < 0.05)
multiplicities[name] = max(1, m)
iso = check_four_state_isolation(
conv_lat_mat, (na, nb, nc), iso_couplings,
target_species=conv_target_species, # enumerate ALL bond directions
dimer_vectors=iso_dimer_vecs,
max_multiplicity=max_mult_cap,
dimer_multiplicity=multiplicities,
)
print("")
print(f" Four-state isolation check — supercell axes: "
f"a={iso['axes'][0]:.4f}, b={iso['axes'][1]:.4f}, "
f"c={iso['axes'][2]:.4f} Å")
d_max = iso.get("d_max", 0.0)
by_name = {s["name"]: s for s in iso["per_shell"]}
for shell in iso["per_shell"]:
status = "OK" if shell["ok"] else "FAIL"
u = shell["dimer_dir"]
u_str = (f"u=({u[0]:+.3f}, {u[1]:+.3f}, {u[2]:+.3f})"
if u is not None else "u=(direction unknown)")
if shell["ok"]:
mult = shell["multiplicity"]
extra = f", multiplicity {mult}" if (mult and mult > 1) else ""
why = f" only Jk between the dimer atoms{extra}"
elif shell["cross"]:
why = (" ✗ cross-contamination by " + ", ".join(shell["cross"])
+ " (a different J links the dimer atoms)")
else:
# multiplicity cap: clean but the dimer wraps onto its own image.
why = (f" ✗ multiplicity {shell['multiplicity']} > {max_mult_cap} "
f"(dimer wraps onto its own image; exceeds {cap_flag})")
print(f" {shell['name']}: d = {shell['d']:.4f} Å, dimer {u_str} "
f"[{status}]{why}")
if not iso["ok"]:
# The four-state method is only valid for an isolated dimer, so an
# under-isolated cell is a hard error — never silently generated.
rec = iso["recommended"]
found = tuple(rec) != (na, nb, nc)
bad = [s for s in iso["per_shell"] if not s["ok"]]
contaminated = [s for s in bad if s["cross"]]
msg_lines = [
"",
"!" * 78,
f" FOUR-STATE ISOLATION ERROR — supercell {na}×{nb}×{nc} does not "
f"isolate every J shell; no four-state inputs were written",
"!" * 78,
]
if contaminated:
msg_lines.append(
" Cross-contamination — a DIFFERENT coupling links the dimer "
"atoms, so the four-state combination would mix couplings:")
for s in bad:
if s["cross"]:
msg_lines.append(
f" - {s['name']} (d = {s['d']:.4f} Å): also linked by "
+ ", ".join(s["cross"]))
else:
# The multiplicity cap rejected an otherwise-clean shell.
msg_lines.append(
f" - {s['name']} (d = {s['d']:.4f} Å): multiplicity "
f"{s['multiplicity']} > {max_mult_cap} (dimer wraps onto its "
f"own image; exceeds {cap_flag})")
if found:
msg_lines += [
"",
" Smallest supercell that isolates every shell:",
f" --supercell {rec[0]} {rec[1]} {rec[2]}",
"",
" Or lower --cutoff so the long-range shells are dropped.",
]
else:
msg_lines += [
"",
" No isolating supercell was found near the requested one. "
"Raise --max-supercell to search larger cells, or lower "
"--cutoff so the long-range shells are dropped.",
]
msg_lines += ["!" * 78, ""]
log.error("\n".join(msg_lines))
log.error("four-state isolation criterion violated for %d shell(s): %s",
len(iso["failing"]), ", ".join(f[0] for f in iso["failing"]))
return 2
log.info("Generated %d configurations across %d couplings.",
len(all_configs), len(couplings))
# ---- Step 7: unique configurations ---------------------------------------
# --check-degeneracy: group by the UNITARY ops only, so configurations the
# antiunitary part (unitary × time reversal) proves degenerate stay
# separate DFT runs; the recorded pairs let mag4-extract verify |ΔE| ≈ 0.
unique_groups = find_unique_configs(
all_configs, sym_ops, time_reversal=not args.check_degeneracy)
tr_pairs = (time_reversal_pairs(unique_groups, sym_ops)
if args.check_degeneracy else None)
log.info("Unique configurations: %d (out of %d total)%s.",
len(unique_groups), len(all_configs),
(f"; {len(tr_pairs)} time-reversal pair(s) kept for the "
"degeneracy check") if args.check_degeneracy else "")
# Magnetic (Shubnikov) group of each unique configuration — BNS number +
# type, primed operations included. Informational; cheap (spglib per
# config on the working supercell).
from mag4.wien2k import analyze_config_magnetic_group
config_msgs = [analyze_config_magnetic_group(
crystal, group, spin_magnitude=moment_magnitude)
for group in unique_groups]
# ---- Step 8: report (written inside the compound dir) -------------------
# In dry-run mode we still print everything but never create files.
compound = (os.path.splitext(os.path.basename(args.cif))[0]
if args.cif else args.name)
if args.dry_run:
output_file = None
log.info("--dry-run active: no directory or file will be written.")
else:
os.makedirs(compound, exist_ok=True)
output_file = os.path.join(compound, f"{args.name}_report.dat")
# Pre-compute per-config reduced cells when the WIEN2k symmetric path
# is enabled with cell reduction. The metadata travels into the report
# so `mag4-extract` can rescale energies for J extraction.
per_config_cells = None
if (args.code in ("wien2k", "both")
and args.gen_struct and args.symmetrize_struct
and args.reduce_cell):
from mag4.geometry import reduce_to_primitive
from mag4.wien2k import (analyze_config_symmetry,
_wien2k_lattice_letter_from_sg,
_wien2k_centering_factor)
per_config_cells = []
for group in unique_groups:
red = reduce_to_primitive(crystal, group,
spin_magnitude=moment_magnitude,
symprec=args.symprec)
sym = analyze_config_symmetry(red, group,
spin_magnitude=moment_magnitude,
symprec=args.symprec)
n_at_conv = sum(len(g) for g in red.all_species)
# WIEN2k's ``:ENE`` is the total energy of the cell that
# contains only the atoms ACTUALLY LISTED in the struct
# file -- for an F/B/I/C/R-centred lattice letter, that
# is the conv-cell count divided by the centring factor
# (F→4, B/I/CXY/CYZ/CXZ→2, R→3, P/H→1). For NiO config1
# (Fm-3m, F lattice) the conv cell has 8 atoms but WIEN2k
# only integrates over 2 (1 Ni + 1 O listed); recording
# the conv count makes ``mag4-extract``'s LCM scaling
# off by a factor of 4 and the resulting J value blows
# up by ~10^6.
lat_letter = _wien2k_lattice_letter_from_sg(
sym.spacegroup_symbol, red.cell
)
centering = _wien2k_centering_factor(lat_letter)
n_at_wien2k = n_at_conv // centering
from pymatgen.core import Lattice
vol = Lattice.from_parameters(*red.cell).volume
per_config_cells.append(dict(
cell=tuple(red.cell),
n_atoms=n_at_wien2k, # what WIEN2k :ENE refers to
n_atoms_conv=n_at_conv, # the conv-cell count (info)
lattice_letter=lat_letter, # F / B / CXY / P / …
volume_A3=vol,
sg_number=sym.spacegroup_number,
sg_symbol=sym.spacegroup_symbol,
))
# Total atom count in the working supercell (all elements).
# ``mag4-extract`` uses this as the reference cell to which each
# config's WIEN2k :ENE is rescaled before applying the four-state
# formula — the formula is *defined* on the supercell, and DFT
# extensivity lets us scale energies from any reduced cell up to
# this size.
n_supercell_atoms = sum(len(g) for g in crystal.all_species)
# ``multiplicities`` (how many Jk bonds link each chosen pair across images;
# the report's 4·M divisor) was computed above, before the isolation check,
# so the cap and the formula divisor agree.
# Per-config energy equations in terms of the couplings: for each unique
# config, c_k = Σ_{⟨ij⟩∈Jk} s_i s_j, so E[#i] = E0 + Σ_k c_k·(Jk·S²). Lets
# the user CHECK that the extraction formulas isolate each Jk (the other
# couplings' coefficients cancel in the four-state combination).
em_couplings = [(c[0], c[1]) for c in couplings]
config_coefficients = config_energy_coefficients(
atoms, crystal.lattice_matrix, em_couplings,
[g[0].spin_vector for g in unique_groups], tol=safe_tol)
# --jring: also compute each config's four-spin RING coefficient so the
# report can show the Jring contribution to the energies AND how much it
# LEAKS into each bilinear J formula (the four-state method has no ring
# term, so a non-zero leak silently biases J — this just makes it visible).
ring_coefficients = None
ring_edge = None
# The Jring check is a free hint (no extra DFT, no change to the generated
# inputs), so it runs by default; --no-jring sets args.jring to False.
jring_auto = args.jring is None
if args.jring is not False:
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, ring_edge = select_ring_plaquettes(
em_couplings, cell_atoms, image_atoms,
prefer_label=args.jring_coupling, prefer_dist=args.jring_dist)
if not plaq:
if jring_auto:
log.info("Jring check (default-on): no square plaquette in "
"this lattice — nothing to report.")
else:
log.warning("--jring: no square plaquette found for the ring "
"term; skipping the Jring contamination report.")
else:
names = [a[0] for a in cell_atoms]
ring_coefficients = [
count_ring(plaq, {names[i]: g[0].spin_vector[i]
for i in range(len(names))})
for g in unique_groups]
# --jprime: detect BEYOND-CUTOFF bilinear shells J' and compute each config's
# coefficient for them, so the report can show their contribution to the
# configN energies AND how much each leaks into the extracted Jk (the
# four-state combination only cancels couplings that do NOT fold onto a
# self-image; a J' with d > --cutoff can survive and bias Jk in a small cell).
spectator_shells = None
spectator_coefficients = None
# Like the Jring check, the beyond-cutoff J' check is a free hint and runs
# by default (it needs --cif; --no-jprime disables it, an explicit
# --jprime [RADIUS] keeps forcing it).
jprime_auto = args.jprime is None and not args.no_jprime
if args.jprime is not None or jprime_auto:
if not args.cif:
if jprime_auto:
log.info("J' check (default-on) skipped: needs --cif to "
"detect beyond-cutoff shells.")
else:
log.warning("--jprime needs --cif to auto-detect beyond-cutoff "
"shells; skipping the J' contamination report.")
else:
lat = np.asarray(crystal.lattice_matrix, dtype=float)
# Default scan radius: 2×--cutoff (P. Blaha's suggestion), but at
# least out to the shortest lattice translation — the first shells
# that fold onto a self-image are the ones that can leak, so a
# small cell must still reach them.
radius = (args.jprime if args.jprime and args.jprime > 0.0
else max(shortest_periodic_image(lat), 2.0 * args.cutoff))
results = get_mm_distances(args.cif, element, radius)
shells = group_inequivalent(results, args.equiv_tol)
spec = [(g["label"], g["distance"]) for g in shells
if g["distance"] > args.cutoff + args.equiv_tol]
if not spec:
if jprime_auto:
log.info("J' check (default-on): no bilinear shell between "
"--cutoff (%.4f Å) and %.4f Å — nothing beyond "
"the model to check.", args.cutoff, radius)
else:
log.warning("--jprime: no bilinear shell between --cutoff "
"(%.4f Å) and %.4f Å; nothing to check.",
args.cutoff, radius)
else:
spectator_shells = spec
spectator_coefficients = config_energy_coefficients(
atoms, crystal.lattice_matrix, spec,
[g[0].spin_vector for g in unique_groups], tol=safe_tol)
print_report(all_configs, unique_groups, atoms, couplings, ref_bath,
len(sym_ops), (na, nb, nc), output_file,
per_config_cells=per_config_cells,
n_supercell_atoms=n_supercell_atoms,
multiplicities=multiplicities,
config_coefficients=config_coefficients,
ring_coefficients=ring_coefficients,
ring_edge=ring_edge,
spectator_shells=spectator_shells,
spectator_coefficients=spectator_coefficients,
tr_pairs=tr_pairs,
config_msgs=config_msgs)
# ---- Step 9: MAGMOM appendix ---------------------------------------------
magmom_lines = build_magmom_lines(unique_groups, crystal, moment_magnitude)
spin_header = (f" MAGMOM (S = {args.spin:g}, |MAGMOM| = 2·S = "
f"{moment_magnitude:g} μ_B)")
if output_file is not None:
with open(output_file, "a") as f:
f.write("\n")
f.write("=" * 80 + "\n")
f.write(spin_header + "\n")
f.write("=" * 80 + "\n")
for uid, (group, mm) in enumerate(zip(unique_groups, magmom_lines), 1):
member_strs = [f"{m.coupling_name}({m.state_label})" for m in group]
f.write(f" Config #{uid} [{', '.join(member_strs)}]\n")
f.write(f" MAGMOM = {mm}\n\n")
print()
print("=" * 80)
print(spin_header)
print("=" * 80)
for uid, (group, mm) in enumerate(zip(unique_groups, magmom_lines), 1):
member_strs = [f"{m.coupling_name}({m.state_label})" for m in group]
print(f" Config #{uid} [{', '.join(member_strs)}]")
print(f" MAGMOM = {mm}")
print()
# Reference-state tip — printed in BOTH dry-run and write modes (before the
# dry-run early return below), so tuning --cutoff with --dry-run still shows
# the larger commensurate cell that hosts the requested order.
if ref_suggestion is not None and tuple(ref_suggestion) != (na, nb, nc):
s = ref_suggestion
print("")
print("-" * 78)
print(f" Reference-state tip: the requested order is frustrated in this "
f"{na}×{nb}×{nc} cell (isolation-minimal).")
if ref_sugg_frac is not None and ref_sugg_frac <= 1e-9:
print(f" This larger commensurate supercell hosts it FULLY (every "
f"reference bond satisfied) — consider:")
else:
print(f" A larger commensurate supercell hosts it with fewer broken "
f"bonds — consider:")
print(f" --supercell {s[0]} {s[1]} {s[2]}")
print("-" * 78)
# ---- Step 10: per-config directories ------------------------------------
if args.dry_run:
print("=" * 80)
print(" --dry-run : skipping directory creation.")
print(" Re-run without --dry-run (and without changing flags) to write")
print(f" {compound}/four_state_report.dat and the configN/ tree.")
print("=" * 80)
return 0
write_vasp = args.code in ("vasp", "both")
write_wien = args.code in ("wien2k", "both")
# POSCAR is always written (also needed by WIEN2k workflow via xyz2struct).
# INCAR + KPOINTS only if a VASP run is wanted. ``--functional PBE+U``
# / ``--ldaul`` / ``--ldauu`` are validated up at the top of main()
# so a bad combination fails before any heavy work.
create_config_dirs(unique_groups, crystal, magmom_lines,
args.name, args.kdens, compound,
write_vasp_extras=write_vasp,
kspacing=args.kspacing,
functional=args.functional,
magnetic_element=args.element,
ldaul=args.ldaul,
ldauu=args.ldauu,
ldauj=args.ldauj,
aexx=args.aexx,
encut=args.encut,
metagga=args.metagga,
lmaxtau=args.lmaxtau)
if write_wien:
use_symmetric = args.gen_struct and args.symmetrize_struct
if use_symmetric:
try:
write_symmetric_case_files(
unique_groups, crystal,
parent_dir=compound,
instgen_path=args.instgen,
spin_magnitude=moment_magnitude,
symprec=args.symprec,
reduce_cell=args.reduce_cell,
)
except RuntimeError as exc:
log.error("WIEN2k symmetric generation failed: %s", exc)
return 1
else:
try:
write_case_inst_files(unique_groups, crystal,
parent_dir=compound,
instgen_path=args.instgen,
spin_magnitude=moment_magnitude)
except RuntimeError as exc:
log.error("WIEN2k case.inst generation failed: %s", exc)
return 1
if args.gen_struct:
try:
write_case_struct_files(crystal, len(unique_groups),
compound)
except (ValueError, OSError) as exc:
log.error("WIEN2k configN.struct generation failed: %s", exc)
return 1
if args.gen_job_scripts:
try:
write_job_scripts(compound, functional=args.functional,
kspacing=args.kspacing,
element=args.element, ldaul=args.ldaul,
ldauu=args.ldauu, ldauj=args.ldauj)
except OSError as exc:
log.error("WIEN2k job script generation failed: %s", exc)
return 1
return 0
if __name__ == "__main__": # pragma: no cover
sys.exit(main())