defectpl.utils

Displacement and configuration coordinates

calc_delR

calc_delR(dR)

Compute the unweighted Cartesian displacement norm ΔR.

Parameters:

Name	Type	Description	Default
dR	ndarray	Atomic displacement matrix, shape (natoms, 3), in Å.	required

Returns:

Type	Description
float	:math:\Delta R = \sqrt{\sum_{a,i} \Delta R_{a,i}^2} in Å.

Examples:

>>> import numpy as np
>>> calc_delR(np.array([[1.0, 2.0, -2.0], [0.0, 0.0, 0.0]]))
3.0

Source code in defectpl/utils.py

def calc_delR(dR: np.ndarray) -> float:
    """
    Compute the unweighted Cartesian displacement norm ΔR.

    Parameters
    ----------
    dR : np.ndarray
        Atomic displacement matrix, shape ``(natoms, 3)``, in Å.

    Returns
    -------
    float
        :math:`\\Delta R = \\sqrt{\\sum_{a,i} \\Delta R_{a,i}^2}` in Å.

    Examples
    --------
    >>> import numpy as np
    >>> calc_delR(np.array([[1.0, 2.0, -2.0], [0.0, 0.0, 0.0]]))
    3.0
    """
    return float(np.sqrt(np.sum(dR**2)))

calc_delQ

calc_delQ(masses, dR)

Compute the mass-weighted configuration coordinate difference ΔQ.

Parameters:

Name	Type	Description	Default
masses	ndarray	Atomic masses in AMU, shape (natoms,).	required
dR	ndarray	Atomic displacement matrix in Å, shape (natoms, 3).	required

Returns:

Type	Description
float	:math:\Delta Q = \sqrt{\sum_a m_a |\Delta\mathbf{R}_a|^2} in :math:\sqrt{\text{amu}} \cdot \text{Å}.

Notes

The closure relation :math:\Delta Q^2 = \sum_k q_k^2 holds when the phonon eigenvectors span the full displacement space.

Examples:

>>> import numpy as np
>>> masses = np.array([12.011, 15.999])
>>> dR = np.array([[0.1, 0.0, 0.0], [0.0, 0.2, 0.0]])
>>> calc_delQ(masses, dR)  # sqrt(12.011*0.01 + 15.999*0.04)
0.8289...

Source code in defectpl/utils.py

def calc_delQ(masses: np.ndarray, dR: np.ndarray) -> float:
    """
    Compute the mass-weighted configuration coordinate difference ΔQ.

    Parameters
    ----------
    masses : np.ndarray
        Atomic masses in AMU, shape ``(natoms,)``.
    dR : np.ndarray
        Atomic displacement matrix in Å, shape ``(natoms, 3)``.

    Returns
    -------
    float
        :math:`\\Delta Q = \\sqrt{\\sum_a m_a |\\Delta\\mathbf{R}_a|^2}` in
        :math:`\\sqrt{\\text{amu}} \\cdot \\text{Å}`.

    Notes
    -----
    The closure relation :math:`\\Delta Q^2 = \\sum_k q_k^2` holds when the
    phonon eigenvectors span the full displacement space.

    Examples
    --------
    >>> import numpy as np
    >>> masses = np.array([12.011, 15.999])
    >>> dR = np.array([[0.1, 0.0, 0.0], [0.0, 0.2, 0.0]])
    >>> calc_delQ(masses, dR)  # sqrt(12.011*0.01 + 15.999*0.04)
    0.8289...
    """
    return float(np.sqrt(np.sum(masses * np.sum(dR**2, axis=1))))

calc_delta_Q

calc_delta_Q(struct1, struct2)

Calculate the mass-weighted configuration coordinate displacement delta Q between two structures.

Parameters:

Name	Type	Description	Default
struct1	Structure	Initial reference equilibrium configuration.	required
struct2	Structure	Final reference equilibrium configuration.	required

Returns:

Type	Description
float	Mass-weighted structural displacement delta Q in amu^(½) * Å.

Raises:

Type	Description
ValueError	If structures have unequal numbers of atoms.

Source code in defectpl/utils.py

def calc_delta_Q(struct1, struct2) -> float:
    """
    Calculate the mass-weighted configuration coordinate displacement delta Q between two structures.

    Parameters
    ----------
    struct1 : pymatgen.core.Structure
        Initial reference equilibrium configuration.
    struct2 : pymatgen.core.Structure
        Final reference equilibrium configuration.

    Returns
    -------
    float
        Mass-weighted structural displacement delta Q in amu^(1/2) * Å.

    Raises
    ------
    ValueError
        If structures have unequal numbers of atoms.
    """
    if len(struct1) != len(struct2):
        raise ValueError("Structures must have the same number of atoms.")

    from pymatgen.util.coord import pbc_shortest_vectors

    masses = np.array([site.specie.atomic_mass for site in struct1.sites])
    lattice = struct1.lattice
    n = len(struct1)
    dR = np.vstack(
        [
            pbc_shortest_vectors(
                lattice, struct1.frac_coords[i], struct2.frac_coords[i]
            )
            for i in range(n)
        ]
    ).reshape(n, 3)
    return calc_delQ(masses, dR)

calc_qks

calc_qks(masses, dR, eigenvectors)

Project atomic displacements onto phonon normal modes to find mode coordinates q_k.

This function uses a loop-based approach to calculate the dimensioned configuration coordinates from real-space displacement vectors.

Parameters:

Name	Type	Description	Default
masses	ndarray	1D array of atomic masses for each atom in the system. Shape: (N_atoms,). Unit: AMU.	required
dR	ndarray	2D array of atomic displacement vectors between two configurations. Shape: (N_atoms, 3). Unit: Angstrom (Å).	required
eigenvectors	ndarray	3D array representing the normal mode eigenvectors matrix. Shape: (N_modes, N_atoms, 3). Dimensionless (normalized).	required

Returns:

Type	Description
ndarray	1D array of projected configuration coordinates for each mode k. Shape: (N_modes,). Unit: SI units (kg^{½} * m).

See Also

calc_qks_vectorized : Faster, loopless alternative for this calculation.

Source code in defectpl/utils.py

def calc_qks(
    masses: np.ndarray, dR: np.ndarray, eigenvectors: np.ndarray
) -> np.ndarray:
    """
    Project atomic displacements onto phonon normal modes to find mode coordinates q_k.

    This function uses a loop-based approach to calculate the dimensioned
    configuration coordinates from real-space displacement vectors.

    Parameters
    ----------
    masses : np.ndarray
        1D array of atomic masses for each atom in the system.
        Shape: (N_atoms,). Unit: AMU.
    dR : np.ndarray
        2D array of atomic displacement vectors between two configurations.
        Shape: (N_atoms, 3). Unit: Angstrom (Å).
    eigenvectors : np.ndarray
        3D array representing the normal mode eigenvectors matrix.
        Shape: (N_modes, N_atoms, 3). Dimensionless (normalized).

    Returns
    -------
    np.ndarray
        1D array of projected configuration coordinates for each mode k.
        Shape: (N_modes,). Unit: SI units (kg^{1/2} * m).

    See Also
    --------
    calc_qks_vectorized : Faster, loopless alternative for this calculation.
    """
    qks = []
    sqrt_m = np.sqrt(masses)
    for k in range(len(eigenvectors)):
        qk = np.sum(
            [sqrt_m[i] * np.dot(dR[i], eigenvectors[k][i]) for i in range(len(masses))]
        )
        qk = qk * ANG2M * np.sqrt(AMU2KG)
        qks.append(qk)
    return np.array(qks)

calc_qks_force_mode

calc_qks_force_mode(masses, forces, eigenvectors, frequencies_eV)

Project atomic force differences onto phonon normal modes using energies in eV.

This function calculates the dimensioned configuration coordinates using a loop-based approach, derived from the excited state forces acting on the ground state geometry. Acoustic or zero-frequency modes are protected to prevent division-by-zero errors.

Parameters:

Name	Type	Description	Default
masses	ndarray	1D array of atomic masses for each atom in the system. Shape: (N_atoms,). Unit: AMU.	required
forces	ndarray	2D array containing the difference in force vectors between the excited and ground states. Shape: (N_atoms, 3). Unit: eV/Angstrom (eV/Å).	required
eigenvectors	ndarray	3D array representing the normal mode eigenvectors matrix. Shape: (N_modes, N_atoms, 3). Dimensionless (normalized).	required
frequencies_eV	ndarray	1D array of Gamma-point phonon mode energies. Shape: (N_modes,). Unit: Electron-volts (eV).	required

Returns:

Type	Description
ndarray	1D array of projected configuration coordinates for each mode k. Shape: (N_modes,). Unit: SI units (kg^{½} * m).

See Also

calc_qks_force_vectorized : Faster, loopless alternative for this calculation.

Source code in defectpl/utils.py

def calc_qks_force_mode(
    masses: np.ndarray,
    forces: np.ndarray,
    eigenvectors: np.ndarray,
    frequencies_eV: np.ndarray,
) -> np.ndarray:
    """
    Project atomic force differences onto phonon normal modes using energies in eV.

    This function calculates the dimensioned configuration coordinates using a
    loop-based approach, derived from the excited state forces acting on the
    ground state geometry. Acoustic or zero-frequency modes are protected to
    prevent division-by-zero errors.

    Parameters
    ----------
    masses : np.ndarray
        1D array of atomic masses for each atom in the system.
        Shape: (N_atoms,). Unit: AMU.
    forces : np.ndarray
        2D array containing the difference in force vectors between the excited
        and ground states. Shape: (N_atoms, 3). Unit: eV/Angstrom (eV/Å).
    eigenvectors : np.ndarray
        3D array representing the normal mode eigenvectors matrix.
        Shape: (N_modes, N_atoms, 3). Dimensionless (normalized).
    frequencies_eV : np.ndarray
        1D array of Gamma-point phonon mode energies.
        Shape: (N_modes,). Unit: Electron-volts (eV).

    Returns
    -------
    np.ndarray
        1D array of projected configuration coordinates for each mode k.
        Shape: (N_modes,). Unit: SI units (kg^{1/2} * m).

    See Also
    --------
    calc_qks_force_vectorized : Faster, loopless alternative for this calculation.
    """
    qks = []
    sqrt_m = np.sqrt(masses)

    for k in range(len(eigenvectors)):
        freq_ev = frequencies_eV[k]

        # Handle the acoustic/zero-frequency modes safely
        if np.isclose(freq_ev, 0.0, atol=1e-5):
            qks.append(0.0)
            continue

        # 1. Raw projection: sum_i (1 / sqrt(M_i)) * dot(F_i, e_{k;i}) -> (eV / Angstrom) / sqrt(AMU)
        proj_sum = np.sum(
            [
                (1.0 / sqrt_m[i]) * np.dot(forces[i], eigenvectors[k][i])
                for i in range(len(masses))
            ]
        )

        # 2. Convert raw projection numerator to SI units -> (Joules / meter) / sqrt(kg)
        proj_sum_SI = proj_sum * (EV2J / ANG2M) / np.sqrt(AMU2KG)

        # 3. Convert energy in eV to angular frequency omega_k (rad/s)
        omega_k = (freq_ev * EV2J) / HBAR_JS

        # 4. Divide by omega_k^2 to find configuration coordinate
        qk_SI = (1.0 / (omega_k**2)) * proj_sum_SI

        qks.append(qk_SI)

    return np.array(qks)

Huang–Rhys factors and spectral function

calc_Sks

calc_Sks(qks, frequencies)

Compute partial Huang–Rhys factors for each Gamma-point phonon mode.

Parameters:

Name	Type	Description	Default
qks	ndarray	Mode-projected configuration coordinates in SI units (kg^(½)·m), shape (nmodes,). Obtained from :func:calc_qks or :func:calc_qks_force_mode.	required
frequencies	ndarray	Phonon mode energies in eV, shape (nmodes,).	required

Returns:

Type	Description
ndarray	Partial Huang–Rhys factors :math:S_k, dimensionless, shape (nmodes,).

Notes

The partial HR factor is [Alkauskas 2014a, Eq. 4]:

.. math::

S_k = \frac{\omega_k q_k^2}{2\hbar}

where :math:q_k is in SI units (kg^(½)·m) and :math:\omega_k is converted from eV via :math:\omega_k = E_k / \hbar. The total HR factor is :math:S = \sum_k S_k.

Source code in defectpl/utils.py

def calc_Sks(qks: np.ndarray, frequencies: np.ndarray) -> np.ndarray:
    """
    Compute partial Huang–Rhys factors for each Gamma-point phonon mode.

    Parameters
    ----------
    qks : np.ndarray
        Mode-projected configuration coordinates in SI units (kg^(1/2)·m),
        shape ``(nmodes,)``.  Obtained from :func:`calc_qks` or
        :func:`calc_qks_force_mode`.
    frequencies : np.ndarray
        Phonon mode energies in eV, shape ``(nmodes,)``.

    Returns
    -------
    np.ndarray
        Partial Huang–Rhys factors :math:`S_k`, dimensionless, shape ``(nmodes,)``.

    Notes
    -----
    The partial HR factor is [Alkauskas 2014a, Eq. 4]:

    .. math::

        S_k = \\frac{\\omega_k q_k^2}{2\\hbar}

    where :math:`q_k` is in SI units (kg^(1/2)·m) and :math:`\\omega_k`
    is converted from eV via :math:`\\omega_k = E_k / \\hbar`.
    The total HR factor is :math:`S = \\sum_k S_k`.
    """
    return frequencies * qks**2 / (2 * HBAR_JS * HBAR_EVS)

calc_S_omega

calc_S_omega(frequencies, Sks, omega_range, sigma=0.006)

Compute the continuous electron–phonon spectral density function S(ω).

Each discrete mode contribution :math:S_k \delta(\omega - \omega_k) is replaced by a normalized Gaussian of width sigma, giving:

.. math::

S(\omega) = \sum_k S_k \, g(\omega - \omega_k, \sigma)

Parameters:

Name	Type	Description	Default
frequencies	ndarray	Phonon mode energies in eV, shape (nmodes,).	required
Sks	ndarray	Partial Huang–Rhys factors, shape (nmodes,).	required
omega_range	list of float	Energy grid specification [start, stop, npoints] in eV.	required
sigma	float	Gaussian broadening width in eV. The default 6 meV is suitable for defect calculations; increase for broad sidebands.	6e-3

Returns:

Type	Description
ndarray	Spectral density :math:S(\omega) evaluated on the energy grid, shape (npoints,).

See Also

calc_St : Transforms S(ω) to the time domain via inverse FFT.

Source code in defectpl/utils.py

def calc_S_omega(
    frequencies: np.ndarray,
    Sks: np.ndarray,
    omega_range: List[float],
    sigma: float = 6e-3,
) -> np.ndarray:
    """
    Compute the continuous electron–phonon spectral density function S(ω).

    Each discrete mode contribution :math:`S_k \\delta(\\omega - \\omega_k)` is
    replaced by a normalized Gaussian of width *sigma*, giving:

    .. math::

        S(\\omega) = \\sum_k S_k \\, g(\\omega - \\omega_k, \\sigma)

    Parameters
    ----------
    frequencies : np.ndarray
        Phonon mode energies in eV, shape ``(nmodes,)``.
    Sks : np.ndarray
        Partial Huang–Rhys factors, shape ``(nmodes,)``.
    omega_range : list of float
        Energy grid specification ``[start, stop, npoints]`` in eV.
    sigma : float, default 6e-3
        Gaussian broadening width in eV.  The default 6 meV is suitable for
        defect calculations; increase for broad sidebands.

    Returns
    -------
    np.ndarray
        Spectral density :math:`S(\\omega)` evaluated on the energy grid,
        shape ``(npoints,)``.

    See Also
    --------
    calc_St : Transforms S(ω) to the time domain via inverse FFT.
    """
    npts = int(omega_range[2])
    omega_start, omega_stop = omega_range[0], omega_range[1]
    dw = (omega_stop - omega_start) / (npts - 1)

    # Pad working grid by 5σ on each side so Gaussian tails from near-edge
    # modes are captured correctly.
    half_w = int(np.ceil(5.0 * sigma / dw))
    npts_pad = npts + 2 * half_w
    omega_start_pad = omega_start - half_w * dw

    # Distribute each S_k between its two nearest grid points (linear
    # interpolation) to eliminate per-mode quantisation error.
    pos = (frequencies - omega_start_pad) / dw
    idx_lo = np.floor(pos).astype(int)
    frac = pos - idx_lo
    S_discrete = np.zeros(npts_pad)
    in_lo = (idx_lo >= 0) & (idx_lo < npts_pad)
    in_hi = (idx_lo + 1 >= 0) & (idx_lo + 1 < npts_pad)
    np.add.at(S_discrete, idx_lo[in_lo], Sks[in_lo] * (1.0 - frac[in_lo]))
    np.add.at(S_discrete, (idx_lo + 1)[in_hi], Sks[in_hi] * frac[in_hi])

    # Gaussian convolution kernel spanning ±5σ
    kx = np.arange(-half_w, half_w + 1) * dw
    kernel = np.exp(-0.5 * (kx / sigma) ** 2) / (sigma * np.sqrt(2.0 * np.pi))

    # Linear (non-circular) FFT convolution
    n_fft = npts_pad + len(kernel) - 1
    S_conv = np.fft.irfft(
        np.fft.rfft(S_discrete, n=n_fft) * np.fft.rfft(kernel, n=n_fft),
        n=n_fft,
    )
    # Peak for a mode at original grid index i appears in S_conv at 2*half_w+i
    return S_conv[2 * half_w : 2 * half_w + npts]

gaussian_broadening

gaussian_broadening(omega, omega_k, sigma)

Evaluate a normalized Gaussian centred at omega_k with width sigma.

Used by :func:calc_S_omega to replace each discrete delta-function mode contribution with a Gaussian of unit area.

Parameters:

Name	Type	Description	Default
omega	float or ndarray	Evaluation energy in eV (can be a grid).	required
omega_k	float or ndarray	Centre energy in eV (phonon mode frequency or broadcast array).	required
sigma	float	Standard deviation (broadening width) in eV.	required

Returns:

Type	Description
float or ndarray	Gaussian value(s) at omega in eV^{-1}.

Source code in defectpl/utils.py

def gaussian_broadening(
    omega: Union[float, np.ndarray], omega_k: Union[float, np.ndarray], sigma: float
) -> Union[float, np.ndarray]:
    """
    Evaluate a normalized Gaussian centred at *omega_k* with width *sigma*.

    Used by :func:`calc_S_omega` to replace each discrete delta-function
    mode contribution with a Gaussian of unit area.

    Parameters
    ----------
    omega : float or np.ndarray
        Evaluation energy in eV (can be a grid).
    omega_k : float or np.ndarray
        Centre energy in eV (phonon mode frequency or broadcast array).
    sigma : float
        Standard deviation (broadening width) in eV.

    Returns
    -------
    float or np.ndarray
        Gaussian value(s) at *omega* in eV^{-1}.
    """
    return (1.0 / (sigma * np.sqrt(2.0 * np.pi))) * np.exp(
        -0.5 * ((omega - omega_k) / sigma) ** 2
    )

Phonon IPR

calc_IPR

calc_IPR(eigenvectors)

Calculate the site-projected Inverse Participation Ratio (IPR) for phonon modes.

For each normal mode i the per-atom weight is:

p_{i,a} = sum_{xyz} e_{i,a,xyz}^2

and the IPR is defined as:

IPR_i = sum_a p_{i,a}^2 / (sum_a p_{i,a})^2

Range: 1/N (fully delocalized over N atoms) to 1 (fully localized on one atom).

Parameters:

Name	Type	Description	Default
eigenvectors	ndarray	Phonon normal mode eigenvectors, shape (nmodes, natoms, 3). Works for both normalized (phonopy) and un-normalized eigenvectors.	required

Returns:

Type	Description
ndarray	IPR value for each normal mode, shape (nmodes,).

Source code in defectpl/utils.py

def calc_IPR(eigenvectors: np.ndarray) -> np.ndarray:
    """
    Calculate the site-projected Inverse Participation Ratio (IPR) for phonon modes.

    For each normal mode i the per-atom weight is:

        p_{i,a} = sum_{xyz} e_{i,a,xyz}^2

    and the IPR is defined as:

        IPR_i = sum_a p_{i,a}^2 / (sum_a p_{i,a})^2

    Range: 1/N (fully delocalized over N atoms) to 1 (fully localized on one atom).

    Parameters
    ----------
    eigenvectors : np.ndarray
        Phonon normal mode eigenvectors, shape (nmodes, natoms, 3).
        Works for both normalized (phonopy) and un-normalized eigenvectors.

    Returns
    -------
    np.ndarray
        IPR value for each normal mode, shape (nmodes,).
    """
    participations = np.sum(eigenvectors * eigenvectors, axis=2)
    return np.sum(participations**2, axis=1) / np.sum(participations, axis=1) ** 2

calc_IPR_alkauskas

calc_IPR_alkauskas(eigenvectors)

Calculate the phonon IPR following Alkauskas et al. (New J. Phys. 16, 073026, 2014), eq. 12.

For each normal mode k the per-atom weight is:

p_{k,α} = Σ_i e_{k,α,i}²

and the Alkauskas IPR is:

IPR_k = (Σ_α p_{k,α})² / Σ_α p_{k,α}²

For phonopy-normalized eigenvectors (Σ_α p_{k,α} = 1) this reduces to 1 / Σ_α p_{k,α}².

Range: 1 (fully localized on one atom) to N (fully delocalized over N atoms).

This is the reciprocal of the traditional IPR returned by :func:calc_IPR, and equals the paper's equation 12 for normalized eigenvectors. The localization ratio β_k = N × IPR_k⁻¹ (traditional) = N / IPR_k (Alkauskas) — both conventions give the same β_k via natoms * calc_IPR(evecs) or natoms / calc_IPR_alkauskas(evecs).

Parameters:

Name	Type	Description	Default
eigenvectors	ndarray	Phonon normal mode eigenvectors, shape (nmodes, natoms, 3). Works for both normalized (phonopy) and un-normalized eigenvectors.	required

Returns:

Type	Description
ndarray	Alkauskas IPR value for each normal mode, shape (nmodes,).

Source code in defectpl/utils.py

def calc_IPR_alkauskas(eigenvectors: np.ndarray) -> np.ndarray:
    """
    Calculate the phonon IPR following Alkauskas et al. (New J. Phys. 16, 073026, 2014), eq. 12.

    For each normal mode k the per-atom weight is:

        p_{k,α} = Σ_i e_{k,α,i}²

    and the Alkauskas IPR is:

        IPR_k = (Σ_α p_{k,α})² / Σ_α p_{k,α}²

    For phonopy-normalized eigenvectors (Σ_α p_{k,α} = 1) this reduces to 1 / Σ_α p_{k,α}².

    Range: 1 (fully localized on one atom) to N (fully delocalized over N atoms).

    This is the reciprocal of the traditional IPR returned by :func:`calc_IPR`, and equals the
    paper's equation 12 for normalized eigenvectors.  The localization ratio β_k = N × IPR_k⁻¹
    (traditional) = N / IPR_k (Alkauskas) — both conventions give the same β_k via
    ``natoms * calc_IPR(evecs)`` or ``natoms / calc_IPR_alkauskas(evecs)``.

    Parameters
    ----------
    eigenvectors : np.ndarray
        Phonon normal mode eigenvectors, shape (nmodes, natoms, 3).
        Works for both normalized (phonopy) and un-normalized eigenvectors.

    Returns
    -------
    np.ndarray
        Alkauskas IPR value for each normal mode, shape (nmodes,).
    """
    participations = np.sum(eigenvectors * eigenvectors, axis=2)
    return np.sum(participations, axis=1) ** 2 / np.sum(participations**2, axis=1)

Temperature-dependent PL and absorption

calc_phonon_occupation

calc_phonon_occupation(frequencies, temperature)

Compute the Bose-Einstein phonon occupation number for each mode.

Parameters:

Name	Type	Description	Default
frequencies	ndarray	Phonon frequencies in eV, shape (nmodes,).	required
temperature	float	Temperature in Kelvin. Pass 0 to get an all-zero array (T = 0 limit).	required

Returns:

Type	Description
ndarray	Occupation numbers :math:\bar{n}_k(T), shape (nmodes,).

Notes

Implements Jin et al. (2021), Eq. 9 / Alkauskas et al. (2014), Eq. 9:

.. math::

\bar{n}_k(T) = \frac{1}{e^{\hbar\omega_k / k_B T} - 1}

At T = 0 all occupation numbers are zero (no thermally excited phonons).

Source code in defectpl/utils.py

def calc_phonon_occupation(frequencies: np.ndarray, temperature: float) -> np.ndarray:
    """
    Compute the Bose-Einstein phonon occupation number for each mode.

    Parameters
    ----------
    frequencies : np.ndarray
        Phonon frequencies in **eV**, shape ``(nmodes,)``.
    temperature : float
        Temperature in **Kelvin**.  Pass 0 to get an all-zero array (T = 0 limit).

    Returns
    -------
    np.ndarray
        Occupation numbers :math:`\\bar{n}_k(T)`, shape ``(nmodes,)``.

    Notes
    -----
    Implements Jin *et al.* (2021), Eq. 9 / Alkauskas *et al.* (2014), Eq. 9:

    .. math::

        \\bar{n}_k(T) = \\frac{1}{e^{\\hbar\\omega_k / k_B T} - 1}

    At T = 0 all occupation numbers are zero (no thermally excited phonons).
    """
    if temperature <= 0.0:
        return np.zeros(len(frequencies))
    x = np.asarray(frequencies, dtype=float) / (KB_EV * temperature)
    return 1.0 / (np.expm1(x))

calc_C_omega

calc_C_omega(frequencies, Sks, nks, omega_range, sigma=0.006)

Build the thermal electron–phonon spectral density C(ℏω, T) on a uniform grid.

Parameters:

Name	Type	Description	Default
frequencies	ndarray	Phonon frequencies in eV, shape (nmodes,).	required
Sks	ndarray	Partial Huang–Rhys factors, shape (nmodes,).	required
nks	ndarray	Bose-Einstein occupation numbers :math:\bar{n}_k(T), shape (nmodes,). Pass an all-zero array for T = 0 (returns a zero grid immediately).	required
omega_range	list[ω_min, ω_max, n_points]	Energy grid parameters in eV.	required
sigma	float or (float, float)	Gaussian broadening in eV. Scalar → uniform; 2-tuple → linearly interpolated from lowest to highest frequency mode. Default 6 meV.	0.006

Returns:

Type	Description
ndarray	:math:C(\hbar\omega, T) in eV\ :sup:-1, shape (n_points,).

Notes

Implements Jin et al. (2021), Eq. 12:

.. math::

C(\hbar\omega, T) = \sum_k \bar{n}_k(T)\, S_k\,
\delta(\hbar\omega - \hbar\omega_k)

with δ-functions replaced by Gaussians (with optional frequency-dependent width σ(ω_k)). At T = 0 the function returns a zero array without allocating the FFT buffers.

Source code in defectpl/utils.py

def calc_C_omega(
    frequencies: np.ndarray,
    Sks: np.ndarray,
    nks: np.ndarray,
    omega_range: List[Union[float, int]],
    sigma: Union[float, Tuple[float, float]] = 6e-3,
) -> np.ndarray:
    """
    Build the thermal electron–phonon spectral density C(ℏω, T) on a uniform grid.

    Parameters
    ----------
    frequencies : np.ndarray
        Phonon frequencies in **eV**, shape ``(nmodes,)``.
    Sks : np.ndarray
        Partial Huang–Rhys factors, shape ``(nmodes,)``.
    nks : np.ndarray
        Bose-Einstein occupation numbers :math:`\\bar{n}_k(T)`, shape ``(nmodes,)``.
        Pass an all-zero array for T = 0 (returns a zero grid immediately).
    omega_range : list [ω_min, ω_max, n_points]
        Energy grid parameters in **eV**.
    sigma : float or (float, float), optional
        Gaussian broadening in **eV**.  Scalar → uniform; 2-tuple → linearly
        interpolated from lowest to highest frequency mode.  Default 6 meV.

    Returns
    -------
    np.ndarray
        :math:`C(\\hbar\\omega, T)` in eV\ :sup:`-1`, shape ``(n_points,)``.

    Notes
    -----
    Implements Jin *et al.* (2021), Eq. 12:

    .. math::

        C(\\hbar\\omega, T) = \\sum_k \\bar{n}_k(T)\\, S_k\\,
        \\delta(\\hbar\\omega - \\hbar\\omega_k)

    with δ-functions replaced by Gaussians (with optional frequency-dependent
    width σ(ω_k)).  At T = 0 the function returns a zero array without
    allocating the FFT buffers.
    """
    if not np.any(nks):
        return np.zeros(int(omega_range[2]))

    weights = nks * Sks
    sigmas = _sigma_per_mode(frequencies, sigma)

    npts = int(omega_range[2])
    omega_start, omega_stop = omega_range[0], omega_range[1]
    dw = (omega_stop - omega_start) / (npts - 1)

    # Per-mode sigma requires individual Gaussian placement rather than a
    # single FFT-convolution kernel.  We build a padded grid, place each
    # mode's Gaussian explicitly, then trim.
    sigma_max = float(np.max(sigmas))
    half_w = int(np.ceil(5.0 * sigma_max / dw))
    npts_pad = npts + 2 * half_w
    omega_start_pad = omega_start - half_w * dw
    omega_pad = omega_start_pad + np.arange(npts_pad) * dw

    C_grid = np.zeros(npts_pad)
    for wk, ck, sk in zip(frequencies, weights, sigmas):
        if ck == 0.0:
            continue
        gauss = np.exp(-0.5 * ((omega_pad - wk) / sk) ** 2) / (
            sk * np.sqrt(2.0 * np.pi)
        )
        C_grid += ck * gauss

    return C_grid[half_w : half_w + npts]

calc_Ct

calc_Ct(C_omega)

Transform C(ℏω, T) to the time domain, returning the real cosine component.

Parameters:

Name	Type	Description	Default
C_omega	ndarray	Thermal spectral density :math:C(\hbar\omega, T) on a uniform grid, shape (npoints,).	required

Returns:

Type	Description
ndarray	Real-valued :math:C(t, T) = \sum_k \bar{n}_k S_k \cos(\omega_k t), shape (npoints,).

Notes

Implements Jin et al. (2021), using the relation:

.. math::

C(t, T) = \operatorname{Re}\!\left[
    \int_0^\infty C(\hbar\omega, T)\, e^{-i\omega t}\, d(\hbar\omega)
\right]

which is the real part of the inverse Fourier transform of C(ω, T), identical in structure to :func:calc_St but restricted to its real part.

Source code in defectpl/utils.py

def calc_Ct(C_omega: np.ndarray) -> np.ndarray:
    """
    Transform C(ℏω, T) to the time domain, returning the real cosine component.

    Parameters
    ----------
    C_omega : np.ndarray
        Thermal spectral density :math:`C(\\hbar\\omega, T)` on a uniform grid,
        shape ``(npoints,)``.

    Returns
    -------
    np.ndarray
        Real-valued :math:`C(t, T) = \\sum_k \\bar{n}_k S_k \\cos(\\omega_k t)`,
        shape ``(npoints,)``.

    Notes
    -----
    Implements Jin *et al.* (2021), using the relation:

    .. math::

        C(t, T) = \\operatorname{Re}\\!\\left[
            \\int_0^\\infty C(\\hbar\\omega, T)\\, e^{-i\\omega t}\\, d(\\hbar\\omega)
        \\right]

    which is the real part of the inverse Fourier transform of C(ω, T),
    identical in structure to :func:`calc_St` but restricted to its real part.
    """
    return np.real(calc_St(C_omega))

calc_C_total

calc_C_total(nks, Sks)

Compute the zero-time thermal correction C(0, T) = Σ_k n̄_k S_k.

Parameters:

Name	Type	Description	Default
nks	ndarray	Bose-Einstein occupation numbers, shape (nmodes,).	required
Sks	ndarray	Partial Huang–Rhys factors, shape (nmodes,).	required

Returns:

Type	Description
float	:math:C(0, T) = \sum_k \bar{n}_k(T)\, S_k.

Notes

This is the thermal analogue of the total Huang–Rhys factor S = Σ_k S_k. At T = 0 it evaluates to zero, keeping the generating function unchanged.

Source code in defectpl/utils.py

def calc_C_total(nks: np.ndarray, Sks: np.ndarray) -> float:
    """
    Compute the zero-time thermal correction C(0, T) = Σ_k n̄_k S_k.

    Parameters
    ----------
    nks : np.ndarray
        Bose-Einstein occupation numbers, shape ``(nmodes,)``.
    Sks : np.ndarray
        Partial Huang–Rhys factors, shape ``(nmodes,)``.

    Returns
    -------
    float
        :math:`C(0, T) = \\sum_k \\bar{n}_k(T)\\, S_k`.

    Notes
    -----
    This is the thermal analogue of the total Huang–Rhys factor S = Σ_k S_k.
    At T = 0 it evaluates to zero, keeping the generating function unchanged.
    """
    return float(np.dot(nks, Sks))

calc_Absorption_Intensity

calc_Absorption_Intensity(Gts, EZPL, resolution)

Compute the absorption spectral function and absorption intensity from G(t, T).

Parameters:

Name	Type	Description	Default
Gts	ndarray	Complex PL generating function G(t, T) from :func:calc_Gts, shape (npoints,).	required
EZPL	float	Zero-phonon line energy in eV.	required
resolution	float	Spectral grid density in points per eV.	required

Returns:

Name	Type	Description
A_abs	ndarray	Absorption spectral function :math:A_{abs}(\hbar\omega), shape (npoints,).
intensity_abs	ndarray	Absorption intensity :math:\alpha(\hbar\omega) \propto \omega\, A_{abs}, shape (npoints,).

Notes

The absorption generating function is the complex conjugate of the PL generating function (Jin et al. 2021; derivation from Eq. 8):

.. math::

G_{abs}(t, T) = G^*(t, T)

because replacing :math:S(t) = \sum_k S_k e^{i\omega_k t} with :math:S^*(t) = \sum_k S_k e^{-i\omega_k t} while leaving the real thermal correction :math:C(t,T) unchanged yields exactly :math:\overline{G(t,T)}.

The sideband appears on the high-energy side of the ZPL (energy axis :math:E_{ZPL} + E), opposite to PL emission. The prefactor is :math:\omega (linear) rather than :math:\omega^3.

Source code in defectpl/utils.py

def calc_Absorption_Intensity(
    Gts: np.ndarray, EZPL: float, resolution: float
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Compute the absorption spectral function and absorption intensity from G(t, T).

    Parameters
    ----------
    Gts : np.ndarray
        Complex PL generating function G(t, T) from :func:`calc_Gts`,
        shape ``(npoints,)``.
    EZPL : float
        Zero-phonon line energy in eV.
    resolution : float
        Spectral grid density in points per eV.

    Returns
    -------
    A_abs : np.ndarray
        Absorption spectral function :math:`A_{abs}(\\hbar\\omega)`,
        shape ``(npoints,)``.
    intensity_abs : np.ndarray
        Absorption intensity :math:`\\alpha(\\hbar\\omega) \\propto \\omega\\, A_{abs}`,
        shape ``(npoints,)``.

    Notes
    -----
    The absorption generating function is the complex conjugate of the PL
    generating function (Jin *et al.* 2021; derivation from Eq. 8):

    .. math::

        G_{abs}(t, T) = G^*(t, T)

    because replacing :math:`S(t) = \\sum_k S_k e^{i\\omega_k t}` with
    :math:`S^*(t) = \\sum_k S_k e^{-i\\omega_k t}` while leaving the real
    thermal correction :math:`C(t,T)` unchanged yields exactly
    :math:`\\overline{G(t,T)}`.

    The sideband appears on the **high-energy side** of the ZPL (energy axis
    :math:`E_{ZPL} + E`), opposite to PL emission.  The prefactor is
    :math:`\\omega` (linear) rather than :math:`\\omega^3`.
    """
    G_abs = np.conj(Gts)
    A1_abs = np.fft.fft(G_abs)
    n = len(A1_abs)
    shift_idx = int(EZPL * resolution)
    j = np.arange(n)
    # Shift in the opposite direction vs PL: sideband at j > shift_idx (higher E)
    A_abs = A1_abs[(j - shift_idx) % n]
    omega_1 = j / resolution
    return A_abs, A_abs * omega_1

calc_effective_phonon_frequency

calc_effective_phonon_frequency(frequencies, delta_Q_k)

Compute the effective phonon frequency for the 1D configurational coordinate model.

Parameters:

Name	Type	Description	Default
frequencies	ndarray	Phonon frequencies in eV, shape (nmodes,).	required
delta_Q_k	ndarray	Mode-projected configurational displacements q_k in amu\ :sup:1/2·Å, shape (nmodes,).	required

Returns:

Type	Description
float	Effective phonon frequency :math:\Omega in eV.

Notes

Implements Jin et al. (2021), Eq. 16:

.. math::

\Omega = \frac{\sum_k \omega_k^2\, \Delta Q_k^2}{\sum_k \Delta Q_k^2}

This weighted average collapses the full phonon spectrum onto a single representative mode for use in the 1D harmonic oscillator lineshape model (:class:~defectpl.defectpl.VibrationalSpectra1D).

Source code in defectpl/utils.py

def calc_effective_phonon_frequency(
    frequencies: np.ndarray, delta_Q_k: np.ndarray
) -> float:
    """
    Compute the effective phonon frequency for the 1D configurational coordinate model.

    Parameters
    ----------
    frequencies : np.ndarray
        Phonon frequencies in **eV**, shape ``(nmodes,)``.
    delta_Q_k : np.ndarray
        Mode-projected configurational displacements q_k in amu\\ :sup:`1/2`·Å,
        shape ``(nmodes,)``.

    Returns
    -------
    float
        Effective phonon frequency :math:`\\Omega` in **eV**.

    Notes
    -----
    Implements Jin *et al.* (2021), Eq. 16:

    .. math::

        \\Omega = \\frac{\\sum_k \\omega_k^2\\, \\Delta Q_k^2}{\\sum_k \\Delta Q_k^2}

    This weighted average collapses the full phonon spectrum onto a single
    representative mode for use in the 1D harmonic oscillator lineshape model
    (:class:`~defectpl.defectpl.VibrationalSpectra1D`).
    """
    dq2 = delta_Q_k**2
    denom = float(np.sum(dq2))
    if denom == 0.0:
        raise ValueError(
            "All delta_Q_k values are zero; effective frequency is undefined."
        )
    return float(np.sum(frequencies**2 * dq2) / denom)

Generating function and spectrum

calc_St

calc_St(S_omega)

Transform the electron–phonon spectral density S(ω) to the time domain S(t).

Parameters:

Name	Type	Description	Default
S_omega	ndarray	Spectral density on a uniform energy grid, shape (npoints,).	required

Returns:

Type	Description
ndarray	Complex-valued time-domain array :math:S(t), shape (npoints,).

Notes

Implements [Alkauskas 2014a, Eq. 9]:

.. math::

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

via an inverse FFT with phase centering.

See Also

calc_Gts : Constructs G(t) from S(t).

Source code in defectpl/utils.py

def calc_St(S_omega: np.ndarray) -> np.ndarray:
    """
    Transform the electron–phonon spectral density S(ω) to the time domain S(t).

    Parameters
    ----------
    S_omega : np.ndarray
        Spectral density on a uniform energy grid, shape ``(npoints,)``.

    Returns
    -------
    np.ndarray
        Complex-valued time-domain array :math:`S(t)`, shape ``(npoints,)``.

    Notes
    -----
    Implements [Alkauskas 2014a, Eq. 9]:

    .. math::

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

    via an inverse FFT with phase centering.

    See Also
    --------
    calc_Gts : Constructs G(t) from S(t).
    """
    Sts = np.fft.ifft(S_omega)
    return 2.0 * np.pi * np.fft.ifftshift(Sts)

calc_Gts

calc_Gts(Sts, total_HR, gamma, resolution, Cts=None, C_total=0.0)

Compute the generating function G(t, T) from the time-domain spectral function.

Parameters:

Name	Type	Description	Default
Sts	ndarray	Time-domain spectral function S(t), shape (npoints,).	required
total_HR	float	Total Huang–Rhys factor :math:S = \sum_k S_k (dimensionless).	required
gamma	float	ZPL broadening parameter in meV; applied as a Lorentzian decay :math:e^{-\gamma|t|} to reproduce finite ZPL linewidth.	required
resolution	float	Spectral grid density in points per eV (resolution = npoints / max_energy).	required
Cts	ndarray	Real-valued time-domain thermal correction C(t, T) from :func:calc_Ct, shape (npoints,). Pass None (default) for the T = 0 limit.	None
C_total	float	Zero-time thermal correction C(0, T) = Σ n̄_k S_k from :func:calc_C_total. Ignored when Cts is None. Default 0.0.	0.0

Returns:

Type	Description
ndarray	Complex generating function :math:G(t, T), shape (npoints,).

Notes

Implements Jin et al. (2021), Eq. 8 / Eq. 10:

.. math::

G(t, T) = e^{S(t) - S + 2C(t,T) - 2C(0,T)}\, e^{-\gamma|t|}

At T = 0 (Cts = None or all-zero) the thermal correction vanishes and the expression reduces to the Alkauskas (2014) formula :math:G(t) = e^{S(t)-S}.

Source code in defectpl/utils.py

def calc_Gts(
    Sts: np.ndarray,
    total_HR: float,
    gamma: float,
    resolution: float,
    Cts: np.ndarray = None,
    C_total: float = 0.0,
) -> np.ndarray:
    """
    Compute the generating function G(t, T) from the time-domain spectral function.

    Parameters
    ----------
    Sts : np.ndarray
        Time-domain spectral function S(t), shape ``(npoints,)``.
    total_HR : float
        Total Huang–Rhys factor :math:`S = \\sum_k S_k` (dimensionless).
    gamma : float
        ZPL broadening parameter in meV; applied as a Lorentzian decay
        :math:`e^{-\\gamma|t|}` to reproduce finite ZPL linewidth.
    resolution : float
        Spectral grid density in points per eV (``resolution = npoints / max_energy``).
    Cts : np.ndarray, optional
        Real-valued time-domain thermal correction C(t, T) from :func:`calc_Ct`,
        shape ``(npoints,)``.  Pass ``None`` (default) for the T = 0 limit.
    C_total : float, optional
        Zero-time thermal correction C(0, T) = Σ n̄_k S_k from :func:`calc_C_total`.
        Ignored when ``Cts`` is ``None``.  Default 0.0.

    Returns
    -------
    np.ndarray
        Complex generating function :math:`G(t, T)`, shape ``(npoints,)``.

    Notes
    -----
    Implements Jin *et al.* (2021), Eq. 8 / Eq. 10:

    .. math::

        G(t, T) = e^{S(t) - S + 2C(t,T) - 2C(0,T)}\\, e^{-\\gamma|t|}

    At T = 0 (``Cts = None`` or all-zero) the thermal correction vanishes and
    the expression reduces to the Alkauskas (2014) formula :math:`G(t) = e^{S(t)-S}`.
    """
    n = len(Sts)
    t = (1.0 / resolution) * (np.arange(n) - n / 2)
    correction = (2.0 * Cts - 2.0 * C_total) if Cts is not None else 0.0
    return np.exp(Sts - total_HR + correction) * np.exp(-gamma * np.abs(t))

calc_Spectrum_Intensity

calc_Spectrum_Intensity(Gts, EZPL, resolution)

Compute the optical spectral function A(ℏω) and PL intensity L(ℏω) from G(t).

Parameters:

Name	Type	Description	Default
Gts	ndarray	Complex generating function G(t), shape (npoints,).	required
EZPL	float	Zero-phonon line energy in eV. Used to shift the spectral axis so that the ZPL appears at the correct absolute photon energy.	required
resolution	float	Spectral grid density in points per eV (npoints / max_energy).	required

Returns:

Name	Type	Description
A	ndarray	Optical spectral function :math:A(\hbar\omega), shape (npoints,). Obtained as the FFT of G(t) with the ZPL shifted to EZPL.
intensity	ndarray	Normalized PL intensity :math:L(\hbar\omega) \propto \omega^3 A(\hbar\omega), shape (npoints,).

Notes

Implements [Alkauskas 2014a, Eq. 7]:

.. math::

A(E_{\text{ZPL}} - \hbar\omega) =
\frac{1}{2\pi} \int_{-\infty}^{\infty} G(t)\, e^{i\omega t}\, dt

The :math:\omega^3 prefactor originates from the photon density of states.

Source code in defectpl/utils.py

def calc_Spectrum_Intensity(
    Gts: np.ndarray, EZPL: float, resolution: float
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Compute the optical spectral function A(ℏω) and PL intensity L(ℏω) from G(t).

    Parameters
    ----------
    Gts : np.ndarray
        Complex generating function G(t), shape ``(npoints,)``.
    EZPL : float
        Zero-phonon line energy in eV.  Used to shift the spectral axis so that
        the ZPL appears at the correct absolute photon energy.
    resolution : float
        Spectral grid density in points per eV (``npoints / max_energy``).

    Returns
    -------
    A : np.ndarray
        Optical spectral function :math:`A(\\hbar\\omega)`, shape ``(npoints,)``.
        Obtained as the FFT of G(t) with the ZPL shifted to *EZPL*.
    intensity : np.ndarray
        Normalized PL intensity :math:`L(\\hbar\\omega) \\propto \\omega^3 A(\\hbar\\omega)`,
        shape ``(npoints,)``.

    Notes
    -----
    Implements [Alkauskas 2014a, Eq. 7]:

    .. math::

        A(E_{\\text{ZPL}} - \\hbar\\omega) =
        \\frac{1}{2\\pi} \\int_{-\\infty}^{\\infty} G(t)\\, e^{i\\omega t}\\, dt

    The :math:`\\omega^3` prefactor originates from the photon density of states.
    """
    A1 = np.fft.fft(Gts)
    n = len(A1)
    shift_idx = int(EZPL * resolution)
    j = np.arange(n)
    A = A1[(shift_idx - j) % n]
    omega_3 = (j / resolution) ** 3
    return A, A * omega_3

1D overlap integrals

calculate_hermite

calculate_hermite(n, x)

Compute the physicist's Hermite polynomial H_n(x) using recurrence relation stability.

Parameters:

Name	Type	Description	Default
n	int	The order index of the target Hermite polynomial.	required
x	float	The specific coordinate position coordinate value to evaluate.	required

Returns:

Type	Description
float	The evaluated numerical result of the physicist's Hermite polynomial.

Source code in defectpl/utils.py

def calculate_hermite(n: int, x: float) -> float:
    """
    Compute the physicist's Hermite polynomial H_n(x) using recurrence relation stability.

    Parameters
    ----------
    n : int
        The order index of the target Hermite polynomial.
    x : float
        The specific coordinate position coordinate value to evaluate.

    Returns
    -------
    float
        The evaluated numerical result of the physicist's Hermite polynomial.
    """
    if n == 0:
        return 1.0
    elif n == 1:
        return 2.0 * x
    y0, y1 = 1.0, 2.0 * x
    for k in range(2, n + 1):
        y0, y1 = y1, 2.0 * x * y1 - 2.0 * (k - 1) * y0
    return y1

calculate_overlap_element

calculate_overlap_element(m, n, rho, cosfi, sinfi)

Calculate the 1D analytical transition overlap integral between two harmonic oscillator states.

This implements analytical Franck-Condon factor calculations supporting models where ground and excited configurations can possess different vibrational frequencies.

Parameters:

Name	Type	Description	Default
m	int	Vibrational quantum number index mapping the excited state.	required
n	int	Vibrational quantum number index mapping the ground state.	required
rho	float	Dimensionless offset displacement parameter between harmonic minima.	required
cosfi	float	Mixing angle cosine property resolving frequency transformation shifts scaling.	required
sinfi	float	Mixing angle sine property resolving frequency transformation shifts scaling.	required

Returns:

Type	Description
float	The computed structural overlap integral element value between states (m, n).

Source code in defectpl/utils.py

def calculate_overlap_element(
    m: int, n: int, rho: float, cosfi: float, sinfi: float
) -> float:
    """
    Calculate the 1D analytical transition overlap integral between two harmonic oscillator states.

    This implements analytical Franck-Condon factor calculations supporting models
    where ground and excited configurations can possess different vibrational frequencies.

    Parameters
    ----------
    m : int
        Vibrational quantum number index mapping the excited state.
    n : int
        Vibrational quantum number index mapping the ground state.
    rho : float
        Dimensionless offset displacement parameter between harmonic minima.
    cosfi : float
        Mixing angle cosine property resolving frequency transformation shifts scaling.
    sinfi : float
        Mixing angle sine property resolving frequency transformation shifts scaling.

    Returns
    -------
    float
        The computed structural overlap integral element value between states (m, n).
    """
    pr1 = (-1) ** n * np.sqrt(2 * cosfi * sinfi) * np.exp(-(rho**2))
    ix = 0.0
    k1, k2 = divmod(m, 2)
    l1, l2 = divmod(n, 2)

    for kx in range(k1 + 1):
        for lx in range(l1 + 1):
            k, lq = 2 * kx + k2, 2 * lx + l2
            try:
                pr2 = (
                    np.sqrt(float(math.factorial(n) * math.factorial(m)))
                    / (
                        math.factorial(k)
                        * math.factorial(lq)
                        * math.factorial(k1 - kx)
                        * math.factorial(l1 - lx)
                    )
                    * 2.0 ** ((k + lq - m - n) // 2)
                )
            except ValueError:
                pr2 = 0
            pr3 = (sinfi**k) * (cosfi**lq)
            f = calculate_hermite(k + lq, rho)
            ix += pr1 * pr2 * pr3 * f

    return ix

Configuration Coordinate Diagram fitting

get_omega_from_pes

get_omega_from_pes(q_values, energies, ax=None, eval_grid=None)

Fit a harmonic parabola E(Q) = ½·M·ω²·Q² to PES data points.

A second-order polynomial is fit to (Q, E) data. The effective angular frequency ω is extracted from the curvature coefficient and returned in eV (as an effective phonon energy ℏω).

Parameters:

Name	Type	Description	Default
q_values	ndarray	Configuration coordinate values Q in amu^(½)·Å.	required
energies	ndarray	Total energies (or energy differences) in eV, zero-referenced to the minimum.	required
ax	Axes	If provided, the fitted parabola is plotted on this axes object.	None
eval_grid	ndarray	Q grid for plotting the fit curve. Required when ax is given.	None

Returns:

Type	Description
float	Effective phonon energy ℏω in eV derived from the parabola curvature.

Notes

The fit uses :func:numpy.polyfit (degree 2). The curvature coefficient a₂ (eV/amu·Å²) is related to the angular frequency by:

.. math::

\hbar\omega = \sqrt{2 a_2 / M_{\text{ref}}}

where an effective reference mass of 1 amu is assumed (since Q is already mass-weighted). The units are handled via the AMU2KG, ANG2M, and EV2J conversion factors.

Source code in defectpl/utils.py

def get_omega_from_pes(
    q_values: np.ndarray,
    energies: np.ndarray,
    ax=None,
    eval_grid: np.ndarray = None,
) -> float:
    """
    Fit a harmonic parabola E(Q) = ½·M·ω²·Q² to PES data points.

    A second-order polynomial is fit to *(Q, E)* data.  The effective angular
    frequency ω is extracted from the curvature coefficient and returned in
    **eV** (as an effective phonon energy ℏω).

    Parameters
    ----------
    q_values : np.ndarray
        Configuration coordinate values Q in amu^(1/2)·Å.
    energies : np.ndarray
        Total energies (or energy differences) in eV, zero-referenced to the minimum.
    ax : matplotlib.axes.Axes, optional
        If provided, the fitted parabola is plotted on this axes object.
    eval_grid : np.ndarray, optional
        Q grid for plotting the fit curve.  Required when *ax* is given.

    Returns
    -------
    float
        Effective phonon energy ℏω in **eV** derived from the parabola curvature.

    Notes
    -----
    The fit uses :func:`numpy.polyfit` (degree 2).  The curvature coefficient
    a₂ (eV/amu·Å²) is related to the angular frequency by:

    .. math::

        \\hbar\\omega = \\sqrt{2 a_2 / M_{\\text{ref}}}

    where an effective reference mass of 1 amu is assumed (since Q is already
    mass-weighted).  The units are handled via the AMU2KG, ANG2M, and EV2J
    conversion factors.
    """
    coeffs = np.polyfit(q_values, energies, 2)
    a2 = coeffs[0]  # coefficient of Q^2 in eV / (amu^1/2 * Å)^2

    # Convert curvature to SI: [eV / (amu * Å^2)] → [J / (kg * m^2)] = [rad^2/s^2]
    a2_SI = a2 * EV2J / (AMU2KG * ANG2M**2)
    omega_SI = np.sqrt(max(2.0 * a2_SI, 0.0))  # rad/s
    omega_eV = omega_SI * HBAR_JS / EV2J

    if ax is not None and eval_grid is not None:
        poly_fn = np.poly1d(coeffs)
        ax.plot(eval_grid, poly_fn(eval_grid), "--", lw=1)

    return float(omega_eV)

extract_important_properties

extract_important_properties(pl_engine, filename='important_properties.txt')

Write a structured summary of scalar PL properties to a text file.

Extracts S, DW factor, ΔQ, ΔR, and configuration metadata from a :class:~defectpl.defectpl.Photoluminescence instance. Force-mode runs omit ΔQ/ΔR since they are not directly computed in that mode.

Parameters:

Name	Type	Description	Default
pl_engine	Photoluminescence	A fully computed Photoluminescence instance.	required
filename	str	Output file path. Default "important_properties.txt".	'important_properties.txt'

Source code in defectpl/utils.py

def extract_important_properties(
    pl_engine, filename: str = "important_properties.txt"
) -> None:
    """
    Write a structured summary of scalar PL properties to a text file.

    Extracts S, DW factor, ΔQ, ΔR, and configuration metadata from a
    :class:`~defectpl.defectpl.Photoluminescence` instance.  Force-mode
    runs omit ΔQ/ΔR since they are not directly computed in that mode.

    Parameters
    ----------
    pl_engine : Photoluminescence
        A fully computed Photoluminescence instance.
    filename : str, optional
        Output file path.  Default ``"important_properties.txt"``.
    """
    import numpy as np

    # 1. Determine the calculation mode based on input state flags
    is_force_mode = (
        hasattr(pl_engine, "dF") and pl_engine.dF is not None and np.any(pl_engine.dF)
    )
    calc_mode = "Force Mode" if is_force_mode else "Displacement Mode"

    # 2. Build the structural layout header
    lines = [
        "===================================================",
        "                PROPERTIES SUMMARY                 ",
        "===================================================",
        f"Calculation Run Mode         : {calc_mode:>12}",
        f"Zero-Phonon Line (ZPL) Energy : {getattr(pl_engine, 'EZPL', 'N/A'):>12} eV",
        f"Total Huang-Rhys (HR) Factor : {getattr(pl_engine, 'HR_factor', 'N/A'):>12.6f}",
        f"Debye-Waller (DW) Factor     : {getattr(pl_engine, 'DW_factor', 'N/A'):>12.6f}",
    ]

    # 3. Conditional Step: Only inject spatial metrics if we are NOT in Force Mode
    if not is_force_mode:
        lines.extend(
            [
                f"Mass-Weighted Delta Q (delQ)  : {getattr(pl_engine, 'delQ', 'N/A'):>12.6f} amu^(1/2)*Å",
                f"Structural Delta R (delR)     : {getattr(pl_engine, 'delR', 'N/A'):>12.6f} Å",
            ]
        )

    # 4. Append operational configuration tracking metadata
    lines.extend(
        [
            "---------------------------------------------------",
            f"Total Number of Atoms (natoms): {getattr(pl_engine, 'natoms', 'N/A'):>12}",
            f"ZPL Broadening Factor (gamma) : {getattr(pl_engine, 'gamma', 'N/A'):>12.2f} meV",
            f"Gaussian Broadening (sigma)   : {getattr(pl_engine, 'sigma', 'N/A'):>12.6f} eV",
            f"Energy Mesh Resolution        : {getattr(pl_engine, 'resolution', 'N/A'):>12} points/eV",
            "===================================================",
        ]
    )

    # Write cleanly to file out
    with open(filename, "w", encoding="utf-8") as f:
        f.write("\n".join(lines) + "\n")

    print(f"Important parameters successfully exported to text summary: {filename}")