1# Copyright 2011-2019 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 

9from keyword import iskeyword 

10import functools 

11import operator 

12import collections 

13 

14from math import sqrt 

15import tinyarray as ta 

16import numpy as np 

17import sympy 

18 

19import kwant.lattice 

20import kwant.builder 

21import kwant.system 

22 

23import kwant.continuum 

24import kwant.continuum._common 

25import kwant.continuum.discretizer 

26from kwant.continuum import momentum_operators, position_operators 

27 

28__all__ = ["to_landau_basis", "discretize_landau"] 

29 

30coordinate_vectors = dict(zip("xyz", np.eye(3))) 

31 

32ladder_lower, ladder_raise = sympy.symbols(r"a a^\dagger", commutative=False) 

33 

34 

35def to_landau_basis(hamiltonian, momenta=None): 

36 r"""Replace two momenta by Landau level ladder operators. 

37 

38 Replaces: 

39 

40 k_0 -> sqrt(B/2) * (a + a^\dagger) 

41 k_1 -> 1j * sqrt(B/2) * (a - a^\dagger) 

42 

43 Parameters 

44 ---------- 

45 hamiltonian : str or sympy expression 

46 The Hamiltonian to which to apply the Landau level transformation. 

47 momenta : sequence of str (optional) 

48 The momenta to replace with Landau level ladder operators. If not 

49 provided, 'k_x' and 'k_y' are used 

50 

51 Returns 

52 ------- 

53 hamiltonian : sympy expression 

54 momenta : sequence of sympy atoms 

55 The momentum operators that have been replaced by ladder operators. 

56 normal_coordinate : sympy atom 

57 The remaining position coordinate. May or may not be present 

58 in 'hamiltonian'. 

59 """ 

60 hamiltonian = kwant.continuum.sympify(hamiltonian) 

61 momenta = _normalize_momenta(momenta) 

62 normal_coordinate = _find_normal_coordinate(hamiltonian, momenta) 

63 

64 # Substitute ladder operators for Landau level momenta 

65 B = sympy.symbols("B") 

66 hamiltonian = hamiltonian.subs( 

67 { 

68 momenta[0]: sympy.sqrt(abs(B) / 2) * (ladder_raise + ladder_lower), 

69 momenta[1]: sympy.I 

70 * sympy.sqrt(abs(B) / 2) 

71 * (ladder_lower - ladder_raise), 

72 } 

73 ) 

74 

75 return hamiltonian, momenta, normal_coordinate 

76 

77 

78def discretize_landau(hamiltonian, N, momenta=None, grid_spacing=1): 

79 """Discretize a Hamiltonian in a basis of Landau levels. 

80 

81 Parameters 

82 ---------- 

83 hamiltonian : str or sympy expression 

84 N : positive integer 

85 The number of Landau levels in the basis. 

86 momenta : sequence of str (optional) 

87 The momenta defining the plane perpendicular to the magnetic field. 

88 If not provided, "k_x" and "k_y" are used. 

89 grid_spacing : float, default: 1 

90 The grid spacing to use when discretizing the normal direction 

91 (parallel to the magnetic field). 

92 

93 Returns 

94 ------- 

95 builder : `~kwant.builder.Builder` 

96 'hamiltonian' discretized in a basis of Landau levels in the plane 

97 defined by 'momenta'. If a third coordinate is present in 'hamiltonian', 

98 then this also has a translational symmetry in that coordinate direction. 

99 The builder has a parameter 'B' in addition to any other parameters 

100 present in the provided 'hamiltonian'. 

101 

102 Notes 

103 ----- 

104 The units of magnetic field are :math:`ϕ₀ / 2 π a²` with :math:`ϕ₀ = h/e` 

105 the magnetic flux quantum and :math:`a` the unit length. 

106 """ 

107 

108 if N <= 0: 

109 raise ValueError("N must be positive") 

110 

111 hamiltonian, momenta, normal_coordinate = to_landau_basis(hamiltonian, momenta) 

112 

113 # Discretize in normal direction and split terms for onsites/hoppings into 

114 # monomials in ladder operators. 

115 tb_hamiltonian, _ = kwant.continuum.discretize_symbolic( 

116 hamiltonian, coords=[normal_coordinate.name] 

117 ) 

118 tb_hamiltonian = { 

119 key: kwant.continuum._common.monomials(value, gens=(ladder_lower, ladder_raise)) 

120 for key, value in tb_hamiltonian.items() 

121 } 

122 

123 # Replace ladder operator monomials by tuple of integers: 

124 # e.g. a^\dagger a^2 -> (+1, -2). 

125 tb_hamiltonian = { 

126 outer_key: { 

127 _ladder_term(inner_key): inner_value 

128 for inner_key, inner_value in outer_value.items() 

129 } 

130 for outer_key, outer_value in tb_hamiltonian.items() 

131 } 

132 

133 # Construct map from LandauLattice HoppingKinds to a sequence of pairs 

134 # that encode the ladder operators and multiplying expression. 

135 tb_hoppings = collections.defaultdict(list) 

136 for outer_key, outer_value in tb_hamiltonian.items(): 

137 for inner_key, inner_value in outer_value.items(): 

138 tb_hoppings[(*outer_key, sum(inner_key))] += [(inner_key, inner_value)] 

139 # Extract the number of orbitals on each site/hopping 

140 random_element = next(iter(tb_hoppings.values()))[0][1] 

141 norbs = 1 if isinstance(random_element, sympy.Expr) else random_element.shape[0] 

142 tb_onsite = tb_hoppings.pop((0, 0), None) 

143 

144 # Construct Builder 

145 if _has_coordinate(normal_coordinate, hamiltonian): 

146 sym = kwant.lattice.TranslationalSymmetry([grid_spacing, 0]) 

147 else: 

148 sym = kwant.system.NoSymmetry() 

149 lat = LandauLattice(grid_spacing, norbs=norbs) 

150 syst = kwant.Builder(sym) 

151 

152 # Add onsites 

153 landau_sites = (lat(0, j) for j in range(N)) 

154 if tb_onsite is None: 

155 syst[landau_sites] = ta.zeros((norbs, norbs)) 

156 else: 

157 syst[landau_sites] = _builder_value( 

158 tb_onsite, normal_coordinate.name, grid_spacing, is_onsite=True 

159 ) 

160 

161 # Add zero hoppings between adjacent Landau levels. 

162 # Necessary to be able to use the Landau level builder 

163 # to populate another builder using builder.fill(). 

164 syst[kwant.builder.HoppingKind((0, 1), lat)] = ta.zeros((norbs, norbs)) 

165 

166 # Add the hoppings from the Hamiltonian 

167 for hopping, parts in tb_hoppings.items(): 

168 syst[kwant.builder.HoppingKind(hopping, lat)] = _builder_value( 

169 parts, normal_coordinate.name, grid_spacing, is_onsite=False 

170 ) 

171 

172 return syst 

173 

174 

175# This has to subclass lattice so that it will work with TranslationalSymmetry. 

176class LandauLattice(kwant.lattice.Monatomic): 

177 """ 

178 A `~kwant.lattice.Monatomic` lattice with a Landau level index per site. 

179 

180 Site tags (see `~kwant.system.SiteFamily`) are pairs of integers, where 

181 the first integer describes the real space position and the second the 

182 Landau level index. 

183 

184 The real space Bravais lattice is one dimensional, oriented parallel 

185 to the magnetic field. The Landau level index represents the harmonic 

186 oscillator basis states used for the Landau quantization in the plane 

187 normal to the field. 

188 

189 Parameters 

190 ---------- 

191 grid_spacing : float 

192 Real space lattice spacing (parallel to the magnetic field). 

193 offset : float (optional) 

194 Displacement of the lattice origin from the real space 

195 coordinates origin. 

196 """ 

197 

198 def __init__(self, grid_spacing, offset=None, name="", norbs=None): 

199 if offset is not None: 199 ↛ 200line 199 didn't jump to line 200, because the condition on line 199 was never true

200 offset = [offset, 0] 

201 # The second vector and second coordinate do not matter (they are 

202 # not used in pos()) 

203 super().__init__([[grid_spacing, 0], [0, 1]], offset, name, norbs) 

204 

205 def pos(self, tag): 

206 return ta.array((self.prim_vecs[0, 0] * tag[0] + self.offset[0],)) 

207 

208 def landau_index(self, tag): 

209 return tag[-1] 

210 

211 

212def _builder_value(terms, normal_coordinate, grid_spacing, is_onsite): 

213 """Construct an onsite/hopping value function from a list of terms 

214 

215 Parameters 

216 ---------- 

217 terms : list 

218 Each element is a pair (ladder_term, sympy expression/matrix). 

219 ladder_term is a tuple of integers that encodes a string of 

220 Landau raising/lowering operators and the sympy expression 

221 is the rest 

222 normal_coordinate : str 

223 grid_spacing : float 

224 The grid spacing in the normal direction 

225 is_onsite : bool 

226 True if we are constructing an onsite value function 

227 """ 

228 ladder_term_symbols = [sympy.Symbol(_ladder_term_name(lt)) for lt, _ in terms] 

229 # Construct a single expression from the terms, multiplying each part 

230 # by the placeholder that represents the prefactor from the ladder operator term. 

231 (ladder, (_, part)), *rest = zip(ladder_term_symbols, terms) 

232 expr = ladder * part 

233 for ladder, (_, part) in rest: 

234 expr += ladder * part 

235 expr = expr.subs( 

236 {kwant.continuum.discretizer._displacements[normal_coordinate]: grid_spacing} 

237 ) 

238 # Construct the return string and temporary variable names 

239 # for function calls. 

240 return_string, map_func_calls, const_symbols, _cache = kwant.continuum.discretizer._return_string( 

241 expr, coords=[normal_coordinate] 

242 ) 

243 

244 # Remove the ladder term placeholders, as these are not parameters 

245 const_symbols = set(const_symbols) 

246 for ladder_term in ladder_term_symbols: 

247 const_symbols.discard(ladder_term) 

248 

249 # Construct the argument part of signature. Arguments 

250 # consist of any constants and functions in the return string. 

251 arg_names = set.union( 

252 {s.name for s in const_symbols}, {str(k.func) for k in map_func_calls} 

253 ) 

254 arg_names.discard("Abs") # Abs function is not a parameter 

255 for arg_name in arg_names: 

256 if not (arg_name.isidentifier() and not iskeyword(arg_name)): 256 ↛ 257line 256 didn't jump to line 257, because the condition on line 256 was never true

257 raise ValueError( 

258 "Invalid name in used symbols: {}\n" 

259 "Names of symbols used in Hamiltonian " 

260 "must be valid Python identifiers and " 

261 "may not be keywords".format(arg_name) 

262 ) 

263 arg_names = ", ".join(sorted(arg_names)) 

264 # Construct site part of the function signature 

265 site_string = "from_site" if is_onsite else "to_site, from_site" 

266 # Construct function signature 

267 if arg_names: 

268 function_header = "def _({}, {}):".format(site_string, arg_names) 

269 else: 

270 function_header = "def _({}):".format(site_string) 

271 # Construct function body 

272 function_body = [] 

273 if "B" not in arg_names: 

274 # B is not a parameter for terms with no ladder operators but we still 

275 # need something to pass to _evaluate_ladder_term 

276 function_body.append("B = +1") 

277 function_body.extend( 

278 [ 

279 "{}, = from_site.pos".format(normal_coordinate), 

280 "_ll_index = from_site.family.landau_index(from_site.tag)", 

281 ] 

282 ) 

283 # To get the correct hopping if B < 0, we need to set the Hermitian 

284 # conjugate of the ladder operator matrix element, which swaps the 

285 # from_site and to_site Landau level indices. 

286 if not is_onsite: 

287 function_body.extend( 

288 ["if B < 0:", 

289 " _ll_index = to_site.family.landau_index(to_site.tag)" 

290 ] 

291 ) 

292 function_body.extend( 

293 "{} = _evaluate_ladder_term({}, _ll_index, B)".format(symb.name, lt) 

294 for symb, (lt, _) in zip(ladder_term_symbols, terms) 

295 ) 

296 function_body.extend( 

297 "{} = {}".format(v, kwant.continuum.discretizer._print_sympy(k)) 

298 for k, v in map_func_calls.items() 

299 ) 

300 function_body.append(return_string) 

301 func_code = "\n ".join([function_header] + function_body) 

302 # Add "I" to namespace just in case sympy again would miss to replace it 

303 # with Python's 1j as it was the case with SymPy 1.2 when "I" was argument 

304 # of some function. 

305 namespace = { 

306 "pi": np.pi, 

307 "I": 1j, 

308 "_evaluate_ladder_term": _evaluate_ladder_term, 

309 "Abs": abs, 

310 } 

311 namespace.update(_cache) 

312 # Construct full source, including cached arrays 

313 source = [] 

314 for k, v in _cache.items(): 

315 source.append("{} = (\n{})\n".format(k, repr(np.array(v)))) 

316 source.append(func_code) 

317 exec(func_code, namespace) 

318 f = namespace["_"] 

319 f._source = "".join(source) 

320 

321 return f 

322 

323 

324def _ladder_term(operator_string): 

325 r"""Return a tuple of integers representing a string of ladder operators 

326 

327 Parameters 

328 ---------- 

329 operator_string : Sympy expression 

330 Monomial in ladder operators, e.g. a^\dagger a^2 a^\dagger. 

331 

332 Returns 

333 ------- 

334 ladder_term : tuple of int 

335 e.g. a^\dagger a^2 -> (+1, -2) 

336 """ 

337 ret = [] 

338 for factor in operator_string.as_ordered_factors(): 

339 ladder_op, exponent = factor.as_base_exp() 

340 if ladder_op == ladder_lower: 

341 sign = -1 

342 elif ladder_op == ladder_raise: 

343 sign = +1 

344 else: 

345 sign = 0 

346 ret.append(sign * int(exponent)) 

347 return tuple(ret) 

348 

349 

350def _ladder_term_name(ladder_term): 

351 """ 

352 Parameters 

353 ---------- 

354 ladder_term : tuple of int 

355 

356 Returns 

357 ------- 

358 ladder_term_name : str 

359 """ 

360 

361 def ns(i): 

362 if i >= 0: 

363 return str(i) 

364 else: 

365 return "_" + str(-i) 

366 

367 return "_ladder_{}".format("_".join(map(ns, ladder_term))) 

368 

369 

370def _evaluate_ladder_term(ladder_term, n, B): 

371 r"""Evaluates the prefactor for a ladder operator on a landau level. 

372 

373 Example: a^\dagger a^2 -> (n - 1) * sqrt(n) 

374 

375 Parameters 

376 ---------- 

377 ladder_term : tuple of int 

378 Represents a string of ladder operators. Positive 

379 integers represent powers of the raising operator, 

380 negative integers powers of the lowering operator. 

381 n : non-negative int 

382 Landau level index on which to act with ladder_term. 

383 B : float 

384 Magnetic field with sign 

385 

386 Returns 

387 ------- 

388 ladder_term_prefactor : float 

389 """ 

390 assert n >= 0 

391 # For negative B we swap a and a^\dagger. 

392 if B < 0: 

393 ladder_term = tuple(-i for i in ladder_term) 

394 ret = 1 

395 for m in reversed(ladder_term): 

396 if m > 0: 

397 factors = range(n + 1, n + m + 1) 

398 elif m < 0: 

399 factors = range(n + m + 1, n + 1) 

400 if n == 0: 

401 return 0 # a|0> = 0 

402 else: 

403 factors = (1,) 

404 ret *= sqrt(functools.reduce(operator.mul, factors)) 

405 n += m 

406 return ret 

407 

408 

409def _normalize_momenta(momenta=None): 

410 """Return Landau level momenta as Sympy atoms 

411 

412 Parameters 

413 ---------- 

414 momenta : None or pair of int or pair of str 

415 The momenta to choose. If None then 'k_x' and 'k_y' 

416 are chosen. If integers, then these are the indices 

417 of the momenta: 0 → k_x, 1 → k_y, 2 → k_z. If strings, 

418 then these name the momenta. 

419 

420 Returns 

421 ------- 

422 The specified momenta as sympy atoms. 

423 """ 

424 

425 # Specify which momenta to substitute for the Landau level basis. 

426 if momenta is None: 

427 # Use k_x and k_y by default 

428 momenta = momentum_operators[:2] 

429 else: 

430 if len(momenta) != 2: 430 ↛ 431line 430 didn't jump to line 431, because the condition on line 430 was never true

431 raise ValueError("Two momenta must be specified.") 

432 

433 k_names = [k.name for k in momentum_operators] 

434 

435 if all([type(i) is int for i in momenta]) and all( 435 ↛ 438line 435 didn't jump to line 438, because the condition on line 435 was never true

436 [i >= 0 and i < 3 for i in momenta] 

437 ): 

438 momenta = [momentum_operators[i] for i in momenta] 

439 elif all([isinstance(momentum, str) for momentum in momenta]) and all( 439 ↛ 446line 439 didn't jump to line 446, because the condition on line 439 was never false

440 [momentum in k_names for momentum in momenta] 

441 ): 

442 momenta = [ 

443 momentum_operators[k_names.index(momentum)] for momentum in momenta 

444 ] 

445 else: 

446 raise ValueError("Momenta must all be integers or strings.") 

447 

448 return tuple(momenta) 

449 

450 

451def _find_normal_coordinate(hamiltonian, momenta): 

452 discrete_momentum = next( 452 ↛ exitline 452 didn't finish the generator expression on line 452

453 momentum for momentum in momentum_operators if momentum not in momenta 

454 ) 

455 normal_coordinate = position_operators[momentum_operators.index(discrete_momentum)] 

456 return normal_coordinate 

457 

458 

459def _has_coordinate(coord, expr): 

460 momentum = momentum_operators[position_operators.index(coord)] 

461 atoms = set(expr.atoms()) 

462 return coord in atoms or momentum in atoms