defectpl.defectpl

Photoluminescence dataclass

Photoluminescence(
    frequencies,
    eigenvectors,
    masses,
    EZPL,
    dR=None,
    dF=None,
    resolution=1000,
    max_energy=5.0,
    sigma=0.006,
    gamma=2.0,
)

Bases: MSONable

Core engine for first-principles photoluminescence lineshape calculations.

Implements the generating-function formalism of Alkauskas et al. (New J. Phys. 16, 073026, 2014) to compute the full multi-phonon PL spectrum of a point defect from DFT inputs.

Two input modes are supported:

• Displacement mode — supply dR (atomic displacements between ground and excited equilibrium geometries).
• Force mode — supply dF (forces on the ground-state atoms when the system is in the excited-state charge state).

Parameters:

Name	Type	Description	Default
frequencies	(ndarray, shape(nmodes))	Phonon mode frequencies at the Γ point in eV.	required
eigenvectors	(ndarray, shape(nmodes, natoms, 3))	Mass-normalised phonon displacement eigenvectors (real part only).	required
masses	(ndarray, shape(natoms))	Atomic masses in atomic mass units (amu).	required
EZPL	float	Zero-phonon line (ZPL) energy in eV.	required
dR	(ndarray, shape(natoms, 3))	Atomic displacement vector (excited − ground equilibrium) in Å. Used in displacement mode.	None
dF	(ndarray, shape(natoms, 3))	Atomic force difference (excited charge − ground charge) in eV/Å. Used in force mode. dR takes priority if both are given.	None
resolution	int	Number of energy grid points per eV; sets spectral resolution. Default 1000.	1000
max_energy	float	Upper bound of the energy axis in eV. Default 5.0.	5.0
sigma	float	Gaussian broadening width applied to the spectral function S(ω) in eV. Default 6 meV.	0.006
gamma	float	Lorentzian (homogeneous) broadening of the ZPL in meV. Default 2.0 meV.	2.0

Attributes:

Name	Type	Description
natoms	int	Number of atoms in the supercell.
delR	float	Root-mean-square atomic displacement ΔR in Å.
delQ	float	Mass-weighted configuration coordinate displacement ΔQ in amu\ :sup:1/2·Å.
qks	(ndarray, shape(nmodes))	Mode-projected configurational displacement q\ :sub:k in amu\ :sup:1/2·Å.
Sks	(ndarray, shape(nmodes))	Partial (mode-resolved) Huang–Rhys factors S\ :sub:k.
HR_factor	float	Total Huang–Rhys factor S = Σ S\ :sub:k.
DW_factor	float	Debye–Waller factor exp(−S).
iprs	(ndarray, shape(nmodes))	Phonon inverse participation ratio (IPR) per mode; range [1/N, 1].
localization_ratio	(ndarray, shape(nmodes))	Localization ratio β\ :sub:k = N · IPR\ :sub:k; range [1, N].
S_omega	ndarray	Continuous Huang–Rhys spectral density S(ω) (eV\ :sup:-1).
Sts	ndarray	Fourier transform of S(ω) used in the generating function.
Gts	ndarray	Generating function G(t) = exp[S(t) − S] · exp(−γ|t|).
A_line	ndarray	Photon energy axis for the lineshape in eV.
intensity	ndarray	Normalised PL intensity as a function of photon energy.

Notes

The PL lineshape is obtained via the Fourier transform of the generating function G(t):

.. math::

L(\hbar\omega) = \int_{-\infty}^{\infty} G(t)\,
e^{i(E_{\mathrm{ZPL}} - \hbar\omega)t/\hbar}\, dt

where

.. math::

G(t) = e^{S(t) - S} \cdot e^{-\gamma|t|}

and :math:S(t) = \int_0^{\infty} S(\hbar\omega)\, e^{-i\omega t}\, d(\hbar\omega).

References

Alkauskas, Buckley, Awschalom & Van de Walle, New J. Phys. 16, 073026 (2014).

Examples:

>>> from defectpl.phonon import read_band_yaml
>>> import numpy as np
>>> freqs, evecs, masses = read_band_yaml("band.yaml")
>>> dR = np.load("dR.npy")          # shape (natoms, 3), in Å
>>> pl = Photoluminescence(
...     frequencies=freqs,
...     eigenvectors=evecs,
...     masses=masses,
...     EZPL=1.945,
...     dR=dR,
... )
>>> print(f"S = {pl.HR_factor:.3f},  DW = {pl.DW_factor:.4f}")

compute_properties

compute_properties()

Run the full calculation pipeline and populate all derived attributes.

Source code in defectpl\defectpl.py

def compute_properties(self):
    """Run the full calculation pipeline and populate all derived attributes."""
    self.delR = utils.calc_delR(self.dR) if self.dR is not None else 0.0
    self.delQ = (
        utils.calc_delQ(self.masses, self.dR) if self.dR is not None else 0.0
    )

    if self.dF is not None and np.any(self.dF):
        self.qks = utils.calc_qks_force_mode(
            self.masses, self.dF, self.eigenvectors, self.frequencies
        )
    elif self.dR is not None and np.any(self.dR):
        self.qks = utils.calc_qks(self.masses, self.dR, self.eigenvectors)
    else:
        raise ValueError(
            "Either dR or dF must be provided and non-zero to compute qks."
        )

    self.Sks = utils.calc_Sks(self.qks, self.frequencies)
    self.HR_factor = float(np.sum(self.Sks))
    self.DW_factor = float(np.exp(-self.HR_factor))

    self.iprs = utils.calc_IPR(self.eigenvectors)
    self.localization_ratio = self.natoms * self.iprs

    self.S_omega = utils.calc_S_omega(
        self.frequencies, self.Sks, self.omega_range, self.sigma
    )
    self.Sts = utils.calc_St(self.S_omega)
    self.Gts = utils.calc_Gts(self.Sts, self.HR_factor, self.gamma, self.resolution)
    self.A_line, self.intensity = utils.calc_Spectrum_Intensity(
        self.Gts, self.EZPL, self.resolution
    )

as_dict

as_dict()

Serialize to a JSON-compatible dictionary.

Complex-valued arrays (Sts, Gts) and derived spectral arrays (A_line, intensity) are intentionally omitted; they are cheaply recomputed by :meth:from_dict.

Source code in defectpl\defectpl.py

def as_dict(self) -> dict:
    """
    Serialize to a JSON-compatible dictionary.

    Complex-valued arrays (``Sts``, ``Gts``) and derived spectral arrays
    (``A_line``, ``intensity``) are intentionally omitted; they are
    cheaply recomputed by :meth:`from_dict`.
    """
    return {
        "@module": self.__class__.__module__,
        "@class": self.__class__.__name__,
        # Core Inputs
        "frequencies": self.frequencies.tolist(),
        "eigenvectors": self.eigenvectors.tolist(),
        "masses": self.masses.tolist(),
        "dR": self.dR.tolist() if self.dR is not None else None,
        "dF": self.dF.tolist() if self.dF is not None else None,
        "EZPL": self.EZPL,
        "gamma": self.gamma,
        "resolution": self.resolution,
        "max_energy": self.max_energy,
        "sigma": self.sigma,
        # Safe Real-Valued Computed Properties
        "natoms": self.natoms,
        "delR": float(self.delR) if hasattr(self.delR, "__float__") else self.delR,
        "delQ": float(self.delQ) if hasattr(self.delQ, "__float__") else self.delQ,
        "qks": self.qks.tolist() if self.qks is not None else None,
        "Sks": self.Sks.tolist() if self.Sks is not None else None,
        "HR_factor": (
            float(self.HR_factor)
            if hasattr(self.HR_factor, "__float__")
            else self.HR_factor
        ),
        "DW_factor": (
            float(self.DW_factor)
            if hasattr(self.DW_factor, "__float__")
            else self.DW_factor
        ),
        "iprs": self.iprs.tolist() if self.iprs is not None else None,
        "localization_ratio": (
            self.localization_ratio.tolist()
            if self.localization_ratio is not None
            else None
        ),
        "omega_range": self.omega_range,
        "S_omega": self.S_omega.tolist() if self.S_omega is not None else None,
        # Explicitly drop complex/spectral properties to avoid serialization bugs
        "Sts": None,
        "Gts": None,
        "A_line": None,
        "intensity": None,
    }

from_dict classmethod

from_dict(d)

Deserialize from a dictionary produced by :meth:as_dict.

Core inputs and real-valued computed properties are loaded directly; complex arrays (Sts, Gts) and the intensity spectrum are recomputed on the fly from the stored S(ω).

Source code in defectpl\defectpl.py

@classmethod
def from_dict(cls, d: dict):
    """
    Deserialize from a dictionary produced by :meth:`as_dict`.

    Core inputs and real-valued computed properties are loaded directly;
    complex arrays (``Sts``, ``Gts``) and the intensity spectrum are
    recomputed on the fly from the stored S(ω).
    """
    # Create an uninitialized instance to handle assignment without tripping standard __post_init__ loops
    obj = cls.__new__(cls)

    # Load Core Inputs
    obj.frequencies = np.array(d["frequencies"])
    obj.eigenvectors = np.array(d["eigenvectors"])
    obj.masses = np.array(d["masses"])
    obj.dR = np.array(d["dR"]) if d.get("dR") is not None else None
    obj.dF = np.array(d["dF"]) if d.get("dF") is not None else None
    obj.EZPL = d["EZPL"]
    obj.gamma = d["gamma"]
    obj.resolution = d.get("resolution", 1000)
    obj.max_energy = d.get("max_energy", 5.0)
    obj.sigma = d.get("sigma", 6e-3)

    # Load Stored Real-Valued Properties
    obj.natoms = d.get("natoms", len(obj.masses))
    obj.delR = d.get("delR")
    obj.delQ = d.get("delQ")
    obj.qks = np.array(d["qks"]) if d.get("qks") is not None else None
    obj.Sks = np.array(d["Sks"]) if d.get("Sks") is not None else None
    obj.HR_factor = d.get("HR_factor")
    obj.DW_factor = d.get("DW_factor")
    obj.iprs = np.array(d["iprs"]) if d.get("iprs") is not None else None
    obj.localization_ratio = (
        np.array(d["localization_ratio"])
        if d.get("localization_ratio") is not None
        else None
    )
    obj.omega_range = d.get(
        "omega_range", [0.0, obj.max_energy, int(obj.max_energy * obj.resolution)]
    )
    obj.S_omega = np.array(d["S_omega"]) if d.get("S_omega") is not None else None

    # Set up placeholders for the missing complex/spectral parameters
    obj.Sts = None
    obj.Gts = None
    obj.A_line = None
    obj.intensity = None

    # Fallback Trigger: Recompute the complex-dependent parts of the pipeline on-the-fly
    if obj.intensity is None:
        # Recompute the remaining downstream properties (Sts, Gts, A_line, intensity)
        obj.Sts = utils.calc_St(obj.S_omega)
        obj.Gts = utils.calc_Gts(obj.Sts, obj.HR_factor, obj.gamma, obj.resolution)
        obj.A_line, obj.intensity = utils.calc_Spectrum_Intensity(
            obj.Gts, obj.EZPL, obj.resolution
        )

    return obj

from_dict_expensive classmethod

from_dict_expensive(d)

Reconstruct by replaying the full pipeline from primary inputs only.

Slower than :meth:from_dict because __post_init__ recomputes everything from scratch; useful when stored arrays may be stale.

Source code in defectpl\defectpl.py

@classmethod
def from_dict_expensive(cls, d: dict):
    """
    Reconstruct by replaying the full pipeline from primary inputs only.

    Slower than :meth:`from_dict` because ``__post_init__`` recomputes
    everything from scratch; useful when stored arrays may be stale.
    """
    return cls(
        frequencies=np.array(d["frequencies"]),
        eigenvectors=np.array(d["eigenvectors"]),
        masses=np.array(d["masses"]),
        dR=np.array(d["dR"]) if d.get("dR") is not None else None,
        dF=np.array(d["dF"]) if d.get("dF") is not None else None,
        EZPL=d["EZPL"],
        gamma=d["gamma"],
        resolution=d.get("resolution", 1000),
        max_energy=d.get("max_energy", 5.0),
        sigma=d.get("sigma", 6e-3),
    )

generate_plots

generate_plots(
    out_dir, max_freq=None, iylim=None, fig_format="pdf"
)

Generate all standard diagnostic plots and save them to out_dir.

Produces ten figures: phonon energy vs mode index, IPR vs energy, localization ratio vs energy, q\ :sub:k vs energy, S\ :sub:k vs energy, S(ω) alone, S(ω)+S\ :sub:k, S(ω)+S\ :sub:k coloured by localization ratio, S(ω)+S\ :sub:k coloured by IPR, and the final PL intensity spectrum.

Parameters:

Name	Type	Description	Default
out_dir	str or Path	Directory where figure files are written (created if absent).	required
max_freq	float	Upper frequency cut-off for S(ω) plots in meV. If None the maximum phonon frequency is used.	None
iylim	tuple of float	(ymin, ymax) for the intensity plot y-axis.	None
fig_format	(pdf, png, svg, jpg)	Output figure format. Default 'pdf'.	'pdf'

Source code in defectpl\defectpl.py

def generate_plots(
    self,
    out_dir: Union[str, Path],
    max_freq: Optional[float] = None,
    iylim=None,
    fig_format="pdf",
):
    r"""
    Generate all standard diagnostic plots and save them to *out_dir*.

    Produces ten figures: phonon energy vs mode index, IPR vs energy,
    localization ratio vs energy, q\ :sub:`k` vs energy,
    S\ :sub:`k` vs energy, S(ω) alone, S(ω)+S\ :sub:`k`, S(ω)+S\ :sub:`k`
    coloured by localization ratio, S(ω)+S\ :sub:`k` coloured by IPR,
    and the final PL intensity spectrum.

    Parameters
    ----------
    out_dir : str or Path
        Directory where figure files are written (created if absent).
    max_freq : float, optional
        Upper frequency cut-off for S(ω) plots in **meV**.  If ``None``
        the maximum phonon frequency is used.
    iylim : tuple of float, optional
        ``(ymin, ymax)`` for the intensity plot y-axis.
    fig_format : {'pdf', 'png', 'svg', 'jpg'}, optional
        Output figure format.  Default ``'pdf'``.
    """
    plotter = Plotter()
    iplot_xlim = (max(0.0, self.EZPL - 2.0), self.EZPL + 1.0)
    freq_limit = (max_freq / 1000.0) if max_freq else None

    plotter.plot_penergy_vs_pmode(
        frequencies=self.frequencies,
        plot=False,
        out_dir=out_dir,
        fig_format=fig_format,
    )
    plotter.plot_ipr_vs_penergy(
        self.frequencies,
        self.iprs,
        plot=False,
        out_dir=out_dir,
        fig_format=fig_format,
    )
    plotter.plot_loc_rat_vs_penergy(
        self.frequencies,
        self.localization_ratio,
        plot=False,
        out_dir=out_dir,
        fig_format=fig_format,
    )
    plotter.plot_qk_vs_penergy(
        self.frequencies,
        self.qks,
        plot=False,
        out_dir=out_dir,
        fig_format=fig_format,
    )
    plotter.plot_HR_factor_vs_penergy(
        self.frequencies,
        self.Sks,
        plot=False,
        out_dir=out_dir,
        fig_format=fig_format,
    )

    plotter.plot_S_omega_vs_penergy(
        self.frequencies,
        self.S_omega,
        self.omega_range,
        plot=False,
        out_dir=out_dir,
        max_freq=freq_limit,
        fig_format=fig_format,
    )
    plotter.plot_S_omega_Sks_vs_penergy(
        self.frequencies,
        self.S_omega,
        self.omega_range,
        self.Sks,
        plot=False,
        out_dir=out_dir,
        max_freq=freq_limit,
        fig_format=fig_format,
    )
    plotter.plot_S_omega_Sks_Loc_rat_vs_penergy(
        self.frequencies,
        self.S_omega,
        self.omega_range,
        self.Sks,
        self.localization_ratio,
        plot=False,
        out_dir=out_dir,
        max_freq=freq_limit,
        fig_format=fig_format,
    )
    plotter.plot_S_omega_Sks_ipr_vs_penergy(
        self.frequencies,
        self.S_omega,
        self.omega_range,
        self.Sks,
        self.iprs,
        plot=False,
        out_dir=out_dir,
        max_freq=freq_limit,
        fig_format=fig_format,
    )

    plotter.plot_intensity_vs_penergy(
        frequencies=self.frequencies,
        I=self.intensity,
        resolution=self.resolution,
        xlim=iplot_xlim,
        plot=False,
        out_dir=out_dir,
        iylim=iylim,
        fig_format=fig_format,
    )
    print("All static visualization plots generated successfully.")

VibrationalSpectra1D dataclass

VibrationalSpectra1D(
    EZPL,
    w1_meV,
    w2_meV,
    DQ,
    T,
    E0,
    dE,
    M,
    NN1=22,
    NN2=52,
    K2EV=8.617333262e-05,
)

Bases: MSONable

1D harmonic-oscillator model for the vibrational lineshape of a luminescence band.

Computes Franck–Condon overlaps between vibrational levels of two displaced harmonic potentials (ground and excited state) with potentially different effective frequencies, following Alkauskas, Yan & Van de Walle (Phys. Rev. B 90, 075202, 2014).

This model is complementary to :class:Photoluminescence — it replaces the multi-phonon generating function with an explicit sum over vibrational quantum numbers, which is exact within the harmonic approximation and more suitable when ω\ :sub:1 ≠ ω\ :sub:2.

Parameters:

Name	Type	Description	Default
EZPL	float	Zero-phonon line energy in eV.	required
w1_meV	float	Ground-state effective phonon frequency in meV.	required
w2_meV	float	Excited-state effective phonon frequency in meV.	required
DQ	float	Configuration coordinate displacement ΔQ in amu\ :sup:1/2·Å.	required
T	float	Temperature in K (sets Boltzmann weights for ground-state levels).	required
E0	float	Starting photon energy for the output lineshape grid in eV.	required
dE	float	Energy step of the output lineshape grid in eV.	required
M	int	Number of energy grid points.	required
NN1	int	Maximum vibrational quantum number included for the ground state. Default 22.	22
NN2	int	Maximum vibrational quantum number included for the excited state. Default 52.	52

Attributes:

Name	Type	Description
overlap_matrix	(ndarray, shape(NN1 + 1, NN2 + 1))	Franck–Condon overlap integrals ⟨i|j⟩.
overlap_data	dict	Keys "contributions" and "energies" — per-transition Boltzmann-weighted FC squared integrals and their photon energies.
spectral_data	dict	Keys "energies", "dos", "dosw3" — the energy grid, Gaussian-broadened transition density, and the ω³-weighted (intensity-corrected) lineshape.

Notes

The overlap integral between vibrational levels :math:i (ground) and :math:j (excited) is evaluated recursively using the displaced harmonic oscillator recurrence relation. The full lineshape is then obtained by Gaussian broadening with width σ = 0.70·ω\ :sub:2.

The ω³-weighted density dosw3 approximates the spontaneous emission rate and is normalized to unit area.

Examples:

>>> spec = VibrationalSpectra1D(
...     EZPL=1.945, w1_meV=65.0, w2_meV=62.0,
...     DQ=1.8, T=300.0, E0=1.3, dE=0.001, M=700
... )
>>> spec.compute_lineshape()
>>> e_peak, _ = spec.get_peak_position()
>>> print(f"Peak at {e_peak:.3f} eV,  FWHM = {spec.get_fwhm():.3f} eV")

compute_overlap_matrix

compute_overlap_matrix()

Fill overlap_matrix with Franck–Condon overlaps ⟨i|j⟩ for all (i, j) pairs.

Source code in defectpl\defectpl.py

def compute_overlap_matrix(self) -> None:
    """Fill ``overlap_matrix`` with Franck–Condon overlaps ⟨i|j⟩ for all (i, j) pairs."""
    for i in range(self.NN1 + 1):
        for j in range(self.NN2 + 1):
            self.overlap_matrix[i, j] = utils.calculate_overlap_element(
                i, j, self.rho, self.cosfi, self.sinfi
            )

compute_spectrum

compute_spectrum()

Compute Boltzmann-weighted Franck–Condon transition intensities.

Populates overlap_data with per-transition contributions and energies, and prints the closure-relation sum (should be ≈ 1).

Source code in defectpl\defectpl.py

def compute_spectrum(self) -> None:
    """
    Compute Boltzmann-weighted Franck–Condon transition intensities.

    Populates ``overlap_data`` with per-transition contributions and
    energies, and prints the closure-relation sum (should be ≈ 1).
    """
    self.compute_overlap_matrix()

    # Guard against zero temperature limits
    if self.TE <= 0:
        Z = 1.0
        weights = np.zeros(self.NN1 + 1)
        weights[0] = 1.0
    else:
        Z = 1.0 / (1.0 - np.exp(-self.w1 / self.TE))
        weights = np.exp(-np.arange(self.NN1 + 1) * self.w1 / self.TE) / Z

    contr, en = [], []
    for i in range(self.NN1 + 1):
        weight = weights[i]
        for j in range(self.NN2 + 1):
            val = self.overlap_matrix[i, j]
            contrib = weight * (val**2)
            contr.append(contrib)
            en.append(self.EZPL - j * self.w2 + i * self.w1)

    print(f"Closure relation sum (should be ~1.0): {sum(contr):.6f}")
    self.overlap_data = {"contributions": contr, "energies": en}

compute_lineshape

compute_lineshape()

Gaussian-broaden the FC transitions to produce the PL lineshape.

Calls :meth:compute_spectrum if not already done, then builds the energy-resolved transition density dos and the ω³-weighted normalised intensity dosw3 on the grid [E0, E0 + M·dE]. Results are stored in spectral_data.

Source code in defectpl\defectpl.py

def compute_lineshape(self) -> None:
    """
    Gaussian-broaden the FC transitions to produce the PL lineshape.

    Calls :meth:`compute_spectrum` if not already done, then builds the
    energy-resolved transition density ``dos`` and the ω³-weighted
    normalised intensity ``dosw3`` on the grid [E0, E0 + M·dE].
    Results are stored in ``spectral_data``.
    """
    if not self.overlap_data:
        self.compute_spectrum()

    energies_grid = self.E0 + np.arange(self.M) * self.dE
    contr = np.array(self.overlap_data["contributions"])
    en = np.array(self.overlap_data["energies"])

    # Vectorized calculation pipeline across grid positions to maximize performance
    delta_E = en[:, np.newaxis] - energies_grid[np.newaxis, :]
    gaussian = np.exp(-(delta_E**2) / (2 * self.sigma**2)) / (
        self.sigma * np.sqrt(2 * np.pi)
    )
    dos = np.dot(contr, gaussian)

    dosw3 = dos * (energies_grid**3)
    norm_factor = np.sum(dosw3) * self.dE
    if norm_factor > 0:
        dosw3 /= norm_factor

    self.spectral_data = {
        "energies": energies_grid.tolist(),
        "dos": dos.tolist(),
        "dosw3": dosw3.tolist(),
    }

save_results

save_results(
    overlap_file="overlap.json",
    lineshape_file="lineshape.json",
)

Write results to JSON files.

Parameters:

Name	Type	Description	Default
overlap_file	str	Destination for the MSONable object JSON (includes all inputs). Default "overlap.json".	'overlap.json'
lineshape_file	str	Destination for the spectral_data dict. Default "lineshape.json".	'lineshape.json'

Source code in defectpl\defectpl.py

def save_results(
    self, overlap_file: str = "overlap.json", lineshape_file: str = "lineshape.json"
) -> None:
    """
    Write results to JSON files.

    Parameters
    ----------
    overlap_file : str, optional
        Destination for the MSONable object JSON (includes all inputs).
        Default ``"overlap.json"``.
    lineshape_file : str, optional
        Destination for the ``spectral_data`` dict.
        Default ``"lineshape.json"``.
    """
    Path(overlap_file).write_text(self.to_json(), encoding="utf-8")
    with open(lineshape_file, "w", encoding="utf-8") as f:
        json.dump(self.spectral_data, f, indent=4)

plot_lineshape

plot_lineshape(save_file=None, figsize=(4.4, 4.4))

Plot the normalised PL lineshape (ω³-weighted intensity vs energy).

Parameters:

Name	Type	Description	Default
save_file	str	If given, save to this path (dpi 600); otherwise plt.show().	None
figsize	tuple of float	Matplotlib figure size in inches. Default (4.4, 4.4).	(4.4, 4.4)

Raises:

Type	Description
ValueError	If :meth:compute_lineshape has not been called yet.

Source code in defectpl\defectpl.py

def plot_lineshape(
    self, save_file: Optional[str] = None, figsize: Tuple[float, float] = (4.4, 4.4)
) -> None:
    """
    Plot the normalised PL lineshape (ω³-weighted intensity vs energy).

    Parameters
    ----------
    save_file : str, optional
        If given, save to this path (dpi 600); otherwise ``plt.show()``.
    figsize : tuple of float, optional
        Matplotlib figure size in inches.  Default ``(4.4, 4.4)``.

    Raises
    ------
    ValueError
        If :meth:`compute_lineshape` has not been called yet.
    """
    if not self.spectral_data:
        raise ValueError("Run compute_lineshape() before plotting.")

    fig, ax = plt.subplots(figsize=figsize)
    ax.plot(self.spectral_data["energies"], self.spectral_data["dosw3"])
    ax.set_xlabel("Energy (eV)")
    ax.set_ylabel("Intensity (arb. u.)")
    ax.set_yticks([])

    if save_file:
        plt.savefig(save_file, dpi=600, bbox_inches="tight")
        plt.close(fig)
        print(f"Lineshape plot saved to {save_file}")
    else:
        plt.show()
        plt.close(fig)

get_peak_position

get_peak_position()

Return the energy and intensity at the lineshape peak.

Returns:

Type	Description
(energy, intensity) : tuple of float	Peak photon energy in eV and the corresponding normalised intensity.

Raises:

Type	Description
ValueError	If :meth:compute_lineshape has not been called yet.

Source code in defectpl\defectpl.py

def get_peak_position(self) -> Tuple[float, float]:
    """
    Return the energy and intensity at the lineshape peak.

    Returns
    -------
    (energy, intensity) : tuple of float
        Peak photon energy in eV and the corresponding normalised intensity.

    Raises
    ------
    ValueError
        If :meth:`compute_lineshape` has not been called yet.
    """
    if not self.spectral_data:
        raise ValueError("Run compute_lineshape() before accessing metrics.")
    dosw3 = np.array(self.spectral_data["dosw3"])
    energies = np.array(self.spectral_data["energies"])
    idx_max = np.argmax(dosw3)
    print(f"Peak position: {energies[idx_max]:.3f} eV at {self.T} K.")
    return float(energies[idx_max]), float(dosw3[idx_max])

get_fwhm

get_fwhm()

Compute the Full Width at Half Maximum (FWHM) of the PL lineshape.

Returns:

Type	Description
float	FWHM in eV. Returns 0.0 if fewer than two points exceed half-maximum.

Raises:

Type	Description
ValueError	If :meth:compute_lineshape has not been called yet.

Source code in defectpl\defectpl.py

def get_fwhm(self) -> float:
    """
    Compute the Full Width at Half Maximum (FWHM) of the PL lineshape.

    Returns
    -------
    float
        FWHM in eV.  Returns 0.0 if fewer than two points exceed half-maximum.

    Raises
    ------
    ValueError
        If :meth:`compute_lineshape` has not been called yet.
    """
    if not self.spectral_data:
        raise ValueError("Run compute_lineshape() before accessing metrics.")
    dosw3 = np.array(self.spectral_data["dosw3"])
    energies = np.array(self.spectral_data["energies"])

    half_max = np.max(dosw3) / 2.0
    indices = np.where(dosw3 >= half_max)[0]

    if len(indices) < 2:
        return 0.0

    fwhm = energies[indices[-1]] - energies[indices[0]]
    print(f"FWHM: {fwhm:.3f} eV at {self.T} K.")
    return float(fwhm)

ConfigurationCoordinateDiagram dataclass

ConfigurationCoordinateDiagram(
    ground_struct, excited_struct
)

Bases: MSONable

Configuration coordinate diagram for a two-state defect transition.

Manages interpolation of structures between the ground and excited equilibrium geometries, DFT input-file generation, potential energy surface (PES) extraction from completed Vasprun files, and harmonic fitting of the resulting parabolic wells.

Parameters:

Name	Type	Description	Default
ground_struct	Structure	DFT-relaxed ground-state supercell geometry.	required
excited_struct	Structure	DFT-relaxed excited-state supercell geometry (same species ordering).	required

Attributes:

Name	Type	Description
dQ	float	Mass-weighted configuration coordinate displacement ΔQ in amu\ :sup:1/2·Å, computed automatically from the two structures.

Notes

The configuration coordinate is defined as

.. math::

\Delta Q = \sqrt{\sum_i m_i\,|\Delta\mathbf{R}_i|^2}

where the sum runs over all atoms and :math:\Delta\mathbf{R}_i is the displacement of atom i between the two equilibrium geometries (minimum-image convention with periodic boundary conditions).

The zero-phonon line energy and Stokes shift can be read from the fitted parabola minima returned by :meth:analyze_ccd.

Examples:

>>> from pymatgen.core import Structure
>>> gs = Structure.from_file("CONTCAR_gs")
>>> es = Structure.from_file("CONTCAR_es")
>>> ccd = ConfigurationCoordinateDiagram(gs, es)
>>> print(f"ΔQ = {ccd.dQ:.3f} amu^1/2 Å")
>>> ccd.setup_calculations(
...     displacements=[-0.5, -0.25, 0.25, 0.5, 0.75],
...     output_dir="ccd_calcs",
...     ground_input_dir="inputs/ground",
...     excited_input_dir="inputs/excited",
... )

generate_structures

generate_structures(displacements, remove_zero=True)

Linearly interpolate structures at fractional displacements along ΔQ.

Parameters:

Name	Type	Description	Default
displacements	array-like of float	Fractional displacements along the ground→excited path. A value of 0.0 is the ground minimum; 1.0 is the excited minimum. Values outside [0, 1] extrapolate.	required
remove_zero	bool	Drop any entries equal to 0.0 (ground minimum already exists). Default True.	True

Returns:

Name	Type	Description
ground_structs	list of Structure	Structures displaced from ground toward excited (for ground-state PES).
excited_structs	list of Structure	Structures displaced from ground + 1.0 (for excited-state PES).

Source code in defectpl\defectpl.py

def generate_structures(
    self, displacements: Union[List[float], np.ndarray], remove_zero: bool = True
) -> Tuple[List[Structure], List[Structure]]:
    """
    Linearly interpolate structures at fractional displacements along ΔQ.

    Parameters
    ----------
    displacements : array-like of float
        Fractional displacements along the ground→excited path.
        A value of 0.0 is the ground minimum; 1.0 is the excited minimum.
        Values outside [0, 1] extrapolate.
    remove_zero : bool, optional
        Drop any entries equal to 0.0 (ground minimum already exists).
        Default True.

    Returns
    -------
    ground_structs : list of Structure
        Structures displaced from ground toward excited (for ground-state PES).
    excited_structs : list of Structure
        Structures displaced from ground + 1.0 (for excited-state PES).
    """
    disp_arr = np.atleast_1d(displacements)
    if remove_zero:
        disp_arr = disp_arr[disp_arr != 0.0]

    ground_structs = self.ground_struct.interpolate(
        self.excited_struct, nimages=disp_arr
    )
    excited_structs = self.ground_struct.interpolate(
        self.excited_struct, nimages=(disp_arr + 1.0)
    )
    return ground_structs, excited_structs

setup_calculations

setup_calculations(
    displacements,
    output_dir,
    ground_input_dir,
    excited_input_dir,
    input_files=None,
)

Write interpolated POSCAR files and copy VASP input files.

Creates the directory tree::

output_dir/
  ground/0/POSCAR  KPOINTS  POTCAR  INCAR
  ground/1/...
  excited/0/...
  ...

Parameters:

Name	Type	Description	Default
displacements	array-like of float	Fractional displacements (see :meth:generate_structures).	required
output_dir	str or Path	Root directory for the calculation tree (created if absent).	required
ground_input_dir	str or Path	Directory containing the VASP input files for ground-state runs.	required
excited_input_dir	str or Path	Directory containing the VASP input files for excited-state runs.	required
input_files	list of str	File names to copy from each input directory. Default ["KPOINTS", "POTCAR", "INCAR"].	None

Source code in defectpl\defectpl.py

def setup_calculations(
    self,
    displacements: Union[List[float], np.ndarray],
    output_dir: Union[str, Path],
    ground_input_dir: Union[str, Path],
    excited_input_dir: Union[str, Path],
    input_files: Optional[List[str]] = None,
) -> None:
    """
    Write interpolated POSCAR files and copy VASP input files.

    Creates the directory tree::

        output_dir/
          ground/0/POSCAR  KPOINTS  POTCAR  INCAR
          ground/1/...
          excited/0/...
          ...

    Parameters
    ----------
    displacements : array-like of float
        Fractional displacements (see :meth:`generate_structures`).
    output_dir : str or Path
        Root directory for the calculation tree (created if absent).
    ground_input_dir : str or Path
        Directory containing the VASP input files for ground-state runs.
    excited_input_dir : str or Path
        Directory containing the VASP input files for excited-state runs.
    input_files : list of str, optional
        File names to copy from each input directory.
        Default ``["KPOINTS", "POTCAR", "INCAR"]``.
    """
    if input_files is None:
        input_files = ["KPOINTS", "POTCAR", "INCAR"]

    out_path = Path(output_dir)
    g_in_path = Path(ground_input_dir)
    e_in_path = Path(excited_input_dir)

    g_structs, e_structs = self.generate_structures(displacements)

    for idx, struct in enumerate(g_structs):
        target_dir = out_path / "ground" / str(idx)
        target_dir.mkdir(parents=True, exist_ok=True)
        struct.to(filename=str(target_dir / "POSCAR"), fmt="poscar")
        for filename in input_files:
            copyfile(g_in_path / filename, target_dir / filename)

    for idx, struct in enumerate(e_structs):
        target_dir = out_path / "excited" / str(idx)
        target_dir.mkdir(parents=True, exist_ok=True)
        struct.to(filename=str(target_dir / "POSCAR"), fmt="poscar")
        for filename in input_files:
            copyfile(e_in_path / filename, target_dir / filename)

extract_pes_profile

extract_pes_profile(vasprun_paths, tol=0.001)

Extract (Q, E) data points from a list of completed Vasprun files.

Parameters:

Name	Type	Description	Default
vasprun_paths	list of str or Path	Paths to vasprun.xml files, one per interpolated image.	required
tol	float	Fractional-coordinate tolerance for projecting each structure onto the configuration coordinate axis. Default 0.001.	0.001

Returns:

Name	Type	Description
q_values	ndarray	Configuration coordinate values Q in amu\ :sup:1/2·Å.
energies	ndarray	Total DFT energies shifted so the minimum is zero (eV).

Source code in defectpl\defectpl.py

def extract_pes_profile(
    self, vasprun_paths: List[Union[str, Path]], tol: float = 0.001
) -> Tuple[np.ndarray, np.ndarray]:
    r"""
    Extract (Q, E) data points from a list of completed Vasprun files.

    Parameters
    ----------
    vasprun_paths : list of str or Path
        Paths to ``vasprun.xml`` files, one per interpolated image.
    tol : float, optional
        Fractional-coordinate tolerance for projecting each structure
        onto the configuration coordinate axis.  Default 0.001.

    Returns
    -------
    q_values : numpy.ndarray
        Configuration coordinate values Q in amu\ :sup:`1/2`·Å.
    energies : numpy.ndarray
        Total DFT energies shifted so the minimum is zero (eV).
    """
    from pymatgen.io.vasp.outputs import Vasprun

    total_runs = len(vasprun_paths)
    q_values = np.zeros(total_runs)
    energies = np.zeros(total_runs)

    for idx, path in enumerate(vasprun_paths):
        vr = Vasprun(str(path), parse_dos=False, parse_eigen=False)
        q_values[idx] = get_q_from_structure(
            self.ground_struct, self.excited_struct, vr.structures[-1], tol=tol
        )
        energies[idx] = vr.final_energy

    return q_values, (energies - np.min(energies))

analyze_ccd

analyze_ccd(
    ground_vaspruns,
    excited_vaspruns,
    dE=0.0,
    plot=True,
    figsize=(3.3, 3.3),
    xlim=(-3.0, 10.0),
    ylim=(-0.5, 4.0),
    save_plot=None,
)

Fit harmonic wells to the ground- and excited-state PES data.

Parameters:

Name	Type	Description	Default
ground_vaspruns	list of str or Path	vasprun.xml paths for the ground-state PES images.	required
excited_vaspruns	list of str or Path	vasprun.xml paths for the excited-state PES images.	required
dE	float	Rigid energy offset applied to the excited-state PES in eV (e.g. to set the excited minimum to the ZPL energy). Default 0.0.	0.0
plot	bool	If True, show a matplotlib figure. Default True.	True
figsize	tuple of float	Figure size in inches. Default (3.3, 3.3).	(3.3, 3.3)
xlim	tuple of float	x-axis limits for the CCD plot in amu\ :sup:1/2·Å.	(-3.0, 10.0)
ylim	tuple of float	y-axis limits for the CCD plot in eV.	(-0.5, 4.0)
save_plot	str	If given, save the figure to this path instead of showing.	None

Returns:

Type	Description
(omega_ground, omega_excited) : tuple of float	Fitted effective angular frequencies in eV for the ground and excited harmonic potentials.

Source code in defectpl\defectpl.py

def analyze_ccd(
    self,
    ground_vaspruns: List[Union[str, Path]],
    excited_vaspruns: List[Union[str, Path]],
    dE: float = 0.0,
    plot: bool = True,
    figsize: Tuple[float, float] = (3.3, 3.3),
    xlim: Tuple[float, float] = (-3.0, 10.0),
    ylim: Tuple[float, float] = (-0.5, 4.0),
    save_plot: Optional[str] = None,
) -> Tuple[float, float]:
    r"""
    Fit harmonic wells to the ground- and excited-state PES data.

    Parameters
    ----------
    ground_vaspruns : list of str or Path
        ``vasprun.xml`` paths for the ground-state PES images.
    excited_vaspruns : list of str or Path
        ``vasprun.xml`` paths for the excited-state PES images.
    dE : float, optional
        Rigid energy offset applied to the excited-state PES in eV
        (e.g. to set the excited minimum to the ZPL energy).
        Default 0.0.
    plot : bool, optional
        If True, show a matplotlib figure.  Default True.
    figsize : tuple of float, optional
        Figure size in inches.  Default ``(3.3, 3.3)``.
    xlim : tuple of float, optional
        x-axis limits for the CCD plot in amu\ :sup:`1/2`·Å.
    ylim : tuple of float, optional
        y-axis limits for the CCD plot in eV.
    save_plot : str, optional
        If given, save the figure to this path instead of showing.

    Returns
    -------
    (omega_ground, omega_excited) : tuple of float
        Fitted effective angular frequencies in eV for the ground and
        excited harmonic potentials.
    """
    q_ground, e_ground = self.extract_pes_profile(ground_vaspruns)
    q_excited, e_excited = self.extract_pes_profile(excited_vaspruns)
    e_excited += dE

    if plot:
        fig, ax = plt.subplots(figsize=figsize)
        ax.scatter(q_ground, e_ground, label="Ground State")
        ax.scatter(q_excited, e_excited, label="Excited State")

        grid_line = np.linspace(xlim[0], xlim[1], 100)
        ground_omega = utils.get_omega_from_pes(
            q_ground, e_ground, ax=ax, eval_grid=grid_line
        )
        excited_omega = utils.get_omega_from_pes(
            q_excited, e_excited, ax=ax, eval_grid=grid_line
        )

        ax.set_xlabel(r"$Q$ ($\mathrm{amu}^{1/2}\cdot\mathrm{\AA}$)")
        ax.set_ylabel("Energy (eV)")
        ax.set_xlim(xlim)
        ax.set_ylim(ylim)
        ax.xaxis.set_minor_locator(MultipleLocator(1))
        ax.yaxis.set_minor_locator(MultipleLocator(1))
        ax.legend()

        if save_plot:
            plt.savefig(save_plot, bbox_inches="tight", dpi=600)
            plt.close(fig)
        else:
            plt.show()
            plt.close(fig)
    else:
        ground_omega = utils.get_omega_from_pes(q_ground, e_ground)
        excited_omega = utils.get_omega_from_pes(q_excited, e_excited)

    return ground_omega, excited_omega

estimate_vertical_transitions

estimate_vertical_transitions(
    ground_omega, excited_omega, dE, eps=1e-06
)

Calculate vertical absorption and emission energy transitions.

Source code in defectpl\defectpl.py

def estimate_vertical_transitions(
    self, ground_omega: float, excited_omega: float, dE: float, eps: float = 1e-6
) -> Tuple[float, float, float, float]:
    """Calculate vertical absorption and emission energy transitions."""
    if ground_omega < eps or excited_omega < eps:
        raise ValueError("Phonon frequencies must be strictly positive.")

    conversion_factor = (self.dQ**2) * AMU2KG * (ANG2M**2) / (HBAR_EVS**2) / EV2J

    fc_e = 0.5 * (excited_omega**2) * conversion_factor
    fc_g = 0.5 * (ground_omega**2) * conversion_factor

    e_abs = fc_e + dE
    e_em = dE - fc_g

    print(f"Absorption Energy: {e_abs:.6f} eV")
    print(f"Emission Energy:   {e_em:.6f} eV")
    print(f"Franck-Condon shift (excited): {fc_e:.6f} eV")
    print(f"Franck-Condon shift (ground):  {fc_g:.6f} eV")

    return e_abs, e_em, fc_e, fc_g