Skip to main content
`Ctrl`+`K`

#

# Source code for ase.dft.bandgap

    # fmt: off

    import warnings
    from dataclasses import dataclass

    import numpy as np

    spin_error = (
        'The spin keyword is no longer supported.  Please call the function '
        'with the energies corresponding to the desired spins.')
    _deprecated = object()

    def get_band_gap(calc, direct=False, spin=_deprecated):
        warnings.warn('Please use ase.dft.bandgap.bandgap() instead!')
        gap, (s1, k1, _n1), (s2, k2, _n2) = bandgap(calc, direct, spin=spin)
        ns = calc.get_number_of_spins()
        if ns == 2:
            return gap, (s1, k1), (s2, k2)
        return gap, k1, k2

    @dataclass
    class GapInfo:
        eigenvalues: np.ndarray

        def __post_init__(self):
            self._gapinfo = _bandgap(self.eigenvalues, direct=False)
            self._direct_gapinfo = _bandgap(self.eigenvalues, direct=True)

        @classmethod
        def fromcalc(cls, calc):
            kpts = calc.get_ibz_k_points()
            nk = len(kpts)
            ns = calc.get_number_of_spins()
            eigenvalues = np.array([[calc.get_eigenvalues(kpt=k, spin=s)
                                     for k in range(nk)]
                                    for s in range(ns)])

            efermi = calc.get_fermi_level()
            return cls(eigenvalues - efermi)

        def gap(self):
            return self._gapinfo

        def direct_gap(self):
            return self._direct_gapinfo

        @property
        def is_metallic(self) -> bool:
            return self._gapinfo[0] == 0.0

        @property
        def gap_is_direct(self) -> bool:
            """Whether the direct and indirect gaps are the same transition."""
            return self._gapinfo[1:] == self._direct_gapinfo[1:]

        def description(self, *, ibz_kpoints=None) -> str:
            """Return human-friendly description of direct/indirect gap.

            If ibz_k_points are given, coordinates are printed as well."""

            lines: list[str] = []
            add = lines.append

            def skn(skn):
                """Convert k-point indices (s, k, n) to string."""
                description = 's={}, k={}, n={}'.format(*skn)
                if ibz_kpoints is not None:
                    coordtxt = '[{:.2f}, {:.2f}, {:.2f}]'.format(
                        *ibz_kpoints[skn[1]])
                    description = f'{description}, {coordtxt}'
                return f'({description})'

            gap, skn1, skn2 = self.gap()
            direct_gap, skn_direct1, skn_direct2 = self.direct_gap()

            if self.is_metallic:
                add('No gap')
            else:
                add(f'Gap: {gap:.3f} eV')
                add('Transition (v -> c):')
                add(f'  {skn(skn1)} -> {skn(skn2)}')

            if self.gap_is_direct:
                add('No difference between direct/indirect transitions')
            else:
                add('Direct/indirect transitions are different')
                add(f'Direct gap: {direct_gap:.3f} eV')
                if skn_direct1[0] == skn_direct2[0]:
                    add(f'Transition at: {skn(skn_direct1)}')
                else:
                    transition = skn((f'{skn_direct1[0]}->{skn_direct2[0]}',
                                      *skn_direct1[1:]))
                    add(f'Transition at: {transition}')

            return '\n'.join(lines)

    def bandgap(calc=None, direct=False, spin=_deprecated,
                eigenvalues=None, efermi=None, output=None, kpts=None):
        """Calculates the band-gap.

        Parameters:

        calc: Calculator object
            Electronic structure calculator object.
        direct: bool
            Calculate direct band-gap.
        eigenvalues: ndarray of shape (nspin, nkpt, nband) or (nkpt, nband)
            Eigenvalues.
        efermi: float
            Fermi level (defaults to 0.0).

        Returns a (gap, p1, p2) tuple where p1 and p2 are tuples of indices of the
        valence and conduction points (s, k, n).

        Example:

        >>> gap, p1, p2 = bandgap(silicon.calc)
        >>> print(gap, p1, p2)
        1.2 (0, 0, 3), (0, 5, 4)
        >>> gap, p1, p2 = bandgap(silicon.calc, direct=True)
        >>> print(gap, p1, p2)
        3.4 (0, 0, 3), (0, 0, 4)
        """

        if spin is not _deprecated:
            raise RuntimeError(spin_error)

        if calc:
            kpts = calc.get_ibz_k_points()
            nk = len(kpts)
            ns = calc.get_number_of_spins()
            eigenvalues = np.array([[calc.get_eigenvalues(kpt=k, spin=s)
                                     for k in range(nk)]
                                    for s in range(ns)])
            if efermi is None:
                efermi = calc.get_fermi_level()

        efermi = efermi or 0.0

        gapinfo = GapInfo(eigenvalues - efermi)

        e_skn = gapinfo.eigenvalues
        if eigenvalues.ndim == 2:
            e_skn = e_skn[np.newaxis]  # spinors

        if not np.isfinite(e_skn).all():
            raise ValueError('Bad eigenvalues!')

        gap, (s1, k1, n1), (s2, k2, n2) = _bandgap(e_skn, direct)

        if eigenvalues.ndim != 3:
            p1 = (k1, n1)
            p2 = (k2, n2)
        else:
            p1 = (s1, k1, n1)
            p2 = (s2, k2, n2)

        return gap, p1, p2

    def _bandgap(e_skn, direct):
        """Helper function."""
        ns, nk, nb = e_skn.shape
        s1 = s2 = k1 = k2 = n1 = n2 = None

        N_sk = (e_skn < 0.0).sum(2)  # number of occupied bands

        # Check for bands crossing the fermi-level
        if ns == 1:
            if np.ptp(N_sk[0]) > 0:
                return 0.0, (None, None, None), (None, None, None)
        else:
            if (np.ptp(N_sk, axis=1) > 0).any():
                return 0.0, (None, None, None), (None, None, None)

        if (N_sk == 0).any() or (N_sk == nb).any():
            raise ValueError('Too few bands!')

        e_skn = np.array([[e_skn[s, k, N_sk[s, k] - 1:N_sk[s, k] + 1]
                           for k in range(nk)]
                          for s in range(ns)])
        ev_sk = e_skn[:, :, 0]  # valence band
        ec_sk = e_skn[:, :, 1]  # conduction band

        if ns == 1:
            s1 = 0
            s2 = 0
            gap, k1, k2 = find_gap(ev_sk[0], ec_sk[0], direct)
            n1 = N_sk[0, 0] - 1
            n2 = n1 + 1
            return gap, (0, k1, n1), (0, k2, n2)

        gap, k1, k2 = find_gap(ev_sk.ravel(), ec_sk.ravel(), direct)
        if direct:
            # Check also spin flips:
            for s in [0, 1]:
                g, k, _ = find_gap(ev_sk[s], ec_sk[1 - s], direct)
                if g < gap:
                    gap = g
                    k1 = k + nk * s
                    k2 = k + nk * (1 - s)

        if gap > 0.0:
            s1, k1 = divmod(k1, nk)
            s2, k2 = divmod(k2, nk)
            n1 = N_sk[s1, k1] - 1
            n2 = N_sk[s2, k2]
            return gap, (s1, k1, n1), (s2, k2, n2)
        return 0.0, (None, None, None), (None, None, None)

    def find_gap(ev_k, ec_k, direct):
        """Helper function."""
        if direct:
            gap_k = ec_k - ev_k
            k = gap_k.argmin()
            return gap_k[k], k, k
        kv = ev_k.argmax()
        kc = ec_k.argmin()
        return ec_k[kc] - ev_k[kv], kv, kc