diff --git a/kwant/physics/dispersion.py b/kwant/physics/dispersion.py
index 9934ad71a8f2ea236c09e31b239094adf5a994ce..5c98bf2618536fea4262df0ef8476dcc9181e721 100644
--- a/kwant/physics/dispersion.py
+++ b/kwant/physics/dispersion.py
@@ -6,6 +6,8 @@
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
+"""Band structure calculation for the leads"""
+
import math
import numpy as np
from .. import system
@@ -35,7 +37,8 @@ class Bands:
An instance of this class can be called like a function. Given a momentum
(currently this must be a scalar as all infinite systems are quasi-1-d), it
returns a NumPy array containing the eigenenergies of all modes at this
- momentum
+ momentum. If `derivative_order > 0` or `return_eigenvectors = True`,
+ additional arrays are returned.
Examples
--------
@@ -45,6 +48,7 @@ class Bands:
>>> pyplot.plot(momenta, energies)
>>> pyplot.show()
"""
+ _crossover_size = 8
def __init__(self, sys, args=(), *, params=None):
syst = sys
@@ -57,9 +61,100 @@ class Bands:
self.hop[:, : hop.shape[1]] = hop
self.hop[:, hop.shape[1]:] = 0
- def __call__(self, k):
- # Note: Equation to solve is
- # (V^\dagger e^{ik} + H + V e^{-ik}) \psi = E \psi
+ def __call__(self, k, derivative_order=0, return_eigenvectors=False):
+ r"""Calculate band energies at a given momentum.
+
+ :math:`E_n, \quad E'_n = dE_n / dk, \quad E''_n = d^2E_n / dk^2,
+ \quad n \in \{0, nbands - 1\}`
+
+ :math:`nbands` is the number of open modes.
+ Eigenvectors are orthonormal.
+
+ Parameters
+ ----------
+ k : float
+ momentum
+ derivative_order : {0, 1, 2}, optional
+ Maximal derivative order to calculate. Default is zero.
+ return_eigenvectors : bool, optional
+ if set to `True` return the eigenvectors as last tuple element.
+ By default, no eigenvectors are returned.
+
+ Returns
+ -------
+ energies : numpy float array, size ``nbands``
+ energies :math:`E`
+ velocities : numpy float array, size ``nbands``
+ velocities (first order energy derivatives :math:`E'`)
+ curvatures : numpy float array, size ``nbands``
+ curvatures (second energy derivatives :math:`E''`)
+ eigenvectors : numpy float array, size ``nbands x nbands``
+ eigenvectors
+
+ The number of returned elements varies. If ``derivative_order > 0``
+ we return all derivatives up to the order ``derivative_order``,
+ and total the number of returned elements is ``derivative_order + 1``
+ If ``return_eigenvectors`` is True, in addition the eigenvectors are
+ returned as the last element. In that case, the total number of
+ returned elements is ``derivative_order + 2``.
+
+ Notes
+ -----
+ * All the output arrays are sorted from the lowest energy band
+ to the highest.
+ * The curvature `E''` can be only calculated for non-degenerate bands.
+ """
+ # Equation to solve is
+ # (V^\dagger e^{ik} + H + V e^{-ik}) \psi = E \psi
+ # we obtain the derivatives by perturbing around momentum k
+ # The perturbed Hamiltonian H(k) is
+ # H(k+dk) = H_0(k) + H_1(k) dk + H_2(k) dk^2 + O(dk^3)
+ # and we use h0, h1 and h2 in the code to refer to above terms
+
mat = self.hop * complex(math.cos(k), -math.sin(k))
- mat += mat.conjugate().transpose() + self.ham
- return np.sort(np.linalg.eigvalsh(mat).real)
+ h0 = mat + mat.conjugate().transpose() + self.ham
+
+ # compute energies and if required eigenvectors
+ # numpy routines eigh and eigvalsh return eigenvalues in ascending order
+ if return_eigenvectors or derivative_order > 0:
+ energies, eigenvectors = np.linalg.eigh(h0)
+ else:
+ energies = np.linalg.eigvalsh(h0)
+ output = (energies.real,)
+
+ if derivative_order >= 1: # compute velocities
+ h1 = 1j*(- mat + mat.conjugate().transpose())
+ # Switch computation from diag(a^\dagger @ b) to sum(a.T * b, axis=0),
+ # which is computationally more efficient for larger matrices.
+ # a and b are NxN matrices, and we switch for N > _crossover_size.
+ # Both ways are numerically (but not necessarily bitwise) equivalent.
+ if derivative_order == 1 and len(energies) > self._crossover_size:
+ velocities = np.sum(eigenvectors.conjugate() * (h1 @ eigenvectors),
+ axis=0).real
+ else:
+ ph1p = eigenvectors.conjugate().transpose() @ h1 @ eigenvectors
+ velocities = np.diag(ph1p).real
+ output += (velocities,)
+
+ if derivative_order >= 2: # compute curvatures
+ # ediff_{i,j} = 1 / (E_i - E_j) if i != j, else 0
+ ediff = energies.reshape((-1, 1)) - energies.reshape((1, -1))
+ ediff = np.divide(1, ediff, where=(ediff != 0))
+ h2 = - (mat + mat.conjugate().transpose())
+ curvatures = np.sum(eigenvectors.conjugate() * (h2 @ eigenvectors)
+ - 2 * ediff * np.abs(ph1p)**2, axis=0).real
+ # above expression is similar to: curvatures =
+ # np.diag(eigenvectors.conjugate().transpose() @ h2 @ eigenvectors
+ # + 2 * ediff @ np.abs(ph1p)**2).real
+ output += (curvatures,)
+
+ if derivative_order > 2:
+ raise NotImplementedError('Derivatives of the energy dispersion ' +
+ 'only implemented up to second order.')
+ if return_eigenvectors:
+ output += (eigenvectors,)
+ if len(output) == 1:
+ # Backwards compatibility: if returning only energies,
+ # don't make it a length-1 tuple.
+ return output[0]
+ return output
diff --git a/kwant/physics/tests/test_dispersion.py b/kwant/physics/tests/test_dispersion.py
index 414ac17854a8e491c8010900686f938b5f35d5ac..a42d6deadf5d88bb0390318cb5442f88229742dc 100644
--- a/kwant/physics/tests/test_dispersion.py
+++ b/kwant/physics/tests/test_dispersion.py
@@ -6,26 +6,33 @@
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
+import numpy as np
from numpy.testing import assert_array_almost_equal, assert_almost_equal
from pytest import raises
+from numpy import linspace
import kwant
from math import pi, cos, sin
-def test_band_energies(N=5):
+
+def make_lead():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
lat = kwant.lattice.square()
syst[[lat(0, 0), lat(0, 1)]] = 3
syst[lat(0, 1), lat(0, 0)] = -1
syst[((lat(1, y), lat(0, y)) for y in range(2))] = -1
+ return syst
+
- band_energies = kwant.physics.Bands(syst.finalized())
+def test_band_energies(N=5):
+ band_energies = kwant.physics.Bands(make_lead().finalized())
for i in range(-N, N):
k = i * pi / N
energies = band_energies(k)
assert_array_almost_equal(sorted(energies),
sorted([2 - 2 * cos(k), 4 - 2 * cos(k)]))
+
def test_same_as_lead():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
lat = kwant.lattice.chain()
@@ -39,6 +46,7 @@ def test_same_as_lead():
for momentum in momenta:
assert_almost_equal(bands(momentum)[0], 0)
+
def test_raise_nonhermitian():
syst = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
lat = kwant.lattice.chain()
@@ -46,3 +54,97 @@ def test_raise_nonhermitian():
syst[lat(0), lat(1)] = complex(cos(0.2), sin(0.2))
syst = syst.finalized()
raises(ValueError, kwant.physics.Bands, syst)
+
+
+def test_band_velocities():
+ syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+ lat = kwant.lattice.square()
+ syst[lat(0, 0)] = 1
+ syst[lat(0, 1)] = 3
+ syst[lat(1, 0), lat(0, 0)] = -1
+ syst[lat(1, 1), lat(0, 1)] = 2
+ bands = kwant.physics.Bands(syst.finalized())
+ eps = 1E-4
+ for k in linspace(-pi, pi, 200):
+ vel = bands(k, derivative_order=1)[1]
+ # higher order formula for first derivative to get required accuracy
+ num_vel = (- bands(k+2*eps) + bands(k-2*eps) +
+ 8*(bands(k+eps) - bands(k-eps))) / (12 * eps)
+ assert_array_almost_equal(vel, num_vel)
+
+
+def test_band_velocity_derivative():
+ syst = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+ lat = kwant.lattice.square()
+ syst[lat(0, 0)] = 1
+ syst[lat(0, 1)] = 3
+ syst[lat(1, 0), lat(0, 0)] = -1
+ syst[lat(1, 1), lat(0, 1)] = 2
+ bands = kwant.physics.Bands(syst.finalized())
+ eps = 1E-4
+ eps2 = eps * eps
+ c3 = 1 / 90
+ c2 = - 3 / 20
+ c1 = 3 / 2
+ c0 = - 49 / 18
+ for k in linspace(-pi, pi, 200):
+ dvel = bands(k, derivative_order=2)[2]
+ # higher order formula for second derivative to get required accuracy
+ num_dvel = (c3 * (bands(k+3*eps) + bands(k-3*eps)) +
+ c2 * (bands(k+2*eps) + bands(k-2*eps)) +
+ c1 * (bands(k+eps) + bands(k-eps)) +
+ c0 * bands(k)) / eps2
+ assert_array_almost_equal(dvel, num_dvel)
+
+
+def test_eigenvector_calculation():
+ lead = make_lead().finalized()
+ bands = kwant.physics.Bands(lead)
+ k = 1.2 # k-point arbitrary
+
+ # check if eigenvalues are sorted and always sorted in the same way
+ energies = bands(k)
+ assert (np.sort(energies) == energies).all()
+ e_1 = bands(k, derivative_order=0)
+ e_2, *_ = bands(k, derivative_order=0, return_eigenvectors=True)
+ assert_array_almost_equal(energies, e_1)
+ assert_array_almost_equal(energies, e_2)
+ for order in [1, 2]:
+ e_1, *_ = bands(k, derivative_order=order)
+ e_2, *_ = bands(k, derivative_order=order, return_eigenvectors=True)
+ assert_array_almost_equal(energies, e_1)
+ assert_array_almost_equal(energies, e_2)
+
+ # check eigenvector
+ unity_matrix = np.identity(len(energies))
+ zero_matrix = 0 * unity_matrix
+ for order in [0, 1, 2]:
+ eval_evec = bands(k, derivative_order=order, return_eigenvectors=True)
+
+ assert len(eval_evec) == order + 2
+ energies, *_, psi = eval_evec
+
+ # completeness relation
+ outer_product = psi @ psi.conjugate().transpose()
+ assert_array_almost_equal(outer_product.real, unity_matrix)
+ assert_array_almost_equal(outer_product.imag, zero_matrix)
+
+ # eigenvalue equation
+ mat = bands.hop * np.exp(- 1j * k)
+ hamiltonian = mat + mat.conjugate().transpose() + bands.ham
+ energies = np.asarray([energies])
+ assert_array_almost_equal(hamiltonian @ psi, energies * psi)
+
+ # check that algorithm switch is correct
+ bands_switch = kwant.physics.Bands(lead)
+ bands_switch._crossover_size = 0
+ for order in [1, 2]:
+ e_1, v_1, *_ = bands(k, order)
+ e_2, v_2, *_ = bands_switch(k, order)
+ assert_array_almost_equal(v_1, v_2)
+
+
+def test_raise_implemented():
+ k = 1 # k-point arbitrary
+ bands = kwant.physics.Bands(make_lead().finalized())
+ raises(NotImplementedError, bands, k, derivative_order=3)