Source code for mag4.energy
"""DFT total-energy parsers (VASP ``OUTCAR``, WIEN2k ``case.scf``).
The :func:`get_energy` dispatcher returns an :class:`EnergyResult` describing
both the converged total energy (in **eV**) and a structured convergence
status, which the upcoming ``mag4-extract`` CLI uses to decide whether each
configuration is safe to feed into the J-extraction formulas.
Only the fields actually needed for the four-state pipeline are extracted:
* the **last** total energy printed in the file (assumed converged);
* a coarse convergence flag (the calculation finished cleanly);
* for WIEN2k: the last ``:DIS`` (charge distance) value and the verbatim
``:DIS`` / ``:ENE`` lines, so the user can verify the parser picked up
the right data and that the SCF converged tightly enough.
Both parsers are **streaming** (line-by-line) and tolerate large output
files without loading them entirely into memory.
"""
from __future__ import annotations
import os
import re
from dataclasses import dataclass, field
from typing import Dict, List, Optional
#: 1 Rydberg in electron-volts (CODATA 2018).
RY_TO_EV: float = 13.605693122994
[docs]
@dataclass
class EnergyResult:
"""Outcome of parsing one DFT output file.
Attributes
----------
energy_eV : float or None
Last total energy reported in the file, expressed in **eV**. ``None``
when no energy could be parsed (typically: the calculation never
produced one — crashed early, file truncated, etc.).
converged : bool
Coarse convergence flag (see backend-specific notes in
:func:`extract_vasp_energy` / :func:`extract_wien2k_energy`).
source : str
Absolute path of the file that was parsed.
n_iter : int
Number of (electronic) iterations detected. ``0`` if unknown.
diagnostics : list of str
Free-form one-line warnings / notes the caller may want to surface.
"""
energy_eV: Optional[float]
converged: bool
source: str
n_iter: int = 0
diagnostics: List[str] = field(default_factory=list)
#: Last ``:DIS`` (charge distance) numerical value, in WIEN2k units.
#: Populated only for WIEN2k; ``None`` for VASP or when no ``:DIS`` line
#: was matched. Useful as a tight convergence check the user can eyeball.
last_dis: Optional[float] = None
#: Verbatim ``:DIS`` line from the SCF file (last occurrence).
last_dis_line: Optional[str] = None
#: Verbatim ``:ENE`` line whose value was used as ``energy_eV``.
last_ene_line: Optional[str] = None
# ---------------------------------------------------------------------------
# VASP — OUTCAR
# ---------------------------------------------------------------------------
# Matches lines like: " free energy TOTEN = -1234.56789012 eV"
_VASP_TOTEN_RE = re.compile(
r"free\s+energy\s+TOTEN\s*=\s*([+-]?\d+\.\d+)\s*eV", re.IGNORECASE
)
# Marks the end of an electronic step: " RMM: 3 -123.45... ..."
# We use these to count iterations and to capture the last electronic step.
_VASP_NELM_HIT_RE = re.compile(
r"WARNING:\s+The number of electronic SCF steps\s+NELM\s*=\s*(\d+)",
re.IGNORECASE,
)
_VASP_END_OF_RUN = "Voluntary context switches:"
_VASP_REQUIRED_ACCURACY = "reached required accuracy"
#: VASP's electronic-SCF convergence marker — "aborting loop because EDIFF is
#: reached". This is the authoritative *converged* signal, written when the SCF
#: meets EDIFF, BEFORE the final energy block and long before the process's
#: ``Voluntary context switches:`` tail. A job OOM-killed / timed-out during the
#: final clean-up loses the tail but is still physically converged, so we accept
#: it on this marker.
_VASP_EDIFF_REACHED = "EDIFF is reached"
[docs]
def extract_vasp_energy(outcar_path: str) -> EnergyResult:
"""Parse a VASP ``OUTCAR`` and return its converged total energy.
Convergence rule (simple, code-agnostic):
``converged`` is ``True`` when **all** of the following hold:
* the file exists and is readable,
* at least one ``free energy TOTEN`` line was parsed,
* **no** ``WARNING: ... NELM = N`` line was emitted,
* the SCF actually converged — either VASP's electronic marker
``aborting loop because EDIFF is reached`` or the ionic
``reached required accuracy`` appears — **or** the run finished
cleanly (``Voluntary context switches:``).
The clean-finish marker is sufficient but **not** required: a job
OOM-killed or timed out during the final clean-up drops that tail
without invalidating the already-converged energy, so the EDIFF
marker rescues it (see ``diagnostics``).
Tighter convergence checks (per-element forces, stress, EDIFF /
EDIFFG numerical thresholds) are intentionally NOT enforced here —
they vary per project; see ``diagnostics`` for the raw flags so the
caller can decide what to do.
Returns
-------
EnergyResult
Always non-``None``; check ``.energy_eV`` and ``.converged`` to
interpret the outcome.
"""
if not os.path.isfile(outcar_path):
return EnergyResult(None, False, outcar_path,
diagnostics=[f"file not found: {outcar_path}"])
last_energy: Optional[float] = None
n_iter = 0
nelm_hit = False
end_marker_seen = False
relax_converged = False
ediff_reached = False
with open(outcar_path, encoding="utf-8", errors="replace") as fh:
for line in fh:
m = _VASP_TOTEN_RE.search(line)
if m:
last_energy = float(m.group(1))
n_iter += 1
continue
if _VASP_NELM_HIT_RE.search(line):
nelm_hit = True
continue
if _VASP_END_OF_RUN in line:
end_marker_seen = True
continue
if _VASP_EDIFF_REACHED in line:
ediff_reached = True
continue
if _VASP_REQUIRED_ACCURACY in line:
relax_converged = True
continue
# The SCF physically converged when VASP reached EDIFF (electronic) or
# the required accuracy (ionic). A clean process finish is a *separate*
# signal — sufficient but not necessary, since an OOM/timeout kill during
# the final clean-up drops the tail without touching the converged energy.
scf_converged = ediff_reached or relax_converged
diags: List[str] = []
if last_energy is None:
diags.append("no 'free energy TOTEN' line found")
if not end_marker_seen:
if scf_converged:
diags.append("run did not finish cleanly (no 'Voluntary context "
"switches:') but the SCF converged (EDIFF reached) — "
"accepting the converged energy (likely OOM/timeout "
"killed during final clean-up)")
else:
diags.append("OUTCAR does not end with 'Voluntary context switches:' "
"and no convergence marker found (run may have been "
"killed before the SCF converged)")
if nelm_hit:
diags.append("NELM was reached on at least one electronic step "
"(SCF may not have converged)")
if relax_converged:
diags.append("ionic relaxation reached required accuracy")
converged = ((last_energy is not None) and (not nelm_hit)
and (scf_converged or end_marker_seen))
return EnergyResult(
energy_eV=last_energy,
converged=converged,
source=outcar_path,
n_iter=n_iter,
diagnostics=diags,
)
# ---------------------------------------------------------------------------
# WIEN2k — case.scf
# ---------------------------------------------------------------------------
# Matches: ":ENE : ********** TOTAL ENERGY IN Ry = -1234.56789012"
_WIEN_ENE_RE = re.compile(
r":ENE\s*:.*TOTAL\s+ENERGY\s+IN\s+Ry\s*=\s*([+-]?\d+\.\d+)",
re.IGNORECASE,
)
# Matches any flavour of charge-distance line and grabs the first float.
# WIEN2k variants in the wild include:
# ":DIS : CHARGE DISTANCE ( 0.0000123 for atom 1 spin 1)"
# ":DIS001: CHARGE DISTANCE ( 0.00001234 for atom 1 spin 1)"
# ":DIS : TOTAL CHARGE DISTANCE = 0.0000123"
_WIEN_DIS_RE = re.compile(
r":DIS\s*\d*\s*:.*?([+-]?\d+\.\d+(?:[eE][+-]?\d+)?)",
re.IGNORECASE,
)
# Iteration banner. The literal word "ITERATION" disappeared in modern
# WIEN2k output, so we only require ``:ITE<digits>:``.
_WIEN_ITE_RE = re.compile(r":ITE\s*\d+\s*:", re.IGNORECASE)
[docs]
def extract_wien2k_energy(scf_path: str) -> EnergyResult:
"""Parse a WIEN2k SCF file and return its converged total energy in **eV**.
Convergence rule (simple):
``converged`` is ``True`` when **all** of the following hold:
* the file exists and is readable,
* at least one ``:ENE`` line was parsed,
* at least one ``:ITE…:`` iteration block was detected,
* the energy was emitted by the **last** iteration block of the file
(i.e. the calculation didn't crash mid-iteration).
This does **not** enforce a numeric threshold on ``:DIS`` (charge
distance) — the value is surfaced in :attr:`EnergyResult.last_dis`
and the verbatim ``:DIS`` line in :attr:`EnergyResult.last_dis_line`
so the user can eyeball that the SCF converged tightly enough.
Conversion: ``E [eV] = E [Ry] × 13.605693122994``.
"""
if not os.path.isfile(scf_path):
return EnergyResult(None, False, scf_path,
diagnostics=[f"file not found: {scf_path}"])
last_energy_ry: Optional[float] = None
last_ene_line: Optional[str] = None
last_dis_value: Optional[float] = None
last_dis_line: Optional[str] = None
n_iter = 0
last_ite_idx = -1
last_ene_line_idx = -1
with open(scf_path, encoding="utf-8", errors="replace") as fh:
for idx, line in enumerate(fh):
if _WIEN_ITE_RE.search(line):
n_iter += 1
last_ite_idx = idx
continue
m = _WIEN_ENE_RE.search(line)
if m:
last_energy_ry = float(m.group(1))
last_ene_line = line.rstrip()
last_ene_line_idx = idx
continue
m = _WIEN_DIS_RE.search(line)
if m:
last_dis_value = float(m.group(1))
last_dis_line = line.rstrip()
continue
diags: List[str] = []
if last_energy_ry is None:
diags.append("no ':ENE' line found")
if n_iter == 0:
diags.append("no ':ITE…:' iteration block detected")
energy_ev = last_energy_ry * RY_TO_EV if last_energy_ry is not None else None
# Sanity: the last :ENE we saw must come from the last iteration
# block (not from an earlier one followed by a crash).
energy_in_last_ite = (last_ene_line_idx > last_ite_idx) if n_iter > 0 else False
converged = (last_energy_ry is not None) and (n_iter > 0) and energy_in_last_ite
if not energy_in_last_ite and last_energy_ry is not None and n_iter > 0:
diags.append("last :ENE is NOT inside the last :ITE block — "
"calculation may have been interrupted")
return EnergyResult(
energy_eV=energy_ev,
converged=converged,
source=scf_path,
n_iter=n_iter,
diagnostics=diags,
last_dis=last_dis_value,
last_dis_line=last_dis_line,
last_ene_line=last_ene_line,
)
# ---------------------------------------------------------------------------
# Dispatcher
# ---------------------------------------------------------------------------
[docs]
def get_energy(config_dir: str, code: str,
scf_basename: Optional[str] = None) -> EnergyResult:
"""Convenience dispatcher: parse the right output for the requested code.
Parameters
----------
config_dir : str
Path to a single ``configN/`` directory.
code : {'vasp', 'wien2k'}
Which output file to look for inside ``config_dir``:
* ``vasp`` → ``OUTCAR``
* ``wien2k`` → SCF file (see ``scf_basename``).
scf_basename : str, optional
WIEN2k only. Basename (with or without the ``.scf`` suffix) of the
SCF file to parse inside ``config_dir``. Useful when several SCF
runs are stored side by side under names like ``pbe.scf``,
``r2scan.scf``. When omitted (default) the standard WIEN2k
convention ``<basename(config_dir)>.scf`` is used — i.e.
``config5/config5.scf`` for ``config5/``.
"""
code = code.lower()
if code == "vasp":
return extract_vasp_energy(os.path.join(config_dir, "OUTCAR"))
if code == "wien2k":
if scf_basename:
base = scf_basename
if base.lower().endswith(".scf"):
base = base[:-4]
else:
base = os.path.basename(os.path.normpath(config_dir))
return extract_wien2k_energy(os.path.join(config_dir, f"{base}.scf"))
raise ValueError(f"Unknown DFT code: {code!r}. Choose 'vasp' or 'wien2k'.")
# ---------------------------------------------------------------------------
# Total magnetization (for ``mag4-extract --method energy-mapping --check-mag``)
# ---------------------------------------------------------------------------
#: VASP OSZICAR ionic-step line: "... d E =0.0 mag= 5.0000"
_VASP_OSZ_MAG_RE = re.compile(r"\bmag=\s*([-+]?\d+(?:\.\d+)?(?:[eE][-+]?\d+)?)")
#: VASP OUTCAR: " number of electron 176.000 magnetization 5.000"
_VASP_OUT_MAG_RE = re.compile(
r"number of electron\s+\S+\s+magnetization\s+"
r"([-+]?\d+(?:\.\d+)?(?:[eE][-+]?\d+)?)")
#: WIEN2k case.scf: ":MMTOT: TOTAL MAGNETIC MOMENT IN CELL = 5.00000"
_WIEN_MMTOT_RE = re.compile(r":MMTOT\s*:.*?=\s*([-+]?\d+\.\d+)")
[docs]
def read_vasp_magnetization(config_dir: str) -> Optional[float]:
"""Total magnetization (μ_B) of a converged VASP run, or ``None``.
Prefers the **last** ``mag=`` in ``OSZICAR``; falls back to the ``OUTCAR``
``number of electron … magnetization …`` line. Collinear (``ISPIN=2``)
runs report a single scalar moment.
"""
osz = os.path.join(config_dir, "OSZICAR")
if os.path.isfile(osz):
val = None
with open(osz) as fh:
for line in fh:
m = _VASP_OSZ_MAG_RE.search(line)
if m:
val = float(m.group(1))
if val is not None:
return val
out = os.path.join(config_dir, "OUTCAR")
if os.path.isfile(out):
val = None
with open(out) as fh:
for line in fh:
m = _VASP_OUT_MAG_RE.search(line)
if m:
val = float(m.group(1))
return val
return None
#: VASP OUTCAR per-ion moment table header (written when LORBIT=10/11):
#: magnetization (x)
#: # of ion s p d tot
_VASP_MAG_BLOCK = "magnetization (x)"
_VASP_ION_ROW_RE = re.compile(
r"^\s*(\d+)\s+((?:[-+]?\d+\.\d+\s+)*[-+]?\d+\.\d+)\s*$")
[docs]
def read_vasp_site_moments(config_dir: str,
ion_indices) -> Dict[int, Optional[float]]:
"""Per-ion magnetic moments (μ_B) from a VASP ``OUTCAR``.
Parses the **last** ``magnetization (x)`` block (the converged one; needs
``LORBIT=10/11`` in the INCAR — mag4 sets ``LORBIT=11`` by default) and
returns the ``tot`` column (last field of each ion row) for the requested
1-based ion indices. Ion order follows the POSCAR; mag4 writes the
magnetic atoms first, so they are ions ``1 … n_mag``.
Returns ``{ion: float or None}`` (``None`` for an ion the block didn't
list, or ``{}``-like all-None when no block/OUTCAR was found).
"""
out_path = os.path.join(config_dir, "OUTCAR")
moments: Dict[int, Optional[float]] = {int(i): None for i in ion_indices}
if not os.path.isfile(out_path):
return moments
with open(out_path, errors="replace") as fh:
lines = fh.readlines()
# locate the last "magnetization (x)" block
starts = [i for i, ln in enumerate(lines) if _VASP_MAG_BLOCK in ln]
if not starts:
return moments
i = starts[-1]
# advance to the "# of ion" header, then past the dashed separator
n = len(lines)
while i < n and "# of ion" not in lines[i]:
i += 1
i += 1
if i < n and set(lines[i].strip()) <= {"-"}:
i += 1
# read ion rows until a non-matching / separator / "tot" line
while i < n:
ln = lines[i]
if "tot" in ln or set(ln.strip()) <= {"-"} or not ln.strip():
break
m = _VASP_ION_ROW_RE.match(ln)
if not m:
break
ion = int(m.group(1))
tot = float(m.group(2).split()[-1]) # last column = total
if ion in moments:
moments[ion] = tot
i += 1
return moments
[docs]
def read_site_moments(config_dir: str, code: str, ion_indices,
scf_basename: Optional[str] = None
) -> Dict[int, Optional[float]]:
"""Per-site magnetic moments for the requested code (dispatcher)."""
code = code.lower()
if code == "vasp":
return read_vasp_site_moments(config_dir, ion_indices)
if code == "wien2k":
from mag4.analysis.wien2k import extract_moments_from_scf, _resolve_scf
return extract_moments_from_scf(
_resolve_scf(config_dir, scf_basename), ion_indices)
return {int(i): None for i in ion_indices}
[docs]
def read_wien2k_magnetization(config_dir: str,
scf_basename: Optional[str] = None) -> Optional[float]:
"""Total magnetic moment in the cell (μ_B) from the last ``:MMTOT`` line.
The SCF file is resolved exactly as in :func:`get_energy`
(``<basename(config_dir)>.scf`` unless ``scf_basename`` overrides it).
"""
if scf_basename:
base = (scf_basename[:-4]
if scf_basename.lower().endswith(".scf") else scf_basename)
else:
base = os.path.basename(os.path.normpath(config_dir))
scf = os.path.join(config_dir, f"{base}.scf")
if not os.path.isfile(scf):
return None
val = None
with open(scf) as fh:
for line in fh:
m = _WIEN_MMTOT_RE.search(line)
if m:
val = float(m.group(1))
return val
[docs]
def read_magnetization(config_dir: str, code: str,
scf_basename: Optional[str] = None) -> Optional[float]:
"""Total magnetization (μ_B) for the requested code, or ``None``."""
code = code.lower()
if code == "vasp":
return read_vasp_magnetization(config_dir)
if code == "wien2k":
return read_wien2k_magnetization(config_dir, scf_basename)
return None
# ---------------------------------------------------------------------------
# Band gap / metallicity (catch a config that converged to a metallic — i.e.
# usually wrong — solution for a magnetic insulator)
# ---------------------------------------------------------------------------
[docs]
@dataclass
class BandGap:
"""Band gap of one configuration (eV).
``gap_eV`` is the conduction–valence separation (0.0 when metallic);
``metallic`` is ``True`` when occupied and empty states overlap (no gap),
``None`` when the file carried no usable eigenvalue information.
"""
gap_eV: Optional[float]
metallic: Optional[bool]
vbm_eV: Optional[float] = None
cbm_eV: Optional[float] = None
source: str = ""
note: str = ""
#: VASP OUTCAR eigenvalue line: `` 12 -1.2345 1.00000`` (band, E, occ).
_VASP_EIG_RE = re.compile(
r"^\s*\d+\s+(-?\d+\.\d+(?:[eE][+-]?\d+)?)\s+(\d+\.\d+(?:[eE][+-]?\d+)?)\s*$")
def _gap_from_states(occ_e, unocc_e, gap_tol):
"""``(gap_eV, metallic, vbm, cbm)`` from occupied/empty band energies (eV)."""
if not occ_e or not unocc_e:
return None, None, None, None
vbm, cbm = max(occ_e), min(unocc_e)
gap = cbm - vbm
return max(0.0, gap), (gap <= gap_tol), vbm, cbm
[docs]
def extract_vasp_gap(outcar_path: str, *, occ_tol: float = 0.5,
gap_tol: float = 0.05) -> BandGap:
"""Band gap from the eigenvalue/occupation block at the end of a VASP OUTCAR.
Reads the **last** eigenvalue dump (after the final ``E-fermi`` line, all
k-points and spin channels), splits bands into occupied (occ > ``occ_tol``)
and empty, and reports ``CBM − VBM``. Metallic when they overlap.
"""
if not os.path.isfile(outcar_path):
return BandGap(None, None, note="OUTCAR not found")
with open(outcar_path, errors="replace") as fh:
lines = fh.readlines()
start = None
for i, ln in enumerate(lines):
if "E-fermi" in ln:
start = i
if start is None:
return BandGap(None, None, note="no eigenvalues (no E-fermi) in OUTCAR")
occ_e, unocc_e, frac = [], [], False
for ln in lines[start:]:
m = _VASP_EIG_RE.match(ln)
if not m:
continue
e, occ = float(m.group(1)), float(m.group(2))
(occ_e if occ > occ_tol else unocc_e).append(e)
if 0.01 < occ < 0.99:
frac = True
gap, metallic, vbm, cbm = _gap_from_states(occ_e, unocc_e, gap_tol)
if gap is None:
return BandGap(None, None, note="could not split occupied/empty bands")
# A partially occupied band means E_F crosses it → metallic, even if the
# smearing-split sub-states sit a little apart in energy.
metallic = metallic or frac
return BandGap(0.0 if metallic else gap, metallic, vbm, cbm,
source="VASP OUTCAR eigenvalues",
note=("partial occupations — E_F crosses a band" if frac
else "bands cross E_F" if metallic else ""))
#: WIEN2k ``:GAP`` line (insulator): ``:GAP : 0.0924 Ry = 1.257 eV``.
_WIEN_GAP_EV_RE = re.compile(r":GAP\b[^=]*=\s*(-?\d+\.\d+)\s*eV", re.IGNORECASE)
_WIEN_GAP_RY_RE = re.compile(r":GAP\b\s*:?\s*(-?\d+\.\d+)\s*Ry", re.IGNORECASE)
#: WIEN2k ``:BANnnnnn: band emin emax occ`` (Ry).
_WIEN_BAN_RE = re.compile(
r":BAN\s*\d+\s*:\s*\d+\s+(-?\d+\.\d+)\s+(-?\d+\.\d+)\s+(\d+\.\d+)")
[docs]
def extract_wien2k_gap(scf_path: str, *, occ_tol: float = 0.5,
gap_tol: float = 0.05) -> BandGap:
"""Band gap from a WIEN2k SCF file.
Prefers the ``:GAP`` line (printed when lapw2 finds a gap). Otherwise — e.g.
with temperature smearing — uses the ``Bandranges`` / ``:BANnnnnn`` block:
VBM = highest ``emax`` of an occupied band (occ > ``occ_tol``), CBM = lowest
``emin`` of an empty band; energies are in Ry and converted to eV.
"""
if not os.path.isfile(scf_path):
return BandGap(None, None, note="scf file not found")
with open(scf_path, errors="replace") as fh:
lines = fh.readlines()
# Spin-polarised runs print per-channel ":GAP (this spin)" lines AND the
# physical ":GAP (global)" one (min over both channels; VBM and CBM may
# even sit in different channels). Prefer the global value explicitly —
# a per-spin gap can badly overstate the true gap (a half-metal has a
# gapped channel yet is globally metallic). Fallback order: last
# "(global)" line > last unqualified ":GAP" line > min over the final
# two "(this spin)" lines (the two channels of the last iteration).
gap_global = None
gap_plain = None
gap_spin: List[float] = []
for ln in lines:
m = _WIEN_GAP_EV_RE.search(ln)
if m:
val = float(m.group(1))
else:
m = _WIEN_GAP_RY_RE.search(ln)
if not m:
continue
val = float(m.group(1)) * RY_TO_EV
low = ln.lower()
if "global" in low:
gap_global = val
elif "this spin" in low:
gap_spin.append(val)
else:
gap_plain = val
if gap_global is not None:
gap_ev, src = gap_global, "WIEN2k :GAP (global)"
elif gap_plain is not None:
gap_ev, src = gap_plain, "WIEN2k :GAP"
elif gap_spin:
gap_ev, src = min(gap_spin[-2:]), "WIEN2k :GAP (min over spins)"
else:
gap_ev, src = None, ""
if gap_ev is not None:
note = "bands overlap" if gap_ev <= gap_tol else ""
if src.endswith("(min over spins)"):
note = (note + "; " if note else "") + \
"no :GAP (global) line — per-spin minimum used"
return BandGap(max(0.0, gap_ev), gap_ev <= gap_tol,
source=src, note=note)
start = None
for i, ln in enumerate(lines):
if "Bandranges" in ln:
start = i
region = lines[start:] if start is not None else lines
occ_emax, unocc_emin, frac = [], [], False
for ln in region:
m = _WIEN_BAN_RE.search(ln)
if not m:
continue
emin, emax, occ = float(m.group(1)), float(m.group(2)), float(m.group(3))
(occ_emax.append(emax) if occ > occ_tol else unocc_emin.append(emin))
if 0.01 < occ < 0.99:
frac = True
if occ_emax and unocc_emin:
vbm, cbm = max(occ_emax), min(unocc_emin)
gap = (cbm - vbm) * RY_TO_EV
metallic = (gap <= gap_tol) or frac
return BandGap(0.0 if metallic else gap, metallic,
vbm_eV=vbm * RY_TO_EV, cbm_eV=cbm * RY_TO_EV,
source="WIEN2k :BAN ranges",
note=("partial band occupations — metallic" if frac
else "bands overlap" if metallic else ""))
return BandGap(None, None, note="no :GAP line or :BAN block in scf")
[docs]
def read_band_gap(config_dir: str, code: str,
scf_basename: Optional[str] = None, *,
gap_tol: float = 0.05) -> BandGap:
"""Band gap (eV) + metallicity for one config, or ``BandGap(None, None)``."""
code = code.lower()
if code == "vasp":
return extract_vasp_gap(os.path.join(config_dir, "OUTCAR"), gap_tol=gap_tol)
if code == "wien2k":
if scf_basename:
base = (scf_basename[:-4]
if scf_basename.lower().endswith(".scf") else scf_basename)
else:
base = os.path.basename(os.path.normpath(config_dir))
return extract_wien2k_gap(os.path.join(config_dir, f"{base}.scf"),
gap_tol=gap_tol)
return BandGap(None, None, note="unknown code")