Coverage for kwant/kpm.py : 90%
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
1# -*- coding: utf-8 -*-
2# Copyright 2011-2019 Kwant authors.
3#
4# This file is part of Kwant. It is subject to the license terms in the file
5# LICENSE.rst found in the top-level directory of this distribution and at
6# https://kwant-project.org/license. A list of Kwant authors can be found in
7# the file AUTHORS.rst at the top-level directory of this distribution and at
8# https://kwant-project.org/authors.
9import math
10from operator import add
11from collections.abc import Iterable
12from functools import reduce
13import numpy as np
14from numpy.polynomial.chebyshev import chebval
15from scipy.sparse import coo_matrix, csr_matrix
16from scipy.integrate import simps
17from scipy.sparse.linalg import eigsh, LinearOperator
18import scipy.fftpack as fft
20from . import system
21from ._common import ensure_rng
22from .operator import (_LocalOperator, _get_tot_norbs, _get_all_orbs,
23 _normalize_site_where)
24from .graph.defs import gint_dtype
26__all__ = ['SpectralDensity', 'Correlator', 'conductivity',
27 'RandomVectors', 'LocalVectors', 'jackson_kernel', 'lorentz_kernel',
28 'fermi_distribution']
30SAMPLING = 2 # number of sampling points to number of moments ratio
33class SpectralDensity:
34 r"""Calculate the spectral density of an operator.
36 This class makes use of the kernel polynomial
37 method (KPM), presented in [1]_, to obtain the spectral density
38 :math:`ρ_A(e)`, as a function of the energy :math:`e`, of some
39 operator :math:`A` that acts on a kwant system or a Hamiltonian.
40 In general
42 .. math::
43 ρ_A(E) = ρ(E) A(E),
45 where :math:`ρ(E) = \sum_{k=0}^{D-1} δ(E-E_k)` is the density of
46 states, and :math:`A(E)` is the expectation value of :math:`A` for
47 all the eigenstates with energy :math:`E`.
49 Parameters
50 ----------
51 hamiltonian : `~kwant.system.FiniteSystem` or matrix Hamiltonian
52 If a system is passed, it should contain no leads.
53 params : dict, optional
54 Additional parameters to pass to the Hamiltonian and operator.
55 operator : operator, dense matrix, or sparse matrix, optional
56 Operator for which the spectral density will be evaluated. If
57 it is callable, the ``densities`` at each energy will have the
58 dimension of the result of ``operator(bra, ket)``. If it has a
59 ``dot`` method, such as ``numpy.ndarray`` and
60 ``scipy.sparse.matrices``, the densities will be scalars.
61 num_vectors : positive int, or None, default: 10
62 Number of vectors used in the KPM expansion. If ``None``, the
63 number of vectors used equals the length of the 'vector_factory'.
64 num_moments : positive int, default: 100
65 Number of moments, order of the KPM expansion. Mutually exclusive
66 with ``energy_resolution``.
67 energy_resolution : positive float, optional
68 The resolution in energy of the KPM approximation to the spectral
69 density. Mutually exclusive with ``num_moments``.
70 vector_factory : iterable, optional
71 If provided, it should contain (or yield) vectors of the size of
72 the system. If not provided, random phase vectors are used.
73 The default random vectors are optimal for most cases, see the
74 discussions in [1]_ and [2]_.
75 bounds : pair of floats, optional
76 Lower and upper bounds for the eigenvalue spectrum of the system.
77 If not provided, they are computed.
78 eps : positive float, default: 0.05
79 Parameter to ensure that the rescaled spectrum lies in the
80 interval ``(-1, 1)``; required for stability.
81 rng : seed, or random number generator, optional
82 Random number generator used for the calculation of the spectral
83 bounds, and to generate random vectors (if ``vector_factory`` is
84 not provided). If not provided, numpy's rng will be used; if it
85 is an Integer, it will be used to seed numpy's rng, and if it is
86 a random number generator, this is the one used.
87 kernel : callable, optional
88 Callable that takes moments and returns stabilized moments.
89 By default, the `~kwant.kpm.jackson_kernel` is used.
90 The Lorentz kernel is also accesible by passing
91 `~kwant.kpm.lorentz_kernel`.
92 mean : bool, default: ``True``
93 If ``True``, return the mean spectral density for the vectors
94 used, otherwise return an array of densities for each vector.
95 accumulate_vectors : bool, default: ``True``
96 Whether to save or discard each vector produced by the vector
97 factory. If it is set to ``False``, it is not possible to
98 increase the number of moments, but less memory is used.
100 Notes
101 -----
102 When passing an operator defined in `~kwant.operator`, the
103 result of ``operator(bra, ket)`` depends on the attribute ``sum``
104 of such operator. If ``sum=True``, densities will be scalars, that
105 is, total densities of the system. If ``sum=False`` the densities
106 will be arrays of the length of the system, that is, local
107 densities.
109 .. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
110 <https://arxiv.org/abs/cond-mat/0504627>`_.
111 .. [2] `Phys. Rev. E 69, 057701 (2004)
112 <https://arxiv.org/abs/cond-mat/0401202>`_
114 Examples
115 --------
116 In the following example, we will obtain the density of states of a
117 graphene sheet, defined as a honeycomb lattice with first nearest
118 neighbors coupling.
120 We start by importing kwant and defining a
121 `~kwant.system.FiniteSystem`,
123 >>> import kwant
124 ...
125 >>> def circle(pos):
126 ... x, y = pos
127 ... return x**2 + y**2 < 100
128 ...
129 >>> lat = kwant.lattice.honeycomb()
130 >>> syst = kwant.Builder()
131 >>> syst[lat.shape(circle, (0, 0))] = 0
132 >>> syst[lat.neighbors()] = -1
134 and after finalizing the system, create an instance of
135 `~kwant.kpm.SpectralDensity`
137 >>> fsyst = syst.finalized()
138 >>> rho = kwant.kpm.SpectralDensity(fsyst)
140 The ``energies`` and ``densities`` can be accessed with
142 >>> energies, densities = rho()
144 or
146 >>> energies, densities = rho.energies, rho.densities
148 Attributes
149 ----------
150 energies : array of floats
151 Array of sampling points with length ``2 * num_moments`` in
152 the range of the spectrum.
153 densities : array of floats
154 Spectral density of the ``operator`` evaluated at the energies.
155 """
157 def __init__(self, hamiltonian, params=None, operator=None,
158 num_vectors=10, num_moments=None, energy_resolution=None,
159 vector_factory=None, bounds=None, eps=0.05, rng=None,
160 kernel=None, mean=True, accumulate_vectors=True):
162 if num_moments and energy_resolution:
163 raise TypeError("either 'num_moments' or 'energy_resolution' "
164 "must be provided.")
166 # self.eps ensures that the rescaled Hamiltonian has a
167 # spectrum strictly in the interval (-1,1).
168 self.eps = eps
170 # Normalize the format of 'ham'
171 if isinstance(hamiltonian, system.System):
172 hamiltonian = hamiltonian.hamiltonian_submatrix(params=params,
173 sparse=True)
174 try:
175 hamiltonian = csr_matrix(hamiltonian)
176 except Exception:
177 raise ValueError("'hamiltonian' is neither a matrix "
178 "nor a Kwant system.")
180 # Normalize 'operator' to a common format.
181 if operator is None:
182 self.operator = None
183 elif isinstance(operator, _LocalOperator):
184 self.operator = operator.bind(params=params)
185 elif callable(operator):
186 self.operator = operator
187 elif hasattr(operator, 'dot'):
188 operator = csr_matrix(operator)
189 self.operator = lambda bra, ket: np.vdot(bra, operator.dot(ket))
190 else:
191 raise ValueError("Parameter 'operator' has no '.dot' "
192 "attribute and is not callable.")
194 self.mean = mean
195 rng = ensure_rng(rng)
196 # store this vector for reproducibility
197 self._v0 = np.exp(2j * np.pi * rng.random_sample(hamiltonian.shape[0]))
199 if eps <= 0:
200 raise ValueError("'eps' must be positive")
202 # Hamiltonian rescaled as in Eq. (24)
203 self.hamiltonian, (self._a, self._b) = _rescale(hamiltonian,
204 eps=self.eps,
205 v0=self._v0,
206 bounds=bounds)
207 self.bounds = (self._b - self._a, self._b + self._a)
209 if energy_resolution:
210 num_moments = math.ceil((1.6 * self._a) / energy_resolution)
211 elif num_moments is None:
212 num_moments = 100
214 if num_moments <= 0 or num_moments != int(num_moments):
215 raise ValueError("'num_moments' must be a positive integer")
217 if vector_factory is None:
218 self._vector_factory = _VectorFactory(
219 RandomVectors(hamiltonian, rng=rng),
220 num_vectors=num_vectors,
221 accumulate=accumulate_vectors)
222 else:
223 if not isinstance(vector_factory, Iterable): 223 ↛ 224line 223 didn't jump to line 224, because the condition on line 223 was never true
224 raise TypeError('vector_factory must be iterable')
225 try:
226 len(vector_factory)
227 except TypeError:
228 if num_vectors is None: 228 ↛ 229line 228 didn't jump to line 229, because the condition on line 228 was never true
229 raise ValueError('num_vectors must be provided if'
230 'vector_factory has no length.')
231 self._vector_factory = _VectorFactory(
232 vector_factory,
233 num_vectors=num_vectors,
234 accumulate=accumulate_vectors)
235 num_vectors = self._vector_factory.num_vectors
237 self._last_two_alphas = []
238 self._moments_list = []
240 self.num_moments = num_moments
241 self._update_moments_list(self.num_moments, num_vectors)
243 # set kernel before calling moments
244 self.kernel = kernel if kernel is not None else jackson_kernel
245 moments = self._moments()
246 self.densities, self._gammas = _calc_fft_moments(moments)
248 @property
249 def energies(self):
250 return (self._a * _chebyshev_nodes(SAMPLING * self.num_moments)
251 + self._b)
253 @property
254 def num_vectors(self):
255 return len(self._moments_list)
257 def __call__(self, energy=None):
258 """Return the spectral density evaluated at ``energy``.
260 Parameters
261 ----------
262 energy : float or sequence of floats, optional
264 Returns
265 -------
266 energies : array of floats
267 Drawn from the nodes of the highest Chebyshev polynomial.
268 Not returned if 'energy' was not provided
269 densities : float or array of floats
270 single ``float`` if the ``energy`` parameter is a single
271 ``float``, else an array of ``float``.
273 Notes
274 -----
275 If ``energy`` is not provided, then the densities are obtained
276 by Fast Fourier Transform of the Chebyshev moments.
277 """
278 if energy is None:
279 return self.energies, self.densities
280 else:
281 energy = np.asarray(energy)
282 e = (energy - self._b) / self._a
283 g_e = (np.pi * np.sqrt(1 - e) * np.sqrt(1 + e))
285 moments = self._moments()
286 # factor 2 comes from the norm of the Chebyshev polynomials
287 moments[1:] = 2 * moments[1:]
289 return np.transpose(chebval(e, moments) / g_e)
291 def integrate(self, distribution_function=None):
292 """Returns the total spectral density.
294 Returns the integral over the whole spectrum with an optional
295 distribution function. ``distribution_function`` should be able
296 to take arrays as input. Defined using Gauss-Chebyshev
297 integration.
298 """
299 # This factor divides the sum to normalize the Gauss integral
300 # and rescales the integral back with ``self._a`` to normal
301 # scale.
302 factor = self._a / (2 * self.num_moments)
303 if distribution_function is None:
304 rho = self._gammas
305 else:
306 # The evaluation of the distribution function should be at
307 # the energies without rescaling.
308 distribution_array = distribution_function(self.energies)
309 rho = np.transpose(self._gammas.transpose() * distribution_array)
310 return factor * np.sum(rho, axis=0)
312 def add_moments(self, num_moments=None, *, energy_resolution=None):
313 """Increase the number of Chebyshev moments.
315 Parameters
316 ----------
317 num_moments: positive int
318 The number of Chebyshev moments to add. Mutually
319 exclusive with ``energy_resolution``.
320 energy_resolution: positive float, optional
321 Features wider than this resolution are visible
322 in the spectral density. Mutually exclusive with
323 ``num_moments``.
324 """
325 if not ((num_moments is None) ^ (energy_resolution is None)):
326 raise TypeError("either 'num_moments' or 'energy_resolution' "
327 "must be provided.")
329 if energy_resolution:
330 if energy_resolution <= 0:
331 raise ValueError("'energy_resolution' must be positive")
332 # factor of 1.6 comes from the fact that we use the
333 # Jackson kernel when calculating the FFT, which has
334 # maximal slope π/2. Rounding to 1.6 ensures that the
335 # energy resolution is sufficient.
336 present_resolution = self._a * 1.6 / self.num_moments
337 if present_resolution < energy_resolution:
338 raise ValueError('Energy resolution is already smaller '
339 'than the requested resolution')
340 num_moments = math.ceil((1.6 * self._a) / energy_resolution)
342 if (num_moments is None or num_moments <= 0
343 or num_moments != int(num_moments)):
344 raise ValueError("'num_moments' must be a positive integer")
346 self._update_moments_list(self.num_moments + num_moments,
347 self.num_vectors)
348 self.num_moments += num_moments
350 # recalculate quantities derived from the moments
351 moments = self._moments()
352 self.densities, self._gammas = _calc_fft_moments(moments)
354 def add_vectors(self, num_vectors=None):
355 """Increase the number of vectors
357 Parameters
358 ----------
359 num_vectors: positive int, optional
360 The number of vectors to add.
361 """
362 self._vector_factory.add_vectors(num_vectors)
363 num_vectors = self._vector_factory.num_vectors - self.num_vectors
365 self._update_moments_list(self.num_moments,
366 self.num_vectors + num_vectors)
368 # recalculate quantities derived from the moments
369 moments = self._moments()
370 self.densities, self._gammas = _calc_fft_moments(moments)
372 def _moments(self):
373 moments = np.real_if_close(self._moments_list)
374 # put moments in the first axis, to return an array of densities
375 moments = np.swapaxes(moments, 0, 1)
376 if self.mean:
377 moments = np.mean(moments, axis=1)
378 # divide by scale factor to reflect the integral rescaling
379 moments /= self._a
380 # stabilized moments with a kernel
381 moments = self.kernel(moments)
382 return moments
384 def _update_moments_list(self, n_moments, num_vectors):
385 """Calculate the Chebyshev moments of an operator's spectral
386 density.
388 The algorithm is based on the KPM method as depicted in `Rev.
389 Mod. Phys., Vol. 78, No. 1 (2006)
390 <https://arxiv.org/abs/cond-mat/0504627>`_.
392 Parameters
393 ----------
394 n_moments : integer
395 Number of Chebyshev moments.
396 num_vectors : integer
397 Number of vectors used for sampling.
398 """
400 if self.num_vectors == num_vectors:
401 r_start = 0
402 new_vectors = 0
403 elif self.num_vectors < num_vectors: 403 ↛ 407line 403 didn't jump to line 407, because the condition on line 403 was never false
404 r_start = self.num_vectors
405 new_vectors = num_vectors - self.num_vectors
406 else:
407 raise ValueError('Cannot decrease number of vectors')
408 self._moments_list.extend([0.] * new_vectors)
409 self._last_two_alphas.extend([0.] * new_vectors)
411 if n_moments == self.num_moments:
412 m_start = 2
413 new_moments = 0
414 if new_vectors == 0: 414 ↛ 416line 414 didn't jump to line 416, because the condition on line 414 was never true
415 # nothing new to calculate
416 return
417 else:
418 if not self._vector_factory.accumulate: 418 ↛ 419line 418 didn't jump to line 419, because the condition on line 418 was never true
419 raise ValueError("Cannot increase the number of moments if "
420 "'accumulate_vectors' is 'False'.")
421 new_moments = n_moments - self.num_moments
422 m_start = self.num_moments
423 if new_moments < 0: 423 ↛ 424line 423 didn't jump to line 424, because the condition on line 423 was never true
424 raise ValueError('Cannot decrease number of moments')
426 if new_vectors != 0: 426 ↛ 427line 426 didn't jump to line 427, because the condition on line 426 was never true
427 raise ValueError("Only 'num_moments' *or* 'num_vectors' "
428 "may be updated at a time.")
430 for r in range(r_start, num_vectors):
431 alpha_zero = self._vector_factory[r]
433 one_moment = [0.] * n_moments
434 if new_vectors > 0:
435 alpha = alpha_zero
436 alpha_next = self.hamiltonian.matvec(alpha)
437 if self.operator is None:
438 one_moment[0] = np.vdot(alpha_zero, alpha_zero)
439 one_moment[1] = np.vdot(alpha_zero, alpha_next)
440 else:
441 one_moment[0] = self.operator(alpha_zero, alpha_zero)
442 one_moment[1] = self.operator(alpha_zero, alpha_next)
444 if new_moments > 0:
445 (alpha, alpha_next) = self._last_two_alphas[r]
446 one_moment[0:self.num_moments] = self._moments_list[r]
447 # Iteration over the moments
448 # Two cases can occur, depicted in Eq. (28) and in Eq. (29),
449 # respectively.
450 # ----
451 # In the first case, self.operator is None and we can use
452 # Eqs. (34) and (35) to obtain the density of states, with
453 # two moments ``one_moment`` for every new alpha.
454 # ----
455 # In the second case, the operator is not None and a matrix
456 # multiplication should be used.
457 if self.operator is None:
458 for n in range(m_start//2, n_moments//2):
459 alpha_save = alpha_next
460 alpha_next = (2 * self.hamiltonian.matvec(alpha_next)
461 - alpha)
462 alpha = alpha_save
463 # Following Eqs. (34) and (35)
464 one_moment[2*n] = (2 * np.vdot(alpha, alpha)
465 - one_moment[0])
466 one_moment[2*n+1] = (2 * np.vdot(alpha_next, alpha)
467 - one_moment[1])
468 if n_moments % 2:
469 # odd moment
470 one_moment[n_moments - 1] = (
471 2 * np.vdot(alpha_next, alpha_next) - one_moment[0])
472 else:
473 for n in range(m_start, n_moments):
474 alpha_save = alpha_next
475 alpha_next = (2 * self.hamiltonian.matvec(alpha_next)
476 - alpha)
477 alpha = alpha_save
478 one_moment[n] = self.operator(alpha_zero, alpha_next)
480 if self._vector_factory.accumulate:
481 self._last_two_alphas[r] = (alpha, alpha_next)
482 self._moments_list[r] = one_moment[:]
483 else:
484 self._moments_list[r] = one_moment
487class Correlator:
488 """Calculates the response of the correlation between two operators.
490 The response tensor :math:`χ_{α β}` of an operator :math:`O_α`
491 to a perturbation in an operator :math:`O_β`, is defined here based
492 on [3]_, and [4]_, and takes the form
494 .. math::
495 χ_{α β}(µ, T) =
496 \\int_{-\\infty}^{\\infty}{\\mathrm{d}E} f(µ-E, T)
497 \\left({O_α ρ(E) O_β \\frac{\\mathrm{d}G^{+}}{\\mathrm{d}E}} -
498 {O_α \\frac{\\mathrm{d}G^{-}}{\\mathrm{d}E} O_β ρ(E)}\\right)
500 .. [3] `Phys. Rev. Lett. 114, 116602 (2015)
501 <https://arxiv.org/abs/1410.8140>`_.
502 .. [4] `Phys. Rev. B 92, 184415 (2015)
503 <https://doi.org/10.1103/PhysRevB.92.184415>`_
505 Internally, the correlation is approximated with a
506 two dimensional KPM expansion,
508 .. math::
510 χ_{α β}(µ, T) =
511 \\int_{-1}^1{\\mathrm{d}E} \\frac{f(µ-E,T)}{(1-E^2)^2}
512 \\sum_{m,n}Γ_{n m}(E)µ_{n m}^{α β},
514 with coefficients
516 .. math::
518 Γ_{m n}(E) =
519 (E - i n \\sqrt{1 - E^2}) e^{i n \\arccos(E)} T_m(E)
521 + (E + i m \\sqrt{1 - E^2}) e^{-i m \\arccos(E)} T_n(E),
523 and moments matrix
524 :math:`µ_{n m}^{α β} = \\mathrm{Tr}(O_α T_m(H) O_β T_n(H))`.
526 The trace is calculated creating two instances of
527 `~kwant.kpm.SpectralDensity`, and saving the vectors
528 :math:`Ψ_{n r} = O_β T_n(H)\\rvert r\\rangle`,
529 and :math:`Ω_{m r} = T_m(H) O_α\\rvert r\\rangle` , where
530 :math:`\\rvert r\\rangle` is a vector provided by the
531 ``vector_factory``.
532 The moments matrix is formed with the product
533 :math:`µ_{m n} = \\langle Ω_{m r} \\rvert Ψ_{n r}\\rangle` for
534 every :math:`\\rvert r\\rangle`.
536 Parameters
537 ----------
538 hamiltonian : `~kwant.system.FiniteSystem` or matrix Hamiltonian
539 If a system is passed, it should contain no leads.
540 operator1, operator2 : operators, dense matrix, or sparse matrix, optional
541 Operators to be passed to two different instances of
542 `~kwant.kpm.SpectralDensity`.
543 **kwargs : dict
544 Keyword arguments to pass to `~kwant.kpm.SpectralDensity`.
546 Notes
547 -----
548 The ``operator1`` must act to the right as :math:`O_α\\rvert r\\rangle`.
549 """
551 def __init__(self, hamiltonian, operator1=None, operator2=None, **kwargs):
553 # Normalize 'operator1' and 'operator2' to functions that take
554 # and return a vector.
555 params = kwargs.get('params')
556 self.mean = kwargs.get('mean', True)
557 accumulate_vectors = kwargs.get('accumulate_vectors', False)
558 kwargs['accumulate_vectors'] = True
559 kwargs.pop('operator', None)
560 self.operator1 = _normalize_operator(operator1, params)
561 self.operator2 = _normalize_operator(operator2, params)
563 # initialize `SpectralDensity` to get `T_n(H)|r>` with a fake operator
564 def fake_op(bra, ket): return ket
566 # The vector factory used is the one passed by the user (or rng)
567 # to save the vectors, accumulate_vectors must be 'True'
568 self._spectrum_R = SpectralDensity(hamiltonian, operator=fake_op,
569 **kwargs)
570 self._a = self._spectrum_R._a
571 self._b = self._spectrum_R._b
572 _a = self._a * (1 - self._spectrum_R.eps / 2)
573 bounds = (self._b - _a, self._b + _a)
574 self.num_vectors = self._spectrum_R.num_vectors
575 self.num_moments = self._spectrum_R.num_moments
576 # apply operator2 to obtain `Psi_{n,r} = op2 T_n(H)|r>`
577 self._update_psi()
579 # instantiate the second `SpectralDensity`
580 # `accumulate_vectors` is set to the user defined option
581 # rewrite the bounds to match the rescaled bounds in the next call
582 kwargs['accumulate_vectors'] = accumulate_vectors
583 kwargs['num_vectors'] = self.num_vectors
584 kwargs['num_moments'] = self.num_moments
585 kwargs['energy_resolution'] = None
586 # Now we must take operator1 applied to the initial
587 # vectors to get `op1|r>`
588 # The vector factory used is defined below to ensure applying the
589 # same initial vectors stored in `self._vector_factory.saved_vectors`
590 kwargs['vector_factory'] = self._op_factory()
591 kwargs['bounds'] = bounds
592 self._spectrum_L = SpectralDensity(hamiltonian, operator=fake_op,
593 **kwargs)
594 # and now self._moments_list is `Omega_{m,r} = T_m(H) op1|r>`
595 # The shape of '_omega' is '(num_vecs, num_moments, dim_output)',
596 # where 'dim_output' is the dimension of the output of 'operator1'
597 self._omega = np.array(self._spectrum_L._moments_list)
599 self._calculate_moments_matrix()
600 self._build_integral_factor()
602 def __call__(self, mu=0, temperature=0):
603 """Returns the linear response :math:`χ_{α β}(µ, T)`
605 Parameters
606 ----------
607 mu : float
608 Chemical potential defined in the same units of energy as
609 the Hamiltonian.
610 temperature : float
611 Temperature in units of energy, the same as defined in the
612 Hamiltonian.
613 """
614 e = self.energies
615 e_rescaled = (e - self._b) / self._a
617 # rescale the energy to compare with the chemical potential
618 distribution_array = fermi_distribution(e, mu, temperature)
619 integrand = np.divide(distribution_array, (1 - e_rescaled ** 2) ** 2)
620 integrand = np.multiply(integrand, self._integral_factor)
621 integral = simps(integrand, x=e_rescaled)
622 # gives the linear response in units of volume * e^2/h
623 prefactor = 2 * 4**2 / ((2 * self._a) ** 2)
624 return prefactor * integral
626 @property
627 def energies(self):
628 return self._spectrum_R.energies
630 def add_moments(self, num_moments=None, *, energy_resolution=None):
631 """Increase the number of Chebyshev moments
633 Parameters
634 ----------
635 num_moments: positive int, optional
636 The number of Chebyshev moments to add. Mutually
637 exclusive with 'energy_resolution'.
638 energy_resolution: positive float, optional
639 Features wider than this resolution are visible
640 in the spectral density. Mutually exclusive with
641 ``num_moments``.
642 """
644 self._spectrum_R.add_moments(num_moments=num_moments,
645 energy_resolution=energy_resolution)
646 self.num_moments = self._spectrum_R.num_moments
647 # apply operator2 to obtain `Psi_{n,r} = op2
648 self._update_psi()
650 self._spectrum_L.add_moments(num_moments=num_moments,
651 energy_resolution=energy_resolution)
652 self._omega = np.array(self._spectrum_L._moments_list)
654 self._calculate_moments_matrix()
655 self._build_integral_factor()
657 def add_vectors(self, num_vectors=None):
658 """Increase the number of vectors
660 Parameters
661 ----------
662 num_vectors: positive int, optional
663 The number of vectors to add.
664 """
665 # get `T_n(H)|r>` with a fake operator
666 self._spectrum_R.add_vectors(num_vectors)
667 # apply operator2 to obtain `Psi_{n,r} = op2 T_n(H)|r>`
668 self._update_psi()
670 # _spectrum_L vector_factory is linked to _spectrum_R vector_factory
671 self._spectrum_L.add_vectors(num_vectors)
672 self.num_vectors = self._spectrum_L.num_vectors
673 # and now self._moments_list is `Omega_{m,r} = T_m(H) op1|r>`
674 self._omega = np.array(self._spectrum_L._moments_list)
676 self._calculate_moments_matrix()
677 self._build_integral_factor()
679 def _calculate_moments_matrix(self):
680 """Return the moments matrix, averaged over the vectors used """
681 # The final matrix is ready to be computed as
682 # `µ_{m,n} = <Omega_{m,r} | Psi_{n,r}>`
683 # for every `r` in `num_vectors`.
684 # 'moments_matrix' will be an array of moments matrix for each vector
685 # the shape of `moments_matrix` is
686 # `(num_vecs, num_moments, num_moments)`
687 self.moments_matrix = self._omega.conjugate() @ self._psi
688 if self.mean: 688 ↛ exitline 688 didn't return from function '_calculate_moments_matrix', because the condition on line 688 was never false
689 self.moments_matrix = np.mean(self.moments_matrix, axis=0)
691 def _op_factory(self):
692 """Factory of vectors ``operator1(vec[idx])``.
694 This factory will get updated with more vectors when
695 ``_spectrum_R._vector_factory`` gets updated to include more
696 vectors.
697 """
698 for vector in self._spectrum_R._vector_factory: 698 ↛ 700line 698 didn't jump to line 700, because the loop on line 698 didn't complete
699 yield self.operator1(vector)
700 return
702 def _update_psi(self):
703 """Axes are swapped in the end the get the shape
704 '(num_vecs, dim_output, num_moments)', where 'dim_output'
705 is the dimension of the output of 'operator2'."""
706 self._psi = np.array([
707 [
708 self.operator2(self._spectrum_R._moments_list[r][n])
709 for n in range(self._spectrum_R.num_moments)
710 ]
711 for r in range(self._spectrum_R.num_vectors)
712 ]).swapaxes(1, 2)
714 def _build_integral_factor(self):
715 """ Build the integral factor
717 .. math::
718 Γ_{m n}(E)
719 = (E - i n \\sqrt{1 - E^2}) e^{i n \\arccos(E)} T_m(E)
721 + (E + i m \\sqrt{1 - E^2}) e^{-i m \\arccos(E)} T_n(E),
723 times the moments matrix :math:`µ_{m n}` and sum over :math:`m`
724 and :math:`n`. :math:`E` is the array of the sampling points
725 selected as the Chebyshev nodes.
726 """
728 n_moments = self.num_moments
730 # get kernel array
731 g_kernel = self._spectrum_R.kernel(np.ones(n_moments))
732 g_kernel[0] /= 2
733 mu_kernel = np.outer(g_kernel, g_kernel) * self.moments_matrix
735 e_scaled = (self.energies - self._b) / self._a
737 m_array = np.arange(n_moments)
739 def _integral_factor(e):
740 # arrays for faster calculation
741 sqrt_e = np.sqrt(1 - e ** 2)
742 arccos_e = np.arccos(e)
744 exp_n = np.exp(1j * arccos_e * m_array)
745 t_n = np.real(exp_n)
747 e_plus = (e - 1j * sqrt_e * m_array)
748 e_plus = e_plus * exp_n
750 big_gamma = e_plus[None, :] * t_n[:, None]
751 big_gamma += big_gamma.conj().T
752 return np.tensordot(mu_kernel, big_gamma.T)
753 self._integral_factor = np.array([_integral_factor(e)
754 for e in e_scaled]).T
757def conductivity(hamiltonian, alpha='x', beta='x', positions=None, **kwargs):
758 """Returns a callable object to obtain the elements of the
759 conductivity tensor using the Kubo-Bastin approach.
761 A `~kwant.kpm.Correlator` instance is created to obtain the
762 correlation between two components of the current operator
764 .. math::
766 σ_{α β}(µ, T) =
767 \\frac{1}{V} \\int_{-\\infty}^{\\infty}{\\mathrm{d}E} f(µ-E, T)
768 \\left({j_α ρ(E) j_β \\frac{\\mathrm{d}G^{+}}{\\mathrm{d}E}} -
769 {j_α \\frac{\\mathrm{d}G^{-}}{\\mathrm{d}E} j_β ρ(E)}\\right),
771 where :math:`V` is the volume where the conductivity is sampled.
772 In this implementation it is assumed that the vectors are normalized
773 and :math:`V=1`, otherwise the result of this calculation must be
774 normalized with the corresponding volume.
776 The equations used here are based on [3]_ and [4]_
778 .. [3] `Phys. Rev. Lett. 114, 116602 (2015)
779 <https://arxiv.org/abs/1410.8140>`_.
780 .. [4] `Phys. Rev. B 92, 184415 (2015)
781 <https://doi.org/10.1103/PhysRevB.92.184415>`_
783 Parameters
784 ----------
785 hamiltonian : `~kwant.system.FiniteSystem` or matrix Hamiltonian
786 If a system is passed, it should contain no leads.
787 alpha, beta : str, or operators
788 If ``hamiltonian`` is a kwant system, or if the ``positions``
789 are provided, ``alpha`` and ``beta`` can be the directions of the
790 velocities as strings {'x', 'y', 'z'}.
791 Otherwise ``alpha`` and ``beta`` should be the proper velocity
792 operators, which can be members of `~kwant.operator` or matrices.
793 positions : array of float, optioinal
794 If ``hamiltonian`` is a matrix, the velocities can be calculated
795 internally by passing the positions of each orbital in the system
796 when ``alpha`` or ``beta`` are one of the directions {'x', 'y', 'z'}.
797 **kwargs : dict
798 Keyword arguments to pass to `~kwant.kpm.Correlator`.
800 Examples
801 --------
802 We will obtain the conductivity of the Haldane model, defined as a
803 honeycomb lattice with first nearest neighbors coupling, and
804 imaginary second nearest neighbors coupling.
806 We start by importing kwant and defining a
807 `~kwant.system.FiniteSystem`,
809 >>> import kwant
810 ...
811 >>> def circle(pos):
812 ... x, y = pos
813 ... return x**2 + y**2 < 100
814 ...
815 >>> lat = kwant.lattice.honeycomb()
816 >>> syst = kwant.Builder()
817 >>> syst[lat.shape(circle, (0, 0))] = 0
818 >>> syst[lat.neighbors()] = -1
819 >>> syst[lat.a.neighbors()] = -0.5j
820 >>> syst[lat.b.neighbors()] = 0.5j
821 >>> fsyst = syst.finalized()
823 Now we can call `~kwant.kpm.conductivity` to calculate the transverse
824 conductivity at chemical potential 0 and temperature 0.01.
826 >>> cond = kwant.kpm.conductivity(fsyst, alpha='x', beta='y')
827 >>> cond(mu=0, temperature=0.01)
828 """
830 if positions is None and not isinstance(hamiltonian, system.System): 830 ↛ 831line 830 didn't jump to line 831, because the condition on line 830 was never true
831 raise ValueError("If 'hamiltonian' is a matrix, positions "
832 "must be provided")
834 params = kwargs.get('params')
835 alpha = _velocity(hamiltonian, params, alpha, positions)
836 beta = _velocity(hamiltonian, params, beta, positions)
838 correlator = Correlator(
839 hamiltonian, operator1=alpha, operator2=beta, **kwargs)
841 return correlator
844class _VectorFactory:
845 """Factory for Hilbert space vectors.
847 Parameters
848 ----------
849 vectors : iterable
850 Iterable of Hilbert space vectors.
851 num_vectors : int, optional
852 Total number of vectors. If not specified, will be set to the
853 length of 'vectors'.
854 accumulate : bool, default: True
855 If True, the attribute 'saved_vectors' will store the vectors
856 generated.
857 """
859 def __init__(self, vectors=None, num_vectors=None, accumulate=True):
860 assert isinstance(vectors, Iterable)
861 try:
862 _len = len(vectors)
863 if num_vectors is None:
864 num_vectors = _len
865 except TypeError:
866 _len = np.inf
867 if num_vectors is None: 867 ↛ 868line 867 didn't jump to line 868, because the condition on line 867 was never true
868 raise ValueError("'num_vectors' must be specified when "
869 "'vectors' has no len() method.")
870 self._max_vectors = _len
871 self._iterator = iter(vectors)
873 self.accumulate = accumulate
874 self.saved_vectors = []
876 self.num_vectors = 0
877 self.add_vectors(num_vectors=num_vectors)
879 self._last_idx = -np.inf
880 self._last_vector = None
882 def _fill_in_saved_vectors(self, index):
883 if index < self._last_idx and not self.accumulate: 883 ↛ 884line 883 didn't jump to line 884, because the condition on line 883 was never true
884 raise ValueError("Cannot get previous values if 'accumulate' "
885 "is False")
887 if index >= self.num_vectors: 887 ↛ 888line 887 didn't jump to line 888, because the condition on line 887 was never true
888 raise IndexError('Requested more vectors than available')
890 self._last_idx = index
891 if self.accumulate:
892 if self.saved_vectors[index] is None:
893 self.saved_vectors[index] = next(self._iterator)
894 else:
895 self._last_vector = next(self._iterator)
897 def __getitem__(self, index):
898 self._fill_in_saved_vectors(index)
899 if self.accumulate:
900 return self.saved_vectors[index]
901 return self._last_vector
903 def add_vectors(self, num_vectors=None):
904 """Increase the number of vectors
906 Parameters
907 ----------
908 num_vectors: positive int, optional
909 The number of vectors to add.
910 """
911 if num_vectors is None: 911 ↛ 912line 911 didn't jump to line 912, because the condition on line 911 was never true
912 raise ValueError("'num_vectors' must be specified.")
913 else:
914 if num_vectors <= 0 or num_vectors != int(num_vectors):
915 raise ValueError("'num_vectors' must be a positive integer")
916 elif self.num_vectors + num_vectors > self._max_vectors: 916 ↛ 917line 916 didn't jump to line 917, because the condition on line 916 was never true
917 raise ValueError("'num_vectors' is larger than available "
918 "vectors")
920 self.num_vectors += num_vectors
922 if self.accumulate:
923 self.saved_vectors.extend([None] * num_vectors)
926def RandomVectors(syst, where=None, rng=None):
927 """Returns a random phase vector iterator for the sites in 'where'.
929 Parameters
930 ----------
931 syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
932 If a system is passed, it should contain no leads.
933 where : sequence of `int` or `~kwant.system.Site`, or callable, optional
934 Spatial range of the vectors produced. If ``syst`` is a
935 `~kwant.system.FiniteSystem`, where behaves as in
936 `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
937 must be a list of integers with the indices where column vectors
938 are nonzero.
939 """
940 rng = ensure_rng(rng)
941 tot_norbs, orbs = _normalize_orbs_where(syst, where)
942 while True:
943 vector = np.zeros(tot_norbs, dtype=complex)
944 vector[orbs] = np.exp(2j * np.pi * rng.random_sample(len(orbs)))
945 yield vector
948class LocalVectors:
949 """Generates a local vector iterator for the sites in 'where'.
951 Parameters
952 ----------
953 syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
954 If a system is passed, it should contain no leads.
955 where : sequence of `int` or `~kwant.system.Site`, or callable, optional
956 Spatial range of the vectors produced. If ``syst`` is a
957 `~kwant.system.FiniteSystem`, where behaves as in
958 `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
959 must be a list of integers with the indices where column vectors
960 are nonzero.
961 """
963 def __init__(self, syst, where=None, *args):
964 self.tot_norbs, self.orbs = _normalize_orbs_where(syst, where)
965 self._idx = 0
967 def __len__(self,):
968 return len(self.orbs)
970 def __iter__(self,):
971 return self
973 def __next__(self,):
974 if self._idx < len(self): 974 ↛ 979line 974 didn't jump to line 979, because the condition on line 974 was never false
975 vector = np.zeros(self.tot_norbs)
976 vector[self.orbs[self._idx]] = 1
977 self._idx = self._idx + 1
978 return vector
979 raise StopIteration('Too many vectors requested from this generator')
981# ### Auxiliary functions
984def fermi_distribution(energy, mu, temperature):
985 """Returns the Fermi distribution f(e, µ, T) evaluated at 'e'.
987 Parameters
988 ----------
989 energy : float or sequence of floats
990 Energy array where the Fermi distribution is evaluated.
991 mu : float
992 Chemical potential defined in the same units of energy as
993 the Hamiltonian.
994 temperature : float
995 Temperature in units of energy, the same as defined in the
996 Hamiltonian.
997 """
998 if temperature < 0: 998 ↛ 999line 998 didn't jump to line 999, because the condition on line 998 was never true
999 raise ValueError("temperature must be non-negative")
1000 elif temperature == 0: 1000 ↛ 1003line 1000 didn't jump to line 1003, because the condition on line 1000 was never false
1001 return np.array(np.less(energy - mu, 0), dtype=float)
1002 else:
1003 return 1 / (1 + np.exp((energy - mu) / temperature))
1006def _from_where_to_orbs(syst, where):
1007 """Returns a list of slices of the orbitals in 'where'"""
1008 assert isinstance(syst, system.System)
1009 where = _normalize_site_where(syst, where)
1010 _site_ranges = np.asarray(syst.site_ranges, dtype=gint_dtype)
1011 offsets, norbs = _get_all_orbs(where, _site_ranges)
1012 # concatenate all the orbitals
1013 orbs = [list(range(start, start+orbs))
1014 for start, orbs in zip(offsets[:, 0], norbs[:, 0])]
1015 orbs = reduce(add, orbs)
1016 return orbs
1019def _normalize_orbs_where(syst, where):
1020 """Return total number of orbitals and a list of slices of
1021 orbitals in 'where'"""
1022 if isinstance(syst, system.System):
1023 tot_norbs = _get_tot_norbs(syst)
1024 orbs = _from_where_to_orbs(syst, where)
1025 else:
1026 try:
1027 tot_norbs = csr_matrix(syst).shape[0]
1028 except TypeError:
1029 raise TypeError("'syst' is neither a matrix "
1030 "nor a Kwant system.")
1031 orbs = (range(tot_norbs) if where is None
1032 else np.asarray(where, int))
1033 return tot_norbs, orbs
1036def _velocity(hamiltonian, params, op_type, positions):
1037 """Construct the velocity operator
1039 Parameters
1040 ----------
1041 hamiltonian : ndarray or a Kwant system
1042 System for which the velocity operator is calculated.
1043 params : dict
1044 Parametres of the system
1045 op_type : str, matrix or operator
1046 If ``op_type`` is a 'str' in {'x', 'y', 'z'}, the velocity operator
1047 is calculated using the ``hamiltonian`` and ``positions``, else
1048 if ``op_type`` is an operator or a matrix, this is the velocity
1049 operator.
1050 positions : ndarray of shape ``(N, dim)``
1051 Positions of each orbital. This parameter is not used if
1052 ``hamiltonian`` is a Kwant system.
1053 """
1054 directions = dict(x=0, y=1, z=2)
1056 if isinstance(op_type, _LocalOperator): 1056 ↛ 1057line 1056 didn't jump to line 1057, because the condition on line 1056 was never true
1057 operator = op_type
1058 elif isinstance(op_type, str):
1059 direction = directions[op_type]
1060 if isinstance(hamiltonian, system.System):
1061 operator, norbs, norbs = hamiltonian.hamiltonian_submatrix(
1062 params=params, sparse=True, return_norb=True
1063 )
1064 positions = np.vstack([[hamiltonian.pos(i)] * norb
1065 for i, norb in enumerate(norbs)])
1066 elif positions is not None: 1066 ↛ 1068line 1066 didn't jump to line 1068, because the condition on line 1066 was never false
1067 operator = coo_matrix(hamiltonian, copy=True)
1068 displacements = positions[operator.col] - positions[operator.row]
1069 if direction > displacements.shape[1]: 1069 ↛ 1070line 1069 didn't jump to line 1070, because the condition on line 1069 was never true
1070 raise ValueError("{} is not an allowed direction.".format(op_type))
1071 operator.data *= 1j * displacements[:, direction]
1072 operator = operator.tocsr()
1073 else:
1074 try:
1075 operator = csr_matrix(op_type)
1076 except Exception:
1077 raise ValueError("Velocity operator must be provided as a matrix, "
1078 "a kwant operator, or a direction.")
1079 return operator
1082def _normalize_operator(op, params):
1083 """Normalize 'op' to a function that takes and returns a vector."""
1084 if op is None:
1085 def r_op(v): return v
1086 elif isinstance(op, _LocalOperator): 1086 ↛ 1087line 1086 didn't jump to line 1087, because the condition on line 1086 was never true
1087 op = op.bind(params=params)
1088 r_op = op.act
1089 elif callable(op): 1089 ↛ 1090line 1089 didn't jump to line 1090, because the condition on line 1089 was never true
1090 r_op = op
1091 elif hasattr(op, 'dot'): 1091 ↛ 1095line 1091 didn't jump to line 1095, because the condition on line 1091 was never false
1092 op = csr_matrix(op)
1093 r_op = op.dot
1094 else:
1095 raise TypeError("The operators must have a '.dot' "
1096 "attribute or must be callable.")
1097 return r_op
1100def jackson_kernel(moments):
1101 """Convolutes ``moments`` with the Jackson kernel.
1103 Taken from Eq. (71) of `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
1104 <https://arxiv.org/abs/cond-mat/0504627>`_.
1105 """
1107 n_moments = len(moments)
1108 m = np.arange(n_moments)
1109 kernel_array = ((n_moments - m + 1) *
1110 np.cos(np.pi * m/(n_moments + 1)) +
1111 np.sin(np.pi * m/(n_moments + 1)) /
1112 np.tan(np.pi/(n_moments + 1)))
1113 kernel_array /= n_moments + 1
1115 # transposes handle the case where operators have vector outputs
1116 conv_moments = np.transpose(moments.transpose() * kernel_array)
1117 return conv_moments
1120def lorentz_kernel(moments, l=4):
1121 """Convolutes ``moments`` with the Lorentz kernel.
1123 Taken from Eq. (71) of `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
1124 <https://arxiv.org/abs/cond-mat/0504627>`_.
1126 The additional parameter ``l`` controls the decay of the kernel.
1127 """
1129 n_moments = len(moments)
1131 m = np.arange(n_moments)
1132 kernel_array = np.sinh(l * (1 - m / n_moments)) / np.sinh(l)
1134 # transposes handle the case where operators have vector outputs
1135 conv_moments = np.transpose(moments.transpose() * kernel_array)
1136 return conv_moments
1139def _rescale(hamiltonian, eps, v0, bounds):
1140 """Rescale a Hamiltonian and return a LinearOperator
1142 Parameters
1143 ----------
1144 hamiltonian : 2D array
1145 Hamiltonian of the system.
1146 eps : scalar
1147 Ensures that the bounds are strict.
1148 v0 : random vector, or None
1149 Used as the initial residual vector for the algorithm that
1150 finds the lowest and highest eigenvalues.
1151 bounds : tuple, or None
1152 Boundaries of the spectrum. If not provided the maximum and
1153 minimum eigenvalues are calculated.
1154 """
1155 # Relative tolerance to which to calculate eigenvalues. Because after
1156 # rescaling we will add eps / 2 to the spectral bounds, we don't need
1157 # to know the bounds more accurately than eps / 2.
1158 tol = eps / 2
1160 if bounds:
1161 lmin, lmax = bounds
1162 else:
1163 lmax = float(eigsh(hamiltonian, k=1, which='LA',
1164 return_eigenvectors=False, tol=tol, v0=v0))
1165 lmin = float(eigsh(hamiltonian, k=1, which='SA',
1166 return_eigenvectors=False, tol=tol, v0=v0))
1168 a = np.abs(lmax-lmin) / (2. - eps)
1169 b = (lmax+lmin) / 2.
1171 if lmax - lmin <= abs(lmax + lmin) * tol / 2:
1172 raise ValueError(
1173 'The Hamiltonian has a single eigenvalue, it is not possible to '
1174 'obtain a spectral density.')
1176 def rescaled(v):
1177 return (hamiltonian.dot(v) - b * v) / a
1179 rescaled_ham = LinearOperator(shape=hamiltonian.shape, matvec=rescaled)
1181 return rescaled_ham, (a, b)
1184def _chebyshev_nodes(n_sampling):
1185 """Return an array of 'n_sampling' points in the interval (-1,1)"""
1186 raw, step = np.linspace(np.pi, 0, n_sampling,
1187 endpoint=False, retstep=True)
1188 return np.cos(raw + step / 2)
1191def _calc_fft_moments(moments):
1192 """This function takes the stabilized moments and returns an array
1193 of points and an array of the evaluated function at those points.
1194 """
1195 moments = np.asarray(moments)
1196 # extra_shape handles the case where operators have vector outputs
1197 n_moments, *extra_shape = moments.shape
1198 n_sampling = SAMPLING * n_moments
1199 moments_ext = np.zeros([n_sampling] + extra_shape, dtype=moments.dtype)
1201 # special points at the abscissas of Chebyshev integration
1202 e_rescaled = _chebyshev_nodes(n_sampling)
1204 # transposes handle the case where operators have vector outputs
1205 moments_ext[:n_moments] = moments
1206 # The function evaluated in these special data points is the FFT of
1207 # the moments as in Eq. (83).
1208 # The order of gammas must be reversed to match the energies in
1209 # ascending order
1210 gammas = np.transpose(fft.dct(moments_ext.transpose(), type=3))[::-1]
1212 # Element-wise division of moments_FFT over gk, as in Eq. (83).
1213 gk = np.pi * np.sqrt(1 - e_rescaled ** 2)
1214 rho = np.transpose(np.divide(gammas.transpose(), gk))
1216 return rho, gammas