kwant.kpm – Kernel Polynomial Method

4.3. kwant.kpm – Kernel Polynomial Method#

class kwant.kpm.Correlator(hamiltonian, operator1=None, operator2=None, **kwargs)[source]#

Bases: object

Calculates the response of the correlation between two operators.

The response tensor \(χ_{α β}\) of an operator \(O_α\) to a perturbation in an operator \(O_β\), is defined here based on [3], and [4], and takes the form

\[χ_{α β}(µ, T) = \int_{-\infty}^{\infty}{\mathrm{d}E} f(µ-E, T) \left({O_α ρ(E) O_β \frac{\mathrm{d}G^{+}}{\mathrm{d}E}} - {O_α \frac{\mathrm{d}G^{-}}{\mathrm{d}E} O_β ρ(E)}\right)\]

Internally, the correlation is approximated with a two dimensional KPM expansion,

\[χ_{α β}(µ, T) = \int_{-1}^1{\mathrm{d}E} \frac{f(µ-E,T)}{(1-E^2)^2} \sum_{m,n}Γ_{n m}(E)µ_{n m}^{α β},\]

with coefficients

\[ \begin{align}\begin{aligned}Γ_{m n}(E) = (E - i n \sqrt{1 - E^2}) e^{i n \arccos(E)} T_m(E)\\+ (E + i m \sqrt{1 - E^2}) e^{-i m \arccos(E)} T_n(E),\end{aligned}\end{align} \]

and moments matrix \(µ_{n m}^{α β} = \mathrm{Tr}(O_α T_m(H) O_β T_n(H))\).

The trace is calculated creating two instances of SpectralDensity, and saving the vectors \(Ψ_{n r} = O_β T_n(H)\rvert r\rangle\), and \(Ω_{m r} = T_m(H) O_α\rvert r\rangle\) , where \(\rvert r\rangle\) is a vector provided by the vector_factory. The moments matrix is formed with the product \(µ_{m n} = \langle Ω_{m r} \rvert Ψ_{n r}\rangle\) for every \(\rvert r\rangle\).

Parameters:

  • hamiltonian (FiniteSystem or matrix Hamiltonian) – If a system is passed, it should contain no leads.

  • operator1 (operators, dense matrix, or sparse matrix, optional) – Operators to be passed to two different instances of SpectralDensity.

  • operator2 (operators, dense matrix, or sparse matrix, optional) – Operators to be passed to two different instances of SpectralDensity.

  • **kwargs (dict) – Keyword arguments to pass to SpectralDensity.

Notes

The operator1 must act to the right as \(O_α\rvert r\rangle\).

__call__(mu=0, temperature=0)[source]#

Returns the linear response \(χ_{α β}(µ, T)\)

Parameters:

  • mu (float) – Chemical potential defined in the same units of energy as the Hamiltonian.

  • temperature (float) – Temperature in units of energy, the same as defined in the Hamiltonian.

add_moments(num_moments=None, *, energy_resolution=None)[source]#

Increase the number of Chebyshev moments

Parameters:

  • num_moments (positive int, optional) – The number of Chebyshev moments to add. Mutually exclusive with ‘energy_resolution’.

  • energy_resolution (positive float, optional) – Features wider than this resolution are visible in the spectral density. Mutually exclusive with num_moments.

add_vectors(num_vectors=None)[source]#

Increase the number of vectors

Parameters:

num_vectors (positive int, optional) – The number of vectors to add.

class kwant.kpm.LocalVectors(syst, where=None, *args)[source]#

Bases: object

Generates a local vector iterator for the sites in ‘where’.

Parameters:

  • syst (FiniteSystem or matrix Hamiltonian) – If a system is passed, it should contain no leads.

  • where (sequence of int or Site, or callable, optional) – Spatial range of the vectors produced. If syst is a FiniteSystem, where behaves as in Density. If syst is a matrix, where must be a list of integers with the indices where column vectors are nonzero.

kwant.kpm.RandomVectors(syst, where=None, rng=None)[source]#

Returns a random phase vector iterator for the sites in ‘where’.

Parameters:

  • syst (FiniteSystem or matrix Hamiltonian) – If a system is passed, it should contain no leads.

  • where (sequence of int or Site, or callable, optional) – Spatial range of the vectors produced. If syst is a FiniteSystem, where behaves as in Density. If syst is a matrix, where must be a list of integers with the indices where column vectors are nonzero.

class kwant.kpm.SpectralDensity(hamiltonian, params=None, operator=None, num_vectors=10, num_moments=None, energy_resolution=None, vector_factory=None, bounds=None, eps=0.05, rng=None, kernel=None, mean=True, accumulate_vectors=True)[source]#

Bases: object

Calculate the spectral density of an operator.

This class makes use of the kernel polynomial method (KPM), presented in [1], to obtain the spectral density \(ρ_A(e)\), as a function of the energy \(e\), of some operator \(A\) that acts on a kwant system or a Hamiltonian. In general

\[ρ_A(E) = ρ(E) A(E),\]

where \(ρ(E) = \sum_{k=0}^{D-1} δ(E-E_k)\) is the density of states, and \(A(E)\) is the expectation value of \(A\) for all the eigenstates with energy \(E\).

Parameters:

  • hamiltonian (FiniteSystem or matrix Hamiltonian) – If a system is passed, it should contain no leads.

  • params (dict, optional) – Additional parameters to pass to the Hamiltonian and operator.

  • operator (operator, dense matrix, or sparse matrix, optional) – Operator for which the spectral density will be evaluated. If it is callable, the densities at each energy will have the dimension of the result of operator(bra, ket). If it has a dot method, such as numpy.ndarray and scipy.sparse.matrices, the densities will be scalars.

  • num_vectors (positive int, or None, default: 10) – Number of vectors used in the KPM expansion. If None, the number of vectors used equals the length of the ‘vector_factory’.

  • num_moments (positive int, default: 100) – Number of moments, order of the KPM expansion. Mutually exclusive with energy_resolution.

  • energy_resolution (positive float, optional) – The resolution in energy of the KPM approximation to the spectral density. Mutually exclusive with num_moments.

  • vector_factory (iterable, optional) – If provided, it should contain (or yield) vectors of the size of the system. If not provided, random phase vectors are used. The default random vectors are optimal for most cases, see the discussions in [1] and [2].

  • bounds (pair of floats, optional) – Lower and upper bounds for the eigenvalue spectrum of the system. If not provided, they are computed.

  • eps (positive float, default: 0.05) – Parameter to ensure that the rescaled spectrum lies in the interval (-1, 1); required for stability.

  • rng (seed, or random number generator, optional) – Random number generator used for the calculation of the spectral bounds, and to generate random vectors (if vector_factory is not provided). If not provided, numpy’s rng will be used; if it is an Integer, it will be used to seed numpy’s rng, and if it is a random number generator, this is the one used.

  • kernel (callable, optional) – Callable that takes moments and returns stabilized moments. By default, the jackson_kernel is used. The Lorentz kernel is also accesible by passing lorentz_kernel.

  • mean (bool, default: True) – If True, return the mean spectral density for the vectors used, otherwise return an array of densities for each vector.

  • accumulate_vectors (bool, default: True) – Whether to save or discard each vector produced by the vector factory. If it is set to False, it is not possible to increase the number of moments, but less memory is used.

Notes

When passing an operator defined in operator, the result of operator(bra, ket) depends on the attribute sum of such operator. If sum=True, densities will be scalars, that is, total densities of the system. If sum=False the densities will be arrays of the length of the system, that is, local densities.

Examples

In the following example, we will obtain the density of states of a graphene sheet, defined as a honeycomb lattice with first nearest neighbors coupling.

We start by importing kwant and defining a FiniteSystem,

>>> import kwant
...
>>> def circle(pos):
...     x, y = pos
...     return x**2 + y**2 < 100
...
>>> lat = kwant.lattice.honeycomb()
>>> syst = kwant.Builder()
>>> syst[lat.shape(circle, (0, 0))] = 0
>>> syst[lat.neighbors()] = -1

and after finalizing the system, create an instance of SpectralDensity

>>> fsyst = syst.finalized()
>>> rho = kwant.kpm.SpectralDensity(fsyst)

The energies and densities can be accessed with

>>> energies, densities = rho()

or

>>> energies, densities = rho.energies, rho.densities
energies[source]#

Array of sampling points with length 2 * num_moments in the range of the spectrum.

Type:

array of floats

densities[source]#

Spectral density of the operator evaluated at the energies.

Type:

array of floats

__call__(energy=None)[source]#

Return the spectral density evaluated at energy.

Parameters:

energy (float or sequence of floats, optional)

Returns:

  • energies (array of floats) – Drawn from the nodes of the highest Chebyshev polynomial. Not returned if ‘energy’ was not provided

  • densities (float or array of floats) – single float if the energy parameter is a single float, else an array of float.

Notes

If energy is not provided, then the densities are obtained by Fast Fourier Transform of the Chebyshev moments.

add_moments(num_moments=None, *, energy_resolution=None)[source]#

Increase the number of Chebyshev moments.

Parameters:

  • num_moments (positive int) – The number of Chebyshev moments to add. Mutually exclusive with energy_resolution.

  • energy_resolution (positive float, optional) – Features wider than this resolution are visible in the spectral density. Mutually exclusive with num_moments.

add_vectors(num_vectors=None)[source]#

Increase the number of vectors

Parameters:

num_vectors (positive int, optional) – The number of vectors to add.

integrate(distribution_function=None)[source]#

Returns the total spectral density.

Returns the integral over the whole spectrum with an optional distribution function. distribution_function should be able to take arrays as input. Defined using Gauss-Chebyshev integration.

kwant.kpm.conductivity(hamiltonian, alpha='x', beta='x', positions=None, **kwargs)[source]#

Returns a callable object to obtain the elements of the conductivity tensor using the Kubo-Bastin approach.

A Correlator instance is created to obtain the correlation between two components of the current operator

\[σ_{α β}(µ, T) = \frac{1}{V} \int_{-\infty}^{\infty}{\mathrm{d}E} f(µ-E, T) \left({j_α ρ(E) j_β \frac{\mathrm{d}G^{+}}{\mathrm{d}E}} - {j_α \frac{\mathrm{d}G^{-}}{\mathrm{d}E} j_β ρ(E)}\right),\]

where \(V\) is the volume where the conductivity is sampled. In this implementation it is assumed that the vectors are normalized and \(V=1\), otherwise the result of this calculation must be normalized with the corresponding volume.

The equations used here are based on [5] and [6]

Parameters:

  • hamiltonian (FiniteSystem or matrix Hamiltonian) – If a system is passed, it should contain no leads.

  • alpha (str, or operators) – If hamiltonian is a kwant system, or if the positions are provided, alpha and beta can be the directions of the velocities as strings {‘x’, ‘y’, ‘z’}. Otherwise alpha and beta should be the proper velocity operators, which can be members of operator or matrices.

  • beta (str, or operators) – If hamiltonian is a kwant system, or if the positions are provided, alpha and beta can be the directions of the velocities as strings {‘x’, ‘y’, ‘z’}. Otherwise alpha and beta should be the proper velocity operators, which can be members of operator or matrices.

  • positions (array of float, optioinal) – If hamiltonian is a matrix, the velocities can be calculated internally by passing the positions of each orbital in the system when alpha or beta are one of the directions {‘x’, ‘y’, ‘z’}.

  • **kwargs (dict) – Keyword arguments to pass to Correlator.

Examples

We will obtain the conductivity of the Haldane model, defined as a honeycomb lattice with first nearest neighbors coupling, and imaginary second nearest neighbors coupling.

We start by importing kwant and defining a FiniteSystem,

>>> import kwant
...
>>> def circle(pos):
...     x, y = pos
...     return x**2 + y**2 < 100
...
>>> lat = kwant.lattice.honeycomb()
>>> syst = kwant.Builder()
>>> syst[lat.shape(circle, (0, 0))] = 0
>>> syst[lat.neighbors()] = -1
>>> syst[lat.a.neighbors()] = -0.5j
>>> syst[lat.b.neighbors()] = 0.5j
>>> fsyst = syst.finalized()

Now we can call conductivity to calculate the transverse conductivity at chemical potential 0 and temperature 0.01.

>>> cond = kwant.kpm.conductivity(fsyst, alpha='x', beta='y')
>>> cond(mu=0, temperature=0.01)
kwant.kpm.fermi_distribution(energy, mu, temperature)[source]#

Returns the Fermi distribution f(e, µ, T) evaluated at ‘e’.

Parameters:

  • energy (float or sequence of floats) – Energy array where the Fermi distribution is evaluated.

  • mu (float) – Chemical potential defined in the same units of energy as the Hamiltonian.

  • temperature (float) – Temperature in units of energy, the same as defined in the Hamiltonian.

kwant.kpm.jackson_kernel(moments)[source]#

Convolutes moments with the Jackson kernel.

Taken from Eq. (71) of Rev. Mod. Phys., Vol. 78, No. 1 (2006).

kwant.kpm.lorentz_kernel(moments, l=4)[source]#

Convolutes moments with the Lorentz kernel.

Taken from Eq. (71) of Rev. Mod. Phys., Vol. 78, No. 1 (2006).

The additional parameter l controls the decay of the kernel.

Contents