"""VASP input generation: per-config POSCAR / INCAR / KPOINTS for the four-state pipeline.
Layer: **VASP-specific** (one of the two DFT back-ends; see also ``wien2k.py``).
"""
from __future__ import annotations
import logging
import os
from typing import List, Optional, Tuple, Union
import numpy as np
from mag4.constants import CrystalData, FourStateConfig
log = logging.getLogger(__name__)
def _compress_magmom(vals: list) -> str:
"""Compress consecutive equal values: ``[7, 7, 7, -7]`` → ``'3*7 -7'``."""
result = []
i = 0
while i < len(vals):
v = vals[i]
count = 1
while i + count < len(vals) and vals[i + count] == v:
count += 1
result.append(f"{count}*{v:g}" if count > 1 else f"{v:g}")
i += count
return " ".join(result)
[docs]
def build_magmom_lines(unique_groups: List[List[FourStateConfig]],
crystal: CrystalData,
spin_magnitude: float) -> List[str]:
"""Return one MAGMOM string per unique configuration.
Atom order matches the POSCAR (``crystal.all_species`` order). Magnetic
atoms get ``±spin_magnitude``; all others get 0.
"""
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)}
lines: 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)
lines.append(_compress_magmom(vals))
return lines
#: Recognised values for the ``functional`` argument of :func:`incar_content`.
#:
#: - ``PBE`` — plain GGA, no on-site Hubbard term, no exact exchange.
#: - ``PBE+U`` — Dudarev/Liechtenstein DFT+U on the magnetic element only.
#: - ``PBE0`` — unscreened hybrid, ``AEXX=0.25`` by default.
#: - ``R2SCAN`` — convenience alias forcing ``METAGGA=R2SCAN``.
#: - ``METAGGA`` — meta-GGA family; ``--metagga`` selects the flavour
#: (R2SCAN / SCAN / RSCAN / ...).
FUNCTIONAL_CHOICES = ("PBE", "PBE+U", "PBE0", "R2SCAN", "METAGGA")
#: Backward-compatible alias (the historic flag was ``--xc``).
XC_CHOICES = FUNCTIONAL_CHOICES
# ---------------------------------------------------------------------------
# Element-block helpers (used by the meta-GGA LMAXTAU heuristic).
# ---------------------------------------------------------------------------
def _is_f_element(symbol: str) -> bool:
"""Return ``True`` for lanthanoids/actinoids.
Their standard PAW potentials carry f partial waves (lmax=3) — even La
(empty 4f) and Lu (filled 4f, which pymatgen labels block 'd') — so a
meta-GGA run needs ``LMAXTAU = 8``. ``block == 'f'`` alone misses
Lu/Lr, hence the ``is_lanthanoid or is_actinoid`` check.
"""
try:
from pymatgen.core import Element
e = Element(symbol)
return bool(e.is_lanthanoid or e.is_actinoid)
except Exception:
return False
def _ldau_section(
elements: List[str],
magnetic_element: str,
ldaul: int,
ldauu: float,
ldauj: float,
) -> List[str]:
"""Build the LDAU block: the magnetic element gets the user-supplied
``(L, U, J)`` values, every other element gets ``(-1, 0, 0)``.
POSCAR element order is preserved via ``crystal.elements`` so the
VASP ``LDAU*`` arrays line up with the species block.
"""
l_vals = [ldaul if el == magnetic_element else -1 for el in elements]
u_vals = [ldauu if el == magnetic_element else 0.0 for el in elements]
j_vals = [ldauj if el == magnetic_element else 0.0 for el in elements]
elements_str = " ".join(f"{el:>6s}" for el in elements)
lmaxmix = 6 if ldaul >= 3 else (4 if ldaul >= 2 else 2)
return [
"# ── DFT+U (Liechtenstein, LDAUTYPE=2) ──────────────────────",
" LDAU = .TRUE.",
" LDAUTYPE = 2",
f"# {elements_str}",
f" LDAUL = {' '.join(f'{v:>6d}' for v in l_vals)}"
f" # 0=s, 1=p, 2=d, 3=f, -1=no U",
f" LDAUU = {' '.join(f'{v:>6.2f}' for v in u_vals)} # eV",
f" LDAUJ = {' '.join(f'{v:>6.2f}' for v in j_vals)} # eV",
" LDAUPRINT = 2",
" LASPH = .TRUE.",
f" LMAXMIX = {lmaxmix}",
]
def _hybrid_section(aexx: float = 0.25) -> List[str]:
"""Standard PBE0 settings (Adamo / Barone 1999) — ``aexx`` exact
exchange, no range separation. Use ``HSE06`` (HFSCREEN=0.2) if
screening is desired.
"""
return [
f"# ── Hybrid functional (PBE0, {aexx:g} exact exchange) ──────",
" LHFCALC = .TRUE.",
" GGA = PE",
f" AEXX = {aexx:g} # PBE0 fraction of exact exchange",
" HFSCREEN = 0.0 # 0 = unscreened (HSE06 uses 0.2)",
" PRECFOCK = Normal",
" LASPH = .TRUE.",
" TIME = 0.4",
]
def _metagga_section(
metagga: str = "R2SCAN",
*,
elements: Optional[List[str]] = None,
lmaxtau: Union[int, str] = "auto",
) -> List[str]:
"""Build the meta-GGA block.
``metagga`` is the value written for the ``METAGGA`` tag (e.g.
``R2SCAN``, ``SCAN``, ``RSCAN``, ``R2SCANL``). ``LMAXTAU`` is the
maximum *l* used in the PAW one-centre expansion of the
kinetic-energy density. VASP's default is 6 (safe for s/p/d); f
valences (lanthanides / actinides) need 8.
``lmaxtau`` accepts ``"auto"`` (default: 8 if any f-element is
present, otherwise VASP's default 6 is silently used) or an explicit
integer in ``{0, 2, 4, 6, 8}`` (s→2, p→4, d→6, f→8). Force 6 for
f-in-core POTCARs (RE_3, La_s) where the 4f is not in the valence.
"""
out = [
f"# ── Meta-GGA functional ({metagga}) ────────────────────────",
f" METAGGA = {metagga} # requires LASPH=.TRUE. (set below)",
" LASPH = .TRUE.",
"# GGA is intentionally NOT set — METAGGA selects the functional",
]
if str(lmaxtau).lower() != "auto":
lv = int(lmaxtau)
note = ""
if elements:
f_elems = [el for el in elements if _is_f_element(el)]
if f_elems and lv <= 6:
note = (" # f-in-core POTCAR (RE_3/_s): 4f not in valence"
if elements else " # user override")
out.append(f" LMAXTAU = {lv}{note}")
elif elements:
f_elems = [el for el in elements if _is_f_element(el)]
if f_elems:
out.append(
f" LMAXTAU = 8 # f-element(s) present: "
f"{', '.join(f_elems)} (default 6 unsafe for f valence)"
)
else:
out.append(
"# LMAXTAU defaults to 6 (safe for s/p/d; only f-elements need 8)"
)
else:
out.append(
"# LMAXTAU defaults to 6 (safe for s/p/d; only f-elements need 8)"
)
return out
[docs]
def incar_content(
magmom: str,
cfg_label: str,
*,
functional: Optional[str] = None,
xc: Optional[str] = None,
elements: Optional[List[str]] = None,
magnetic_element: Optional[str] = None,
ldaul: Optional[int] = None,
ldauu: Optional[float] = None,
ldauj: float = 0.0,
aexx: float = 0.25,
encut: float = 500.0,
metagga: str = "R2SCAN",
lmaxtau: Union[int, str] = "auto",
# Backward-compatible aliases (hyphen-style kwargs from earlier mag4):
ldau_l: Optional[int] = None,
ldau_u: Optional[float] = None,
ldau_j: Optional[float] = None,
) -> str:
"""Return the INCAR file content as a string for a given config label.
Parameters
----------
magmom
Pre-built MAGMOM string for this configuration.
cfg_label
Human-readable tag printed in ``SYSTEM``.
functional
Exchange-correlation choice. One of :data:`FUNCTIONAL_CHOICES`
(``PBE``, ``PBE+U``, ``PBE0``, ``R2SCAN``, ``METAGGA``).
Case-insensitive.
xc
Backward-compatible alias for ``functional`` (the historic flag
was ``--xc``). ``functional`` takes precedence when both are
supplied.
elements
POSCAR element order — needed when ``functional == "PBE+U"`` so
the LDAU arrays carry one slot per element, and used by the
meta-GGA LMAXTAU heuristic.
magnetic_element
Symbol of the magnetic species (Cu, Ni, Gd, …). The LDAU
values apply only to this element; every other element gets
``-1 / 0.0 / 0.0``.
ldaul, ldauu, ldauj
Hubbard parameters for the magnetic element. ``ldaul`` is
the orbital quantum number (``2`` for d, ``3`` for f);
``ldauu`` and ``ldauj`` are in eV. Required when
``functional == "PBE+U"``.
aexx
Fraction of exact exchange in the PBE0 hybrid (default 0.25).
Ignored for non-hybrid functionals.
encut
Plane-wave cut-off in eV written to ``ENCUT`` (default 500).
metagga
Meta-GGA flavour written to the ``METAGGA`` tag when
``functional == "METAGGA"`` (default R2SCAN, but typical
choices include SCAN, RSCAN, R2SCAN, R2SCANL). When
``functional == "R2SCAN"`` the value is forced to R2SCAN.
lmaxtau
Meta-GGA ``LMAXTAU`` tag (s→2, p→4, d→6, f→8). ``"auto"``
(default): 8 if a lanthanide/actinide is present, else VASP's
default 6 is used (line commented out). Force 6 explicitly
for f-in-core POTCARs (RE_3, La_s) where the 4f is not in the
valence.
ldau_l, ldau_u, ldau_j
Deprecated hyphen-style aliases for ``ldaul``/``ldauu``/``ldauj``
— kept so older callers (and tests) keep working.
"""
# Honor old hyphen-style kwargs as fallback aliases.
if ldaul is None and ldau_l is not None:
ldaul = ldau_l
if ldauu is None and ldau_u is not None:
ldauu = ldau_u
if ldau_j is not None and ldauj == 0.0:
ldauj = ldau_j
# ``functional`` is the canonical kwarg; ``xc`` is the old alias.
func = functional if functional is not None else xc
if func is None:
func = "PBE"
func = func.upper()
if func not in FUNCTIONAL_CHOICES:
raise ValueError(
f"unknown functional={func!r}; choose from {FUNCTIONAL_CHOICES} "
f"(historic alias: ``xc=...``)"
)
# ALGO depends on the functional family.
if func == "PBE0":
algo_line = " ALGO = All # required for hybrids"
elif func in ("R2SCAN", "METAGGA"):
algo_line = " ALGO = All # CG: stable for meta-GGA"
else:
algo_line = " ALGO = Normal"
lines = [
f"# INCAR -- {cfg_label} [{func}]",
f"# Review ALL tags before submitting.",
f"# ── General ────────────────────────────────────────────────",
f" SYSTEM = {cfg_label}",
f" PREC = Accurate",
f" ENCUT = {encut:g} # adjust to your POTCAR",
f" EDIFF = 1E-6",
f" NSW = 0 # single-point (set > 0 for relaxation)",
f"# ── Electronic ─────────────────────────────────────────────",
algo_line,
f" NELM = 100",
f" NELMIN = 6",
f" ISMEAR = 0 # Gaussian smearing ",
f" SIGMA = 0.01",
f"# ── Spin polarisation (no spin-orbit) ───────────────────────",
f" ISPIN = 2",
f" MAGMOM = {magmom}",
f" LORBIT = 11 # projected DOS + magnetic moments",
]
if func == "PBE":
lines.append("# ── Exchange-correlation : PBE (GGA, no on-site term) ─────")
lines.append(" GGA = PE")
elif func == "PBE+U":
if elements is None or magnetic_element is None:
raise ValueError(
"incar_content(functional='PBE+U', ...) requires `elements` "
"and `magnetic_element`"
)
if ldaul is None or ldauu is None:
raise ValueError(
"incar_content(functional='PBE+U', ...) requires `ldaul` "
"and `ldauu` (alias `ldau_l`/`ldau_u` also accepted)"
)
lines.append("# ── Exchange-correlation : PBE+U (GGA + Hubbard U) ─────────")
lines.append(" GGA = PE")
lines.extend(_ldau_section(
elements, magnetic_element, ldaul, ldauu, ldauj,
))
elif func == "PBE0":
lines.extend(_hybrid_section(aexx=aexx))
elif func in ("R2SCAN", "METAGGA"):
mgga = "R2SCAN" if func == "R2SCAN" else metagga
lines.extend(_metagga_section(
mgga, elements=elements, lmaxtau=lmaxtau,
))
lines.extend([
"# ── Output ──────────────────────────────────────────────────",
" LWAVE = .FALSE.",
" LCHARG = .FALSE.",
" LVTOT = .FALSE.",
" NWRITE = 2",
])
return "\n".join(lines) + "\n"
[docs]
def write_incar(magmom: str, name: str, total_atoms: int, **kwargs) -> str:
"""Write an INCAR for Config #1. Returns the file path.
Extra keyword arguments are forwarded to :func:`incar_content`
(``functional`` / ``xc``, ``elements``, ``magnetic_element``,
``ldaul``/``ldauu``/``ldauj``, ``aexx``, ``encut``, ``metagga``,
``lmaxtau``).
"""
incar_path = f"INCAR_{name}"
with open(incar_path, "w") as f:
f.write(incar_content(magmom, f"{name}_config1", **kwargs))
log.info("INCAR written: %s (Config #1 MAGMOM)", incar_path)
return incar_path
#: Default reciprocal-space k-point spacing (Å⁻¹), VASP-style. Denser than
#: VASP's own 0.5 Å⁻¹ default — well converged for the magnetic insulators /
#: semiconductors this is usually applied to. METALS need a finer mesh (smaller
#: KSPACING, ~0.15–0.2 Å⁻¹, with appropriate smearing): pass --kspacing
#: explicitly for those.
DEFAULT_KSPACING = 0.3
def _kmesh_from_kspacing(lat_mat: np.ndarray, kspacing: float) -> Tuple[int, int, int]:
"""Γ-mesh from VASP's KSPACING rule: ``n_i = max(1, ⌈|b_i| / KSPACING⌉)``.
``b_i`` are the reciprocal-lattice vectors **including the 2π factor**, taken
from the full reciprocal matrix (so oblique/hexagonal cells are handled
correctly, unlike a per-axis ``2π/a_i``). See
https://www.vasp.at/wiki/KSPACING.
"""
recip = 2.0 * np.pi * np.linalg.inv(lat_mat).T # columns are b_i
lens = np.linalg.norm(recip, axis=0)
return tuple(max(1, int(np.ceil(L / kspacing))) for L in lens)
def _kmesh_from_kdens(lat_mat: np.ndarray, kdens: float) -> Tuple[int, int, int]:
"""Legacy real-space density rule: ``n_i = max(1, round(kdens / a_i))``."""
a = [float(np.linalg.norm(lat_mat[:, i])) for i in range(3)]
return tuple(max(1, round(kdens / x)) for x in a)
[docs]
def resolve_kmesh(lat_mat: np.ndarray, *,
kspacing: Optional[float] = None,
kdens: Optional[float] = None
) -> Tuple[Tuple[int, int, int], str, float]:
"""Resolve the Γ-mesh and how it was chosen.
``kspacing`` (VASP-style, Å⁻¹) takes precedence; else ``kdens`` (legacy,
Å); else :data:`DEFAULT_KSPACING`. Returns ``(mesh, mode, value)`` with
``mode`` in ``{"kspacing", "kdens"}``.
"""
if kspacing is not None:
return _kmesh_from_kspacing(lat_mat, kspacing), "kspacing", kspacing
if kdens is not None:
return _kmesh_from_kdens(lat_mat, kdens), "kdens", kdens
return (_kmesh_from_kspacing(lat_mat, DEFAULT_KSPACING),
"kspacing", DEFAULT_KSPACING)
[docs]
def write_kpoints(lat_mat: np.ndarray, path: str, *,
kspacing: Optional[float] = None,
kdens: Optional[float] = None) -> str:
"""Write a Γ-centred automatic KPOINTS file.
The mesh follows VASP's KSPACING rule by default (denser than VASP's 0.5
Å⁻¹); pass ``kdens`` to use the legacy real-space density instead.
"""
(na, nb, nc), mode, val = resolve_kmesh(lat_mat, kspacing=kspacing, kdens=kdens)
tag = (f"KSPACING={val:g} A^-1" if mode == "kspacing"
else f"kdens={val:g} A")
with open(path, "w") as f:
f.write(f"Automatic k-mesh ({tag})\n")
f.write("0\n")
f.write("Gamma\n")
f.write(f" {na} {nb} {nc}\n")
f.write(" 0 0 0\n")
return path
[docs]
def create_config_dirs(unique_groups: List[List[FourStateConfig]],
crystal: CrystalData,
magmom_lines: List[str],
name: str,
kdens: Optional[float] = None,
compound: str = "",
write_vasp_extras: bool = True,
*,
kspacing: Optional[float] = None,
functional: Optional[str] = None,
xc: Optional[str] = None,
magnetic_element: Optional[str] = None,
ldaul: Optional[int] = None,
ldauu: Optional[float] = None,
ldauj: float = 0.0,
aexx: float = 0.25,
encut: float = 500.0,
metagga: str = "R2SCAN",
lmaxtau: Union[int, str] = "auto",
# Backward-compat hyphen-style aliases.
ldau_l: Optional[int] = None,
ldau_u: Optional[float] = None,
ldau_j: Optional[float] = None) -> None:
"""Create one directory per unique configuration with POSCAR (+ optional INCAR + KPOINTS).
Layout::
<compound>/config1/POSCAR
<compound>/config1/INCAR (only if write_vasp_extras)
<compound>/config1/KPOINTS (only if write_vasp_extras)
<compound>/config2/...
Parameters
----------
write_vasp_extras : bool, default ``True``
When ``False``, only POSCAR is written (useful for the WIEN2k workflow,
where POSCAR is later converted to ``case.struct`` via WIEN2k's
``xyz2struct`` tool, and the magnetic information lives in ``case.inst``
instead).
"""
# Honor old hyphen-style kwargs as fallback aliases.
if ldaul is None and ldau_l is not None:
ldaul = ldau_l
if ldauu is None and ldau_u is not None:
ldauu = ldau_u
if ldau_j is not None and ldauj == 0.0:
ldauj = ldau_j
if functional is None and xc is not None:
functional = xc
parent = compound if compound else name
os.makedirs(parent, exist_ok=True)
lat = crystal.lattice_matrix
all_atoms = [atom for group in crystal.all_species for atom in group]
def poscar_content(cfg_label: str) -> str:
lines = [cfg_label, "1.0"]
for k in range(3):
v = lat[:, k]
lines.append(f" {v[0]:20.10f} {v[1]:20.10f} {v[2]:20.10f}")
lines.append(" " + " ".join(crystal.elements))
lines.append(" " + " ".join(str(n) for n in crystal.nb_species))
lines.append("Direct")
for atom in all_atoms:
lines.append(f" {atom[1]:16.10f} {atom[2]:16.10f} {atom[3]:16.10f}")
return "\n".join(lines) + "\n"
print()
print("=" * 80)
if write_vasp_extras:
# Identical wording to the original four_state_magnetic.py for non-regression.
print(" VASP input directories")
else:
print(" Geometry directories (POSCAR only — for WIEN2k xyz2struct)")
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]
cfg_label = f"{name}_config{uid}"
dirpath = os.path.join(parent, f"config{uid}")
os.makedirs(dirpath, exist_ok=True)
with open(os.path.join(dirpath, "POSCAR"), "w") as f:
f.write(poscar_content(cfg_label))
if not write_vasp_extras:
print(f" {dirpath}/ [{', '.join(member_strs)}]")
print(f" POSCAR written (on a WIEN2k machine: cp POSCAR xyz2struct.xyz; xyz2struct)")
print()
continue
with open(os.path.join(dirpath, "INCAR"), "w") as f:
f.write(incar_content(
mm, cfg_label,
functional=functional,
elements=list(crystal.elements),
magnetic_element=magnetic_element,
ldaul=ldaul,
ldauu=ldauu,
ldauj=ldauj,
aexx=aexx,
encut=encut,
metagga=metagga,
lmaxtau=lmaxtau,
))
write_kpoints(lat, os.path.join(dirpath, "KPOINTS"),
kspacing=kspacing, kdens=kdens)
mesh, mode, val = resolve_kmesh(lat, kspacing=kspacing, kdens=kdens)
unit = "A^-1" if mode == "kspacing" else "A"
print(f" {dirpath}/ [{', '.join(member_strs)}]")
print(f" MAGMOM = {mm}")
print(f" KPOINTS = {mesh[0]}x{mesh[1]}x{mesh[2]} "
f"({mode}={val:g} {unit})")
print()
print(f" {len(unique_groups)} director(y/ies) created under {parent}/")
print()