[atomate2] parsing of f orbitals in dos

jshen · May 18, 2022, 10:15pm

It looks like a call the obtain the complete DOS:

materialsproject/atomate2/blob/c3c5a0b7bef64db4e9c0ba0dc7ece0e370764df7/src/atomate2/vasp/schemas/calculation.py#L854

      
        
                OrbitalType.f,
            ]:
                orb_name = orb_type.name
                if (
                    (el.block == "s" and orb_name in ["p", "d", "f"])
                    or (el.block == "p" and orb_name in ["d", "f"])
                    or (el.block == "d" and orb_name == "f")
                ):
                    continue
                dosprops = {
                    "filling": complete_dos.get_band_filling(band=orb_type, elements=[el]),
                    "center": complete_dos.get_band_center(band=orb_type, elements=[el]),
                    "bandwidth": complete_dos.get_band_width(band=orb_type, elements=[el]),
                    "skewness": complete_dos.get_band_skewness(
                        band=orb_type, elements=[el]
                    ),
                    "kurtosis": complete_dos.get_band_kurtosis(
                        band=orb_type, elements=[el]
                    ),
                    "upper_edge": complete_dos.get_upper_band_edge(
                        band=orb_type, elements=[el]

Fails for f-orbitals since there is a call to get_element_spd_dos in get_band_filling which I think is incompatible with f-orbitals.

github.com

materialsproject/pymatgen/blob/ec9464528fbe3b95d86d30c4361722746196b854/pymatgen/electronic_structure/dos.py#L843

      
        
            
            
Returns:
                band filling in eV, often denoted f_d for the d-band
            """
            # Get the projected DOS
            if elements and sites:
                raise ValueError("Both element and site cannot be specified.")
            
            
if elements:
                for i, el in enumerate(elements):
                    spd_dos = self.get_element_spd_dos(el)[band]
                    if i == 0:
                        densities = spd_dos.densities
                    else:
                        densities = add_densities(densities, spd_dos.densities)
                dos = Dos(self.efermi, self.energies, densities)
            elif sites:
                for i, site in enumerate(sites):
                    spd_dos = self.get_site_spd_dos(site)[band]
                    if i == 0:
                        densities = spd_dos.densities

alex · May 23, 2022, 6:19pm

Thanks for raising this. @arosen provided a fix which is now available in the master branch on GitHub.

jshen · May 23, 2022, 7:24pm

Thanks @alex and @arosen