Projected density of states for a list of sites

[wladerer]

Hi everyone,

I am wondering how use the CompleteDos.get_site_dos() method. I am succesful in getting the complete dos for each site, but then adding them together is quite a pain, especially if they are decomposed by momenta.

Is there a way to provide a range of sites and get their summed DOS?

I am currently doing this

bare_dos_calcs = [ Vasprun(f"{bare_parent}/{i}/dos/vasprun.xml", parse_potcar_file=False) for i in range(2,19) ]
total_pdos_list = []
pdos_p_list = []
pdos_d_list = []

for dos_calc in bare_dos_calcs:
    dos = dos_calc.complete_dos
    structure = dos_calc.final_structure
    sites = sorted(structure, key=lambda x: x.frac_coords[2], reverse=True)
    efermi = dos.efermi


    pdos_list = [ dos.get_site_spd_dos(site) for site in sites[:16] ]

    # energy is same window
    energy_range = pdos_list[0][OrbitalType.s].energies

    # A dict of spin densities, e.g., {Spin.up: […], Spin.down: […]}. is pdos[OrbitalType.s].densities
    pdos_s_densities = sum([ pdos[OrbitalType.s].densities[Spin.up] for pdos in pdos_list ])
    pdos_p_densities = sum([ pdos[OrbitalType.p].densities[Spin.up] for pdos in pdos_list ])
    pdos_d_densities = sum([ pdos[OrbitalType.d].densities[Spin.up] for pdos in pdos_list ])
    pdos_total_densities = pdos_s_densities + pdos_p_densities + pdos_d_densities

    # create Dos object class Dos(efermi: float, energies: ArrayLike, densities: Mapping[Spin, ArrayLike], norm_vol: float | None = None)
    pdos_s = Dos(efermi, energy_range, {Spin.up: pdos_s_densities})
    pdos_p = Dos(efermi, energy_range, {Spin.up: pdos_p_densities})
    pdos_d = Dos(efermi, energy_range, {Spin.up: pdos_d_densities})
    pdos_total = Dos(efermi, energy_range, {Spin.up: pdos_total_densities})

    pdos_p_list.append(pdos_p)
    total_pdos_list.append(pdos_total)

Where the above code iterates over all of the sites I am interested in (I am sorting by height in this case). It then sums the densities of each l-decomposed DOS and adds them together, to reconstruct the density of states of the sites I am interested in. Then, I create a Dos object from the energies, densities, and fermi energy.

Is there an easier way of doing this?

This may also be related to the following open discussion:
#3742 (comment)

[fraricci]

Hi there!

I think currently there is not a function to sum doses over specific sites (contrarily to the projected bands).
I came up with this trick that allows to make use of the get_element_dos() and get_element_spd_dos().
The idea is to convert all the sites you are interested into fake sites that have the same element.

Here is the code that should do that:

from pymatgen.core.sites import PeriodicSite
vrun = Vasprun('./gaas/vasprun.xml',parse_projected_eigen=True)
cdos = vrun.complete_dos
sites_selected = cdos.structure.sites[:16]

for site in sites_selected:
    fake_site = PeriodicSite("H",site.frac_coords,site.lattice)
    cdos.pdos[fake_site] = cdos.pdos.pop(site) #or add this additional site to the pdos to be less destructive

tot_elems_dos = get_element_dos() # get the tot dos for "H"
elems_spd_dos = get_element_spd_dos("H")

Does it make sense? I haven't fully tested. Let me know if it works.

[wladerer]

Hahah this is a fun workaround, thank you! I think this should work for my purposes.