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.