Coverage for kwant/kpm.py : 97%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*- # Copyright 2011-2016 Kwant authors. # # This file is part of Kwant. It is subject to the license terms in the file # LICENSE.rst found in the top-level directory of this distribution and at # http://kwant-project.org/license. A list of Kwant authors can be found in # the file AUTHORS.rst at the top-level directory of this distribution and at # http://kwant-project.org/authors.
r"""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 :math:`ρ_A(e)`, as a function of the energy :math:`e`, of some operator :math:`A` that acts on a kwant system or a Hamiltonian. In general
.. math:: ρ_A(e) = ρ(e) A(e),
where :math:`ρ(e) = \sum_{k=0}^{D-1} δ(E-E_k)` is the density of states, and :math:`A(e)` is the expectation value of :math:`A` for all the eigenstates with energy :math:`e`.
Parameters ---------- hamiltonian : `~kwant.system.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. If no operator is provided, the density of states is calculated with a faster algorithm. num_vectors : positive int, default: 10 Number of random vectors for the KPM method. 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 : function, optional The user defined function ``f(n)`` generates random vectors of length ``n`` that will be used in the algorithm. 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 by ``vector_factory``. 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.
Notes ----- When passing an ``operator`` defined in `~kwant.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.
.. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006) <https://arxiv.org/abs/cond-mat/0504627>`_. .. [2] `Phys. Rev. E 69, 057701 (2004) <https://arxiv.org/abs/cond-mat/0401202>`_
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 `~kwant.system.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 `~kwant.kpm.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
Attributes ---------- energies : array of floats Array of sampling points with length ``2 * num_moments`` in the range of the spectrum. densities : array of floats Spectral density of the ``operator`` evaluated at the energies. """
num_vectors=10, num_moments=None, energy_resolution=None, vector_factory=None, bounds=None, eps=0.05, rng=None):
"must be provided.")
# self.eps ensures that the rescaled Hamiltonian has a # spectrum strictly in the interval (-1,1).
# Normalize the format of 'ham' sparse=True) "nor a Kwant system.")
# Normalize 'operator' to a common format. else: 'attribute and is not callable.')
(lambda n: np.exp(2j * np.pi * rng.random_sample(n))) # store this vector for reproducibility # Hamiltonian rescaled as in Eq. (24) eps=self.eps, v0=self._v0, bounds=bounds)
self._vector_factory(self.hamiltonian.shape[0]))
moments, 2 * self.num_moments)
"""Return the spectral density evaluated at ``energy``.
Parameters ---------- energy : float or sequence of float, optional
Returns ------- float, if ``energy`` is float, array of float if ``energy`` is a sequence, a tuple (energies, densities) if ``energy`` was not provided.
Notes ----- If ``energy`` is not provided, then the densities are obtained by Fast Fourier Transform of the Chebyshev moments. """ else: * np.sqrt(1 + rescaled_energy))
np.cos(np.pi * m/(self.num_moments + 1)) + np.sin(np.pi * m/(self.num_moments + 1)) / np.tan(np.pi/(self.num_moments + 1)))
# transposes handle the case where operators have vector outputs
rescaled_energy, coef_cheb) / g_e).real
"""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. """ # This factor divides the sum to normalize the Gauss integral # and rescales the integral back with ``self._a`` to normal # scale. else: # The evaluation of the distribution function should be at # the energies without rescaling.
"""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'. """ "must be provided.")
.format(energy_resolution)) # factor of 1.6 comes from the fact that we use the # Jackson kernel when calculating the FFT, which has # maximal slope π/2. Rounding to 1.6 ensures that the # energy resolution is sufficient. 'than the requested resolution')
or num_moments != int(num_moments)):
self.num_vectors)
# recalculate quantities derived from the moments moments, 2 * self.num_moments)
"""Increase the number of random vectors.
Parameters ---------- num_vectors: positive int The number of random vectors to add. """ self._vector_factory(self.hamiltonian.shape[0])) self.num_vectors + num_vectors)
# recalculate quantities derived from the moments moments, 2 * self.num_moments)
# sum moments of all random vectors # divide by the number of random vectors # divide by scale factor to reflect the integral rescaling
"""Calculate the Chebyshev moments of an operator's spectral density.
The algorithm is based on the KPM method as depicted in `Rev. Mod. Phys., Vol. 78, No. 1 (2006) <https://arxiv.org/abs/cond-mat/0504627>`_.
Parameters ---------- n_moments : integer Number of Chebyshev moments. n_rand : integer Number of random vectors used for sampling. """
else: raise ValueError('Cannot decrease number of random vectors')
# nothing new to calculate return else: raise ValueError('Cannot decrease number of moments')
raise ValueError("Only 'num_moments' *or* 'num_vectors' " "may be updated at a time.")
else:
# Iteration over the moments # Two cases can occur, depicted in Eq. (28) and in Eq. (29), # respectively. # ---- # In the first case, self.operator is None and we can use # Eqs. (34) and (35) to obtain the density of states, with # two moments ``one_moment`` for every new alpha. # ---- # In the second case, the operator is not None and a matrix # multiplication should be used. - alpha) # Following Eqs. (34) and (35) - one_moment[0]) - one_moment[1]) # odd moment 2 * np.vdot(alpha_next, alpha_next) - one_moment[0]) else: - alpha)
# ### Auxiliary functions
"""Rescale a Hamiltonian and return a LinearOperator
Parameters ---------- hamiltonian : 2D array Hamiltonian of the system. eps : scalar Ensures that the bounds 'a' and 'b' are strict. v0 : random vector, or None Used as the initial residual vector for the algorithm that finds the lowest and highest eigenvalues. bounds : tuple, or None Boundaries of the spectrum. If not provided the maximum and minimum eigenvalues are calculated. """ # Relative tolerance to which to calculate eigenvalues. Because after # rescaling we will add eps / 2 to the spectral bounds, we don't need # to know the bounds more accurately than eps / 2.
else: return_eigenvectors=False, tol=tol, v0=v0)) return_eigenvectors=False, tol=tol, v0=v0))
'The Hamiltonian has a single eigenvalue, it is not possible to ' 'obtain a spectral density.')
"""This function takes the normalised moments and returns an array of points and an array of the evaluated function at those points. """ # extra_shape handles the case where operators have vector outputs
# Jackson kernel, as in Eq. (71), and kernel improved moments, # as in Eq. (81). np.sin(np.pi * m / (n_moments + 1)) / np.tan(np.pi / (n_moments + 1))) / (n_moments + 1)
# special points at the abscissas of Chebyshev integration # prefactor in Eq. (80)
# transposes handle the case where operators have vector outputs # The function evaluated in these special data points is the FFT of # the moments as in Eq. (83).
# Element-wise division of moments_FFT over gk, as in Eq. (83).
# Reverse energies and densities to set ascending order. |