def save_dos(calc):
    k, spin = 0, 0
    sigma   = 0.1
    e_min   = -10.0
    e_max   = +10.0
    e_step  = 0.01

    def gauss_ij(energy_i, energy_j, sigma):
        denergy_ij = energy_i[:, np.newaxis] - energy_j[np.newaxis, :]
        norm = 1.0 / (sigma * np.sqrt(2 * np.pi))
        return norm * np.exp(-0.5 * denergy_ij**2 / sigma**2)

    energy_g = np.arange(e_min, e_max, e_step) # Energy grid
    eig_n    = calc.get_eigenvalues(k, spin) - calc.get_fermi_level()
    dos      = np.sum(gauss_ij(energy_g, eig_n, sigma), axis = 1)

    parprint('Saving DOS {}'.format(data_obj.dos.filename))
    if mpi.world.rank == 0:
        np.save(data_obj.dos.filename, (energy_g, dos))