Hide keyboard shortcuts

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. 

8 

9"""Band structure calculation for the leads""" 

10 

11import math 

12import numpy as np 

13from .. import system 

14from .._common import ensure_isinstance, deprecate_args 

15 

16__all__ = ['Bands'] 

17 

18 

19class Bands: 

20 """ 

21 Class of callable objects for the computation of energy bands. 

22 

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'. 

34 

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. 

42 

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 

52 

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 

65 

66 def __call__(self, k, derivative_order=0, return_eigenvectors=False): 

67 r"""Calculate band energies at a given momentum. 

68 

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\}` 

71 

72 :math:`nbands` is the number of open modes. 

73 Eigenvectors are orthonormal. 

74 

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. 

84 

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 

95 

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``. 

102 

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 

115 

116 mat = self.hop * complex(math.cos(k), -math.sin(k)) 

117 h0 = mat + mat.conjugate().transpose() + self.ham 

118 

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,) 

126 

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,) 

140 

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,) 

152 

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