#!/usr/bin/env python3
"""Interactive assignment of VASP core-level XAS peaks.

The script combines files from one converged CH_LSPEC calculation:

* OUTCAR   : XAS tensor and core-state eigenenergies
* EIGENVAL : k-point/band energies and occupations
* PROCAR   : local orbital character of the excited (core-hole) atom
* DOSCAR   : total DOS and the excited-atom projected DOS

Only Gamma-point transitions are drawn as vertical assignment guides.  Hovering
over a guide prints and displays its final-state information.  Clicking the XAS
curve prints the strongest candidates from *all* k points in an energy window.
All k-point data, including reciprocal-lattice fractional k coordinates, are
written to an Excel workbook.

By default, the lowest allowed EIGENVAL/PROCAR assignment guide is aligned with
the first XAS peak.  This removes the arbitrary offset between the OUTCAR core
eigenenergy and the CH_LSPEC photon-energy axis while preserving every final-
state spacing.  Use --align none to inspect raw E-Ecore values.

The PROCAR projection is an assignment guide, not the dipole matrix element
used internally by VASP.  In particular, off-diagonal tensor components cannot
be reconstructed from ordinary PROCAR weights because phase information is
missing.

Example
-------
    python vasp_assign_xas.py avg ./xch --emin 283 --emax 315 --sigma 0.3
    python vasp_assign_xas.py xx  ./xch --core-atom 1 --core-label 1s

Required Python packages: numpy, matplotlib, openpyxl
"""

from __future__ import annotations

import argparse
import math
import os
import re
import sys
from fractions import Fraction
from pathlib import Path
from typing import Any, Iterable, Sequence

import numpy as np


TENSOR_COLUMNS = ("xx", "yy", "zz", "xy", "yz", "zx")
SCRIPT_VERSION = "2026.09.01-4"
PREFERRED_ORBITALS = (
    "s",
    "px",
    "py",
    "pz",
    "p_total",
    "dxy",
    "dyz",
    "dz2",
    "dxz",
    "dx2-y2",
    "d_total",
    "f_total",
)


class VaspFormatError(ValueError):
    """Raised when a required VASP output block cannot be parsed safely."""


def vasp_float(text: str) -> float:
    """Convert ordinary, Fortran-D, and VASP exponent-without-E numbers."""
    value = text.strip().replace("D", "E").replace("d", "e")
    if "e" not in value.lower():
        value = re.sub(r"(?<=\d)([+-]\d{2,3})$", r"e\1", value)
    return float(value)


def numeric_fields(line: str) -> list[float] | None:
    fields = line.split()
    if not fields:
        return None
    try:
        return [vasp_float(item) for item in fields]
    except ValueError:
        return None


def read_text(path: Path) -> str:
    try:
        return path.read_text(encoding="utf-8", errors="replace")
    except OSError as exc:
        raise VaspFormatError(f"Cannot read {path}: {exc}") from exc


def read_incar(path: Path) -> dict[str, str]:
    values: dict[str, str] = {}
    if not path.is_file():
        return values
    for raw_line in read_text(path).splitlines():
        line = re.split(r"[!#]", raw_line, maxsplit=1)[0]
        for statement in line.split(";"):
            if "=" not in statement:
                continue
            key, value = statement.split("=", 1)
            key = key.strip().upper()
            if key:
                values[key] = value.strip()
    return values


def incar_int(values: dict[str, str], key: str, default: int | None = None) -> int | None:
    value = values.get(key)
    if value is None:
        return default
    try:
        return int(float(value.split()[0]))
    except (ValueError, IndexError) as exc:
        raise VaspFormatError(f"INCAR: {key}={value!r} is not an integer") from exc


def incar_float(
    values: dict[str, str], key: str, default: float | None = None
) -> float | None:
    value = values.get(key)
    if value is None:
        return default
    try:
        return vasp_float(value.split()[0])
    except (ValueError, IndexError) as exc:
        raise VaspFormatError(f"INCAR: {key}={value!r} is not numeric") from exc


def incar_bool(values: dict[str, str], key: str, default: bool = False) -> bool:
    value = values.get(key)
    if value is None:
        return default
    token = value.split()[0].strip().strip(".").upper()
    if token in ("TRUE", "T", "1"):
        return True
    if token in ("FALSE", "F", "0"):
        return False
    raise VaspFormatError(f"INCAR: {key}={value!r} is not a Boolean value")


def parse_outcar_scalar(text: str, name: str, default: float | int | None = None):
    matches = re.findall(rf"\b{re.escape(name)}\s*=\s*([-+0-9.EeDd]+)", text)
    if not matches:
        return default
    value = vasp_float(matches[-1])
    return int(value) if isinstance(default, int) else value


def parse_outcar_fermi(outcar_text: str) -> float | None:
    matches = re.findall(r"\bE-fermi\s*:\s*([-+0-9.EeDd]+)", outcar_text)
    return vasp_float(matches[-1]) if matches else None


def clean_fractional_coordinate(value: float, tolerance: float = 1.0e-10) -> float:
    """Remove reciprocal-coordinate roundoff while preserving real mesh values."""
    if abs(value) <= tolerance:
        return 0.0
    nearest_integer = round(value)
    if abs(value - nearest_integer) <= tolerance:
        return float(nearest_integer)
    return float(value)


def fractional_component_label(value: float, max_denominator: int = 96) -> str:
    value = clean_fractional_coordinate(value)
    fraction = Fraction(value).limit_denominator(max_denominator)
    if abs(float(fraction) - value) <= 1.0e-8:
        if fraction.denominator == 1:
            return str(fraction.numerator)
        return f"{fraction.numerator}/{fraction.denominator}"
    return f"{value:.8g}"


def fractional_k_label(kpoint: Sequence[float]) -> str:
    return "(" + ", ".join(fractional_component_label(value) for value in kpoint) + ")"


def read_outcar_fractional_kpoints(
    outcar_text: str, nkpts: int
) -> list[dict[str, Any]] | None:
    """Read VASP's explicitly labelled reciprocal-lattice k-point table."""
    lines = outcar_text.splitlines()
    starts = [
        i
        for i, line in enumerate(lines)
        if "k-points in reciprocal lattice and weights" in line.lower()
    ]
    for start in reversed(starts):
        rows: list[dict[str, float]] = []
        for line in lines[start + 1 :]:
            values = numeric_fields(line)
            if values is None or len(values) < 4:
                if rows:
                    break
                continue
            kx, ky, kz, weight = values[:4]
            rows.append(
                {
                    "k_frac_x": clean_fractional_coordinate(kx),
                    "k_frac_y": clean_fractional_coordinate(ky),
                    "k_frac_z": clean_fractional_coordinate(kz),
                    "k_weight_outcar": float(weight),
                }
            )
            if len(rows) == nkpts:
                for row in rows:
                    row["k_frac_label"] = fractional_k_label(
                        (row["k_frac_x"], row["k_frac_y"], row["k_frac_z"])
                    )
                return rows
    return None


def apply_fractional_kpoints(
    levels: list[dict[str, Any]], kpoints: Sequence[dict[str, Any]] | None
) -> str:
    """Attach authoritative fractional k coordinates to every band row."""
    if kpoints is None:
        for row in levels:
            row["k_frac_x"] = clean_fractional_coordinate(row["k_eigen_x"])
            row["k_frac_y"] = clean_fractional_coordinate(row["k_eigen_y"])
            row["k_frac_z"] = clean_fractional_coordinate(row["k_eigen_z"])
            row["k_frac_label"] = fractional_k_label(
                (row["k_frac_x"], row["k_frac_y"], row["k_frac_z"])
            )
        return "EIGENVAL (fallback; interpreted as reciprocal fractional)"
    for row in levels:
        ik = int(row["ik"])
        if not 1 <= ik <= len(kpoints):
            raise VaspFormatError(f"EIGENVAL k-point index {ik} is outside OUTCAR table.")
        row.update(kpoints[ik - 1])
    return "OUTCAR: k-points in reciprocal lattice and weights"


def parse_poscar_elements(path: Path) -> list[str]:
    """Return one element label per atom for an ordinary VASP-5 POSCAR."""
    if not path.is_file():
        return []
    lines = read_text(path).splitlines()
    if len(lines) < 7:
        return []
    symbols = lines[5].split()
    counts_line = lines[6].split()
    try:
        counts = [int(item) for item in counts_line]
    except ValueError:
        # VASP-4 POSCAR: line 6 contains counts and has no species names.
        try:
            counts = [int(item) for item in symbols]
        except ValueError:
            return []
        symbols = [f"X{i + 1}" for i in range(len(counts))]
    if len(symbols) != len(counts):
        return []
    result: list[str] = []
    for symbol, count in zip(symbols, counts):
        result.extend([symbol] * count)
    return result


def read_xas_tensor(outcar_text: str) -> tuple[np.ndarray, dict[str, np.ndarray]]:
    """Read the last density-density imaginary dielectric-function block."""
    lines = outcar_text.splitlines()
    starts = [
        i
        for i, line in enumerate(lines)
        if "IMAGINARY DIELECTRIC FUNCTION" in line.upper()
        and "DENSITY" in line.upper()
    ]
    if not starts:
        raise VaspFormatError(
            "OUTCAR has no density-density IMAGINARY DIELECTRIC FUNCTION block. "
            "Check CH_LSPEC=.TRUE. and calculation completion."
        )

    # The last complete block corresponds to the final result if VASP printed
    # more than one dielectric block.
    for start in reversed(starts):
        rows: list[list[float]] = []
        data_started = False
        for line in lines[start + 1 :]:
            numbers = numeric_fields(line)
            if numbers is not None and len(numbers) >= 7:
                rows.append(numbers[:7])
                data_started = True
            elif data_started:
                break
        if rows:
            array = np.asarray(rows, dtype=float)
            tensor = {
                name: array[:, i + 1].copy() for i, name in enumerate(TENSOR_COLUMNS)
            }
            tensor["avg"] = (tensor["xx"] + tensor["yy"] + tensor["zz"]) / 3.0
            return array[:, 0].copy(), tensor
    raise VaspFormatError("The XAS header was found, but its numeric table is empty.")


def gaussian_broaden(
    energy: np.ndarray, values: np.ndarray, sigma: float
) -> np.ndarray:
    """Apply Gaussian convolution with sigma in the units of *energy*.

    VASP XAS and DOS tables use uniform meshes.  Explicit zero padding keeps
    the output length fixed even when sigma is large relative to the table.
    """
    energy = np.asarray(energy, dtype=float)
    values = np.asarray(values, dtype=float)
    if energy.ndim != 1 or values.ndim != 1 or len(energy) != len(values):
        raise VaspFormatError("Gaussian broadening requires equal-length 1-D arrays.")
    if sigma <= 0.0 or len(energy) < 2:
        return values.copy()
    spacing = np.diff(energy)
    if np.any(spacing <= 0.0):
        raise VaspFormatError("Energy mesh must be strictly increasing for broadening.")
    step = float(np.median(spacing))
    if not np.allclose(spacing, step, rtol=2.0e-4, atol=max(1.0e-10, step * 1.0e-6)):
        raise VaspFormatError("Energy mesh is not uniform enough for Gaussian broadening.")
    radius = max(1, int(math.ceil(6.0 * sigma / step)))
    offsets = np.arange(-radius, radius + 1, dtype=float) * step
    kernel = np.exp(-0.5 * (offsets / sigma) ** 2)
    kernel /= kernel.sum()
    padded = np.pad(values, (radius, radius), mode="constant")
    return np.convolve(padded, kernel, mode="valid")


def broaden_tensor(
    energy: np.ndarray, tensor: dict[str, np.ndarray], sigma: float
) -> dict[str, np.ndarray]:
    broadened = {
        name: gaussian_broaden(energy, values, sigma)
        for name, values in tensor.items()
        if name != "avg"
    }
    broadened["avg"] = (
        broadened["xx"] + broadened["yy"] + broadened["zz"]
    ) / 3.0
    return broadened


def broaden_dos(dos: dict[str, Any], sigma: float) -> dict[str, Any]:
    result = dict(dos)
    energy = np.asarray(dos["E"], dtype=float)
    total_energy = np.asarray(dos["E_total"], dtype=float)
    result["tDOS"] = gaussian_broaden(total_energy, np.asarray(dos["tDOS"]), sigma)
    for label in ("s", "p", "d", "f"):
        result[label] = gaussian_broaden(energy, np.asarray(dos[label]), sigma)
    return result


def find_xas_anchor(
    energy: np.ndarray, values: np.ndarray, method: str, threshold: float
) -> float:
    """Return the onset or first local peak of the selected XAS component."""
    energy = np.asarray(energy, dtype=float)
    magnitude = np.abs(np.asarray(values, dtype=float))
    if len(energy) != len(magnitude) or not len(energy):
        raise VaspFormatError("Cannot align an empty XAS spectrum.")
    maximum = float(np.nanmax(magnitude))
    if not np.isfinite(maximum) or maximum <= 0.0:
        raise VaspFormatError("Cannot align an XAS component with zero intensity.")
    above = np.flatnonzero(magnitude >= threshold * maximum)
    if not len(above):
        raise VaspFormatError("No XAS point passes the alignment threshold.")
    onset = int(above[0])
    if method == "onset":
        return float(energy[onset])
    for i in range(max(onset, 1), len(magnitude) - 1):
        if (
            magnitude[i] >= threshold * maximum
            and magnitude[i] >= magnitude[i - 1]
            and magnitude[i] > magnitude[i + 1]
        ):
            return float(energy[i])
    return float(energy[int(np.nanargmax(magnitude))])


def read_core_states(outcar_text: str, elements: Sequence[str]) -> list[dict[str, Any]]:
    """Read the final 'core state eigenenergies' table from OUTCAR."""
    lines = outcar_text.splitlines()
    markers = [i for i, line in enumerate(lines) if "core state eigenenergies" in line.lower()]
    if not markers:
        raise VaspFormatError("OUTCAR has no 'the core state eigenenergies are' block.")

    result: list[dict[str, Any]] = []
    atom_line = re.compile(r"^\s*(\d+)\s*-\s*(.*)$")
    pair = re.compile(r"(\d+[spdfgh])\s+([-+0-9.EeDd]+)", re.IGNORECASE)
    for line in lines[markers[-1] + 1 :]:
        match = atom_line.match(line)
        if not match:
            if result and not line.strip():
                break
            continue
        atom = int(match.group(1))
        for label, energy in pair.findall(match.group(2)):
            result.append(
                {
                    "atom": atom,
                    "element": elements[atom - 1] if 0 < atom <= len(elements) else "",
                    "core_state": label.lower(),
                    "E_core": vasp_float(energy),
                }
            )
    if not result:
        raise VaspFormatError("The core-state header was found, but no energies were parsed.")
    return result


def choose_core_state(
    core_states: Sequence[dict[str, Any]],
    incar: dict[str, str],
    atom_override: int | None,
    label_override: str | None,
    energy_override: float | None,
) -> dict[str, Any]:
    atom = atom_override if atom_override is not None else incar_int(incar, "CLNT", 1)
    if label_override:
        label = label_override.lower()
    else:
        n = incar_int(incar, "CLN", None)
        ell = incar_int(incar, "CLL", None)
        letters = "spdfgh"
        label = f"{n}{letters[ell]}" if n is not None and ell is not None and 0 <= ell < len(letters) else "1s"

    candidates = [row for row in core_states if row["atom"] == atom and row["core_state"] == label]
    if energy_override is not None:
        element = candidates[0]["element"] if candidates else ""
        return {
            "atom": atom,
            "element": element,
            "core_state": label,
            "E_core": float(energy_override),
            "source": "--core-energy",
        }
    if len(candidates) != 1:
        available = ", ".join(
            f"atom {row['atom']} {row['core_state']}={row['E_core']:.6f} eV"
            for row in core_states
        )
        raise VaspFormatError(
            f"Cannot uniquely select atom {atom}, core state {label}. Available: {available}. "
            "Use --core-atom/--core-label/--core-energy if needed."
        )
    selected = dict(candidates[0])
    selected["source"] = "OUTCAR"
    return selected


def _next_nonempty(fp, description: str) -> str:
    for line in fp:
        if line.strip():
            return line
    raise VaspFormatError(f"EIGENVAL ended while reading {description}.")


def read_eigenval(path: Path, ispin: int) -> tuple[list[dict[str, Any]], dict[str, int]]:
    if not path.is_file():
        raise VaspFormatError(f"EIGENVAL not found: {path}")
    levels: list[dict[str, Any]] = []
    with path.open("r", encoding="utf-8", errors="replace") as fp:
        for _ in range(5):
            if fp.readline() == "":
                raise VaspFormatError("EIGENVAL is shorter than its six-line header.")
        header = _next_nonempty(fp, "NELECT/NKPTS/NBANDS header").split()
        if len(header) < 3:
            raise VaspFormatError("EIGENVAL NELECT/NKPTS/NBANDS header is invalid.")
        nelect, nkpts, nbands = (int(float(header[i])) for i in range(3))

        for ik in range(1, nkpts + 1):
            k_values = numeric_fields(_next_nonempty(fp, f"k point {ik}"))
            if k_values is None or len(k_values) < 4:
                raise VaspFormatError(f"EIGENVAL k-point header {ik} is invalid.")
            k_eigen_x, k_eigen_y, k_eigen_z, weight = k_values[:4]
            for expected_band in range(1, nbands + 1):
                values = numeric_fields(_next_nonempty(fp, f"k point {ik}, band {expected_band}"))
                if values is None:
                    raise VaspFormatError(f"EIGENVAL band row {expected_band} is invalid.")
                band = int(values[0])
                if ispin == 2:
                    if len(values) < 5:
                        raise VaspFormatError(
                            "ISPIN=2 but EIGENVAL band rows have fewer than five columns."
                        )
                    energies = values[1:3]
                    occupations = values[3:5]
                else:
                    if len(values) < 3:
                        raise VaspFormatError("EIGENVAL band rows have fewer than three columns.")
                    energies = [values[1]]
                    occupations = [values[2]]
                for spin, (energy, occupation) in enumerate(zip(energies, occupations), 1):
                    levels.append(
                        {
                            "ik": ik,
                            "k_eigen_x": k_eigen_x,
                            "k_eigen_y": k_eigen_y,
                            "k_eigen_z": k_eigen_z,
                            "k_weight": weight,
                            "spin": spin,
                            "ilevel": band,
                            "E": energy,
                            "occupation": occupation,
                        }
                    )
    return levels, {"nelect": nelect, "nkpts": nkpts, "nbands": nbands}


def normalize_orbital(name: str) -> str:
    value = name.strip().lower().replace("_", "")
    aliases = {
        "dx2": "dx2-y2",
        "x2-y2": "dx2-y2",
        "dx2-y2": "dx2-y2",
        "tot": "tot",
    }
    return aliases.get(value, value)


def read_procar(path: Path, core_atom: int) -> tuple[dict[tuple[int, int, int], dict[str, float]], set[str]]:
    """Read only the selected atom's first projection table per band."""
    if not path.is_file():
        raise VaspFormatError(f"PROCAR not found: {path}. Run VASP with LORBIT=11.")

    projections: dict[tuple[int, int, int], dict[str, float]] = {}
    orbitals_seen: set[str] = set()
    spin = 1
    ik: int | None = None
    band: int | None = None
    headers: list[str] | None = None
    spin_re = re.compile(r"spin\s+component\s+(\d+)", re.IGNORECASE)
    k_re = re.compile(r"k-point\s+(\d+)\s*:", re.IGNORECASE)
    band_re = re.compile(r"^\s*band\s+(\d+)\b", re.IGNORECASE)

    with path.open("r", encoding="utf-8", errors="replace") as fp:
        for line in fp:
            match = spin_re.search(line)
            if match:
                spin = int(match.group(1))
                continue
            match = k_re.search(line)
            if match:
                ik = int(match.group(1))
                continue
            match = band_re.match(line)
            if match:
                band = int(match.group(1))
                headers = None
                continue

            fields = line.split()
            if fields and fields[0].lower() == "ion":
                headers = [normalize_orbital(item) for item in fields[1:] if item.lower() != "tot"]
                orbitals_seen.update(headers)
                continue
            if headers is None or ik is None or band is None or not fields:
                continue
            if not fields[0].isdigit() or int(fields[0]) != core_atom:
                continue
            key = (spin, ik, band)
            if key in projections:
                # With phase output/non-collinear output, later tables can repeat
                # an ion header.  The first table is the scalar orbital weight.
                continue
            try:
                numbers = [vasp_float(item) for item in fields[1 : 1 + len(headers)]]
            except ValueError:
                continue
            if len(numbers) != len(headers):
                continue
            projections[key] = dict(zip(headers, numbers))

    if not projections:
        raise VaspFormatError(
            f"No projection for atom {core_atom} was read from PROCAR. "
            "Check atom numbering and use LORBIT=11."
        )
    return projections, orbitals_seen


PDOS_LABELS = {
    1: ["s"],
    3: ["s", "p", "d"],
    4: ["s", "p", "d", "f"],
    9: ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2"],
    16: [
        "s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2",
        "fy(3x2-y2)", "fxyz", "fyz2", "fz3", "fxz2", "fz(x2-y2)", "fx(x2-3y2)",
    ],
}


def l_sum(values: dict[str, np.ndarray], letter: str, length: int) -> np.ndarray:
    if letter in values:
        return values[letter].copy()
    matches = [array for name, array in values.items() if name.startswith(letter)]
    return np.sum(matches, axis=0) if matches else np.zeros(length, dtype=float)


def read_doscar(path: Path, core_atom: int, ispin: int) -> dict[str, np.ndarray | float | int]:
    if not path.is_file():
        raise VaspFormatError(f"DOSCAR not found: {path}")
    with path.open("r", encoding="utf-8", errors="replace") as fp:
        first = fp.readline().split()
        if not first:
            raise VaspFormatError("DOSCAR first line is empty.")
        nions = int(float(first[0]))
        for _ in range(4):
            if fp.readline() == "":
                raise VaspFormatError("DOSCAR ended inside its header.")
        total_header_line = fp.readline()
        total_header = total_header_line.split()
        try:
            if len(total_header) < 4:
                raise ValueError
            nedos_value = vasp_float(total_header[2])
            nedos = int(nedos_value)
            if nedos_value != nedos:
                raise ValueError
            efermi = vasp_float(total_header[3])
        except (ValueError, IndexError):
            # Older/fixed-width VASP output can concatenate Emin and NEDOS.
            # Match VASP's 16+16+5+16 field layout in that case.
            if len(total_header_line) < 53:
                raise VaspFormatError("DOSCAR total-DOS header is invalid.")
            try:
                nedos_value = vasp_float(total_header_line[32:37])
                nedos = int(nedos_value)
                if nedos_value != nedos:
                    raise ValueError
                efermi = vasp_float(total_header_line[37:53])
            except ValueError as exc:
                raise VaspFormatError("DOSCAR fixed-width total-DOS header is invalid.") from exc

        total_rows: list[list[float]] = []
        for _ in range(nedos):
            values = numeric_fields(fp.readline())
            if values is None or len(values) < 3:
                raise VaspFormatError("DOSCAR total-DOS block is incomplete.")
            total_rows.append(values)
        total_array = np.asarray(total_rows, dtype=float)
        if ispin == 2 and total_array.shape[1] >= 5:
            total_dos = total_array[:, 1] + total_array[:, 2]
        else:
            total_dos = total_array[:, 1]

        core_rows: list[list[float]] | None = None
        for ion in range(1, nions + 1):
            block_header = fp.readline()
            if not block_header:
                break
            rows: list[list[float]] = []
            for _ in range(nedos):
                values = numeric_fields(fp.readline())
                if values is None or len(values) < 2:
                    raise VaspFormatError(f"DOSCAR projected-DOS block for ion {ion} is incomplete.")
                rows.append(values)
            if ion == core_atom:
                core_rows = rows

    if core_rows is None:
        raise VaspFormatError(
            f"DOSCAR has no projected-DOS block for atom {core_atom}. "
            "Check LORBIT and DOSCAR completeness."
        )
    core = np.asarray(core_rows, dtype=float)
    nvalues = core.shape[1] - 1
    if ispin == 2:
        if nvalues % 2:
            raise VaspFormatError(
                f"Spin-polarized PDOS has an odd number ({nvalues}) of value columns."
            )
        norbitals = nvalues // 2
    else:
        norbitals = nvalues
    labels = PDOS_LABELS.get(norbitals)
    if labels is None:
        raise VaspFormatError(
            f"Unsupported DOSCAR projection layout: {norbitals} orbitals "
            f"({nvalues} value columns, ISPIN={ispin})."
        )

    orbital_arrays: dict[str, np.ndarray] = {}
    for i, label in enumerate(labels):
        if ispin == 2:
            orbital_arrays[label] = core[:, 1 + 2 * i] + core[:, 2 + 2 * i]
        else:
            orbital_arrays[label] = core[:, 1 + i]
    length = len(core)
    return {
        "E": core[:, 0].copy(),
        "E_total": total_array[:, 0].copy(),
        "tDOS": total_dos.copy(),
        "s": l_sum(orbital_arrays, "s", length),
        "p": l_sum(orbital_arrays, "p", length),
        "d": l_sum(orbital_arrays, "d", length),
        "f": l_sum(orbital_arrays, "f", length),
        "Efermi": efermi,
        "NEDOS": nedos,
    }


def finite_sum(values: Iterable[float]) -> float:
    finite = [value for value in values if np.isfinite(value)]
    return float(sum(finite)) if finite else math.nan


def projection_columns(projection: dict[str, float]) -> dict[str, float]:
    def get(name: str) -> float:
        value = projection.get(name, math.nan)
        return float(value)

    px, py, pz = get("px"), get("py"), get("pz")
    p_total = get("p") if np.isfinite(get("p")) else finite_sum((px, py, pz))
    d_names = ("dxy", "dyz", "dz2", "dxz", "dx2-y2")
    d_values = [get(name) for name in d_names]
    d_total = get("d") if np.isfinite(get("d")) else finite_sum(d_values)
    f_components = [value for name, value in projection.items() if name.startswith("f")]
    f_total = get("f") if np.isfinite(get("f")) else finite_sum(f_components)
    result = {
        "s": get("s"),
        "px": px,
        "py": py,
        "pz": pz,
        "p_total": p_total,
        "dxy": get("dxy"),
        "dyz": get("dyz"),
        "dz2": get("dz2"),
        "dxz": get("dxz"),
        "dx2-y2": get("dx2-y2"),
        "d_total": d_total,
        "f_total": f_total,
    }
    return result


def directional_projection(row: dict[str, Any], mode: str) -> float:
    px, py, pz = (row.get(name, math.nan) for name in ("px", "py", "pz"))
    p_total = row.get("p_total", math.nan)

    # LORBIT=10 only provides p_total.  An equal one-third partition is a
    # neutral visualization fallback, not a polarization assignment.
    if not all(np.isfinite(value) for value in (px, py, pz)):
        return max(float(p_total) / 3.0, 0.0) if np.isfinite(p_total) else 0.0
    if mode == "xx":
        return max(float(px), 0.0)
    if mode == "yy":
        return max(float(py), 0.0)
    if mode == "zz":
        return max(float(pz), 0.0)
    if mode == "xy":
        return math.sqrt(max(float(px), 0.0) * max(float(py), 0.0))
    if mode == "yz":
        return math.sqrt(max(float(py), 0.0) * max(float(pz), 0.0))
    if mode == "zx":
        return math.sqrt(max(float(pz), 0.0) * max(float(px), 0.0))
    return max(float(p_total), 0.0) / 3.0 if np.isfinite(p_total) else 0.0


def infer_full_occupation(eigen_levels: Sequence[dict[str, Any]]) -> float:
    """Infer whether a fully occupied EIGENVAL state carries 1 or 2 electrons."""
    occupations = [
        float(row["occupation"])
        for row in eigen_levels
        if np.isfinite(float(row["occupation"]))
    ]
    if not occupations:
        raise VaspFormatError("EIGENVAL contains no finite occupations.")
    maximum = max(occupations)
    if maximum <= 1.05:
        return 1.0
    if maximum <= 2.05:
        return 2.0
    raise VaspFormatError(
        f"Cannot infer a full EIGENVAL occupation from maximum occupation {maximum:.6g}."
    )


def is_gamma(row: dict[str, Any], tolerance: float) -> bool:
    return all(
        abs(float(row[name]) - round(float(row[name]))) <= tolerance
        for name in ("k_frac_x", "k_frac_y", "k_frac_z")
    )


def build_level_rows(
    eigen_levels: list[dict[str, Any]],
    projections: dict[tuple[int, int, int], dict[str, float]],
    core_energy: float,
    mode: str,
    occ_full: float,
    gamma_tolerance: float,
) -> tuple[list[dict[str, Any]], int]:
    rows: list[dict[str, Any]] = []
    missing = 0
    for level in eigen_levels:
        row = dict(level)
        projection = projections.get((row["spin"], row["ik"], row["ilevel"]))
        if projection is None:
            missing += 1
            projection = {}
        row.update(projection_columns(projection))
        row["E_core"] = core_energy
        row["E-Ecore"] = row["E"] - core_energy
        row["E_transition_plot"] = row["E-Ecore"]
        row["unoccupied_fraction"] = float(np.clip(1.0 - row["occupation"] / occ_full, 0.0, 1.0))
        row["projection_guide"] = directional_projection(row, mode)
        row["weighted_guide"] = (
            row["k_weight"] * row["unoccupied_fraction"] * row["projection_guide"]
        )
        row["is_gamma"] = is_gamma(row, gamma_tolerance)
        rows.append(row)
    return rows, missing


def align_transition_energies(
    transitions: Sequence[dict[str, Any]],
    xas_energy: np.ndarray,
    xas_values: np.ndarray,
    method: str,
    threshold: float,
    min_local_p: float,
    manual_shift: float,
) -> dict[str, float | str | None]:
    """Align the lowest allowed guide to an XAS onset/first peak."""
    allowed = [
        row
        for row in transitions
        if row["E-Ecore"] > 0.0
        and np.isfinite(row["p_total"])
        and row["p_total"] >= min_local_p
        and row["weighted_guide"] > 0.0
    ]
    if not allowed:
        raise VaspFormatError(
            "No allowed EIGENVAL/PROCAR transition is available. "
            "Lower --min-local-p."
        )
    transition_anchor = min(float(row["E-Ecore"]) for row in allowed)
    if method == "none":
        xas_anchor = None
        auto_shift = 0.0
    else:
        xas_anchor = find_xas_anchor(xas_energy, xas_values, method, threshold)
        auto_shift = xas_anchor - transition_anchor
    total_shift = auto_shift + manual_shift
    for row in transitions:
        row["alignment_shift_eV"] = total_shift
        row["E_XAS_axis"] = row["E-Ecore"] + total_shift
        row["E_transition_plot"] = row["E-Ecore"] + total_shift
    return {
        "method": method,
        "xas_anchor_eV": xas_anchor,
        "transition_anchor_raw_eV": transition_anchor,
        "auto_shift_eV": auto_shift,
        "manual_shift_eV": manual_shift,
        "total_shift_eV": total_shift,
    }


LEVEL_COLUMNS = [
    "ik", "k_frac_label", "k_frac_x", "k_frac_y", "k_frac_z",
    "k_eigen_x", "k_eigen_y", "k_eigen_z", "k_weight", "k_weight_outcar",
    "spin", "ilevel", "occupation",
    "unoccupied_fraction", "E", "E_core", "E-Ecore", "alignment_shift_eV",
    "E_XAS_axis", "E_transition_plot",
    *PREFERRED_ORBITALS, "projection_guide", "weighted_guide", "is_gamma",
]


def choose_gamma_lines(
    transitions: list[dict[str, Any]],
    emin: float,
    emax: float,
    min_local_p: float,
    max_lines: int,
) -> tuple[list[dict[str, Any]], int]:
    candidates = [
        row
        for row in transitions
        if row["is_gamma"]
        and emin <= row["E_transition_plot"] <= emax
        and np.isfinite(row["p_total"])
        and row["p_total"] >= min_local_p
    ]
    omitted = max(0, len(candidates) - max_lines)
    if omitted:
        candidates = sorted(candidates, key=lambda row: row["weighted_guide"], reverse=True)[:max_lines]
    return sorted(candidates, key=lambda row: row["E_transition_plot"]), omitted


def excel_value(value: Any) -> Any:
    if isinstance(value, np.generic):
        value = value.item()
    if isinstance(value, float) and not math.isfinite(value):
        return None
    if isinstance(value, Path):
        return str(value)
    return value


def write_dict_sheets(
    workbook,
    base_name: str,
    rows: Sequence[dict[str, Any]],
    columns: Sequence[str],
    max_data_rows: int = 1_000_000,
) -> None:
    if not rows:
        sheet = workbook.create_sheet(base_name[:31])
        sheet.append(list(columns))
        return
    for part, start in enumerate(range(0, len(rows), max_data_rows), 1):
        name = base_name if part == 1 else f"{base_name}_{part}"
        sheet = workbook.create_sheet(name[:31])
        sheet.freeze_panes = "A2"
        sheet.append(list(columns))
        for row in rows[start : start + max_data_rows]:
            sheet.append([excel_value(row.get(column)) for column in columns])


def write_excel(
    path: Path,
    config_rows: list[dict[str, Any]],
    core_states: list[dict[str, Any]],
    transitions: list[dict[str, Any]],
    gamma_rows: list[dict[str, Any]],
    xas_energy: np.ndarray,
    tensor_raw: dict[str, np.ndarray],
    tensor_display: dict[str, np.ndarray],
    selected_mode: str,
    dos_raw: dict[str, Any],
    dos_display: dict[str, Any],
    core_energy: float,
    total_shift: float,
    dos_reference_shift: float,
    dos_cutoff_plot: float,
) -> None:
    try:
        from openpyxl import Workbook
    except ImportError as exc:
        raise RuntimeError("Excel output requires openpyxl: python -m pip install openpyxl") from exc

    workbook = Workbook(write_only=True)
    write_dict_sheets(workbook, "config", config_rows, ("parameter", "value"))
    write_dict_sheets(
        workbook,
        "core_states",
        core_states,
        ("atom", "element", "core_state", "E_core", "selected"),
    )
    write_dict_sheets(workbook, "transitions_all", transitions, LEVEL_COLUMNS)
    write_dict_sheets(workbook, "transitions_selected", gamma_rows, LEVEL_COLUMNS)

    spectrum_rows = []
    for i, energy in enumerate(xas_energy):
        row = {"E_XAS": energy}
        for component in (*TENSOR_COLUMNS, "avg"):
            row[f"{component}_raw"] = tensor_raw[component][i]
            row[component] = tensor_display[component][i]
        row["selected_raw"] = tensor_raw[selected_mode][i]
        row["selected"] = tensor_display[selected_mode][i]
        spectrum_rows.append(row)
    spectrum_columns = ["E_XAS"]
    for component in (*TENSOR_COLUMNS, "avg"):
        spectrum_columns.extend((f"{component}_raw", component))
    spectrum_columns.extend(("selected_raw", "selected"))
    write_dict_sheets(
        workbook,
        "spectrum",
        spectrum_rows,
        spectrum_columns,
    )

    dos_rows = []
    energy = np.asarray(dos_raw["E"])
    total_energy = np.asarray(dos_raw["E_total"])
    # Standard DOSCAR uses the same energy mesh for total and projected DOS.
    if len(total_energy) != len(energy) or not np.allclose(total_energy, energy, atol=1e-7):
        raise VaspFormatError("DOSCAR total and projected DOS energy meshes differ.")
    for i, raw_energy in enumerate(energy):
        eigen_reference_energy = raw_energy + dos_reference_shift
        plot_energy = eigen_reference_energy - core_energy + total_shift
        dos_rows.append(
            {
                "E_DOSCAR": raw_energy,
                "E_eigen_reference": eigen_reference_energy,
                "E-Ecore": eigen_reference_energy - core_energy,
                "DOS_reference_shift_eV": dos_reference_shift,
                "alignment_shift_eV": total_shift,
                "E_XAS_axis": plot_energy,
                "E_plot": plot_energy,
                "is_unoccupied_region": plot_energy >= dos_cutoff_plot - 1.0e-9,
                "tDOS_raw": dos_raw["tDOS"][i],
                "tDOS": dos_display["tDOS"][i],
                "s_raw": dos_raw["s"][i],
                "s": dos_display["s"][i],
                "p_raw": dos_raw["p"][i],
                "p": dos_display["p"][i],
                "d_raw": dos_raw["d"][i],
                "d": dos_display["d"][i],
                "f_raw": dos_raw["f"][i],
                "f": dos_display["f"][i],
            }
        )
    write_dict_sheets(
        workbook,
        "dos",
        dos_rows,
        (
            "E_DOSCAR", "DOS_reference_shift_eV", "E_eigen_reference",
            "E-Ecore", "alignment_shift_eV", "E_XAS_axis", "E_plot",
            "is_unoccupied_region", "tDOS_raw", "tDOS",
            "s_raw", "s", "p_raw", "p", "d_raw", "d", "f_raw", "f",
        ),
    )
    path.parent.mkdir(parents=True, exist_ok=True)
    workbook.save(path)


def short_transition(row: dict[str, Any]) -> str:
    orbital_text = (
        f"s={row['s']:.4g}, px={row['px']:.4g}, py={row['py']:.4g}, "
        f"pz={row['pz']:.4g}, d={row['d_total']:.4g}"
    )
    return (
        f"k#{row['ik']} fractional={row['k_frac_label']}, "
        f"band={row['ilevel']}, spin={row['spin']}\n"
        f"E_band={row['E']:.6f} eV, raw E-Ecore={row['E-Ecore']:.6f} eV\n"
        f"alignment shift={row['alignment_shift_eV']:+.6f} eV, "
        f"XAS-axis E={row['E_XAS_axis']:.6f} eV\n"
        f"occ={row['occupation']:.5g}, k-weight={row['k_weight']:.6g}, "
        f"guide={row['weighted_guide']:.5g}\n{orbital_text}"
    )


def print_transition(row: dict[str, Any], prefix: str = "") -> None:
    text = short_transition(row).replace("\n", " | ")
    print(f"{prefix}{text}")


def plot_results(
    xas_energy: np.ndarray,
    xas_y: np.ndarray,
    mode: str,
    dos: dict[str, Any],
    core_energy: float,
    transition_shift: float,
    gamma_rows: list[dict[str, Any]],
    all_transitions: list[dict[str, Any]],
    emin: float,
    emax: float,
    click_window: float,
    top_n: int,
    title: str,
    sigma: float,
    dos_reference_shift: float,
    dos_cutoff_plot: float,
    show_occupied_dos: bool,
    figure_path: Path | None,
    show: bool,
) -> None:
    import matplotlib.pyplot as plt

    fig, (ax_xas, ax_dos) = plt.subplots(
        2,
        1,
        figsize=(11.5, 8.5),
        sharex=True,
        gridspec_kw={"height_ratios": (2.1, 1.35)},
    )
    sigma_text = f", sigma={sigma:g} eV" if sigma > 0.0 else ""
    ax_xas.plot(
        xas_energy, xas_y, color="black", linewidth=1.35,
        label=f"XAS {mode}{sigma_text}",
    )
    ax_xas.set_ylabel("XAS intensity (arb. units)")
    ax_xas.set_title(title)
    ax_xas.grid(alpha=0.18)

    line_artists = []
    max_guide = max((row["weighted_guide"] for row in gamma_rows), default=0.0)
    for row in gamma_rows:
        relative = row["weighted_guide"] / max_guide if max_guide > 0 else 0.0
        height = 0.18 + 0.70 * math.sqrt(max(relative, 0.0))
        alpha = 0.28 + 0.62 * math.sqrt(max(relative, 0.0))
        artist = ax_xas.axvline(
            row["E_transition_plot"], ymin=0.0, ymax=height,
            color="tab:purple", linewidth=0.85, alpha=alpha,
        )
        line_artists.append((artist, row))
    if gamma_rows:
        ax_xas.plot([], [], color="tab:purple", linewidth=1.0, label="Gamma assignment guide")
    ax_xas.legend(loc="best")
    ax_xas.set_xlim(emin, emax)

    dos_energy = (
        np.asarray(dos["E"]) + dos_reference_shift - core_energy + transition_shift
    )
    total_dos = np.asarray(dos["tDOS"])
    dos_mask = np.ones(len(dos_energy), dtype=bool)
    if not show_occupied_dos:
        dos_mask = dos_energy >= dos_cutoff_plot - 1.0e-9
    total_dos_plot = np.where(dos_mask, total_dos, np.nan)
    dos_scope = "all states" if show_occupied_dos else "unoccupied"
    ax_dos.plot(
        dos_energy, total_dos_plot, color="0.35", linewidth=1.0,
        label=f"tDOS ({dos_scope}, cell){sigma_text}",
    )
    ax_dos.fill_between(
        dos_energy, 0.0, np.nan_to_num(total_dos_plot), where=dos_mask,
        color="0.80", alpha=0.45,
    )
    ax_dos.set_ylabel("tDOS (states/eV)", color="0.30")
    ax_dos.tick_params(axis="y", labelcolor="0.30")
    ax_dos.grid(alpha=0.18)

    ax_pdos = ax_dos.twinx()
    pdos_colors = {"s": "tab:blue", "p": "tab:red", "d": "tab:green", "f": "tab:orange"}
    for label in ("s", "p", "d", "f"):
        values = np.where(dos_mask, np.asarray(dos[label]), np.nan)
        if np.any(np.abs(values) > 0):
            ax_pdos.plot(
                dos_energy, values, color=pdos_colors[label], linewidth=1.15,
                label=f"{label}-PDOS (core atom)",
            )
    if not show_occupied_dos:
        ax_dos.axvline(
            dos_cutoff_plot, color="0.45", linewidth=0.8, linestyle="--",
            label="lowest allowed final state",
        )
    ax_pdos.set_ylabel("core-atom PDOS (states/eV)")
    handles1, labels1 = ax_dos.get_legend_handles_labels()
    handles2, labels2 = ax_pdos.get_legend_handles_labels()
    ax_pdos.legend(handles1 + handles2, labels1 + labels2, loc="best", fontsize=9)
    if abs(transition_shift) > 1.0e-12:
        ax_dos.set_xlabel("Aligned final-state energy (eV)")
    else:
        ax_dos.set_xlabel(r"Final-state energy $E-E_{core}$ (eV)")

    hover_annotation = ax_xas.annotate(
        "", xy=(0, 0), xytext=(12, 15), textcoords="offset points",
        bbox={"boxstyle": "round", "fc": "white", "alpha": 0.94},
        arrowprops={"arrowstyle": "->", "color": "0.3"}, fontsize=8,
    )
    hover_annotation.set_visible(False)
    click_annotation = ax_xas.annotate(
        "", xy=(0, 0), xytext=(12, -55), textcoords="offset points",
        bbox={"boxstyle": "round", "fc": "#fff8dc", "alpha": 0.94},
        arrowprops={"arrowstyle": "->", "color": "0.3"}, fontsize=8,
    )
    click_annotation.set_visible(False)
    hover_state = {"index": None}

    def on_motion(event) -> None:
        if event.inaxes is not ax_xas or event.x is None or not line_artists:
            if hover_annotation.get_visible():
                hover_annotation.set_visible(False)
                hover_state["index"] = None
                fig.canvas.draw_idle()
            return
        display_x = np.asarray(
            [ax_xas.transData.transform((row["E_transition_plot"], 0.0))[0] for _, row in line_artists]
        )
        index = int(np.argmin(np.abs(display_x - event.x)))
        if abs(display_x[index] - event.x) > 6.0:
            if hover_annotation.get_visible():
                hover_annotation.set_visible(False)
                hover_state["index"] = None
                fig.canvas.draw_idle()
            return
        row = line_artists[index][1]
        hover_annotation.xy = (row["E_transition_plot"], event.ydata or 0.0)
        hover_annotation.set_text(short_transition(row))
        hover_annotation.set_visible(True)
        if hover_state["index"] != index:
            print_transition(row, prefix="[hover] ")
            hover_state["index"] = index
        fig.canvas.draw_idle()

    def on_click(event) -> None:
        if event.inaxes is not ax_xas or event.xdata is None or event.button != 1:
            return
        energy = float(event.xdata)
        nearby = [
            row for row in all_transitions
            if abs(row["E_transition_plot"] - energy) <= click_window
        ]
        nearby.sort(key=lambda row: row["weighted_guide"], reverse=True)
        selected = nearby[:top_n]
        print(
            f"\n[click] E={energy:.6f} eV; candidates within +/-{click_window:.3f} eV: "
            f"{len(nearby)} (showing {len(selected)})"
        )
        for i, row in enumerate(selected, 1):
            print_transition(row, prefix=f"  {i:2d}. ")
        if selected:
            lines = [
                f"all-k candidates near {energy:.3f} eV",
                *[
                    f"#{row['ik']} b{row['ilevel']} s{row['spin']}: "
                    f"XAS-axis={row['E_XAS_axis']:.3f} eV, "
                    f"guide={row['weighted_guide']:.3g}"
                    for row in selected[:5]
                ],
            ]
        else:
            lines = [f"No candidate within +/-{click_window:.3f} eV"]
        y = float(np.interp(energy, xas_energy, xas_y))
        click_annotation.xy = (energy, y)
        click_annotation.set_text("\n".join(lines))
        click_annotation.set_visible(True)
        fig.canvas.draw_idle()

    fig.canvas.mpl_connect("motion_notify_event", on_motion)
    fig.canvas.mpl_connect("button_press_event", on_click)
    fig.tight_layout()
    if figure_path is not None:
        figure_path.parent.mkdir(parents=True, exist_ok=True)
        fig.savefig(figure_path, dpi=180, bbox_inches="tight")
    if show:
        print("\nHover over a purple Gamma line for its final state.")
        print(f"Click the XAS curve to list all-k candidates within +/-{click_window:.3f} eV.")
        plt.show()
    else:
        plt.close(fig)


def normalize_mode(mode: str) -> str:
    aliases = {"iso": "avg", "x": "xx", "y": "yy", "z": "zz"}
    return aliases.get(mode.lower(), mode.lower())


def output_path(base: Path, user_value: str | None, default_name: str) -> Path:
    if user_value is None:
        return base / default_name
    path = Path(user_value).expanduser()
    return path if path.is_absolute() else base / path


def build_parser() -> argparse.ArgumentParser:
    parser = argparse.ArgumentParser(
        description="Assign VASP XAS features using EIGENVAL/PROCAR/DOSCAR final states.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument("--version", action="version", version=f"%(prog)s {SCRIPT_VERSION}")
    parser.add_argument(
        "mode", nargs="?", default="avg",
        choices=("avg", "iso", "xx", "yy", "zz", "xy", "yz", "zx", "x", "y", "z"),
        help="XAS tensor component; avg=(xx+yy+zz)/3",
    )
    parser.add_argument("directory", nargs="?", default=".", help="VASP calculation directory")
    parser.add_argument("--emin", type=float, help="minimum photon energy to plot")
    parser.add_argument("--emax", type=float, help="maximum photon energy to plot")
    parser.add_argument("--core-atom", type=int, help="1-based excited atom; default: INCAR CLNT")
    parser.add_argument("--core-label", help="core state such as 1s or 2p; default: INCAR CLN/CLL")
    parser.add_argument("--core-energy", type=float, help="override OUTCAR core-state eigenenergy (eV)")
    parser.add_argument(
        "--sigma", type=float, default=0.0,
        help="additional Gaussian sigma applied to both XAS and DOS (eV; 0 disables)",
    )
    parser.add_argument(
        "--align", choices=("first-peak", "onset", "none"), default="first-peak",
        help="align the lowest allowed guide to the selected XAS feature",
    )
    parser.add_argument(
        "--align-threshold", type=float, default=0.02,
        help="fraction of maximum XAS magnitude used to find onset/first peak",
    )
    parser.add_argument(
        "--transition-shift", type=float, default=0.0,
        help="additional manual shift after automatic alignment (eV)",
    )
    parser.add_argument(
        "--unoccupied-dos-only", action="store_true",
        help="hide occupied valence DOS and plot only the XAS final-state region",
    )
    parser.add_argument(
        "--show-occupied-dos", action="store_true", help=argparse.SUPPRESS,
    )
    parser.add_argument(
        "--min-unoccupied", type=float, default=0.01,
        help="minimum unoccupied fraction for transitions_all and interactive search",
    )
    parser.add_argument(
        "--min-local-p", type=float, default=1.0e-4,
        help="minimum core-atom p_total for a Gamma vertical line",
    )
    parser.add_argument("--max-lines", type=int, default=300, help="maximum Gamma lines displayed")
    parser.add_argument("--gamma-tol", type=float, default=1.0e-6, help="Gamma-coordinate tolerance")
    parser.add_argument(
        "--click-window", type=float, default=0.30,
        help="half-width for all-k candidates printed after an XAS click (eV)",
    )
    parser.add_argument("--top", type=int, default=10, help="number of all-k click candidates printed")
    parser.add_argument("--excel", help="Excel filename; relative paths are placed in calculation directory")
    parser.add_argument("--figure", help="figure filename; relative paths are placed in calculation directory")
    parser.add_argument("--no-excel", action="store_true", help="do not write the Excel workbook")
    parser.add_argument("--no-save", action="store_true", help="do not save the PNG figure")
    parser.add_argument("--no-show", action="store_true", help="do not open the interactive Matplotlib window")
    return parser


def validate_args(args: argparse.Namespace) -> None:
    if args.emin is not None and args.emax is not None and args.emin >= args.emax:
        raise VaspFormatError("--emin must be smaller than --emax.")
    if args.core_atom is not None and args.core_atom < 1:
        raise VaspFormatError("--core-atom must be at least 1.")
    if not 0.0 <= args.min_unoccupied <= 1.0:
        raise VaspFormatError("--min-unoccupied must be between 0 and 1.")
    if args.sigma < 0.0:
        raise VaspFormatError("--sigma must be non-negative.")
    if not 0.0 < args.align_threshold <= 1.0:
        raise VaspFormatError("--align-threshold must be greater than 0 and at most 1.")
    if args.min_local_p < 0 or args.max_lines < 1 or args.gamma_tol < 0:
        raise VaspFormatError("Projection/line/tolerance arguments must be non-negative.")
    if args.click_window < 0 or args.top < 1:
        raise VaspFormatError("--click-window must be non-negative and --top must be positive.")


def main(argv: Sequence[str] | None = None) -> int:
    args = build_parser().parse_args(argv)
    validate_args(args)
    show_occupied_dos = not args.unoccupied_dos_only
    mode = normalize_mode(args.mode)
    base = Path(args.directory).expanduser().resolve()
    if base.is_file():
        base = base.parent
    if not base.is_dir():
        raise VaspFormatError(f"Calculation directory not found: {base}")

    paths = {name: base / name for name in ("INCAR", "POSCAR", "OUTCAR", "EIGENVAL", "PROCAR", "DOSCAR")}
    for required in ("OUTCAR", "EIGENVAL", "PROCAR", "DOSCAR"):
        if not paths[required].is_file():
            raise VaspFormatError(f"Required file not found: {paths[required]}")

    print(f"Calculation directory: {base}")
    print(f"XAS component: {mode}")
    incar = read_incar(paths["INCAR"])
    if incar_bool(incar, "LNONCOLLINEAR") or incar_bool(incar, "LSORBIT"):
        raise VaspFormatError(
            "Non-collinear/SOC PROCAR and DOSCAR layouts are not supported by this version. "
            "Use a collinear ISPIN=1 or ISPIN=2 calculation."
        )
    outcar_text = read_text(paths["OUTCAR"])
    elements = parse_poscar_elements(paths["POSCAR"])
    ispin = int(parse_outcar_scalar(outcar_text, "ISPIN", incar_int(incar, "ISPIN", 1) or 1))
    xas_energy, tensor_raw = read_xas_tensor(outcar_text)
    tensor_display = broaden_tensor(xas_energy, tensor_raw, args.sigma)
    core_states = read_core_states(outcar_text, elements)
    selected_core = choose_core_state(
        core_states, incar, args.core_atom, args.core_label, args.core_energy
    )
    core_atom = int(selected_core["atom"])
    core_energy = float(selected_core["E_core"])
    print(
        f"Core state: atom {core_atom} {selected_core.get('element', '')} "
        f"{selected_core['core_state']}, Ecore={core_energy:.6f} eV"
    )

    eigen_levels, eigen_meta = read_eigenval(paths["EIGENVAL"], ispin)
    outcar_kpoints = read_outcar_fractional_kpoints(outcar_text, eigen_meta["nkpts"])
    k_coordinate_source = apply_fractional_kpoints(eigen_levels, outcar_kpoints)
    print(f"Fractional k-point source: {k_coordinate_source}")
    projections, orbitals_seen = read_procar(paths["PROCAR"], core_atom)
    dos_raw = read_doscar(paths["DOSCAR"], core_atom, ispin)
    dos_display = broaden_dos(dos_raw, args.sigma)
    outcar_fermi = parse_outcar_fermi(outcar_text)
    dos_fermi = float(dos_raw["Efermi"])
    dos_reference_shift = 0.0 if outcar_fermi is None else outcar_fermi - dos_fermi
    if abs(dos_reference_shift) > 1.0e-7:
        print(
            f"DOS reference correction from E-fermi: {dos_reference_shift:+.6f} eV "
            f"(OUTCAR={outcar_fermi:.6f}, DOSCAR={dos_fermi:.6f})"
        )
    occ_full = infer_full_occupation(eigen_levels)
    print(f"Full EIGENVAL occupation inferred from data: {occ_full:g}")
    levels, missing = build_level_rows(
        eigen_levels, projections, core_energy, mode, occ_full,
        args.gamma_tol,
    )
    if missing:
        print(
            f"Warning: PROCAR projections are missing for {missing}/{len(levels)} "
            "EIGENVAL spin/k-point/band rows. Missing projections are blank in Excel."
        )
    if not {"px", "py", "pz"}.issubset(orbitals_seen):
        print(
            "Warning: PROCAR does not contain separate px/py/pz columns. "
            "Directional guides use p_total/3; use LORBIT=11 for component assignment."
        )
    if mode in ("xy", "yz", "zx"):
        print(
            "Warning: off-diagonal PROCAR guides use sqrt(Pi*Pj). They cannot reproduce "
            "the sign or phase of the true off-diagonal transition matrix element."
        )

    transitions = [
        row for row in levels
        if row["E-Ecore"] > 0.0 and row["unoccupied_fraction"] >= args.min_unoccupied
    ]
    alignment = align_transition_energies(
        transitions, xas_energy, tensor_display[mode], args.align,
        args.align_threshold, args.min_local_p, args.transition_shift,
    )
    total_shift = float(alignment["total_shift_eV"])
    dos_cutoff_plot = float(alignment["transition_anchor_raw_eV"]) + total_shift
    if args.align == "none":
        print(
            f"Energy alignment: none; manual shift={args.transition_shift:+.6f} eV, "
            f"total shift={total_shift:+.6f} eV"
        )
    else:
        print(
            f"Energy alignment: {args.align}; XAS anchor={alignment['xas_anchor_eV']:.6f} eV, "
            f"raw transition anchor={alignment['transition_anchor_raw_eV']:.6f} eV, "
            f"auto shift={alignment['auto_shift_eV']:+.6f} eV, "
            f"manual shift={args.transition_shift:+.6f} eV"
        )
    if args.sigma > 0.0:
        print(f"Additional Gaussian broadening: sigma={args.sigma:g} eV (XAS and DOS)")
        ch_sigma = incar_float(incar, "CH_SIGMA")
        if ch_sigma is not None:
            print(
                f"Note: this is post-processing broadening on top of INCAR "
                f"CH_SIGMA={ch_sigma:g} eV. It cannot restore intensity lost to an "
                "under-resolved VASP output mesh."
            )
    emin = float(xas_energy.min()) if args.emin is None else max(float(args.emin), float(xas_energy.min()))
    emax = float(xas_energy.max()) if args.emax is None else min(float(args.emax), float(xas_energy.max()))
    if emin >= emax:
        raise VaspFormatError("Requested plot range does not overlap the XAS energy grid.")
    gamma_rows, omitted = choose_gamma_lines(
        transitions, emin, emax, args.min_local_p, args.max_lines
    )
    if omitted:
        print(
            f"Gamma lines: displaying the {len(gamma_rows)} strongest guides; "
            f"{omitted} additional candidates remain in Excel."
        )
    else:
        print(f"Gamma lines in plot range: {len(gamma_rows)}")
    if not any(row["is_gamma"] for row in levels):
        print("Warning: EIGENVAL contains no Gamma-equivalent k point; no vertical lines are drawn.")

    in_range = [
        row["E_transition_plot"] for row in transitions
        if emin <= row["E_transition_plot"] <= emax
    ]
    if not in_range:
        print(
            "Warning: no EIGENVAL transition energy overlaps the XAS plot range. "
            "Verify that OUTCAR/EIGENVAL belong to the same calculation and inspect "
            "Ecore; inspect --align and the recorded alignment values in Excel."
        )

    for row in core_states:
        row["selected"] = (
            row["atom"] == core_atom
            and row["core_state"] == selected_core["core_state"]
        )
    excel_path = None if args.no_excel else output_path(base, args.excel, f"xas_assignment_{mode}.xlsx")
    figure_path = None if args.no_save else output_path(base, args.figure, f"xas_assignment_{mode}.png")
    config_rows = [
        {"parameter": "calculation_directory", "value": str(base)},
        {"parameter": "mode", "value": mode},
        {"parameter": "core_atom", "value": core_atom},
        {"parameter": "core_element", "value": selected_core.get("element", "")},
        {"parameter": "core_state", "value": selected_core["core_state"]},
        {"parameter": "E_core_eV", "value": core_energy},
        {"parameter": "energy_alignment", "value": alignment["method"]},
        {"parameter": "align_threshold", "value": args.align_threshold},
        {"parameter": "XAS_anchor_eV", "value": alignment["xas_anchor_eV"]},
        {"parameter": "raw_transition_anchor_eV", "value": alignment["transition_anchor_raw_eV"]},
        {"parameter": "automatic_shift_eV", "value": alignment["auto_shift_eV"]},
        {"parameter": "manual_transition_shift_eV", "value": args.transition_shift},
        {"parameter": "total_shift_eV", "value": total_shift},
        {"parameter": "DOSCAR_Efermi_eV", "value": dos_fermi},
        {"parameter": "OUTCAR_Efermi_eV", "value": outcar_fermi},
        {"parameter": "DOS_reference_correction_eV", "value": dos_reference_shift},
        {"parameter": "DOS_unoccupied_cutoff_plot_eV", "value": dos_cutoff_plot},
        {"parameter": "show_occupied_DOS", "value": show_occupied_dos},
        {
            "parameter": "DOS_display_scope",
            "value": "occupied + unoccupied" if show_occupied_dos else "unoccupied only",
        },
        {"parameter": "postprocess_sigma_eV", "value": args.sigma},
        {"parameter": "INCAR_CH_SIGMA_eV", "value": incar_float(incar, "CH_SIGMA")},
        {"parameter": "ISPIN", "value": ispin},
        {"parameter": "full_occupation_inferred", "value": occ_full},
        {"parameter": "NKPTS", "value": eigen_meta["nkpts"]},
        {"parameter": "NBANDS", "value": eigen_meta["nbands"]},
        {"parameter": "min_unoccupied", "value": args.min_unoccupied},
        {"parameter": "min_local_p", "value": args.min_local_p},
        {"parameter": "gamma_tolerance", "value": args.gamma_tol},
        {"parameter": "k_coordinate_basis", "value": "reciprocal-lattice fractional"},
        {"parameter": "k_coordinate_source", "value": k_coordinate_source},
        {"parameter": "Gamma_lines_plotted", "value": len(gamma_rows)},
        {"parameter": "assignment_note", "value": "PROCAR projection guide; not VASP dipole matrix element"},
    ]
    if excel_path is not None:
        print(f"Writing Excel workbook: {excel_path}")
        write_excel(
            excel_path, config_rows, core_states, transitions, gamma_rows,
            xas_energy, tensor_raw, tensor_display, mode, dos_raw, dos_display,
            core_energy, total_shift, dos_reference_shift, dos_cutoff_plot,
        )

    title = (
        f"VASP XAS assignment: {selected_core.get('element', '')} atom {core_atom} "
        f"{selected_core['core_state']} ({mode})"
    )
    plot_results(
        xas_energy, tensor_display[mode], mode, dos_display, core_energy, total_shift,
        gamma_rows, transitions, emin, emax, args.click_window, args.top,
        title, args.sigma, dos_reference_shift, dos_cutoff_plot,
        show_occupied_dos, figure_path, not args.no_show,
    )
    if figure_path is not None:
        print(f"Saved figure: {figure_path}")
    print(
        "Assignment reminder: line height/opacity follows k-weight x unoccupied fraction "
        "x local p guide, not the actual CH_LSPEC transition matrix element."
    )
    return 0


if __name__ == "__main__":
    try:
        raise SystemExit(main())
    except (VaspFormatError, RuntimeError, OSError, ValueError) as exc:
        print(f"Error: {exc}", file=sys.stderr)
        raise SystemExit(2)