Hi Sergey,
For method 2, I'm assuming you do a shift-invert around each energy.
If you need the LDOS for every site, then this is a good method.
If you want just a few elements, you can calculate the Green's function as
a continued fraction from the Lanczos tridiagonalization.
https://www.cond-mat.de/events/correl11/manuscript/Koch.pdf
This is an iterative method like KPM but I find that this method converges
faster sometimes.
If you're not happy with the energy resolution of KPM, you may have the
same issue here.
I sometimes need many thousands of iterations.
Here is my implementation.
I don't have any leads in my system.
def lanczos(H, i, N):
"""Calculate Lanczos tridiagonalization for the starting vector |i>.
Algorithm from
https://www.cond-mat.de/events/correl11/manuscript/Koch.pdf
Parameters
----------
H : Sparse array
Hamiltonian.
i : int
Starting element.
N : int
Number of iterations.
Returns
-------
a, b : array(N)
Tridiagonal elements.
"""
b = np.zeros(N)
a = np.zeros(N)
v = np.zeros(H.shape[0], H.dtype)
v[i] = 1
Hv = H.dot(v)
a[0] = Hv[i].real # a_0 = v.H.v
w = Hv - a[0]*v # w = H |v_0> - a_0 |v_0>
for n in range(1, N):
b[n] = np.linalg.norm(w)
w /= b[n] # w = |v_n>
v *= -b[n] # v = -b_n |v_{n-1}>
v, w = w, v # swap(v, w)
w += H.dot(v) # w = H |v_n> - b_n |v_{n-1}>
a[n] = np.vdot(v, w).real # a_n = <v_n|H|v_n> - b_n <v_n|v_{n-1}>
w -= a[n]*v # w = b_{n+1} |v_{n+1}>
return a, b
def greens_function(z, a, b):
"""Calculate the Green's function G = (z-H)^{-1} from Lanczos
tridiagonalization.
Algorithm from
https://www.cond-mat.de/events/correl11/manuscript/Koch.pdf
Parameters
----------
z : array(nz)
Energies to evaluate.
For retarded Green's function, z = ω+iδ.
a, b : array(N)
Tridiagonal elements.
Returns
-------
G : array(nz)
Green's function.
"""
G = z - a[-1]
for n in range(len(a)-2, -1, -1):
G = z - a[n] - b[n+1]**2/G
return 1 / G
a, b = lanczos(H, 0, 500)
omega = np.linspace(-4, 4, 401)
G = greens_function(omega+0.1j, a, b)
ldos = -G.imag/np.pi
Thanks,
Eric
On Thu, Aug 20, 2020 at 10:48 AM Sergey Slizovskiy <ser...@gmail.com> wrote:
> 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.
>
>
>
>
