"""WIEN2k input generation: ``configN.inst`` + ``configN.struct`` per four-state
spin configuration.
Layer: **WIEN2k-specific** (one of the two DFT back-ends; see also ``vasp.py``).
In each ``configN/`` directory mag4 writes:
* ``POSCAR`` – same VASP file used by ``--code vasp`` (kept for
cross-checking; also accepted as input by WIEN2k's ``xyz2struct``).
* ``configN.inst`` – produced by piping the ``u``/``d``/``n`` spin pattern
into WIEN2k's ``instgen -ask`` script. Named after the directory so the
``case == basename(pwd)`` convention works out of the box.
* ``configN.struct`` – written directly by Python (no ``xyz2struct`` needed).
Every atom is emitted as inequivalent (``MULT=1``) with a unique label of
the form ``<element><n>`` (``"Cs1"``, ``"V 9"``, ``"I17"`` …) so WIEN2k
cannot merge atoms via symmetry analysis. Only the identity symmetry op
is written. This guarantees the atom order in ``configN.struct`` matches
the order assumed by ``configN.inst``.
Mapping convention
------------------
Each entry of the spin pattern maps as follows::
+1 (or any positive value) → ``u``
-1 (or any negative value) → ``d``
0 (non-magnetic atom) → ``n``
So a MAGMOM line ``-7 2*7 -7 3*7 4*0`` becomes ``d u u d u u u n n n n``.
"""
from __future__ import annotations
import logging
import os
import re
import shutil
import subprocess
import tempfile
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Sequence, Union
import numpy as np
from mag4.constants import CrystalData, FourStateConfig
from mag4.energy import RY_TO_EV
#: Ångström → Bohr conversion factor used by WIEN2k struct files.
ANG_TO_BOHR: float = 1.0 / 0.5291772083
#: The 73 symmorphic space groups (no glide planes / no screw axes).
#: Used by the ``--no-reduce`` branch in
#: :func:`write_symmetric_case_files` to decide whether to refine
#: orbit groupings (`_cell_compatible_symmetry` + `_relabel_…`) or
#: to fall back to "every atom MULT=1" (`_ungrouped_symmetry`).
#:
#: For symmorphic parents, the refined point group at origin is closed
#: under composition modulo the supercell lattice, so the union-find
#: in `_cell_compatible_symmetry` produces correct
#: supercell-as-P-Bravais orbits. For non-symmorphic parents the
#: parent's intrinsic glide / screw translations alias with
#: supercell-induced internal translations and the union-find
#: over-merges (e.g. MULT=16 in a Pmna 2×2×1 supercell where
#: Pmna's max Wyckoff multiplicity is 8).
SYMMORPHIC_SG_NUMBERS: frozenset = frozenset({
# Triclinic, monoclinic
1, 2, 3, 5, 6, 8, 10, 12,
# Orthorhombic
16, 21, 22, 23, 25, 35, 38, 42, 44, 47, 65, 69, 71,
# Tetragonal
75, 79, 81, 82, 83, 87, 89, 97, 99, 107, 111, 115, 119, 121,
123, 139,
# Trigonal / rhombohedral / hexagonal
143, 146, 147, 148, 149, 150, 155, 156, 157, 160, 162, 164,
166, 168, 174, 175, 177, 183, 187, 189, 191,
# Cubic
195, 196, 197, 200, 202, 204, 207, 209, 211, 215, 216, 217,
221, 225, 229,
})
log = logging.getLogger(__name__)
# ---------------------------------------------------------------------------
# Spin → u/d/n translation
# ---------------------------------------------------------------------------
[docs]
def magmom_to_udn(values: Sequence[float]) -> List[str]:
"""Translate a flat list of MAGMOM values to the ``u/d/n`` letters.
Parameters
----------
values : sequence of float
One value per atom in POSCAR order. Sign carries the spin direction;
zero means non-magnetic.
Returns
-------
list of str
One letter per atom: ``"u"`` for positive, ``"d"`` for negative,
``"n"`` for zero.
Examples
--------
>>> magmom_to_udn([-7, 7, 7, -7, 7, 7, 7, 0, 0, 0, 0])
['d', 'u', 'u', 'd', 'u', 'u', 'u', 'n', 'n', 'n', 'n']
"""
out: List[str] = []
for v in values:
if v > 0:
out.append("u")
elif v < 0:
out.append("d")
else:
out.append("n")
return out
[docs]
def build_udn_lines(unique_groups: List[List[FourStateConfig]],
crystal: CrystalData,
spin_magnitude: float = 1.0) -> List[List[str]]:
"""Per-unique-config list of u/d/n letters in POSCAR order.
Returns a list of length ``len(unique_groups)``; each element is the list
of u/d/n letters for one configuration. Order matches
:func:`mag4.vasp.build_magmom_lines`.
"""
all_atoms = [atom for group in crystal.all_species for atom in group]
mag_label_to_idx = {atom[0]: i for i, atom in enumerate(crystal.target_species)}
out: List[List[str]] = []
for group in unique_groups:
sv = group[0].spin_vector
vals = []
for atom in all_atoms:
lbl = atom[0]
if lbl in mag_label_to_idx:
vals.append(sv[mag_label_to_idx[lbl]] * spin_magnitude)
else:
vals.append(0.0)
out.append(magmom_to_udn(vals))
return out
# ---------------------------------------------------------------------------
# Stub case.struct (only enough for instgen to discover the atom list)
# ---------------------------------------------------------------------------
[docs]
def wien2k_label(element_symbol: str, idx: int) -> str:
"""Format an atom label per WIEN2k case.struct convention.
* Two-letter element: e.g. ``"Cs1"`` → cols 1-2 = ``"Cs"``.
* One-letter element: ``"V 1"`` (letter + space) → cols 1-2 = ``"V "``.
The convention matters because :program:`instgen` reads the element from
the first two characters of each ``RMT=`` line in ``case.struct``.
"""
if len(element_symbol) == 1:
return f"{element_symbol} {idx}"
return f"{element_symbol}{idx}"
[docs]
def write_stub_struct(crystal: CrystalData, path: str) -> None:
"""Write a minimal ``case.struct`` good enough for :program:`instgen`.
Only the per-atom ``RMT=`` line is emitted (one per atom in POSCAR order)
because that is the only thing :program:`instgen` parses. Use
:func:`write_struct_file` to produce a full, ready-to-use struct file.
"""
lines: List[str] = []
for group_idx, group in enumerate(crystal.all_species):
elem = crystal.elements[group_idx]
for k in range(1, len(group) + 1):
label = wien2k_label(elem, k)
lines.append(
f"{label:<10s} NPT= 781 R0=0.00010000 RMT= 2.30000 Z: 0.0"
)
with open(path, "w") as f:
f.write("\n".join(lines) + "\n")
# ---------------------------------------------------------------------------
# Full case.struct writer (replaces the xyz2struct step)
# ---------------------------------------------------------------------------
def _strip_to_element(label: str) -> str:
"""Return just the chemical-symbol prefix of a mag4 atom label.
Examples
--------
>>> _strip_to_element("Cs1")
'Cs'
>>> _strip_to_element("V1_010")
'V'
>>> _strip_to_element("I12_002")
'I'
"""
m = re.match(r"^([A-Za-z]+)", label)
return m.group(1) if m else "X"
def _atomic_number(symbol: str) -> int:
"""Return Z for a chemical symbol via pymatgen (cached)."""
from pymatgen.core import Element
return int(Element(symbol).Z)
def _default_r0(z: int) -> float:
"""A safe per-Z default for the radial mesh first point R0 (bohr).
Heavier elements need a smaller R0 (steeper inner cusp).
"""
if z >= 36: # Kr and heavier
return 0.00001
if z >= 18: # Ar–Kr
return 0.00005
return 0.0001 # H–Ar
def _is_hexagonal_cell(cell: Sequence[float]) -> bool:
"""``True`` for a genuine hexagonal Bravais cell: α=β=90°, γ=120° **and
a=b**.
The ``a=b`` test matters for the four-state supercells: a 1×2×1 expansion of
a hexagonal cell keeps the angles (γ=120°) but doubles ``b``, so ``a≠b`` and
the cell is no longer hexagonal — WIEN2k's ``H`` lattice (which assumes
a=b) would misread the geometry. Such a cell must be ``P`` instead.
"""
a, b, _c, al, be, ga = cell
return (abs(al - 90) < 1e-3 and abs(be - 90) < 1e-3
and abs(ga - 120) < 1e-3
and abs(a - b) < 1e-3 * max(abs(a), abs(b), 1.0))
def _detect_lattice_letter(cell: Sequence[float]) -> str:
"""Pick the WIEN2k single-letter lattice code from the cell geometry.
The four-state pipeline always works in a primitive supercell with
explicit atom positions, so we never need centred lattices. The only
helpful distinction is a genuine hexagonal cell (α=β=90, γ=120, **a=b**) →
``H`` versus everything else → ``P``. A hexagonal cell expanded
anisotropically (e.g. 1×2×1 → a≠b) is therefore ``P``, not ``H``.
"""
return "H" if _is_hexagonal_cell(cell) else "P"
def _wien2k_lattice_letter_from_sg(sg_symbol: str, cell: Sequence[float]) -> str:
"""Map a spglib international SG symbol to the WIEN2k lattice letter.
WIEN2k's naming differs from the International Tables in three places
(WIEN2k user's guide, Table 4.4):
* International "I" (body-centred) → WIEN2k ``B``
* International "A"/"B"/"C" (base-centred, orthorhombic & monoclinic)
→ WIEN2k ``CYZ`` / ``CXZ`` / ``CXY`` (single letter "C" is not used)
* **A genuine hexagonal cell (γ=120° AND a=b) → WIEN2k ``H``**, regardless
of the International symbol's first character. Trigonal SGs such as
``P-3m1`` (#164) and hexagonal SGs such as ``P6_3/mmc`` (#194)
both start with ``"P"`` but their WIEN2k lattice letter is
``H``. Without this override WIEN2k's ``sgroup`` rewrites the
struct, which then trips ``mag4-extract``'s scaling because the
effective SCF cell differs from what mag4 recorded. A hexagonal cell
expanded to ``a≠b`` (e.g. a 1×2×1 supercell) is **not** hexagonal and
maps via the symbol instead (→ ``P``).
Falls back to :func:`_detect_lattice_letter` when the symbol is
empty or carries a character outside the recognised set — so the
legacy non-symmetric code path keeps its current behaviour.
"""
if not sg_symbol:
return _detect_lattice_letter(cell)
# A genuine hexagonal cell (γ=120° AND a=b) ⇒ "H" regardless of the
# SG-symbol prefix. An anisotropically expanded one (a≠b) is NOT
# hexagonal — fall through to the symbol mapping (→ "P" etc.).
if _is_hexagonal_cell(cell):
return "H"
first = sg_symbol[0].upper()
mapping = {
"P": "P",
"F": "F",
"I": "B", # body-centred → WIEN2k 'B'
"R": "R",
"H": "H",
"A": "CYZ", # A-base-centred (yz plane)
"B": "CXZ", # B-base-centred (xz plane) — incl. monoclinic
"C": "CXY", # C-base-centred (xy plane)
}
return mapping.get(first, _detect_lattice_letter(cell))
def _wien2k_centering_factor(lattice_letter: str) -> int:
"""How many lattice points the WIEN2k centering letter generates per
conventional cell — used to fold orbits by centering equivalence.
P, H → 1 (primitive)
F → 4 (face-centred)
B → 2 (body-centred, = international "I")
CXY, CYZ, CXZ → 2 (base-centred)
R → 3 (rhombohedral in hexagonal axes)
"""
return {
"P": 1, "H": 1,
"F": 4,
"B": 2,
"CXY": 2, "CYZ": 2, "CXZ": 2,
"R": 3,
}.get(lattice_letter, 1)
# ---------------------------------------------------------------------------
# Per-configuration symmetry analysis (spglib)
# ---------------------------------------------------------------------------
[docs]
@dataclass
class ConfigSymmetry:
"""Symmetry analysis of a single spin-decorated configuration.
The cell is kept verbatim from the input :class:`CrystalData`; only the
atom ordering (within the cell) and the symmetry-op list are derived.
Attributes
----------
spacegroup_number, spacegroup_symbol
Space group reported by spglib for the spin-decorated structure.
rotations, translations
Symmetry operations in the *input* cell basis, ready to be written
verbatim into a WIEN2k struct file (``rotations`` is integer; the
translations are in fractional coordinates).
origin_shift
Origin shift that spglib would apply to put the structure into the
standard setting. Kept for diagnostics; **not** applied to atom
positions here — we keep input positions so the rotations/translations
in the input basis act directly.
new_to_old
``new_to_old[i]`` is the original POSCAR index of the atom that now
sits at position ``i`` in the symmetrised order. Length = number
of atoms.
orbits
List of Wyckoff orbits in the **new** order: each orbit is a list of
consecutive new-order indices that share a symmetry equivalent class.
``len(orbits)`` is the number of inequivalent sites (== WIEN2k's
``NEQUIV``).
frac_coords_new
Fractional coordinates in the new order (shape ``(n_atoms, 3)``).
"""
spacegroup_number: int
spacegroup_symbol: str
rotations: np.ndarray
translations: np.ndarray
origin_shift: np.ndarray
new_to_old: List[int]
orbits: List[List[int]]
frac_coords_new: np.ndarray
#: For each atom in the **new** order, the index of the corresponding
#: atom in spglib's primitive cell. Used by the WIEN2k struct writer
#: to fold orbits by lattice-centering equivalence (only one rep per
#: primitive index is emitted under an F/B/C-centred letter, matching
#: WIEN2k's "list one of four / one of two" convention).
mapping_to_primitive: List[int] = field(default_factory=list)
@property
def n_atoms(self) -> int:
return len(self.new_to_old)
@property
def n_inequiv(self) -> int:
return len(self.orbits)
def _elements_per_atom(crystal: CrystalData) -> List[str]:
"""Chemical-symbol string per atom in POSCAR (``all_species``) order."""
out: List[str] = []
for grp_idx, grp in enumerate(crystal.all_species):
elem = crystal.elements[grp_idx]
out.extend([elem] * len(grp))
return out
def _frac_coords_array(crystal: CrystalData) -> np.ndarray:
"""``(n_atoms, 3)`` array of fractional coords in POSCAR order."""
rows: List[List[float]] = []
for grp in crystal.all_species:
for atom in grp:
rows.append([float(atom[1]), float(atom[2]), float(atom[3])])
return np.asarray(rows, dtype=float)
def _spin_per_atom(crystal: CrystalData,
group: List[FourStateConfig],
spin_magnitude: float) -> List[float]:
"""Signed spin magnitude per atom in POSCAR order; 0.0 for non-magnetic."""
all_atoms = [a for grp in crystal.all_species for a in grp]
mag_label_to_idx = {a[0]: i for i, a in enumerate(crystal.target_species)}
sv = group[0].spin_vector
out: List[float] = []
for atom in all_atoms:
lbl = atom[0]
if lbl in mag_label_to_idx:
out.append(float(sv[mag_label_to_idx[lbl]]) * spin_magnitude)
else:
out.append(0.0)
return out
def _spglib_atom_types(crystal: CrystalData,
group: List[FourStateConfig],
spin_magnitude: float) -> List[int]:
"""Atom-type code per atom, suitable for spglib's coloring.
Encoding:
* non-magnetic atom of element X → ``Z(X)``
* magnetic atom of element X with spin > 0 → ``100000 + 10*Z(X) + 1``
* magnetic atom of element X with spin < 0 → ``100000 + 10*Z(X) + 2``
The offset guarantees no collision between the magnetic and bare branches
(atomic numbers stop at 118). Two atoms share a type iff they should be
considered identical for symmetry purposes — a magnetic atom is *never*
equivalent to a bare atom of the same element under this coloring, but
a non-magnetic atom with spin = 0 of the same element ``X`` shares the
bare ``Z(X)`` type with any other zero-spin X atom (per the spec).
"""
spins = _spin_per_atom(crystal, group, spin_magnitude)
elems = _elements_per_atom(crystal)
out: List[int] = []
for elem, s in zip(elems, spins):
try:
z = _atomic_number(elem)
except Exception: # pragma: no cover - pymatgen unavailable / unknown sym
z = 0
if s > 0:
out.append(100000 + 10 * z + 1)
elif s < 0:
out.append(100000 + 10 * z + 2)
else:
out.append(z)
return out
def _fallback_p1_symmetry(crystal: CrystalData) -> ConfigSymmetry:
"""Identity-only P1 fallback used when spglib finds no symmetry."""
frac = _frac_coords_array(crystal)
n = frac.shape[0]
return ConfigSymmetry(
spacegroup_number=1,
spacegroup_symbol="P1",
rotations=np.eye(3, dtype=int).reshape(1, 3, 3),
translations=np.zeros((1, 3), dtype=float),
origin_shift=np.zeros(3, dtype=float),
new_to_old=list(range(n)),
orbits=[[i] for i in range(n)],
frac_coords_new=frac,
)
def _ungrouped_symmetry(config_sym: ConfigSymmetry) -> ConfigSymmetry:
"""Each atom is its own ``MULT=1`` orbit.
Used by the ``--no-reduce`` branch of
:func:`write_symmetric_case_files` when the colored parent SG is
**non-symmorphic** (Cmce, Pmna, P4_2/mnm, …). Combined with the
shared-label scheme produced by :func:`_spin_aware_orbit_labels`,
the struct lists every supercell atom on its own ATOM line but
re-uses one label per ``(element, spin)`` class (``Cu1`` for all
UP Cu, ``Cu2`` for all DOWN Cu, ``La1`` for all La, ``O 1`` for
all O, …). WIEN2k's SYMMETRY then groups atoms sharing a label
into proper Wyckoff orbits via ``struct_sgroup``.
"""
n = config_sym.n_atoms
return ConfigSymmetry(
spacegroup_number=config_sym.spacegroup_number,
spacegroup_symbol=config_sym.spacegroup_symbol,
rotations=config_sym.rotations,
translations=config_sym.translations,
origin_shift=config_sym.origin_shift,
new_to_old=list(config_sym.new_to_old),
orbits=[[i] for i in range(n)],
frac_coords_new=config_sym.frac_coords_new,
mapping_to_primitive=list(config_sym.mapping_to_primitive),
)
def _spin_aware_orbit_labels(
config_sym: ConfigSymmetry,
crystal: CrystalData,
group: List[FourStateConfig],
spin_magnitude: float = 1.0,
) -> List[str]:
"""Per-orbit WIEN2k label that encodes ``(element, spin)``.
Returns a list ``labels`` of length ``len(config_sym.orbits)``
where each entry is a label like ``"Cu1"`` (UP magnetic),
``"Cu2"`` (DOWN magnetic) or ``"La1"`` (non-magnetic). The
suffix is fixed by the magnetic class, not by an
orbit-position counter, so multiple ATOM blocks in the resulting
struct share the same label whenever they refer to atoms of the
same element AND the same spin direction.
The writer (:func:`write_struct_file_symmetric`) consumes this
list verbatim instead of auto-generating per-block labels — that
way each non-symmorphic ``--no-reduce`` config produces, in
effect, a "list-each-atom-MULT=1-but-group-by-label" struct.
WIEN2k's ``symmetry`` will then merge same-label atoms at
SG-equivalent positions into the correct Wyckoff orbits in
``struct_sgroup``.
"""
elem_per_atom = _elements_per_atom(crystal)
spin_per_atom = _spin_per_atom(crystal, group, spin_magnitude)
labels: List[str] = []
for orbit in config_sym.orbits:
rep_new = orbit[0]
rep_old = config_sym.new_to_old[rep_new]
elem = elem_per_atom[rep_old]
s = spin_per_atom[rep_old]
if abs(s) <= 1e-6:
suffix = "1" # non-magnetic
elif s > 0:
suffix = "1" # UP magnetic
else:
suffix = "2" # DOWN magnetic
labels.append(f"{elem:<2s}{suffix}")
return labels
def _cell_compatible_symmetry(config_sym: ConfigSymmetry) -> ConfigSymmetry:
"""Re-derive orbits keeping ONLY the symmetry operations compatible
with the working supercell as a primitive Bravais lattice.
spglib reports the *parent* SG of the colored structure (e.g. Fm-3m
for an all-FM NiO 1×1×2 supercell), so the ops list includes
"internal centring" translations like ``(½,½,0)``, ``(0,0,½)``,
``(½,½,½)`` that are NOT lattice vectors of the supercell viewed as
P-Bravais. Under those over-symmetric ops, all 8 Ni atoms of the
NiO supercell collapse into a single MULT=8 orbit — which is
inconsistent with the tetragonal cell shape WIEN2k actually sees.
This helper keeps only the operations whose translation component
is zero modulo the input cell (= "point group at the origin"),
re-derives Wyckoff orbits under those operations via union-find,
and reorders atoms so each orbit's members sit contiguously. The
resulting ``ConfigSymmetry`` reports the SAME parent SG (kept as
metadata) but its ``orbits`` field gives the grouping compatible
with the cell — e.g. tetragonal 4/mmm for an all-FM NiO 1×1×2,
where the 8 Ni split into 5 orbits (sizes 1, 1, 4, 1, 1) instead
of one cubic orbit of size 8.
"""
rots = np.asarray(config_sym.rotations, dtype=int)
trans = np.asarray(config_sym.translations, dtype=float)
n_ops = len(rots)
if n_ops == 0:
return config_sym
# Identify *pure* internal translations of the SG: ops with R = I
# and t ≠ 0 (mod 1). In the standard Bravais types these
# translations are part of the centring (e.g. (½,½,0) for F or
# (½,½,½) for I). But when the working supercell is *not* the
# standardised conventional cell of the colored SG — exactly the
# ``--no-reduce`` case — the SG includes *extra* translations
# generated by the supercell-vs-conv-cell mismatch (e.g. (0,0,½)
# along the doubled c-axis of NiO 1×1×2). These extra translations
# are *not* lattice vectors of the supercell viewed as a P-Bravais
# crystal — but they ARE in spglib's op list, so propagating those
# ops would over-symmetrise the supercell back to the parent SG.
I3 = np.eye(3, dtype=int)
internal_translations: List[np.ndarray] = []
for i in range(n_ops):
if np.array_equal(rots[i], I3):
t_mod = np.mod(trans[i] + 0.5, 1.0) - 0.5
if np.max(np.abs(t_mod)) >= 1e-3:
internal_translations.append(t_mod)
# For each op (R, t), drop duplicates that differ from an earlier
# kept op by an internal translation — this is the right way to
# "factor out" the centring without assuming spglib's origin is
# at mag4's coordinate origin. When ``internal_translations`` is
# empty (primitive Bravais already) every op is kept.
keep_idx: List[int] = []
for i in range(n_ops):
dup = False
for j in keep_idx:
if not np.array_equal(rots[j], rots[i]):
continue
dt = trans[j] - trans[i]
dt_mod = np.mod(dt + 0.5, 1.0) - 0.5
# Either zero (already kept) or an internal translation.
candidates = [np.zeros(3)] + internal_translations
for pt in candidates:
diff = np.mod(dt_mod - pt + 0.5, 1.0) - 0.5
if np.max(np.abs(diff)) < 1e-3:
dup = True
break
if dup:
break
if not dup:
keep_idx.append(i)
pure_rots = rots[keep_idx]
pure_trans = trans[keep_idx]
frac = np.asarray(config_sym.frac_coords_new, dtype=float)
n_atoms = frac.shape[0]
if n_atoms == 0:
return config_sym
# Original orbit each atom belongs to — point-group ops should not
# merge atoms across distinct colored orbits (different element /
# spin types), so we keep this as a sanity gate.
orig_orbit_of = [0] * n_atoms
for orig_idx, orbit in enumerate(config_sym.orbits):
for a in orbit:
orig_orbit_of[a] = orig_idx
# Union-find over atoms: i and j unioned iff some pure point op
# maps i to j (modulo cell) AND they share the same parent-SG
# orbit (= same colored type).
parent = list(range(n_atoms))
def find(x: int) -> int:
while parent[x] != x:
parent[x] = parent[parent[x]]
x = parent[x]
return x
def union(a: int, b: int) -> None:
ra, rb = find(a), find(b)
if ra != rb:
parent[ra] = rb
for R, t in zip(pure_rots, pure_trans):
for i in range(n_atoms):
new_pos = np.mod(R @ frac[i] + t, 1.0)
for j in range(n_atoms):
if orig_orbit_of[j] != orig_orbit_of[i]:
continue
diff = np.mod(frac[j] - new_pos + 0.5, 1.0) - 0.5
if np.max(np.abs(diff)) < 1e-3:
union(i, j)
break
# Group atoms by their root and order orbits by their lowest index
# so the resulting struct is deterministic.
by_root: Dict[int, List[int]] = {}
for i in range(n_atoms):
by_root.setdefault(find(i), []).append(i)
new_orbits_old_order = [sorted(v) for v in by_root.values()]
new_orbits_old_order.sort(key=lambda o: o[0])
# Re-flatten atoms in orbit order so analyze_config_symmetry's
# contract (atoms in the same orbit sit contiguously) holds.
new_order: List[int] = []
for orbit in new_orbits_old_order:
new_order.extend(orbit)
pos_map = {old: new for new, old in enumerate(new_order)}
new_orbits = [[pos_map[i] for i in o] for o in new_orbits_old_order]
new_frac = frac[new_order]
new_to_old_old = config_sym.new_to_old
new_to_old_combined = [new_to_old_old[i] for i in new_order]
map_prim_old = config_sym.mapping_to_primitive
new_map_prim = (
[map_prim_old[i] for i in new_order] if map_prim_old else []
)
return ConfigSymmetry(
spacegroup_number=config_sym.spacegroup_number,
spacegroup_symbol=config_sym.spacegroup_symbol,
rotations=pure_rots,
translations=pure_trans,
origin_shift=config_sym.origin_shift,
new_to_old=new_to_old_combined,
orbits=new_orbits,
frac_coords_new=new_frac,
mapping_to_primitive=new_map_prim,
)
def _relabel_cell_compatible_sg(
crystal: CrystalData,
refined_sym: ConfigSymmetry,
symprec: float = 1e-3,
) -> ConfigSymmetry:
"""Replace ``spacegroup_number`` / ``spacegroup_symbol`` of a refined
``ConfigSymmetry`` with the actual cell-compatible SG.
Strategy: re-color each atom by its **refined** orbit and ask
spglib for the SG. By giving every refined orbit a distinct type
we forbid spglib's internal-centring-related ops (which would
cross orbit boundaries), so the detected SG is exactly the
maximal cell-compatible one — e.g. P4/mmm (#123) for NiO 1×1×2
all-FM where ``_cell_compatible_symmetry`` has already collapsed
the parent Fm-3m's 8 internal translations.
"""
try:
import spglib # noqa: WPS433
except ImportError: # pragma: no cover - spglib ships with pymatgen
return refined_sym
lat = np.asarray(crystal.lattice_vectors, dtype=float)
frac = _frac_coords_array(crystal)
n_total = frac.shape[0]
types = [0] * n_total
for orbit_idx, orbit in enumerate(refined_sym.orbits, start=1):
for new_idx in orbit:
old_idx = refined_sym.new_to_old[new_idx]
types[old_idx] = orbit_idx
try:
ds = spglib.get_symmetry_dataset((lat, frac, types),
symprec=symprec)
except Exception: # pragma: no cover - spglib internal failure
return refined_sym
if ds is None:
return refined_sym
from mag4._spglib_compat import ds_get
return ConfigSymmetry(
spacegroup_number=int(ds_get(ds, "number")),
spacegroup_symbol=str(ds_get(ds, "international")),
rotations=refined_sym.rotations,
translations=refined_sym.translations,
origin_shift=refined_sym.origin_shift,
new_to_old=refined_sym.new_to_old,
orbits=refined_sym.orbits,
frac_coords_new=refined_sym.frac_coords_new,
mapping_to_primitive=refined_sym.mapping_to_primitive,
)
[docs]
@dataclass
class MagneticGroupInfo:
"""Magnetic (Shubnikov) space group of one spin-decorated configuration.
spglib classifies the configuration (collinear moments included, so the
PRIMED operations — unitary op × time reversal — are counted) into one of
the 1651 magnetic space-group types. spglib provides the standard
numbers but not the symbol string; the BNS number is the citable label
(symbol lookup: Bilbao MGENPOS / ISO-MAG).
``msg_type`` is the Shubnikov type: 1 = colorless (no primed ops),
2 = grey (contains pure time reversal — only for zero-moment structures),
3 = black-white (half the point operations primed), 4 = black-white with
anti-translations (a translation × T maps the sublattices).
"""
uni_number: int # 1..1651 (UNI / BNS+OG unified index)
bns_number: str # e.g. "65.490"
og_number: str # e.g. "69.6.610"
msg_type: int # Shubnikov type 1..4
parent_number: int # ordinary space group of the parent (family)
parent_symbol: str # e.g. "Cmmm"
n_unitary: int # unprimed operations
n_primed: int # operations combined with time reversal
mpg_symbol: str = "" # magnetic POINT group, see mpg_kind
mpg_kind: str = "" # "colourless" | "grey" | "black-white"
@property
def type_roman(self) -> str:
return {1: "I", 2: "II", 3: "III", 4: "IV"}.get(self.msg_type, "?")
[docs]
def short(self) -> str:
"""Compact display string, e.g. ``BNS 65.490 (type IV, parent Cmmm)``."""
return (f"BNS {self.bns_number} (type {self.type_roman}, "
f"parent {self.parent_symbol})")
[docs]
def point_group(self) -> str:
"""Magnetic point group with its kind, e.g. ``mmm1' (grey)``.
Dropping the translations from the magnetic space group gives the
magnetic point group exactly: an anti-translation (type IV) becomes
pure time reversal 1', so the point group is GREY ``G1'``; a type-III
group keeps half its point operations primed and is written in the
Shubnikov ``G(H)`` notation (parent, unitary subgroup); a type-I
group is the ordinary (colourless) symbol.
"""
return (f"{self.mpg_symbol} ({self.mpg_kind})"
if self.mpg_symbol else "")
[docs]
def analyze_config_magnetic_group(
crystal: CrystalData,
group: List[FourStateConfig],
*,
spin_magnitude: float = 1.0,
symprec: float = 1e-3,
) -> Optional[MagneticGroupInfo]:
"""Magnetic space group (Shubnikov/BNS) of a spin-decorated configuration.
Unlike :func:`analyze_config_symmetry` — which colours up/down spins as
different species and therefore sees only the UNITARY part of the
magnetic group — this passes the collinear moments to spglib's magnetic
dataset, so operations combined with time reversal (primed) are included.
Returns ``None`` when spglib is too old (< 2.0) or fails; callers should
then simply omit the magnetic-group line.
"""
try:
import spglib # noqa: WPS433 — runtime dep via pymatgen
except ImportError: # pragma: no cover
return None
if not hasattr(spglib, "get_magnetic_symmetry_dataset"):
log.warning("spglib %s has no magnetic dataset support (needs >= 2.0);"
" magnetic-group labels omitted.",
getattr(spglib, "__version__", "?"))
return None
lat = np.asarray(crystal.lattice_vectors, dtype=float)
frac = _frac_coords_array(crystal)
elems = _elements_per_atom(crystal)
try:
numbers = [_atomic_number(e) for e in elems]
except Exception: # pragma: no cover - unknown element symbol
return None
magmoms = np.asarray(_spin_per_atom(crystal, group, spin_magnitude),
dtype=float)
def _get(ds, key):
return getattr(ds, key) if hasattr(ds, key) else ds[key]
try:
ds = spglib.get_magnetic_symmetry_dataset(
(lat, frac, numbers, magmoms), symprec=symprec)
if ds is None:
return None
uni = int(_get(ds, "uni_number"))
trs = np.asarray(_get(ds, "time_reversals")).astype(bool)
t = spglib.get_magnetic_spacegroup_type(uni)
bns = str(_get(t, "bns_number"))
og = str(_get(t, "og_number"))
msg_type = int(_get(t, "type"))
parent = int(_get(t, "number"))
try:
from pymatgen.symmetry.groups import SpaceGroup
parent_symbol = SpaceGroup.from_int_number(parent).symbol
except Exception: # pragma: no cover - pymatgen lookup failure
parent_symbol = f"#{parent}"
# Magnetic POINT group: drop translations, dedupe (rotation, θ) pairs.
# A primed identity (an anti-translation modded out) means pure time
# reversal survives → GREY group G1'. Otherwise: no primed ops →
# colourless G; primed ops but no pure 1' → black-white, written in
# the exact Shubnikov G(H) notation (H = unitary subgroup).
mpg_symbol, mpg_kind = "", ""
try:
rots = np.asarray(_get(ds, "rotations"), dtype=int)
pairs = {}
for W, is_primed in zip(rots, trs):
key = (W.tobytes(), bool(is_primed))
pairs.setdefault(key, W)
all_W = np.array([pairs[k] for k in pairs])
uni_W = np.array([pairs[k] for k in pairs if not k[1]])
has_pure_tr = any(k[1] and np.array_equal(
np.frombuffer(k[0], dtype=int).reshape(3, 3), np.eye(3))
for k in pairs)
g_sym = str(spglib.get_pointgroup(all_W)[0]).strip()
if not trs.any():
mpg_symbol, mpg_kind = g_sym, "colourless"
elif has_pure_tr:
mpg_symbol, mpg_kind = f"{g_sym}1'", "grey"
else:
h_sym = str(spglib.get_pointgroup(uni_W)[0]).strip()
mpg_symbol, mpg_kind = f"{g_sym}({h_sym})", "black-white"
except Exception: # pragma: no cover - point-group naming failure
pass
return MagneticGroupInfo(
uni_number=uni, bns_number=bns, og_number=og, msg_type=msg_type,
parent_number=parent, parent_symbol=parent_symbol,
n_unitary=int((~trs).sum()), n_primed=int(trs.sum()),
mpg_symbol=mpg_symbol, mpg_kind=mpg_kind)
except Exception as exc: # pragma: no cover - spglib internal failure
log.warning("magnetic-group analysis failed (%s); label omitted.", exc)
return None
[docs]
def analyze_config_symmetry(
crystal: CrystalData,
group: List[FourStateConfig],
*,
spin_magnitude: float = 1.0,
symprec: float = 1e-3,
) -> ConfigSymmetry:
"""Run spglib on the spin-decorated configuration and pack the result.
Atoms in :class:`CrystalData` are coloured according to their spin
(see :func:`_spglib_atom_types`), so the detected space group is the
*magnetic* space group of the configuration projected onto the user's
cell — equivalent to "highest symmetry compatible with this cell".
The returned :class:`ConfigSymmetry` carries:
* the rotations/translations in the input cell basis (so we can write
them straight into a WIEN2k struct file);
* the permutation that groups atoms by their Wyckoff orbit, with orbits
ordered by the index of their representative — this becomes the new
atom order on disk.
On any spglib error or empty dataset we fall back to a P1
(identity-only) :class:`ConfigSymmetry` so the rest of the pipeline can
proceed unchanged.
"""
try:
import spglib # noqa: WPS433 — runtime dep via pymatgen
except ImportError: # pragma: no cover - spglib ships with pymatgen
log.warning("spglib unavailable; falling back to P1 for this config")
return _fallback_p1_symmetry(crystal)
lat = np.asarray(crystal.lattice_vectors, dtype=float)
frac = _frac_coords_array(crystal)
types = _spglib_atom_types(crystal, group, spin_magnitude)
try:
ds = spglib.get_symmetry_dataset((lat, frac, types), symprec=symprec)
except Exception as exc: # pragma: no cover - spglib internal failure
log.warning("spglib failed (%s); falling back to P1", exc)
return _fallback_p1_symmetry(crystal)
if ds is None:
log.warning("spglib returned no dataset; falling back to P1")
return _fallback_p1_symmetry(crystal)
from mag4._spglib_compat import ds_get
equiv = np.asarray(ds_get(ds, "equivalent_atoms"), dtype=int)
n = len(equiv)
# Group atoms by their representative index, deterministic order.
orbit_map: Dict[int, List[int]] = {}
for i, rep in enumerate(equiv):
orbit_map.setdefault(int(rep), []).append(i)
orbit_reps = sorted(orbit_map.keys())
new_to_old: List[int] = []
for r in orbit_reps:
new_to_old.extend(orbit_map[r])
old_to_new = [0] * n
for new_i, old_i in enumerate(new_to_old):
old_to_new[old_i] = new_i
orbits_new = [[old_to_new[a] for a in orbit_map[r]] for r in orbit_reps]
frac_new = frac[new_to_old]
# mapping_to_primitive: for each atom in the *input* (= new) order, the
# index of its corresponding atom in spglib's primitive cell. Atoms
# related by lattice-centering translations share the same value; we
# use that to emit only one atom per centering coset in F/B/C-centred
# WIEN2k structs.
mapping_raw = np.asarray(
ds_get(ds, "mapping_to_primitive", default=[]), dtype=int
)
if mapping_raw.size == n:
mapping_new = [int(mapping_raw[old_i]) for old_i in new_to_old]
else:
# spglib didn't return the field for some reason — fall back to
# identity (each atom in its own coset; no centering folding).
mapping_new = list(range(n))
return ConfigSymmetry(
spacegroup_number=int(ds_get(ds, "number")),
spacegroup_symbol=str(ds_get(ds, "international")),
rotations=np.asarray(ds_get(ds, "rotations"), dtype=int),
translations=np.asarray(ds_get(ds, "translations"), dtype=float),
origin_shift=np.asarray(ds_get(ds, "origin_shift"), dtype=float),
new_to_old=new_to_old,
orbits=orbits_new,
frac_coords_new=frac_new,
mapping_to_primitive=mapping_new,
)
[docs]
def write_struct_file(
crystal: CrystalData,
path: str,
*,
title: Optional[str] = None,
rmt: Union[float, Dict[str, float]] = 2.0,
npt: int = 781,
r0: Optional[Union[float, Dict[str, float]]] = None,
calc_mode: str = "RELA",
lattice: Optional[str] = None,
) -> None:
"""Write a complete WIEN2k ``case.struct`` from a :class:`CrystalData`.
Every atom is emitted as a separate inequivalent site (``MULT=1``) with a
unique label of the form ``<element><n>`` (``"Cs1"``, ``"V 9"``,
``"I17"`` …). Only the identity symmetry operation is written, so
WIEN2k will not merge or reorder atoms. This guarantees that the atom
order in ``configN.struct`` matches the order assumed by
``configN.inst``.
See :func:`write_struct_file_symmetric` for the spin-aware variant that
uses real-space symmetry (default in the CLI; this function is kept as
a fallback behind ``--no-symmetrize-struct``).
Parameters
----------
crystal
The (super-)cell as produced by :mod:`mag4.supercell`.
path
Where to write the struct file (full path including filename).
title
Free-text title written on line 1. Defaults to a mag4-tagged note.
rmt
Muffin-tin radius in bohr. Either a single value applied to every
element (default 2.0) or a dict ``{element_symbol: value}``.
npt
Number of radial mesh points (default 781; must be odd).
r0
First radial mesh point in bohr. ``None`` (default) → picked
per-Z via :func:`_default_r0`; a single float applies to all
elements; a dict allows per-element overrides.
calc_mode
``"RELA"`` (relativistic, default) or ``"NREL"`` (non-relativistic).
lattice
WIEN2k single-letter lattice code. Auto-detected from the cell
angles when ``None``.
Notes
-----
The per-atom header line follows the format declared in the WIEN2k user
guide: ``(A10, 5X, I5, 5X, F10.8, 5X, F10.5, 5X, F10.5)`` for
``name, NPT, R0, RMT, Z`` respectively. The 5X gaps carry decorative
text (``" NPT="``, ``" R0="``, ``" RMT="``, ``" Z:"``) for
readability — WIEN2k skips them when parsing.
"""
if npt % 2 == 0:
raise ValueError(f"NPT must be odd, got {npt}")
atoms_flat: list = [a for grp in crystal.all_species for a in grp]
n_atoms = len(atoms_flat)
if n_atoms == 0:
raise ValueError("Cannot write struct file: crystal has no atoms")
if n_atoms > 999:
raise ValueError(
f"This generator uses I3 for NONEQUIV.ATOMS and I4 for atom "
f"indices, so it caps at 999 atoms. Got {n_atoms}."
)
a, b, c, alpha, beta, gamma = crystal.cell
a_b, b_b, c_b = a * ANG_TO_BOHR, b * ANG_TO_BOHR, c * ANG_TO_BOHR
if lattice is None:
lattice = _detect_lattice_letter(crystal.cell)
if title is None:
title = "generated by mag4 (every atom inequivalent, identity sym op)"
lines: List[str] = []
# Line 1: title (pad to 80 chars for cosmetics).
lines.append(f"{title:<80s}")
# Line 2: ``<H> LATTICE,NONEQUIV.ATOMS:<I3> <sg_num> <sg_sym>``
# — lattice letter in cols 1-4 (A4), 'LATTICE,NONEQUIV.ATOMS:' in
# cols 5-27 (A23), NONEQUIV count as I3 in cols 28-30, then the
# space-group number and Hermann-Mauguin symbol as free text.
# The legacy writer doesn't have access to the SG via
# :class:`ConfigSymmetry`, so we leave the SG fields blank --
# ``init_lapw`` re-derives them anyway.
lines.append(
f"{lattice:<4s}{'LATTICE,NONEQUIV.ATOMS:':<23s}{n_atoms:3d}"
)
# Line 3: ``MODE OF CALC=<RELA|NREL> unit=bohr`` — emitting just
# ``RELA`` in col 1-4 (the previous mag4 behaviour) silently makes
# WIEN2k treat the run as non-relativistic, giving huge wrong
# total-energy differences. WIEN2k's parser specifically looks
# for the ``MODE OF CALC=`` prefix and the ``unit=bohr`` suffix.
lines.append(f"MODE OF CALC={calc_mode} unit=bohr")
# Line 4: cell parameters (6F10.6) in bohr; angles stay in degrees.
lines.append(
f"{a_b:10.6f}{b_b:10.6f}{c_b:10.6f}"
f"{alpha:10.6f}{beta:10.6f}{gamma:10.6f}"
)
elem_counter: Dict[str, int] = {}
# We need atoms in the same global order as crystal.all_species,
# but the original ATOM block uses a single running 1..N index.
running_idx = 0
for group_idx, group in enumerate(crystal.all_species):
elem_symbol = crystal.elements[group_idx]
try:
z = _atomic_number(elem_symbol)
except Exception:
z = 0
# Resolve RMT and R0 for this element.
this_rmt = (
float(rmt[elem_symbol]) if isinstance(rmt, dict) and elem_symbol in rmt
else (float(rmt) if not isinstance(rmt, dict) else 2.0)
)
if r0 is None:
this_r0 = _default_r0(z)
elif isinstance(r0, dict):
this_r0 = float(r0.get(elem_symbol, _default_r0(z)))
else:
this_r0 = float(r0)
for atom in group:
running_idx += 1
xf, yf, zf = float(atom[1]), float(atom[2]), float(atom[3])
elem_counter[elem_symbol] = elem_counter.get(elem_symbol, 0) + 1
sub_idx = elem_counter[elem_symbol]
# Label: chemical symbol left-padded to 2 chars, then index.
# Examples: "Cs1", "V 9" (V + space + 9), "I17".
label = f"{elem_symbol:<2s}{sub_idx:<d}"
label_a10 = f"{label:<10s}"
# ATOM block header (negative index per WIEN2k convention).
lines.append(
f"ATOM{-running_idx:4d}: X={xf:10.8f} Y={yf:10.8f} Z={zf:10.8f}"
)
# MULT=1 for every atom so WIEN2k cannot merge sites.
lines.append(" MULT= 1 ISPLIT= 8")
# Atom parameter line — see format note in the docstring.
lines.append(
f"{label_a10} NPT={npt:5d} R0={this_r0:10.8f}"
f" RMT={this_rmt:10.5f} Z:{float(z):10.5f}"
)
# Local rotation matrix (identity).
lines.append("LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000")
lines.append(" 0.0000000 1.0000000 0.0000000")
lines.append(" 0.0000000 0.0000000 1.0000000")
# Per WIEN2k user's guide §4.3: nsym = 0 makes SYMMETRY derive the
# operations itself from the cell + atom positions + lattice letter.
# This legacy writer used to list the identity only (the spin-pattern
# atom labels already prevent WIEN2k from over-symmetrising), but
# zero-nsym is more idiomatic and matches the symmetric writer.
lines.append(f"{0:4d} NUMBER OF SYMMETRY OPERATIONS")
with open(path, "w") as f:
f.write("\n".join(lines) + "\n")
# ---------------------------------------------------------------------------
# Symmetric struct writer (one inequivalent block per Wyckoff orbit)
# ---------------------------------------------------------------------------
def _resolve_rmt(rmt: Union[float, Dict[str, float]], elem: str) -> float:
if isinstance(rmt, dict):
return float(rmt.get(elem, 2.0))
return float(rmt)
def _resolve_r0(r0: Optional[Union[float, Dict[str, float]]],
elem: str, z: int) -> float:
if r0 is None:
return _default_r0(z)
if isinstance(r0, dict):
return float(r0.get(elem, _default_r0(z)))
return float(r0)
def _format_sym_op(rot: np.ndarray, trans: np.ndarray, idx: int) -> List[str]:
"""Format one WIEN2k symmetry operation block (3 rotation rows + index)."""
out: List[str] = []
for row, t in zip(rot, trans):
r1, r2, r3 = int(row[0]), int(row[1]), int(row[2])
out.append(f"{r1:2d}{r2:2d}{r3:2d}{t:11.8f}")
out.append(f"{idx:8d}")
return out
[docs]
def write_struct_file_symmetric(
crystal: CrystalData,
config_sym: ConfigSymmetry,
path: str,
*,
title: Optional[str] = None,
rmt: Union[float, Dict[str, float]] = 2.0,
npt: int = 781,
r0: Optional[Union[float, Dict[str, float]]] = None,
calc_mode: str = "RELA",
lattice: Optional[str] = None,
labels: Optional[List[str]] = None,
) -> None:
"""Write a WIEN2k ``case.struct`` that exposes the configuration's
detected space-group symmetry.
One inequivalent block is emitted per Wyckoff orbit (``MULT`` = orbit
size), with the orbit's representative listed on the leading ``ATOM``
line and the remaining members on continuation lines between the
``MULT=`` line and the per-orbit ``NPT/R0/RMT`` line. Each orbit gets
a **unique** ``<element><n>`` label so WIEN2k cannot merge orbits that
happen to share an element. The trailing symmetry block lists every
rotation/translation reported by spglib in the input cell basis.
Atom positions are written in the order encoded by
``config_sym.frac_coords_new`` — i.e. spglib's equivalence classes
grouped together, with orbits ordered by their representative index.
"""
if npt % 2 == 0:
raise ValueError(f"NPT must be odd, got {npt}")
n_atoms = config_sym.n_atoms
if n_atoms == 0:
raise ValueError("Cannot write struct file: crystal has no atoms")
if n_atoms > 999:
raise ValueError(
f"This generator uses I3 for NONEQUIV.ATOMS and I4 for atom "
f"indices, so it caps at 999 atoms. Got {n_atoms}."
)
a, b, c, alpha, beta, gamma = crystal.cell
a_b, b_b, c_b = a * ANG_TO_BOHR, b * ANG_TO_BOHR, c * ANG_TO_BOHR
if lattice is None:
# Take the WIEN2k lattice letter from the colored SG symbol so
# F/B(=I)/CXY/CYZ/CXZ/R-centred groups are emitted in the cell
# setting WIEN2k's ``init_lapw`` expects.
lattice = _wien2k_lattice_letter_from_sg(
config_sym.spacegroup_symbol, crystal.cell)
centering = _wien2k_centering_factor(lattice)
# For centred lattices WIEN2k expects only one atom per centering
# coset in each Wyckoff orbit (the lattice letter implies the
# centering translations). Fold each orbit by mapping_to_primitive:
# atoms with the same primitive index sit in the same centering
# coset, so we keep one of them as the orbit representative.
folded_orbits: List[List[int]] = []
if centering > 1 and config_sym.mapping_to_primitive:
for orbit in config_sym.orbits:
seen_prim: Dict[int, int] = {}
for a in orbit:
pidx = config_sym.mapping_to_primitive[a]
seen_prim.setdefault(pidx, a)
# Preserve original orbit ordering as much as possible.
folded_orbits.append(
[seen_prim[pidx] for pidx in sorted(seen_prim,
key=lambda p: orbit.index(seen_prim[p]))]
)
else:
folded_orbits = [list(o) for o in config_sym.orbits]
n_inequiv_emitted = len(folded_orbits)
# spglib reports *all* space-group operations (point ops × centering
# translations). WIEN2k's struct file convention only stores the
# point-group ops — the centering translations are implicit in the
# lattice letter (F, B, CXY, …). So the count quoted in the title
# must be divided by the centering factor to match what
# ``init_lapw``'s SYMMETRY program will list after re-running.
n_ops_total = len(config_sym.rotations)
n_ops_wien2k = max(1, n_ops_total // centering)
if title is None:
title = (
f"generated by mag4 - SG #{config_sym.spacegroup_number} "
f"({config_sym.spacegroup_symbol}), "
f"{n_inequiv_emitted} inequiv site(s), "
f"{n_ops_wien2k} sym op(s)"
)
elem_per_atom_orig = _elements_per_atom(crystal)
elem_counter: Dict[str, int] = {}
lines: List[str] = []
lines.append(f"{title:<80s}")
# Line 2: lattice letter (A4) + ``LATTICE,NONEQUIV.ATOMS:`` (A23)
# + NONEQUIV (I3) + ``<sg_num> <sg_sym>``. The trailing SG
# number + Hermann-Mauguin symbol are the standard WIEN2k
# convention; ``init_lapw`` uses them as a hint for ``sgroup``.
lines.append(
f"{lattice:<4s}{'LATTICE,NONEQUIV.ATOMS:':<23s}"
f"{n_inequiv_emitted:3d}"
f" {config_sym.spacegroup_number} {config_sym.spacegroup_symbol}"
)
# Line 3: ``MODE OF CALC=<RELA|NREL> unit=bohr``. Bare ``RELA``
# at col 1-4 (mag4's previous behaviour) trips WIEN2k into running
# the non-relativistic branch, which silently corrupts every
# total energy — the four-state subtraction then produces J
# values that are orders of magnitude off.
lines.append(f"MODE OF CALC={calc_mode} unit=bohr")
lines.append(
f"{a_b:10.6f}{b_b:10.6f}{c_b:10.6f}"
f"{alpha:10.6f}{beta:10.6f}{gamma:10.6f}"
)
for orbit_idx, orbit in enumerate(folded_orbits):
mult = len(orbit)
# Representative element (orbit[0] is the first new-order index;
# map back through new_to_old to find the original POSCAR position
# and thus the chemical element).
rep_old = config_sym.new_to_old[orbit[0]]
elem = elem_per_atom_orig[rep_old]
try:
z = _atomic_number(elem)
except Exception:
z = 0
if labels is not None:
# Caller supplied a per-orbit label (e.g. ``Cu1`` for every
# UP magnetic Cu block) — re-use it verbatim and don't
# increment the per-element counter.
label = labels[orbit_idx]
else:
elem_counter[elem] = elem_counter.get(elem, 0) + 1
sub_idx = elem_counter[elem]
label = f"{elem:<2s}{sub_idx:<d}"
label_a10 = f"{label:<10s}"
this_rmt = _resolve_rmt(rmt, elem)
this_r0 = _resolve_r0(r0, elem, z)
iatnr = -(orbit_idx + 1)
# Leading ATOM line — first atom of the orbit.
x0, y0, z0 = config_sym.frac_coords_new[orbit[0]]
lines.append(
f"ATOM{iatnr:4d}: X={x0:10.8f} Y={y0:10.8f} Z={z0:10.8f}"
)
lines.append(f" MULT={mult:2d} ISPLIT= 8")
# Continuation lines for the remaining atoms in the orbit.
for k in orbit[1:]:
xf, yf, zf = config_sym.frac_coords_new[k]
lines.append(
f" {iatnr:4d}: X={xf:10.8f} Y={yf:10.8f} Z={zf:10.8f}"
)
lines.append(
f"{label_a10} NPT={npt:5d} R0={this_r0:10.8f}"
f" RMT={this_rmt:10.5f} Z:{float(z):10.5f}"
)
lines.append("LOCAL ROT MATRIX: 1.0000000 0.0000000 0.0000000")
lines.append(" 0.0000000 1.0000000 0.0000000")
lines.append(" 0.0000000 0.0000000 1.0000000")
# Per WIEN2k user's guide §4.3: writing "0 NUMBER OF SYMMETRY
# OPERATIONS" tells SYMMETRY to derive the full op list from the
# cell + atom positions + lattice letter. That's the most robust
# choice: mag4's coloured-SG ops are correct in the conventional
# cell we emit, but listing them explicitly has historically tripped
# ``init_lapw`` on F-centred groups (where sgroup rewrites the
# Bravais lattice and the listed ops then sit in a basis that no
# longer matches). Letting WIEN2k regenerate the ops sidesteps
# the whole class of basis-mismatch issues.
lines.append(f"{0:4d} NUMBER OF SYMMETRY OPERATIONS")
with open(path, "w") as f:
f.write("\n".join(lines) + "\n")
[docs]
def write_case_struct_files(
crystal: CrystalData,
n_configs: int,
parent_dir: str,
*,
rmt: Union[float, Dict[str, float]] = 2.0,
npt: int = 781,
r0: Optional[Union[float, Dict[str, float]]] = None,
) -> List[str]:
"""Write ``configN.struct`` into every ``configN/`` subdirectory.
The cell and atom list are identical across configurations; only the
spin pattern (in ``configN.inst``) differs. We therefore write the
same struct content N times, naming each file after its directory.
"""
written: List[str] = []
print()
print("=" * 80)
print(f" WIEN2k configN.struct files (no xyz2struct needed)")
print("=" * 80)
for cfg_idx in range(1, n_configs + 1):
cfg_name = f"config{cfg_idx}"
cfg_dir = os.path.join(parent_dir, cfg_name)
os.makedirs(cfg_dir, exist_ok=True)
dst = os.path.join(cfg_dir, f"{cfg_name}.struct")
write_struct_file(crystal, dst, rmt=rmt, npt=npt, r0=r0)
written.append(dst)
print(f" {dst}")
print(f" {len(written)} configN.struct file(s) written.")
print()
return written
# ---------------------------------------------------------------------------
# Symmetric orchestrator: per-config struct + inst with matching atom order
# ---------------------------------------------------------------------------
def _write_stub_struct_orbits(config_sym: ConfigSymmetry,
elem_per_atom_orig: Sequence[str],
path: str) -> None:
"""Stub ``case.struct`` with one ``RMT=`` line per orbit, in orbit order.
Used by :func:`_run_instgen_for_orbits` to drive instgen so that the
resulting ``case.inst`` has one entry per inequivalent site — matching
the per-orbit blocks emitted by :func:`write_struct_file_symmetric`.
"""
lines: List[str] = []
elem_counter: Dict[str, int] = {}
for orbit in config_sym.orbits:
rep_old = config_sym.new_to_old[orbit[0]]
elem = elem_per_atom_orig[rep_old]
elem_counter[elem] = elem_counter.get(elem, 0) + 1
label = wien2k_label(elem, elem_counter[elem])
lines.append(
f"{label:<10s} NPT= 781 R0=0.00010000 RMT= 2.30000 Z: 0.0"
)
with open(path, "w") as f:
f.write("\n".join(lines) + "\n")
def _orbit_udn_letters(crystal: CrystalData,
group: List[FourStateConfig],
config_sym: ConfigSymmetry,
spin_magnitude: float) -> List[str]:
"""One ``u``/``d``/``n`` letter per orbit (taken from the representative)."""
spins_orig = _spin_per_atom(crystal, group, spin_magnitude)
out: List[str] = []
for orbit in config_sym.orbits:
rep_old = config_sym.new_to_old[orbit[0]]
s = spins_orig[rep_old]
if s > 0:
out.append("u")
elif s < 0:
out.append("d")
else:
out.append("n")
return out
def _run_instgen_for_orbits(udn_per_orbit: Sequence[str],
crystal: CrystalData,
config_sym: ConfigSymmetry,
instgen_path: str) -> str:
"""Run :program:`instgen` once with one answer per orbit."""
tmp = tempfile.mkdtemp(prefix="mag4_instgen_")
try:
stub_struct = os.path.join(tmp, "stub.struct")
elem_per_atom = _elements_per_atom(crystal)
_write_stub_struct_orbits(config_sym, elem_per_atom, stub_struct)
proc = subprocess.run(
[instgen_path, "-f", "stub", "-ask"],
input="\n".join(udn_per_orbit) + "\n",
capture_output=True,
text=True,
cwd=tmp,
check=False,
)
inst_path = os.path.join(tmp, "stub.inst")
if proc.returncode != 0 or not os.path.isfile(inst_path):
raise RuntimeError(
f"instgen failed (exit {proc.returncode}).\n"
f" command: {instgen_path} -f stub -ask\n"
f" cwd : {tmp}\n"
f" stdout : {proc.stdout}\n"
f" stderr : {proc.stderr}"
)
with open(inst_path) as f:
return f.read()
finally:
shutil.rmtree(tmp, ignore_errors=True)
[docs]
def write_symmetric_case_files(
unique_groups: List[List[FourStateConfig]],
crystal: CrystalData,
parent_dir: str,
*,
instgen_path: Optional[str] = None,
spin_magnitude: float = 1.0,
symprec: float = 1e-3,
rmt: Union[float, Dict[str, float]] = 2.0,
npt: int = 781,
r0: Optional[Union[float, Dict[str, float]]] = None,
reduce_cell: bool = False,
) -> Dict[str, List]:
"""Write ``configN.struct`` and ``configN.inst`` per config with a
shared, symmetry-derived atom order.
For each configuration:
1. Run :func:`analyze_config_symmetry` once on the spin-decorated cell.
2. Write ``configN.struct`` via :func:`write_struct_file_symmetric` —
one block per Wyckoff orbit, full symmetry-op tail, unique
per-orbit labels.
3. Drive :program:`instgen` with one ``u/d/n`` answer per orbit
(the spin of the orbit's representative) and save the resulting
``configN.inst``.
Because both files are produced from the same :class:`ConfigSymmetry`,
their atom orders are guaranteed to agree without any post hoc
reordering.
Returns a dict with four keys: ``"struct"`` (list of struct paths),
``"inst"`` (list of inst paths), ``"sym"`` (list of
:class:`ConfigSymmetry` results, one per config — handy for tests),
and ``"crystals"`` (the per-config reduced :class:`CrystalData`
used to write each struct — used by the report writer to record
per-config cell metadata for :program:`mag4-extract`'s cell-size
rescaling).
When ``reduce_cell`` is true (``--reduce``; the default is off, so
every config keeps the one working supercell and shares identical
numerics) each config is rewritten in
its colored-structure primitive cell via
:func:`mag4.geometry.reduce_to_primitive` — multiplicities then
match the parent SG's standard Wyckoff table, and downstream DFT
runs are correspondingly cheaper for high-symmetry configs.
"""
from mag4.geometry import reduce_to_primitive
if instgen_path is None:
instgen_path = find_instgen()
if instgen_path is None:
raise RuntimeError(
"Could not locate 'instgen' on this machine. "
"Pass --instgen PATH or set $WIEN2k_INSTGEN / $WIENROOT."
)
log.info("Using instgen at: %s", instgen_path)
print()
print("=" * 80)
print(f" WIEN2k symmetric configN files (instgen = {instgen_path})")
print(f" symprec = {symprec:g} Å")
print("=" * 80)
structs: List[str] = []
insts: List[str] = []
syms: List[ConfigSymmetry] = []
crystals: List[CrystalData] = []
for cfg_idx, group in enumerate(unique_groups, start=1):
cfg_name = f"config{cfg_idx}"
cfg_dir = os.path.join(parent_dir, cfg_name)
os.makedirs(cfg_dir, exist_ok=True)
cfg_crystal = (
reduce_to_primitive(crystal, group,
spin_magnitude=spin_magnitude,
symprec=symprec)
if reduce_cell else crystal
)
crystals.append(cfg_crystal)
sym = analyze_config_symmetry(cfg_crystal, group,
spin_magnitude=spin_magnitude,
symprec=symprec)
# For ``--no-reduce`` mode, branch on whether the colored
# parent SG is symmorphic.
#
# Symmorphic parents (Fm-3m, Pm-3m, P4/mmm, I4/mmm, …): refine
# orbits via ``_cell_compatible_symmetry`` + relabel. This
# works because the parent's point ops have a ``t=0``
# representative for each rotation, hence form a closed
# point-group at origin in supercell coords. NiO 1×1×2
# all-FM keeps its [1, 4, 1, 1, 1] tetragonal Ni grouping and
# the relabelled SG title (P4/mmm).
#
# Non-symmorphic parents (Cmce, Pmna, P4_2/mnm, …): the
# parent's intrinsic glide / screw translations alias with
# supercell-induced internal translations (e.g. Pmna's
# ``t = (½, 0, 0)`` glide is a half-supercell-a translation
# in a 2×2×1 supercell), so ``_cell_compatible_symmetry``'s
# union-find over-merges orbits — La2CuO4 2×2×1 would emit
# MULT=16 in a struct that claims Pmna (max Wyckoff = 8).
# Fallback: every atom on its own ``MULT=1`` block, but with
# labels SHARED per ``(element, spin)`` class (``Cu1`` for all
# UP, ``Cu2`` for all DOWN, ``La1`` for all La, ``O 1`` for
# all O, …). WIEN2k's SYMMETRY uses positions (not labels)
# to derive orbits, then writes a canonical ``struct_sgroup``
# whose ATOM blocks merge same-label atoms at equivalent
# positions.
explicit_labels: Optional[List[str]] = None
if not reduce_cell:
if sym.spacegroup_number in SYMMORPHIC_SG_NUMBERS:
sym = _cell_compatible_symmetry(sym)
sym = _relabel_cell_compatible_sg(cfg_crystal, sym,
symprec=symprec)
else:
sym = _ungrouped_symmetry(sym)
explicit_labels = _spin_aware_orbit_labels(
sym, cfg_crystal, group, spin_magnitude
)
syms.append(sym)
# Lattice letter for the struct: when reducing we trust
# spglib's SG-symbol mapping (F, B, CXY, …); when keeping the
# supercell we use the angle-based P/H detector so we never
# claim a centring the supercell does not support as a Bravais
# lattice. Both branches also skip orbit-coset folding via
# the ``lattice`` override (centring factor 1 for P/H).
struct_path = os.path.join(cfg_dir, f"{cfg_name}.struct")
if reduce_cell:
write_struct_file_symmetric(cfg_crystal, sym, struct_path,
rmt=rmt, npt=npt, r0=r0)
else:
write_struct_file_symmetric(
cfg_crystal, sym, struct_path,
rmt=rmt, npt=npt, r0=r0,
lattice=_detect_lattice_letter(cfg_crystal.cell),
labels=explicit_labels,
)
structs.append(struct_path)
udn_orbit = _orbit_udn_letters(cfg_crystal, group, sym, spin_magnitude)
inst_text = _run_instgen_for_orbits(udn_orbit, cfg_crystal, sym, instgen_path)
inst_path = os.path.join(cfg_dir, f"{cfg_name}.inst")
with open(inst_path, "w") as f:
f.write(inst_text)
insts.append(inst_path)
n_atoms = sum(len(g) for g in cfg_crystal.all_species)
reduction_note = (
f" | reduced to {n_atoms} atoms"
if reduce_cell and n_atoms < sum(len(g) for g in crystal.all_species)
else f" | {n_atoms} atoms (no reduction)"
)
# WIEN2k counts only point-group ops (centring translations are
# implicit in the lattice letter), so divide spglib's full count
# by the centring factor of the lattice letter we actually emit.
# When ``--no-reduce`` is in effect we force P/H (the working
# supercell is not the conv cell of the colored SG), so the
# centring factor is 1 and we quote the full op count.
if reduce_cell:
lattice_letter = _wien2k_lattice_letter_from_sg(
sym.spacegroup_symbol, cfg_crystal.cell)
else:
lattice_letter = _detect_lattice_letter(cfg_crystal.cell)
n_ops_wien2k = max(
1, len(sym.rotations) // _wien2k_centering_factor(lattice_letter)
)
print(f" {cfg_name}: SG #{sym.spacegroup_number} "
f"({sym.spacegroup_symbol})"
f" | {sym.n_inequiv} inequiv site(s)"
f" | {n_ops_wien2k} sym op(s)"
f"{reduction_note}")
print(f" {struct_path}")
print(f" {inst_path} [udn={' '.join(udn_orbit)}]")
print()
print(f" {len(structs)} configN.struct + {len(insts)} configN.inst written.")
print()
return {"struct": structs, "inst": insts, "sym": syms, "crystals": crystals}
# ---------------------------------------------------------------------------
# WIEN2k driver scripts (job.init / job.run) at the compound directory
# ---------------------------------------------------------------------------
_JOB_INIT_TEMPLATE = """\
#!/bin/tcsh -f
ls
foreach i (`/bin/ls -d config*`)
cd $i
setrmt # -a Cu:2.1,O:1.7 ....
cp $i.struct_setrmt $i.struct
# NOTE: HDLOs can be kept for plain PBE but not for PBE+U
init:
init_lapw -sp -prec 2n -noreduc -nohdlo -nokshift{init_opts}
# sgroup may suggest an origin shift to have inversion symmetry:
if ($status == 3) then
set shift=`grep SHIFTED $i.outputs`
if ( $#shift != 1) then
echo ' '
echo "init_lapw stopped for $i"
echo "sgroup produced $i.struct_sgroup with origin shift"
echo "compare original struct file with sgroup-struct"
echo 'if lattice and #of atoms do not change, accept the new file'
echo "if lattice changes (eg. from P to CXY), do not accept "
echo ' '
echo " $i.struct:"
head -4 $i.struct
echo " $i.struct_sgroup:"
head -4 $i.struct_sgroup
echo ' '
echo "Accept the new struct file ? (Y/n)"
set yn=($<)
if("$yn" == n ) then
echo Check and correct symmetry in $i using shifting with
echo "x supercell, and atom-labelling until x nn and x symmetry does not complain"
exit 9
endif
endif
cp $i.struct_sgroup $i.struct
goto init
endif
# init_mgga
{orb_block}
cd ..
end
echo "initialization of WIEN2k finished, check for errors on the screen"
echo "create .machines for parallel execution and execute job.run"
"""
_JOB_RUN_TEMPLATE = """\
#!/bin/tcsh -f
ls
foreach i (`/bin/ls -d config*`)
cd $i
cp ../.machines .
runsp -p -ec 0.000001 -cc 0.00001{orb_opt}
save pbe_prec2n # modify name to your case
cd ..
end
echo "scf calculations finished. Now you can analyse your results using mag4 -extract/magnon/angles/moments"
"""
#: Bohr radius in Å — converts a VASP-style --kspacing (Å⁻¹) into the
#: bohr⁻¹ spacing that ``init_lapw -numk -1 <spacing>`` expects.
_BOHR_ANGSTROM = 0.529177
_ORB_BLOCK_COMMENTED = """\
# For PBE+U insert element and U-value and uncomment 4 lines
# (init_orb wants U and J in Ry: U[Ry] = U[eV] / 13.6057):
#init_orb -orb -f <<EOF
#element 2 U(Ry) 0.0
#
#EOF
"""
[docs]
def write_job_scripts(parent_dir: str, *,
functional: Optional[str] = None,
kspacing: Optional[float] = None,
element: Optional[str] = None,
ldaul: Optional[int] = None,
ldauu: Optional[float] = None,
ldauj: Optional[float] = None) -> List[str]:
"""Write ``job.init`` and ``job.run`` into ``parent_dir`` and chmod +x them.
Both scripts iterate over the sibling ``config*/`` directories via
``foreach``; they belong at the compound level, not inside each configN/.
Existing files are overwritten so re-running ``mag4-magnetic`` is safe.
The template follows P. Blaha's recommended workflow: ``setrmt`` (so the
muffin-tin radii mag4 writes into the struct files are placeholders),
``init_lapw -sp -prec 2n -noreduc -nohdlo -nokshift`` (``-nokshift`` is
always present, so the k-mesh is unshifted by default), and an
interactive check when
sgroup proposes an origin-shifted struct (accept only if lattice type and
atom count are unchanged — sgroup may e.g. reduce P to CXY, which would
silently change the cell the report's energy formulas refer to).
``kspacing`` (Å⁻¹, VASP convention; defaults to mag4's 0.3) is converted
to bohr⁻¹ and inserted as ``init_lapw -numk -1 <spacing>`` so the WIEN2k
k-mesh matches the VASP one. When ``functional == "PBE+U"`` and ``ldauu``
is given, the ``init_orb -orb`` heredoc is generated *active* with
``<element> <l> <U> <J>`` from ``element``/``ldaul``/``ldauu``/``ldauj``,
and ``job.run`` calls ``runsp`` with ``-orb``; otherwise the block stays
commented out as a template.
``ldauu``/``ldauj`` are taken in **eV** (VASP convention — the very same
values the INCAR's ``LDAUU``/``LDAUJ`` receive) and divided by
:data:`~mag4.energy.RY_TO_EV` on the way out, because ``init_orb`` reads U
and J in **Rydberg**. ``ldaul`` is a quantum number and is passed through
unchanged. A comment above the heredoc echoes both the eV input and the Ry
output so the conversion is checkable in the generated script.
"""
pbeu = (functional or "").upper() == "PBE+U"
ks = 0.3 if kspacing is None else float(kspacing)
init_opts = (f" -numk -1 {ks * _BOHR_ANGSTROM:.5f}"
f" # = --kspacing {ks:g} A^-1 in bohr^-1; or other options")
if pbeu and ldauu is not None:
# init_orb reads U and J in **Rydberg**, while --ldauu/--ldauj are eV
# (VASP convention, shared with the INCAR), so convert here. l is a
# quantum number and is passed through unchanged.
#
# The eV→Ry echo goes ABOVE ``<<EOF``, never on the data line: every
# line between <<EOF and EOF is piped to init_orb's stdin verbatim, so
# a trailing "# ..." there would be read as data (unlike the -numk
# comment, which sits on a tcsh command line and is stripped by tcsh).
#
# NOTE (unverified with P. Blaha): this targets the magnetic element by
# symbol on a single line. mag4 writes every site as its own
# inequivalent atom (Ni1, Ni2, ... — spin-up and spin-down Ni are
# distinct orbits), so if init_orb expects one entry per inequivalent
# atom rather than an element symbol, this line needs expanding.
u_ev = float(ldauu)
j_ev = 0.0 if ldauj is None else float(ldauj)
orb_block = (
"# PBE+U (element l U J) — init_orb wants U and J in Ry:\n"
f"# --ldauu {u_ev:g} eV = {u_ev / RY_TO_EV:.5f} Ry"
f" / --ldauj {j_ev:g} eV = {j_ev / RY_TO_EV:.5f} Ry\n"
"init_orb -orb -f <<EOF\n"
f"{element or 'element'} {ldaul if ldaul is not None else 2:g}"
f" {u_ev / RY_TO_EV:.5f} {j_ev / RY_TO_EV:.5f}\n"
"\n"
"EOF\n")
else:
orb_block = _ORB_BLOCK_COMMENTED
job_init = _JOB_INIT_TEMPLATE.format(init_opts=init_opts,
orb_block=orb_block)
job_run = _JOB_RUN_TEMPLATE.format(orb_opt=" -orb" if pbeu else " # -orb")
os.makedirs(parent_dir, exist_ok=True)
written: List[str] = []
print()
print("=" * 80)
print(" WIEN2k driver scripts")
print("=" * 80)
for name, body in (("job.init", job_init),
("job.run", job_run)):
dst = os.path.join(parent_dir, name)
with open(dst, "w") as f:
f.write(body)
os.chmod(dst, 0o755)
written.append(dst)
print(f" {dst}")
print()
return written
# ---------------------------------------------------------------------------
# instgen discovery and invocation
# ---------------------------------------------------------------------------
[docs]
def find_instgen() -> Optional[str]:
"""Locate the WIEN2k ``instgen`` script.
Lookup order:
1. ``$WIEN2k_INSTGEN`` environment variable, if set.
2. ``instgen`` on ``$PATH``.
3. ``instgen_lapw`` on ``$PATH`` (the canonical WIEN2k name).
4. ``$WIENROOT/instgen_lapw``.
5. ``~/bin/instgen``.
Returns the absolute path or ``None`` if not found.
"""
env = os.environ.get("WIEN2k_INSTGEN")
if env and os.access(env, os.X_OK):
return env
for name in ("instgen", "instgen_lapw"):
path = shutil.which(name)
if path:
return path
wienroot = os.environ.get("WIENROOT")
if wienroot:
candidate = os.path.join(wienroot, "instgen_lapw")
if os.access(candidate, os.X_OK):
return candidate
fallback = os.path.expanduser("~/bin/instgen")
if os.access(fallback, os.X_OK):
return fallback
return None
[docs]
def run_instgen(udn_letters: Sequence[str],
crystal: CrystalData,
instgen_path: str) -> str:
"""Run ``instgen -ask`` once and return the produced ``case.inst`` text.
Internally uses a temporary directory containing a stub ``case.struct``,
pipes the ``u/d/n`` letters as answers, runs :program:`instgen`, reads
back the resulting ``case.inst`` and cleans up.
Raises
------
RuntimeError
If :program:`instgen` exits non-zero or fails to produce ``case.inst``.
"""
tmp = tempfile.mkdtemp(prefix="mag4_instgen_")
try:
stub_struct = os.path.join(tmp, "stub.struct")
write_stub_struct(crystal, stub_struct)
proc = subprocess.run(
[instgen_path, "-f", "stub", "-ask"],
input="\n".join(udn_letters) + "\n",
capture_output=True,
text=True,
cwd=tmp,
check=False,
)
inst_path = os.path.join(tmp, "stub.inst")
if proc.returncode != 0 or not os.path.isfile(inst_path):
raise RuntimeError(
f"instgen failed (exit {proc.returncode}).\n"
f" command: {instgen_path} -f stub -ask\n"
f" cwd : {tmp}\n"
f" stdout : {proc.stdout}\n"
f" stderr : {proc.stderr}"
)
with open(inst_path) as f:
return f.read()
finally:
shutil.rmtree(tmp, ignore_errors=True)
# ---------------------------------------------------------------------------
# Top-level: write configN.inst into each configN/ directory
# ---------------------------------------------------------------------------
[docs]
def write_case_inst_files(unique_groups: List[List[FourStateConfig]],
crystal: CrystalData,
parent_dir: str,
instgen_path: Optional[str] = None,
spin_magnitude: float = 1.0) -> List[str]:
"""For every unique configuration, write a ``configN.inst`` file in ``configN/``.
The output file is named after the directory (``config1.inst`` inside
``config1/``, etc.) so that WIEN2k's standard ``case == basename(pwd)``
convention works directly on the user's HPC machine.
Parameters
----------
unique_groups, crystal, spin_magnitude
Same conventions as :func:`mag4.vasp.build_magmom_lines`.
parent_dir : str
Compound (or run) directory hosting the ``configN/`` subdirectories.
instgen_path : str, optional
Path to the ``instgen`` script. Auto-detected via
:func:`find_instgen` if omitted.
Returns
-------
list of str
Paths of the written ``configN.inst`` files, in config order.
"""
if instgen_path is None:
instgen_path = find_instgen()
if instgen_path is None:
raise RuntimeError(
"Could not locate 'instgen' on this machine. "
"Pass --instgen PATH or set $WIEN2k_INSTGEN / $WIENROOT."
)
log.info("Using instgen at: %s", instgen_path)
udn_per_cfg = build_udn_lines(unique_groups, crystal, spin_magnitude)
written: List[str] = []
print()
print("=" * 80)
print(f" WIEN2k configN.inst files (instgen = {instgen_path})")
print("=" * 80)
for cfg_idx, (group, udn) in enumerate(zip(unique_groups, udn_per_cfg), 1):
cfg_name = f"config{cfg_idx}"
cfg_dir = os.path.join(parent_dir, cfg_name)
os.makedirs(cfg_dir, exist_ok=True)
inst_text = run_instgen(udn, crystal, instgen_path)
# WIEN2k convention: case basename = directory name
dst = os.path.join(cfg_dir, f"{cfg_name}.inst")
with open(dst, "w") as f:
f.write(inst_text)
written.append(dst)
member_strs = [f"{m.coupling_name}({m.state_label})" for m in group]
udn_str = " ".join(udn)
print(f" {dst} [{', '.join(member_strs)}]")
print(f" u/d/n = {udn_str}")
print()
print(f" {len(written)} configN.inst file(s) written.")
print()
return written