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# Copyright 2011-2013 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
9"""Band structure calculation for the leads"""
11import math
12import numpy as np
13from .. import system
14from .._common import ensure_isinstance, deprecate_args
16__all__ = ['Bands']
19class Bands:
20 """
21 Class of callable objects for the computation of energy bands.
23 Parameters
24 ----------
25 sys : `kwant.system.InfiniteSystem`
26 The low level infinite system for which the energies are to be
27 calculated.
28 args : tuple, defaults to empty
29 Positional arguments to pass to the ``hamiltonian`` method.
30 Deprecated in favor or 'params' (and mutually exclusive with it).
31 params : dict, optional
32 Dictionary of parameter names and their values. Mutually exclusive
33 with 'args'.
35 Notes
36 -----
37 An instance of this class can be called like a function. Given a momentum
38 (currently this must be a scalar as all infinite systems are quasi-1-d), it
39 returns a NumPy array containing the eigenenergies of all modes at this
40 momentum. If `derivative_order > 0` or `return_eigenvectors = True`,
41 additional arrays are returned.
43 Examples
44 --------
45 >>> bands = kwant.physics.Bands(some_syst)
46 >>> momenta = numpy.linspace(-numpy.pi, numpy.pi, 101)
47 >>> energies = [bands(k) for k in momenta]
48 >>> pyplot.plot(momenta, energies)
49 >>> pyplot.show()
50 """
51 _crossover_size = 8
53 @deprecate_args
54 def __init__(self, sys, args=(), *, params=None):
55 syst = sys
56 ensure_isinstance(syst, (system.InfiniteSystem,
57 system.InfiniteVectorizedSystem))
58 self.ham = syst.cell_hamiltonian(args, params=params)
59 if not np.allclose(self.ham, self.ham.T.conj()):
60 raise ValueError('The cell Hamiltonian is not Hermitian.')
61 hop = syst.inter_cell_hopping(args, params=params)
62 self.hop = np.empty(self.ham.shape, dtype=complex)
63 self.hop[:, : hop.shape[1]] = hop
64 self.hop[:, hop.shape[1]:] = 0
66 def __call__(self, k, derivative_order=0, return_eigenvectors=False):
67 r"""Calculate band energies at a given momentum.
69 :math:`E_n, \quad E'_n = dE_n / dk, \quad E''_n = d^2E_n / dk^2,
70 \quad n \in \{0, nbands - 1\}`
72 :math:`nbands` is the number of open modes.
73 Eigenvectors are orthonormal.
75 Parameters
76 ----------
77 k : float
78 momentum
79 derivative_order : {0, 1, 2}, optional
80 Maximal derivative order to calculate. Default is zero.
81 return_eigenvectors : bool, optional
82 if set to `True` return the eigenvectors as last tuple element.
83 By default, no eigenvectors are returned.
85 Returns
86 -------
87 energies : numpy float array, size ``nbands``
88 energies :math:`E`
89 velocities : numpy float array, size ``nbands``
90 velocities (first order energy derivatives :math:`E'`)
91 curvatures : numpy float array, size ``nbands``
92 curvatures (second energy derivatives :math:`E''`)
93 eigenvectors : numpy float array, size ``nbands x nbands``
94 eigenvectors
96 The number of returned elements varies. If ``derivative_order > 0``
97 we return all derivatives up to the order ``derivative_order``,
98 and total the number of returned elements is ``derivative_order + 1``
99 If ``return_eigenvectors`` is True, in addition the eigenvectors are
100 returned as the last element. In that case, the total number of
101 returned elements is ``derivative_order + 2``.
103 Notes
104 -----
105 * All the output arrays are sorted from the lowest energy band
106 to the highest.
107 * The curvature `E''` can be only calculated for non-degenerate bands.
108 """
109 # Equation to solve is
110 # (V^\dagger e^{ik} + H + V e^{-ik}) \psi = E \psi
111 # we obtain the derivatives by perturbing around momentum k
112 # The perturbed Hamiltonian H(k) is
113 # H(k+dk) = H_0(k) + H_1(k) dk + H_2(k) dk^2 + O(dk^3)
114 # and we use h0, h1 and h2 in the code to refer to above terms
116 mat = self.hop * complex(math.cos(k), -math.sin(k))
117 h0 = mat + mat.conjugate().transpose() + self.ham
119 # compute energies and if required eigenvectors
120 # numpy routines eigh and eigvalsh return eigenvalues in ascending order
121 if return_eigenvectors or derivative_order > 0:
122 energies, eigenvectors = np.linalg.eigh(h0)
123 else:
124 energies = np.linalg.eigvalsh(h0)
125 output = (energies.real,)
127 if derivative_order >= 1: # compute velocities
128 h1 = 1j*(- mat + mat.conjugate().transpose())
129 # Switch computation from diag(a^\dagger @ b) to sum(a.T * b, axis=0),
130 # which is computationally more efficient for larger matrices.
131 # a and b are NxN matrices, and we switch for N > _crossover_size.
132 # Both ways are numerically (but not necessarily bitwise) equivalent.
133 if derivative_order == 1 and len(energies) > self._crossover_size:
134 velocities = np.sum(eigenvectors.conjugate() * (h1 @ eigenvectors),
135 axis=0).real
136 else:
137 ph1p = eigenvectors.conjugate().transpose() @ h1 @ eigenvectors
138 velocities = np.diag(ph1p).real
139 output += (velocities,)
141 if derivative_order >= 2: # compute curvatures
142 # ediff_{i,j} = 1 / (E_i - E_j) if i != j, else 0
143 ediff = energies.reshape((-1, 1)) - energies.reshape((1, -1))
144 ediff = np.divide(1, ediff, where=(ediff != 0))
145 h2 = - (mat + mat.conjugate().transpose())
146 curvatures = np.sum(eigenvectors.conjugate() * (h2 @ eigenvectors)
147 - 2 * ediff * np.abs(ph1p)**2, axis=0).real
148 # above expression is similar to: curvatures =
149 # np.diag(eigenvectors.conjugate().transpose() @ h2 @ eigenvectors
150 # + 2 * ediff @ np.abs(ph1p)**2).real
151 output += (curvatures,)
153 if derivative_order > 2:
154 raise NotImplementedError('Derivatives of the energy dispersion ' +
155 'only implemented up to second order.')
156 if return_eigenvectors:
157 output += (eigenvectors,)
158 if len(output) == 1:
159 # Backwards compatibility: if returning only energies,
160 # don't make it a length-1 tuple.
161 return output[0]
162 return output