Implement the KPM method for calculating density of states for finite systems
The idea is to implement the algorithms explored here in Kwant.
On a very high-level the idea would be to write a solver that, given a system and an operator, calculates the spectral density of the operator (not sure whether
this is the correct term, what we are trying to describe is the quantity denoted a(E) rho(E)
in the linked repository).
The KPM allows us to calculate not just this quantity at a single energy, but actually over the whole spectrum.
The idea would be to have something like:
# in module "kpm"
def spectral_density(syst, *, args=(), operator=None, num_rand_vecs=10, num_moments=10, energy_bounds=()):
pass # actually calculate the moments here
return SpectralDensity(all_the_moments)
class SpectralDensity:
def __init__(self, all_the_moments):
pass # store all the moments
def __call__(self, energy):
pass # evaluate the spectral density @ energy
def average(self, distribution_function):
pass # integrate with distribution function