Dear Colleagues, (sorry, this is a repost of the same question to a new thread)

I have a taks to compute local density of states at a single point for a large homogeneous sample with narrow STS tip potential and in magnetic field. I have tried several methods in Kwant to do this:

1) Naive diagonalization of H for system with no leads and direct calculation of LDOS. It works, but not for large systems, so, I see strong finite-size effects.

2) Sparse Lanczos method with several eigenvalues/functions extracted near each point in the grid of energy values.

3) KPM method. I am not happy with the energy resolution I can get with it.

4) Take a smaller main sample (with the size of the tip potential), but avoid the finite-size effects by attaching a lot of leads pointing in , e.g. 6 directions. Then, use the kwant ldos function to scan LDOS over a grid of energies. Use magnetic_gauge to implement B field in both sample and scattering region.

I see the last possibility as the quickest, but it gives me something that seems wrong to me.

Does kwant ldos function work correctly if all the leads are insulating for a given energy? Does it correctly catch the localized states in the scattering region?

Thank you for any feedback,

Sergey

Below is a copy of a reply by Pablo Piskunow

Dear Sergey,

In terms of scaling, the most convenient is the KPM method, since it scales linearly with the
inverse energy resolution (that is, the number of moments), and linearly with the system size.

Could you clarify what is the issue with the approach 3)?

I am not sure what is the energy resolution that you want to achieve. You can conveniently set
this value when creating a `kwant.kpm.SpectralDensity` instance via the arguments
`energy_resolution` or (exclusive) `num_moments`. The latter is related to the energy resoution
by `\delta \pi / N`, where `\delta` is the full width of the spectrum, and `N` the size of your system.
In the case of a 2D, that will be similar to `L^{^2}`, where `L` is the linear size.

Furthermore, you could progressively increase the energy resolution by adding moments:
```
spectrum = kwant.kpm.SpectralDensity(...)
spectrum.add_moments(...)
```

If you want to resolve individual energy levels then the approach 1) and 3) will have the same scaling,
however, the KPM method will not give you eigenvalues.

Finally, since you want to get the LDOS at a single site or a small set of sites near the STS tip, you should
use the `kwant.kpm.LocalVectors` vector factory, setting the `where` argument to the desired spot.

```
def center(site):
    return site.tag == (0, 0)
vector_factory = kwant.kpm.LocalVectors(syst, where=center)
spectrum = kwant.kpm.SpectralDensity(syst, num_moments=100, num_vectors=None, vector_factory=vector_factory)
```
Note that `num_vectors=None` will let the `kwant.kpm` module to figure out how many vectors does the vector factory produce.
Otherwise, it should be at most the *total number of orbitals* defined by the `where` function, that is, sites times orbitals per site.

As you can see, the `kwant.kpm` module is the most suited for this problem, specially when the sample is very large.
I hope this helps.

Best regards,
Pablo

n the meantime, you can create a `spectrum` instance, and evaluate the KPM expansion at any energy,

for example by

```

spectrum = kwant.kpm.SpectralDensity(...)

energy_array = np.linspace(0, 1, 200)

density_array = spectrum(energy_array)

```

Another note, if the Hamiltonian is huge, you will save some overhead by passing the bounds of the spectrum explicitly:

`bounds=(e_min, e_max)`. These values must be strictly below and above the edges of the spectrum, otherwise

the KPM expansion diverges.

ah, I forgot.

The automatic energy points (`spectrum.energies`) are fixed once you fix the bounds of the spectrum.

So if for a range of magnetic fields you pass the same bounds, then the energies will be at the same

points and so the densities will be easier to plot or compare.

The width of the spectrum might change for different values of the magnetic field; and even if not, there is

small randomness associated with finding the bounds, that you can get rid of by passing a seed to the

random number generator like `rng=0`. If you choose to provide the bounds explicitly, you can disregard

this random number generator stuff.