This calculates the band filling correction based on excess of electrons/holes in CB/VB...

        Note that the total free holes and electrons may also be used for a "shallow donor/acceptor"
               correction with specified band shifts: +num_elec_cbm * Delta E_CBM (or -num_hole_vbm * Delta E_VBM)
               [this is done in the LevelShiftingCorrection class]

def perform_bandfill_corr(self, eigenvalues, kpoint_weights, potalign, vbm, cbm):
        """
        This calculates the band filling correction based on excess of electrons/holes in CB/VB...

        Note that the total free holes and electrons may also be used for a "shallow donor/acceptor"
               correction with specified band shifts: +num_elec_cbm * Delta E_CBM (or -num_hole_vbm * Delta E_VBM)
               [this is done in the LevelShiftingCorrection class]
        """
        bf_corr = 0.

        self.metadata["potalign"] = potalign
        self.metadata["num_hole_vbm"] = 0.
        self.metadata["num_elec_cbm"] = 0.

        if len(eigenvalues.keys()) == 1:  # needed because occupation of non-spin calcs is still 1... should be 2
            spinfctr = 2.
        elif len(eigenvalues.keys()) == 2:
            spinfctr = 1.
        else:
            raise ValueError("Eigenvalue keys greater than 2")

        # for tracking mid gap states...
        resolution = self.resolution
        shifted_cbm = potalign + cbm  # shift cbm with potential alignment
        shifted_vbm = potalign + vbm  # shift vbm with potential alignment

        occupied_midgap = {en: [] for en in np.arange(shifted_vbm, shifted_cbm + resolution, resolution)}
        occupation = {en: 0. for en in np.arange(shifted_vbm, shifted_cbm + resolution, resolution)}
        unoccupied_midgap = {en: [] for en in np.arange(shifted_vbm, shifted_cbm + resolution, resolution)}

        for spinset in eigenvalues.values():
            for kptset, weight in zip(spinset, kpoint_weights):
                for eig in kptset:  # eig[0] is eigenvalue and eig[1] is occupation
                    if (eig[1] and (eig[0] > shifted_cbm)):  # donor MB correction
                        bf_corr += weight * spinfctr * eig[1] * (eig[0] - shifted_cbm)  # "move the electrons down"
                        self.metadata["num_elec_cbm"] += weight * spinfctr * eig[1]
                    elif (eig[1] != 1.) and (eig[0] <= shifted_vbm):  # acceptor MB correction
                        bf_corr += weight * spinfctr * (1. - eig[1]) * (shifted_vbm - eig[0])  # "move the holes up"
                        self.metadata["num_hole_vbm"] += weight * spinfctr * (1. - eig[1])
                    elif (eig[0] > shifted_vbm) and (eig[0] < shifted_cbm):
                        for en in np.arange(shifted_vbm, shifted_cbm + resolution, resolution):
                            if (eig[0] < en + resolution) and (eig[0] > en):
                                if eig[1]:
                                    occupied_midgap[en].append(eig[0])
                                    occupation[en] += eig[1] * weight * spinfctr
                                else:
                                    unoccupied_midgap[en].append(eig[0])
                                continue

        bf_corr *= -1  # need to take negative of this shift for energetic correction

        # summarize defect level results
        self.metadata["total_occupation_defect_levels"] = 0.
        self.metadata["occupied_def_levels"] = []
        self.metadata["unoccupied_def_levels"] = []
        for en in occupied_midgap.keys():
            if len(occupied_midgap[en]):
                self.metadata["occupied_def_levels"].append([np.mean(occupied_midgap[en]), occupation[en]])
                self.metadata["total_occupation_defect_levels"] += occupation[en]
            elif len(unoccupied_midgap[en]):
                self.metadata["unoccupied_def_levels"].append(np.mean(unoccupied_midgap[en]))

        return bf_corr