defectpl.io.vasp

All VASP-specific I/O lives in defectpl.io.vasp.

Eigenvalue parsing

read_eigenval_file

read_eigenval_file(filename, k_idx=0)

Parse a VASP EIGENVAL file and return spin-resolved eigenvalues at one k-point.

Parameters:

Name	Type	Description	Default
filename	str or Path		required
k_idx	int	Zero-based k-point index (default 0, i.e. Γ-point).	0

Returns:

Type	Description
dict — see :func:`defectpl.vasp.read_eigenval_file` for full key list.

Raises:

Type	Description
ValueError	If the calculation is not spin-polarised (ISPIN ≠ 2).

Source code in defectpl/io/vasp.py

def read_eigenval_file(filename: Union[str, Path], k_idx: int = 0) -> Dict[str, Any]:
    """
    Parse a VASP EIGENVAL file and return spin-resolved eigenvalues at one k-point.

    Parameters
    ----------
    filename : str or Path
    k_idx : int
        Zero-based k-point index (default 0, i.e. Γ-point).

    Returns
    -------
    dict  — see :func:`defectpl.vasp.read_eigenval_file` for full key list.

    Raises
    ------
    ValueError
        If the calculation is not spin-polarised (ISPIN ≠ 2).
    """
    from pymatgen.electronic_structure.core import Spin
    from pymatgen.io.vasp.outputs import Eigenval

    from defectpl.ks_analysis import get_homo_lumo_idx

    data: Dict[str, Any] = {}
    eig = Eigenval(filename, separate_spins=True)
    if eig.ispin != 2:
        raise ValueError("The calculation is not spin polarized.")

    print(f"Selecting the {k_idx}-th k-point from {eig.nkpt} k-points.")
    print(f"Selected k-point: {eig.kpoints[k_idx]}")

    data["up"] = list(eig.eigenvalues[Spin.up][k_idx])
    data["down"] = list(eig.eigenvalues[Spin.down][k_idx])

    data["homo_up_idx"], data["lumo_up_idx"] = get_homo_lumo_idx(data["up"])
    data["homo_down_idx"], data["lumo_down_idx"] = get_homo_lumo_idx(data["down"])

    data["homo_up"] = eig.eigenvalue_band_properties[2][0]
    data["homo_down"] = eig.eigenvalue_band_properties[2][1]
    data["lumo_up"] = eig.eigenvalue_band_properties[1][0]
    data["lumo_down"] = eig.eigenvalue_band_properties[1][1]
    data["hl_gap_up"] = eig.eigenvalue_band_properties[0][0]
    data["hl_gap_down"] = eig.eigenvalue_band_properties[0][1]

    data["nelect"] = eig.nelect
    data["nbands"] = eig.nbands
    data["nkpt"] = eig.nkpt
    data["selected_kpoint"] = [k_idx, list(eig.kpoints[k_idx])]
    data["spin_multiplicity"] = get_spin_multiplicity(
        data["homo_up_idx"], data["homo_down_idx"]
    )
    return data

get_spin_multiplicity

get_spin_multiplicity(homo_up_idx, homo_down_idx)

Calculate the spin multiplicity (2S+1) from HOMO level indices.

Source code in defectpl/io/vasp.py

def get_spin_multiplicity(homo_up_idx: int, homo_down_idx: int) -> float:
    """Calculate the spin multiplicity (2S+1) from HOMO level indices."""
    S = abs(homo_up_idx - homo_down_idx) / 2.0
    return 2.0 * S + 1.0

OUTCAR parsing

OutcarParser

OutcarParser(filename)

Lightweight VASP OUTCAR parser (no pymatgen import at construction time).

Parameters:

Name	Type	Description	Default
filename	str or Path	Path to the VASP OUTCAR file.	required

Source code in defectpl/io/vasp.py

def __init__(self, filename: Union[str, Path]):
    self.filename_path = Path(filename).resolve()
    self.natoms = self.get_natoms()

get_structures_and_forces

get_structures_and_forces(outcar_path, poscar_path=None)

Extract all ionic-step structures and forces from a VASP OUTCAR.

Lattice matrices are updated per ionic step from the OUTCAR itself. Species are read natively from POTCAR entries unless poscar_path is given.

Returns:

Type	Description
(structures, forces)	structures : list of pymatgen.core.Structure forces : list of numpy.ndarray, shape (NIONS, 3), in eV/Å

Source code in defectpl/io/vasp.py

def get_structures_and_forces(
    outcar_path: Union[str, Path],
    poscar_path: Optional[Union[str, Path]] = None,
) -> Tuple[List["Structure"], List[np.ndarray]]:
    """
    Extract all ionic-step structures and forces from a VASP OUTCAR.

    Lattice matrices are updated per ionic step from the OUTCAR itself.
    Species are read natively from POTCAR entries unless *poscar_path* is given.

    Returns
    -------
    (structures, forces)
        structures : list of pymatgen.core.Structure
        forces : list of numpy.ndarray, shape (NIONS, 3), in eV/Å
    """
    from pymatgen.core import Structure

    outcar_path = Path(outcar_path)
    natoms = get_nions(outcar_path)

    if poscar_path:
        from pymatgen.io.vasp import Poscar

        poscar_path = Path(poscar_path)
        if not poscar_path.is_file():
            raise FileNotFoundError(f"POSCAR reference file not found at {poscar_path}")
        species = Poscar.from_file(str(poscar_path)).structure.species
    else:
        species = get_species_and_index_map(outcar_path)

    current_lattice = None
    structures: List[Structure] = []
    forces: List[np.ndarray] = []

    with open(outcar_path, "r", encoding="utf-8", errors="ignore") as f:
        iterator = iter(f)
        for line in iterator:
            try:
                if "VOLUME and BASIS-vectors are now :" in line:
                    for _ in range(4):
                        next(iterator)
                    lattice_matrix = []
                    for _ in range(3):
                        lattice_matrix.append(
                            [float(x) for x in next(iterator).split()[:3]]
                        )
                    current_lattice = np.array(lattice_matrix)

                if "POSITION" in line and "TOTAL-FORCE" in line:
                    if current_lattice is None:
                        raise ValueError(
                            f"Parsed a POSITION block before finding a lattice matrix "
                            f"in {outcar_path}. The file layout might be corrupted."
                        )

                    next(iterator)  # skip dashed separator
                    coords = np.zeros((natoms, 3))
                    step_forces = np.zeros((natoms, 3))

                    for i in range(natoms):
                        data = next(iterator).split()
                        coords[i] = [float(data[0]), float(data[1]), float(data[2])]
                        step_forces[i] = [
                            float(data[3]),
                            float(data[4]),
                            float(data[5]),
                        ]

                    struct = Structure(
                        lattice=current_lattice,
                        species=species,
                        coords=coords,
                        coords_are_cartesian=True,
                    )
                    structures.append(struct)
                    forces.append(step_forces)

            except StopIteration:
                raise ValueError(
                    f"Premature end of file encountered while parsing {outcar_path}"
                )

    return structures, forces

get_final_structure_and_forces_from_outcar

get_final_structure_and_forces_from_outcar(outcar_path, poscar_path=None)

Return only the last structure and forces from an OUTCAR.

Source code in defectpl/io/vasp.py

def get_final_structure_and_forces_from_outcar(
    outcar_path: Union[str, Path],
    poscar_path: Optional[Union[str, Path]] = None,
) -> Tuple["Structure", np.ndarray]:
    """Return only the last structure and forces from an OUTCAR."""
    structures, forces = get_structures_and_forces(outcar_path, poscar_path=poscar_path)
    return structures[-1], forces[-1]

get_first_structure_and_forces_from_outcar

get_first_structure_and_forces_from_outcar(outcar_path, poscar_path=None)

Return only the first structure and forces from an OUTCAR.

Source code in defectpl/io/vasp.py

def get_first_structure_and_forces_from_outcar(
    outcar_path: Union[str, Path],
    poscar_path: Optional[Union[str, Path]] = None,
) -> Tuple["Structure", np.ndarray]:
    """Return only the first structure and forces from an OUTCAR."""
    structures, forces = get_structures_and_forces(outcar_path, poscar_path=poscar_path)
    return structures[0], forces[0]

check_outcar_convergence

check_outcar_convergence(outcar_path)

Check electronic and structural convergence from a VASP OUTCAR file.

Returns:

Type	Description
dict with keys ``structural_converged``, ``electronic_converged``,
``finished_cleanly``.

Source code in defectpl/io/vasp.py

def check_outcar_convergence(outcar_path: Union[str, Path]) -> Dict[str, bool]:
    """
    Check electronic and structural convergence from a VASP OUTCAR file.

    Returns
    -------
    dict with keys ``structural_converged``, ``electronic_converged``,
    ``finished_cleanly``.
    """
    outcar_path = Path(outcar_path)
    if not outcar_path.is_file():
        raise FileNotFoundError(f"OUTCAR file not found at {outcar_path}")

    results = {
        "structural_converged": False,
        "electronic_converged": True,
        "finished_cleanly": False,
    }

    tail_lines: deque = deque(maxlen=100)
    has_content = False

    with open(outcar_path, "r", encoding="utf-8", errors="ignore") as f:
        for line in f:
            has_content = True
            tail_lines.append(line)
            if "ELECTRONIC CONVERGENCE MINIMIZATION" in line and "not achieved" in line:
                results["electronic_converged"] = False

    if not has_content:
        return {k: False for k in results}

    for line in tail_lines:
        if "reached required accuracy" in line:
            results["structural_converged"] = True
        tokens = line.split()
        if len(tokens) >= 2 and ("User" in tokens[0] and "time" in tokens[1]):
            results["finished_cleanly"] = True
        elif "Total CPU time" in line:
            results["finished_cleanly"] = True

    return results

get_nions

get_nions(outcar_path)

Extract the NIONS count from a VASP OUTCAR file.

Source code in defectpl/io/vasp.py

def get_nions(outcar_path: Union[str, Path]) -> int:
    """Extract the NIONS count from a VASP OUTCAR file."""
    outcar_path = Path(outcar_path)
    if not outcar_path.is_file():
        raise FileNotFoundError(f"OUTCAR file not found at {outcar_path}")

    with open(outcar_path, "r", encoding="utf-8", errors="ignore") as f:
        for line in f:
            if "NIONS =" in line:
                return int(line.split("=")[-1])

    raise ValueError(f"Could not find 'NIONS =' token within {outcar_path}")

get_species_and_index_map

get_species_and_index_map(outcar_path)

Build a flat per-atom element list from OUTCAR POTCAR entries.

Handles multi-occurrence POTCAR definitions (e.g. N, C, N, C) by aligning to the 'ions per type' array.

Source code in defectpl/io/vasp.py

def get_species_and_index_map(outcar_path: Union[str, Path]) -> List[str]:
    """
    Build a flat per-atom element list from OUTCAR POTCAR entries.

    Handles multi-occurrence POTCAR definitions (e.g. N, C, N, C) by
    aligning to the 'ions per type' array.
    """
    outcar_path = Path(outcar_path)
    if not outcar_path.is_file():
        raise FileNotFoundError(f"OUTCAR file not found at {outcar_path}")

    species_types: List[str] = []
    ions_per_type: List[int] = []

    with open(outcar_path, "r", encoding="utf-8", errors="ignore") as f:
        for line in f:
            if "POTCAR:" in line and "PAW" in line:
                tokens = line.split()
                if len(tokens) >= 3:
                    species_types.append(tokens[2].split("_")[0])
            if "ions per type =" in line:
                ions_per_type = [int(x) for x in line.split("=")[-1].split()]
                break

    if not species_types or not ions_per_type:
        raise ValueError(
            f"Could not fully parse species maps from {outcar_path}. "
            f"Found species types: {species_types}, Counts: {ions_per_type}"
        )

    if len(species_types) < len(ions_per_type):
        raise ValueError(
            f"OUTCAR processing error: Parsed fewer POTCAR species entries "
            f"({len(species_types)}) than blocks in 'ions per type' ({len(ions_per_type)})."
        )

    species_types = species_types[: len(ions_per_type)]
    species_map: List[str] = []
    for element, count in zip(species_types, ions_per_type):
        species_map.extend([element] * count)

    return species_map

Displacement helpers

calc_dR

calc_dR(contcar_gs, contcar_es)

Compute PBC-safe Cartesian displacement vectors ΔR = R(excited) − R(ground).

Returns:

Type	Description
np.ndarray, shape (natoms, 3) — displacements in Å.

Source code in defectpl/io/vasp.py

def calc_dR(
    contcar_gs: Union["Structure", str, Path],
    contcar_es: Union["Structure", str, Path],
) -> np.ndarray:
    """
    Compute PBC-safe Cartesian displacement vectors ΔR = R(excited) − R(ground).

    Returns
    -------
    np.ndarray, shape (natoms, 3)  — displacements in Å.
    """
    from pymatgen.util.coord import pbc_shortest_vectors

    struct_gs = _to_structure(contcar_gs)
    struct_es = _to_structure(contcar_es)

    if len(struct_gs) != len(struct_es):
        raise ValueError(
            "Lattice structure mismatch: Input configurations have differing site counts."
        )

    n = len(struct_gs)
    lattice = struct_gs.lattice
    dR = np.vstack(
        [
            pbc_shortest_vectors(
                lattice, struct_gs.frac_coords[i], struct_es.frac_coords[i]
            )
            for i in range(n)
        ]
    ).reshape(n, 3)
    return dR

calc_delta_Q

calc_delta_Q(struct1, struct2)

Mass-weighted configuration coordinate displacement ΔQ between two structures.

Returns:

Type	Description
float — ΔQ in amu^(1/2) · Å.

Source code in defectpl/io/vasp.py

def calc_delta_Q(struct1: "Structure", struct2: "Structure") -> float:
    """
    Mass-weighted configuration coordinate displacement ΔQ between two structures.

    Returns
    -------
    float  — ΔQ in amu^(1/2) · Å.
    """
    if len(struct1) != len(struct2):
        raise ValueError("Structures must have the same number of atoms.")

    masses = np.array([site.specie.atomic_mass for site in struct1.sites])
    dR = calc_dR(struct1, struct2)
    return calc_delQ(masses, dR)

get_q_from_structure

get_q_from_structure(ground, excited, struct, tol=0.0001, nround=4)

Project a displaced structure onto the ground→excited configuration coordinate.

Returns:

Type	Description
float — Q in amu^(1/2) · Å.

Source code in defectpl/io/vasp.py

def get_q_from_structure(
    ground: "Structure",
    excited: "Structure",
    struct: Union["Structure", str, Path],
    tol: float = 1e-4,
    nround: int = 4,
) -> float:
    """
    Project a displaced structure onto the ground→excited configuration coordinate.

    Returns
    -------
    float  — Q in amu^(1/2) · Å.
    """
    from pymatgen.util.coord import pbc_shortest_vectors

    if isinstance(struct, (str, Path)):
        from pymatgen.core import Structure

        tstruct = Structure.from_file(str(struct))
    else:
        tstruct = struct

    if len(ground) != len(excited) or len(ground) != len(tstruct):
        raise ValueError(
            "Lattice structure mismatch: Input geometries have differing site counts."
        )

    masses = np.array([site.specie.atomic_mass for site in ground], dtype=float)
    lattice = ground.lattice

    dr_excited_raw = pbc_shortest_vectors(
        lattice, ground.frac_coords, excited.frac_coords
    )
    dr_excited_raw = np.reshape(dr_excited_raw, (len(ground), 3))

    total_dQ = float(np.sqrt(np.sum(masses * np.sum(dr_excited_raw**2, axis=1))))

    dx_struct = pbc_shortest_vectors(lattice, ground.frac_coords, tstruct.frac_coords)
    dx_struct = np.reshape(dx_struct, (len(ground), 3))

    active_mask = np.abs(dr_excited_raw) > tol
    if not np.any(active_mask):
        return 0.0

    ratios = dx_struct[active_mask] / dr_excited_raw[active_mask]
    rounded_ratios = np.round(ratios, decimals=nround)

    values, counts = np.unique(rounded_ratios, return_counts=True)
    scaling_factor = float(values[np.argmax(counts)])
    return scaling_factor * total_dQ

Force-difference helpers

calc_dF

calc_dF(ground_data, excited_data)

Compute the force difference dF = F_excited − F_ground at the same geometry.

Both dicts must have keys "structure" (pymatgen Structure) and "forces" (numpy array, shape (NIONS, 3)).

Source code in defectpl/io/vasp.py

def calc_dF(ground_data: dict, excited_data: dict) -> np.ndarray:
    """
    Compute the force difference dF = F_excited − F_ground at the same geometry.

    Both dicts must have keys ``"structure"`` (pymatgen Structure) and
    ``"forces"`` (numpy array, shape (NIONS, 3)).
    """
    if ground_data["structure"] != excited_data["structure"]:
        raise ValueError(
            "Ground and excited state structures do not match. Cannot calculate dF."
        )
    return excited_data["forces"] - ground_data["forces"]

prepare_dF_files

prepare_dF_files(ground_outcar, excited_outcar, ground_poscar=None, excited_poscar=None)

Extract dF = F_excited − F_ground from two VASP OUTCAR files.

Uses the final structure/forces from each OUTCAR.

Source code in defectpl/io/vasp.py

def prepare_dF_files(
    ground_outcar: Union[str, Path],
    excited_outcar: Union[str, Path],
    ground_poscar: Optional[Union[str, Path]] = None,
    excited_poscar: Optional[Union[str, Path]] = None,
) -> np.ndarray:
    """
    Extract dF = F_excited − F_ground from two VASP OUTCAR files.

    Uses the **final** structure/forces from each OUTCAR.
    """
    ground_structure, ground_forces = get_final_structure_and_forces_from_outcar(
        ground_outcar, poscar_path=ground_poscar
    )
    excited_structure, excited_forces = get_final_structure_and_forces_from_outcar(
        excited_outcar, poscar_path=excited_poscar
    )
    ground_data = {"structure": ground_structure, "forces": ground_forces}
    excited_data = {"structure": excited_structure, "forces": excited_forces}
    return calc_dF(ground_data, excited_data)

VASP Protocol reader

VaspReader

VASP-specific implementation of the :class:~defectpl.io.base.PhononReader and :class:~defectpl.io.base.ElectronicReader protocols.

All DFT-code-specific logic lives here; the physics layer only sees the generic :class:~defectpl.core.structures.PhononData / :class:~defectpl.core.structures.EigenvalData containers.